如何用gprMax实现对MATLAB自定义发射信号的仿真

如何用gprMax实现对MATLAB自定义发射信号的仿真

gprMax支持用户自定义信号的输入，但是在更加复杂的环境下，输入信号的波形可能会经常性地发生改变，而若是使用gprMax的自定义波形功能，则需要在每一次重新定义波形以后，重新进行建模处理，这一操作在B-SCAN或3D建模的情况下会消耗大量的时间。

本文中的大部分内容为作者个人的理解，由于刚刚接触GPR，一些理解有可能是错的，读者请自行取舍。

本人不常用CSDN，如有相关问题，可发送电子邮件至2538703415@qq.com

在不重新执行建模的前提下对施加给模型的信号进行快速修改，理论上来说可以做到，思路大致如下：

基本原理

将线性时不变系统(LTI)求解响应的原理引入

对于输入信号x(t)，若已知系统的冲激响应h(t)，则能够得到对于输入信号x(t)的响应y(t):

$$y(t)=x(t)\ast h(t)$$

在gprMax中，开发者在波形选项中给出了冲激信号这一选项，其相关函数位于gprMax文件夹下的waveform.py中：

.
.
.

 elif self.type == 'impulse':
     # time < dt condition required to do impulsive magnetic dipole
     if time == 0 or time < dt:
         ampvalue = 1
     elif time >= dt:
         ampvalue = 0

 elif self.type == 'user':
     ampvalue = self.userfunc(time)

 ampvalue *= self.amp

 .
 .
 .

A-SCAN模式下在实际求解冲激响应模型时，可以直接将.in文件中原本的waveform类型直接改成impulse，其他地方可以不修改

#waveform: impulse 1 1.5e9 my_ricker

在执行仿真的时候可能会出现警告，但不影响结果（暂时没因为此警告出现问题）

---

仿真步骤

在开始前要留意目前为止本人认为冲激信号.in文件中唯一要关注的一个地方：

#time_window: 6e-9

在求解冲激响应模型时，这个时间窗的大小务必按照你最后希望得到的时间窗进行设置（可以设的更大，但是绝对不能比期望值更小）

假设你想要观测10ns里模型接收信号的表现，那么这个时间窗的长度就不能小于10ns

下面是具体的仿真步骤：

1. 用gprMax执行设置好的xxx.in文件：

   文件如下：

    #title: A-scan from a metal cylinder buried in a dielectric half-space
    #domain: 0.240 0.210 0.002
    #dx_dy_dz: 0.002 0.002 0.002
    #time_window: 6e-9
   
    #material: 6 0 1 0 half_space
   
    #waveform: impulse 1 1.5e9 my_ricker
    #hertzian_dipole: z 0.100 0.170 0 my_ricker
    #rx: 0.140 0.170 0
   
    #box: 0 0 0 0.240 0.170 0.002 half_space
    #cylinder: 0.120 0.080 0 0.120 0.080 0.002 0.010 pec
   
    #geometry_view: 0 0 0 0.240 0.210 0.002 0.002 0.002 0.002 cylinder_half_space n
   

   执行指令如下：

    python -m gprMax xxx.in
   

2. 将冲激响应数据导入MATLAB

   执行之后会得到一个.out文件，这个是hdf5格式的文件，可以用matlab直接读，具体怎么找数据的路径不在此赘述，本人的data通过此种方式获得：

   data = h5read('xxx.out', '/rxs/rx1/Ez');
   % 想看结果的话可以直接对data进行plot操作
   % plot(data);
   

3. 以gprMax官方给出的高斯信号为例进行验证

   官方给出的Gaussian信号定义如下：

    .
    .
    .
    if self.type == 'gaussian' or self.type == 'gaussiandot' or self.type == 'gaussiandotnorm'  or self.type == 'gaussianprime' or self.type == 'gaussiandoubleprime':
        self.chi = 1 / self.freq
        self.zeta = 2 * np.pi**2 * self.freq**2
    .
    .
    .
    # Waveforms
    if self.type == 'gaussian':
        delay = time - self.chi
        ampvalue = np.exp(-self.zeta * delay**2)
    .
    .
    .
   

   根据python代码在MATLAB中还原此高斯信号：

   freq = 1.5e9;
   chi = 1/freq;
   
   % simulate a gaussian pulse (gprMax source code)
   t = 1: length(data);
   t = t*6e-9/length(data);
   a = 2*pi^2*freq^2;
   x = exp(-a*(t-chi).^2);
   

   使用x与data进行卷积操作

   % get convolution result
   res = conv(data,x, 'full');
   

   将各种结果打印出来

   subplot(3,1,1)
   plot(t, data);
   title('impulse response');
   xlabel('time (s)');
   ylabel('field strength (V/m)');
   
   subplot(3,1,2)
   plot(t, x);
   title('gaussian pulse');
   xlabel('time (s)');
   ylabel('field strength (V/m)');
   
   subplot(3,1,3)
   axisx = (1:length(res))*6e-9/length(data);
   plot(axisx(1:length(data)), res(1:length(data)));
   title('convolution result');
   xlabel('time (s)');
   ylabel('field strength (V/m)');
   

   得到最终的输出结果，并和gprMax直接设置Gaussian waveform得到的结果进行比较

   下图是MATLAB做卷积运算得到的结果
   

   下图是gprMax直接仿真得到的结果（Ez部分为左下角）：
   

两者得到的结果一致

注意事项

根据gprMax的官方手册可以知道，模型中dx,dy,dz的大小会被发射信号的频率以及介质的相对介电常数限制，其要求dxdydz的取值要小于模型中最小波长的十分之一，即：
$$dxdydz\le \frac{{\lambda }_{min}}{10}$$
这个值可以通过如下方式进行计算：
$$dxdydz=\frac{c}{10f_{max}\sqrt{{ϵ}_{max}}}$$
在模型建立时，要根据模型中介质的最大相对介电常数以及发射信号的最高频率，合理设置dxdydz的值，使dx，dy，dz的值不大于上面公式计算出的结果。
模型中dxdydz的数值如果设置不合理，同样会导致高频信号在与冲激响应卷积时产生错误结果。

如何使用冲激响应法对SFCW信号进行仿真