基于GprMax 的地下冬笋正演模拟研究

余 铖，周凯琪，于 航，温春雨，樊梦成，杨科成

（浙江农林大学光机电工程学院，浙江 杭州 311300）

0 引言

我国是世界上竹林资源最丰富的国家[1-2]，2022年竹林面积达701万hm2，冬笋[3]是竹林收益的主要来源之一。由于冬笋埋藏于地下，目前多依赖于农民的经验采收，效率低下，盲目挖掘也会破坏植物根茎。

探地雷达[4]（Ground Penetrating Radar, GPR）是指利用收发天线发出并接收高频电磁波，来探测地下物质分布规律的一种无损检测方式，具有精度高、效率高、检测特性平稳、对环境适应性好等优势，被广泛应用于植物根系探测领域。利用探地雷达探测冬笋是一种比较好的方法，但是探地雷达的缺陷就在于，要求研究人员必须具备解译其信息的知识，也因为探地雷达信息的复杂性，解译探地雷达信息时往往要求具有丰富的专业知识和经验[5]。尤其是地下环境复杂，使探地雷达信息的复杂度更高。

通过正演模拟可以高效积累雷达图像解译的经验，减少实测成本。本文基于GprMax[6-7]仿真软件，以地下冬笋为对象，研究了不同干扰物、不同大小冬笋与不同深度冬笋的探地雷达正演模型，总结了冬笋的雷达探测经验。

1 探地雷达及正演原理

1.1 探地雷达原理介绍

探地雷达主要由雷达主机、显示单元、接收发天线组成，是利用10 MHz 以上的高频电磁波束对地下介质成像的无损探测技术。探地雷达的工作原理如图1 所示，工作时通过发射天线（TX）向地下发射高频电磁波，当遇到具有电性差异的界面时发生反射，由接收天线（RX）收集反射波，并经采样电路采集信号传输至雷达主机，计算出反射波的相位信息，拟合成探地雷达的二维剖面图。

图1 探地雷达工作原理

1.2 正演仿真原理及GprMax

在宏观尺度上，所有的电磁现象都可由麦克斯韦方程组描述。以下为麦克斯韦方程组描述基本电磁场能量与其对源的影响情况的一阶偏微分方程：

式中，E、H、B、D分别是电场强度、磁场强度、磁场感应强度和电位移矢量；J是电流密度；q是电荷密度；t为时间。从模型的几何形式和初始条件出发，通过利用麦克斯韦方程组，可以进行探地雷达数据模拟，但是在模拟仿真过程中，麦克斯韦方程连续性的特点，使得它在计算机中永远无法求解。

Y e e[8]在1 9 9 6 年提出了时域有限差分法（FDTD），通过利用离散空间与时间的连续性，并利用在时域的递推模拟了波的扩散行为，由此得到了场分布图。Yee 网络如图2 所示，是离散型FDTD 网格的代表形式。

图2 FDTD Yee网络

此外，在进行电磁场仿真模拟时，若在空间无限远时将场传递至零值，则仍不能利用有限的计算空间来求解无限的仿真问题。因此，在距离发射源足够远的地方应用完全匹配层（PML）边界吸收条件对电磁场计算进行截断，从而限制计算空间，边界吸收条件通过吸收传播过来的场，模拟出一个无边界的空间。

GprMax 是一个可以进行电磁波传输模拟的开源软件，在采用了时域有限差分法（FDTD）和PML[9]（Perfectly Matched Layer）的边界吸收条件下，以Yee元胞为电磁场的离散单元，将Maxwell 旋度微分方程转化为差分方程，以达到正演模拟的效果。因在一定时间上迭代求解，具有很高的可靠性和收敛性。GprMax 设计之初主要用作探地雷达（GPR）的模型，同时也可能用作一些其他的应用模型。GprMax 还具有内置的天线模型库，可以进行粗糙地基的模拟搭建、复杂各向异性结构的模拟等，可以适应大部分环境的搭建要求。

2 实验研究与讨论

冬笋通常和竹子的根系一起生长，周围难免会有石头等介质的干扰，使用探地雷达探测地下冬笋时，需要分辨出这些干扰信息。为了在不破坏生态的情况下，挖掘出合适生长时期的冬笋，还需判断出冬笋的埋藏深度与大小。

本实验通过GprMax 模拟了不同形状的干扰物，包括圆形截面的根系、矩形截面的石头，同时模拟了不同深度、不同大小的冬笋的正演模拟图像。模型参数如表1 所示[10]。

表1 模型参数

2.1 不同形状介质正演模拟

首先，探究不同形状、相同埋藏深度的模型的探地雷达正演模拟，设计模型整体宽1 m，深0.5 m。由于在GprMax 的运行中，发射天线和接收天线都需要占据模型中一定的空间，为了给发射天线和接收天线预留位置，介质层的厚度设置为0.1 m，空间步长设置为0.1 m。探测深度在0.05 m～0.15 m之间，使用Ricker 子波源模拟Gpr 的天线，激发频率选择中心频率为1 500 MHz。

首先对不同介质模型进行正演模拟，如图3 所示，模型a 为圆形截面的根系模型，模型b 为矩形截面的岩石模型，模型c为冬笋模型。

图3 不同介质GprMax模型构建示意图

不同形状介质的GprMax 正演模拟结果如图4 所示。从形状类别上来看，可以根据图像的双曲线特征来判断介质的形状，圆形介质的B 扫描图像平滑、规整且没有交叠的现象；而矩形介质的B 扫描图像顶端出现了交叠和颜色加深现象，在交点处有绕射；冬笋介质的内部有多次反射，且内部的反射波杂乱无章。以此特征可作为判断冬笋是否存在的依据。

图4 不同形状介质的GprMax正演模拟结果

2.2 不同大小冬笋正演模拟

与上述模型相同，构建土壤与空气层模型，然后设置埋藏深度相同、不同大小的冬笋，如图5 所示，从左至右分别为长度15 cm、10 cm 与5 cm 的冬笋回波示意图，埋藏深度均为5 cm。

图5 不同大小冬笋GprMax模型构建示意图

不同大小冬笋的正演模拟结果如图6 所示，从图中可以看出，冬笋的上下界面反射出多道回波，冬笋的长度越长，则上下界面的双曲线间隔距离越大。因此，可以通过两条双曲线出现的时间间隔估算冬笋的纵向尺寸。

图6 不同大小冬笋的正演模拟结果

2.3 不同深度冬笋正演模拟

最后构建大小相同、不同埋藏深度的冬笋模型，如图7 所示，从左至右分别为埋藏深度5 cm、10 cm、15 cm的冬笋，冬笋长度均设置为15 cm。

图7 不同埋藏深度冬笋GprMax模型构建示意图

相同长度、不同埋藏深度冬笋的GprMax 正演模拟结果如图8 所示，从回波数据来看，双曲线的顶点位置反映出冬笋的埋藏深度，通过分析此特征可以判断出介质的埋藏深度。冬笋的埋藏深度越大，双曲线的出现时间则越晚，随着深度增加，双曲线的开口大小也有所增加。因此，埋藏越深的冬笋，双曲线出现时间越晚、开口越大，埋藏越浅的冬笋则相反。

图8 不同埋藏深度冬笋的GprMax正演模拟结果

3 结束语

为给探地雷达在冬笋探测领域的应用提供经验，本文利用GprMax 软件对不同形状干扰介质、不同大小及不同埋藏深度的冬笋等模型进行正演模拟，并进行了模型间的对比。结果表明，正演模拟的图像能较好地反映各不同模型的成像特征，能有效反馈出介质模型的物理特性如大小、形状、深度等信息。截面形状不同的物体的双曲线回波形状不同，双曲线的上下界面距离反映出物体的大小，顶点位置则反映出介质的埋藏深度。通过分析这些特征可以有效探测识别出地下冬笋的存在、生长时期（大小）与埋藏深度。本文所应用的GprMax 正演模拟方法有效揭示了探地雷达的图像规律，也证实了探地雷达在冬笋探测领域的可行性。