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

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

（均以宽角法数据为例子，剖面法数据处理方式都是相同的）假设，我们现在已经获得了一个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.

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mzph.cn/bicheng/28118.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！