探地雷达正演模拟,基于时域有限差分方法,四

        突然发现第三章后半部分已经讲了使用接收记录成像的问题,所以这一章只讲解简单的数据分析。

        (均以宽角法数据为例子,剖面法数据处理方式都是相同的)假设,我们现在已经获得了一个GPR记录,可以是常用的.sgy格式,也可以是其他的格式,只要读出来是个二维矩阵文件就可以,主要要做的内容:提取三瞬:瞬时相位、振幅和频率,这节主要使用Python,因为C++的话,一般是要写相应的函数的,但是Python的numpy库已经为我们集成了这一章所有要用到的函数,用起来会很简单。

        首先,正演一套数据,可以转换成.sgy文件,但是没必要,正演完就是一个二维矩阵,直接操作就好。根据本章参考文献的讲解,三瞬的提取方法是对数据做希尔伯特变换(Hilbert transform),将实信号转换为复信号,再提取三瞬即可。

        简单讲一下希尔伯特变换,这个变换时非常简单的,用一个公式就可以解释清楚:

x_{img}(t)=x_{real}(t)*\frac{1}{\pi t}

signal(t)=x_{real}(t)+jx_{img}(t )

如果相对希尔伯特变换有更加深入的了解,可以参考CSDN中的相应博客和教科书,这里不做过多赘述。

在对数据进行希尔伯特变换之后,就得到了这个信号的的复信号,接下来根据参考文献中给出的公式,计算三瞬即可:

1、瞬时振幅:

A=\sqrt(x_{real}^2+x_{img}^2)

2、瞬时相位:

\theta=tan^{-1}\frac{x_{img}}{x_{real}}

3、瞬时频率:

F=\frac{d}{dt}\theta

现在,对正演资料进行处理:

# editor: WangBoHua
# date:2024.6.16
# A:瞬时振幅;angle:瞬时相位角;F:瞬时频率
import numpy as np
import cmath
import matplotlib.pyplot as plt
# Read Gpr Signal
# dt=0.001 #dt是实际记录或者接收记录的采样率
GprData=np.genfromtxt("C:/Users/12981/Desktop/signal.data")
trace=GprData[0].size
Samplelength=int(GprData.size/trace)
Hibert=np.complex128(GprData)
# calculate Hibert transform of each trace
t=np.linspace(1,Samplelength,Samplelength,dtype='float')#*dt
Trace=np.linspace(0,trace,trace,dtype='float')
# Time,Trace=np.meshgrid(t,Trace)
hibert=1/(np.pi*t)
for i in range(0,trace):
    v=np.convolve(GprData[:,i],hibert,'same')
    Hibert[:,i]=GprData[:,i]+cmath.sqrt(-1)*v
# calculate A:
Ap=np.abs(Hibert)
# calculate angle:
Angle=np.angle(Hibert)
# calculate frequency:
Frequency=np.diff(Angle)
# draw
fig=plt.figure(figsize=(16,9),dpi=800,layout="constrained")
# amplitude
ax1=fig.add_subplot(3,1,1)
a=ax1.pcolormesh(Trace,t,Ap, rasterized=True,cmap='jet',shading="nearest",vmin=0,vmax=80)
ax1.invert_yaxis()
# angle
ax2=fig.add_subplot(3,1,2)
b=ax2.pcolormesh(Trace,t,Angle, rasterized=True,cmap='jet',shading="nearest",vmin=-3,vmax=3)
ax2.invert_yaxis()
# frequency
ax3=fig.add_subplot(3,1,3)
Trace2=np.linspace(0,trace,trace-1,dtype='float')
c=ax3.pcolormesh(Trace2,t,Frequency, rasterized=True,cmap='jet',shading="nearest",vmin=-0.3,vmax=0.3)
ax3.invert_yaxis()
c1 = fig.colorbar(a, ax=ax1, shrink=1.0)
c2 = fig.colorbar(b, ax=ax2, shrink=1.0)
c3 = fig.colorbar(c, ax=ax3, shrink=1.0)
plt.savefig('sanshun.png')
fig2=plt.figure(figsize=(16,9),dpi=800,layout="constrained")
ax4=fig2.add_subplot(3,1,1)
ax4.plot(t,Ap)
ax4.grid(True,linestyle='-.')
ax5=fig2.add_subplot(3,1,2)
ax5.plot(t,Angle)
ax5.grid(True,linestyle='-.')
ax6=fig2.add_subplot(3,1,3)
ax6.plot(t,Frequency)
ax6.grid(True,linestyle='-.')
plt.savefig('sanshuncurve.png')
plt.show()

时间原因,就不放图片了,可以自己正演试一试,看看效果怎么样。

声明:欢迎使用上述代码,但请标注公式和图片来源,创作不易,请多多支持,非常感谢!

参考文献:

杜衍庆, et al."公路路基浅层疏松病害地质雷达正演模拟及复信号分析."北京交通大学学报 47.06(2023):147-153.

相关推荐

  1. 基础算法练习】前缀和与模板

    2024-06-16 13:26:02       34 阅读

最近更新

  1. TCP协议是安全的吗?

    2024-06-16 13:26:02       16 阅读
  2. 阿里云服务器执行yum,一直下载docker-ce-stable失败

    2024-06-16 13:26:02       16 阅读
  3. 【Python教程】压缩PDF文件大小

    2024-06-16 13:26:02       15 阅读
  4. 通过文章id递归查询所有评论(xml)

    2024-06-16 13:26:02       18 阅读

热门阅读

  1. gitlab问题记录

    2024-06-16 13:26:02       6 阅读
  2. C# —— 条件分支语句

    2024-06-16 13:26:02       8 阅读
  3. 洛谷题解 - P1036 [NOIP2002 普及组] 选数

    2024-06-16 13:26:02       7 阅读
  4. 深度神经网络

    2024-06-16 13:26:02       5 阅读
  5. ubantu 计算一个文件夹内的文件数量命令

    2024-06-16 13:26:02       9 阅读
  6. vue2和vue 3 的响应式原理

    2024-06-16 13:26:02       6 阅读
  7. 博客摘录「 YOLOv5模型剪枝压缩」2024年5月11日

    2024-06-16 13:26:02       8 阅读
  8. 主流排序算法——python

    2024-06-16 13:26:02       5 阅读
  9. make menuconfig 分析

    2024-06-16 13:26:02       7 阅读