基于gprMax以及Matlab仿真的SFCW雷达信号处理笔记

基于gprMax以及Matlab仿真的SFCW雷达信号处理笔记

GitHub仓库传送门

本人为初学者，且不常使用CSDN，如遇到问题可发送邮件至2538703415@qq.com讨论，这样我能及时看见，但不一定能帮你解决……

当前使用gprMax进行探地雷达仿真主要存在以下两个问题，它们分别是：

• gprMax无法直接对一段完整的SFCW信号进行仿真，若要通过gprMax直接实现SFCW信号的回波仿真，则要将一段完整的SFCW信号拆分成若干段频率单一的余弦信号分别进行仿真，对于拥有上百个频点的SFCW信号而言，需要进行上百次单独的仿真。
• 使用自定义波形或gprMax内置波形进行仿真时，仿真速度偏低且灵活性差，每次在用户需要更改波形时，都要重新开始仿真。

对于提到的这些问题，包括gprMax开发团队在内的很多研究团体都得出了一些解决方案，这份readme文档将对这些解决方案进行归纳总结。这一仓库中的代码最终实现如下结果：

图一：gprMax模型

图二：Matlab仿真结果

根目录文件	文件作用
main.m	运行此文件，执行matlab仿真
simu_single_step.m	运行此文件，执行matlab仿真（无图）
result.mat	对应bscan_impulse_merged.out的仿真结果
result_10m.mat	对应bscan_impulse_merged_10m.out的仿真结果

impulse_lib文件夹	文件作用
bscan_impulse_10m_merged.out	对应bscan_impulse_10m.in的bscan冲激响应输出结果
bscan_impulse_merged.out	对应bscan_impulse.in的bscan冲激响应输出结果

gprMax_source文件夹	文件作用
sfcw_impulse.in	模型1的输入文件，可用于生成vti，对应bscan子文件夹中的bscan_impulse.in
sfcw_impulse_10m.in	模型2的输入文件，可用于生成vti，对应bscan子文件夹中的bscan_impulse_10m.in
impulse_response_10m.vti	模型2（对应图一）的vti文件

---

冲激响应法求解回波模型

为了解决gprMax仿真灵活性差，耗时长的问题，可以采用信号与系统中的基础知识，通过将模拟场景视为一个线性时不变系统（LTI system），将输入信号（也就是用户自定义的发射信号）与这一LTI系统的冲激响应h(t)相卷积，从而得到对应于输入信号的模拟场景回波信号。

图三：冲激响应计算y(t)

在通过一次仿真求得系统的冲激响应之后，用户只需要在matlab中将输入信号与此冲激响应相卷积，就可以得到对应的回波模型。gprMax提供了冲激信号波形（impulse），但是没有在对应的参考文献中列出，此波形可以直接使用，在使用时，其他波形对应的中心频率项对impulse波形不再生效（毕竟冲激信号只有tn=0时有值）。

具体的操作方法不在此处详细列出，对应的生成系统冲激响应的文件在github中可以查看。

这一部分的具体操作可参考：

• https://blog.csdn.net/l_jsaphsj/article/details/127530336?spm=1001.2014.3001.5501
• https://www.youtube.com/watch?v=_FWNeqTr9nc&t=1387s （gprMax官方）

---

SFCW信号回波的仿真（gprMax & Matlab）

SFCW雷达每次发射和接收的信号为单频连续波信号，其频率变化可表示为：
$$f_n=f_0+(n-1)\Delta f$$
其中 $f_n$ 是当前（第n个）频率， $f_0$ 是起始频率， $\Delta f$ 为频率步长。

发射信号可以被表示为：
$$T_x(n)=sin(2\pi f_{n}t)$$
其中n为第n个频点。

假设信号在均匀介质中的传播速度为 $v$ ，则一个距离目标为 $R$ 的物体所反射的回波信号可以表示为：
$$R_x(n)=A_{n}sin[2\pi f_n(t-\frac{2R}{v})]$$
其中 $A_n$ 为回波幅度， $\frac{2R}{v}$ 是信号的双程走时。

第n个频率的相位变化可对应表示为：
$${\Phi }_n=2\pi f_n\frac{2R}{v}$$
接收到的回波信号由此可以表示为：
$$R_x(n)=A_{n}sin(2\pi f_{n}t-{\Phi }_n)$$

不同的频点对应着不同的相移，这些相位信息对于特定的某个频点而言是一个常数，但是将SFCW信号中的几十至几百甚至更多频点的相位信息结合之后，能够得到一个具有周期性的相位信号，将这个信号进行ifft变换，即可得到对应于距离的时域信息。

图四：雷达系统框图

图四中的I和Q分量可通过如下公式进行计算：

$$I_n=A_{n}sin(2\pi f_{n}t-{\Phi }_n)\ast sin(2\pi f_{n}t)=\frac{An}{2}sin(4\pi f_{n}t-{\Phi }_n)-\frac{A_n}{2}cos({\Phi }_n)$$
$$I_n=A_{n}sin(2\pi f_{n}t-{\Phi }_n)\ast cos(2\pi f_{n}t)=\frac{An}{2}sin(4\pi f_{n}t-{\Phi }_n)-\frac{A_n}{2}sin({\Phi }_n)$$

在通过低通滤波器之后，可以提取出I和Q对应的差频分量为：
$$I_n=-\frac{A_n}{2}cos({\Phi }_n)$$
$$Q_n=-\frac{A_n}{2}sin({\Phi }_n)$$

注意在经过LPF之后需要对两部分差频分量进行组合,即输入ADC的信号应当为如下形式：
$$AD{C}_{in}(n)=I_n+j{Q}_n$$

对于经过IFFT处理后得到的时域信号，其横坐标轴的长度由ADC采样点数决定，其横坐标轴的精度由SFCW信号的带宽决定，设SFCW信号的带宽为 $B$ ，则其横坐标轴可表示为：
$$t=0:\frac{1}{B}:\frac{N_t-1}{B}$$
其中$N_t$是ADC采样点数。

---

gprMax生成模型时的注意事项

本人设计的in文件中，包含一项time_step_stability_factor参数，这一参数的具体用法可以参考官方文件，其作用主要是将dt进行手动调整（比如设置为1e-11s）。设置这一参数的主要目的是使gprMax中的dt和Matlab中的采样时间间隔dt相同，以方便仿真（两个dt不相同会怎样这个我也没试过，应该也行问题不大）。

gprMax中网格步长参数的设置

根据官方给出的文档，gprmax模型中的网格步长不能大于模型介质中最短波长的十分之一,即：
$$cal{c}_{dxyz}\le \frac{1}{10}\frac{c}{f_{max}\sqrt{{ϵ}_{max}}}$$

这项参数的最大值可以通过func文件夹中的calc_dxyz函数进行快速计算。

拆分后的单频率信号波形设置

在使用冲激响应方法，通过卷积获得输出信号波形时，如果直接输入一个完整的单频正弦信号，会导致卷积结果中包含大量高频噪声，因此，在设置输入的单频正弦信号时，需要逐渐提升信号幅值，其信号波形如下所示：

图五：单频余弦信号

图五所示的信号在gprMax中对应的波形名称为contsine，其公式如下所示：

$$T_x(t)=\{\begin{array}{ll}k{f}_{n}tsin(2\pi f_{n}t)& k{f}_{n}t<1\\ & \\ sin(2\pi f_{n}t)& k{f}_{n}t\ge 1\end{array}$$

在gprMax的waveforms.py中， $k=0.25$