孙 伟 徐 青 蓝朝桢 李建胜

(1.信息工程大学测绘学院，河南 郑州 450052；2.苏州市数字城市工程研究中心，江苏 苏州 215021)

引 言

探地雷达(GPR)是一种新型的地下目标探测手段[1]。它利用电磁波对地表的穿透能力，从地表上向地下发射脉冲电磁波，电磁波在地下阻抗特性发生变化的界面上发生反射，通过接收反射回波信号，根据其时延、形状及频谱特性等参数，可以解译出目标深度和介质特性。在地下浅层目标的探测和识别方面，除了进行探测设备的研制和改进外，理论上的模拟计算也是这一领域的重要内容之一。GPR的数值模拟是分析探测问题，研究电磁波在介质中传播规律的有效手段，对于提高探测的效果和解释的准确性具有重要意义。此外，利用正演结果进行反演，可验证反演算法的正确性，提高解释精度。目前，国内外学者开展了大量的探地雷达数值模拟工作，取得了很大的进展[2-8]。时域有限差分方法(FDTD)已经成为主流算法被广泛采用[9]，目前正朝着大规模并行计算方向以及快速算法发展以解决三维GPR数值模拟计算量大的问题[10]。数值模拟的结果需要与实测结果进行对比才能判断数值模拟的正确性和可靠性，两者的对比研究不仅在理论上而且在工程应用上都具有重要意义。然而令人遗憾的是，已经公开发表的文献较少涉及两者的对比研究。如何利用数值模拟结果更好地为提高GPR探测效果以及解释的准确性服务变得越来越重要。下面从数值模拟和实验室实测两个方面入手来进行对比研究。

1. 算法原理

GPR技术的应用领域逐渐由低损耗介质向高损耗介质拓展，如堤坝、水体等，而介质的导电性会造成电磁波能量衰减，这是高频GPR只能进行浅层探测的原因。在早期探地雷达的资料处理和解释中,曾不考虑导电率的作用，随着探测精度和要求的提高，介质的导电性越来越重要。含衰减的GPR电磁波传播基本规律可用经典的Maxwell方程来描述

(1)

其中：E为电场强度，单位为V/m；H为磁场强度，单位为A/m；ε为介电常数，单位为F/m；μ为磁导率，单位为H/m；J是电流密度，单位为A/m2；Jm为磁流密度，单位为V/m2；σ为电导率，单位为S/m；σm为等效磁阻率，单位为Ω/m.

在二维剖面情况下，电磁波有TE和TM两种模式，分别对应电场矢量平行和垂直于电磁波入射平面。其中TM模的Maxwell方程如下

(2)

采用FDTD来数值解上述方程组，在TM模式下，把电场安排在整数网格上，把磁场放置在半网格上，离散化方程为

(3)

上式中：Δx，Δy，Δt分别是x，y方向的网格步长以及时间步长，需满足如下Courant稳定性条件

(4)

其中,Vmax是计算域的最大波速。有限的计算空间需要增加吸收层才能模拟电磁波在无限空间传播，采用8层完全匹配层(PML)[11]作为边界条件来吸收外向电磁波。

电磁波激励源采用超宽带Ricker子波

A(t)=A0[1-2(πfmt)]e-(πfmt)2

(5)

式中：fm是超宽带激励源中心频率；A0是脉冲输入最大值。收发天线间距为0.

2. 实验情况

图1给出了GPR探测的实验装置，在一个57 cm×43 cm×34 cm的箱体内充满水，箱体壁厚约5 mm，金属或聚氯乙烯(PVC)管线插入水中。使用仪器为意大利IDS公司生产的RIS型探地雷达，天线中心频率1.6 GHz.测线布置在箱体的外侧面中心水平方向，GPR天线紧贴着箱壁沿着测线从右向左水平扫描。之所以采取侧面扫描，而不采用水面朝下扫描的方式，主要是为了实验的简便以及GPR天线的定位。管线采用较有代表性的金属管线和PVC管线两种。由于天线频率较高，水中电磁波的衰减较大，因此，探测深度有限，管线放置的位置距离前壁300 mm以内。

图1 GPR探测实验装置

3. 结果分析与讨论

实验装置中PVC管线的壁厚仅仅 1.5 mm，介电常数约为3，内外都充满介电常数为81的水。为了建立薄壁PVC管线模型，必须用很小的网格尺寸。设置网格步长为0.2 mm，计算空间总网格数为2 900×2 100.由于单次计算量较大，为了减少总的计算量，计算过程中仅均匀设置44个测点，每个测点间距1 cm.另外，计算模型包含了箱体前壁，但不包含侧壁和后壁，箱体材料的介电常数约为3.

3.1 淡水中的Ф50 mm金属和Ф50 mm PVC管线

淡水中介质损耗相对较小，电导率设为0.5 mS/m，磁导率为0.金属管线和PVC管线的直径均为50 mm，两种管线距离前壁为100 mm，两者间距170 mm，他们的实验布局和数值模拟模型如图2所示。

图2 实验与数值模拟的管线相对位置

(a) 实测 (b) 数值模拟图3 Ф50 mm金属和Ф50 mm PVC管线成像比较

图3给出了实验探测的雷达剖面和数值模拟成像的对比。从图中可以清晰辨别右边较亮的双曲线属于金属管线成像，而左边较暗的双曲线属于PVC管线成像。左右两张图的双曲线位置、形状参数基本一致。因为水箱前壁会部分反射电磁波，反射下行的电磁波又被管线散射后上行进入雷达天线，所以可以看到对于一个管线的雷达剖面是由许多个上下排列、强弱分布的双曲线共同组成。实测成像和数值模拟结果很好地证实了这个现象。在实际应用当中，类似的情况同样会出现，下行的电磁波在被探测物体与地表间来回反射，形成一系列明暗的双曲线。然而，该现象会对图像的解释产生干扰。由于GPR本身的信噪比以及环境因素，GPR剖面经过一系列后处理以后，多次反射产生的后续双曲线可能要强于第一个双曲线，或者金属的后续双曲线可能强于介质的第一个双曲线，正如图3的PVC管线成像所示。这对图像的解释能力提出了很高的要求，如何从雷达剖面里区分出金属管线和非金属管线以及它们的相对位置需要有丰富的经验。但要从图3左图中识别出管线类型以及每种管线的数量是困难的，数值模拟可以在一定程度上帮助解释人员完成这项任务。

3.2 淡水中的Ф75 mm PVC和Ф50 mm PVC管线

相同材料、不同直径的管线的识别在实际应用中很重要，这需要探地雷达有足够的横向分辨率。图4给出了直径分别为75 mm和50 mm PVC管线的实测和数值模拟雷达剖面。它们与前壁距离以及两者间距与3.1节相同。从图4中两图对比可以看到：对于不同的管径，双曲线的形状参数基本能区分出来。右图的双曲线比左图复杂得多，这主要是边界条件残留的反射造成的。和图3比较，可以看出相同材料的管线比不同材料的管线更容易从雷达剖面中识别出来，金属比非金属容易识别。

图5是将两个PVC管线到前壁的距离增大到200 mm后的雷达剖面。显然，距离越大，管线越不容易被探测到。由于水对电磁波的衰减，实际的GPR信噪比较低，双曲线变得模糊且不完整，已经较难从左图的实测图像中根据公式计算出管径。但是，右图管线的双曲线依然非常清晰，这是因为数值模拟没有噪音以及其他各种影响，尽管信号很微弱，但信噪比足够使得双曲线清晰、完整。从这也可以看出：尽管数值模拟能模拟介质衰减等各种情况，但缺乏GPR实际工作条件下的背景噪音，使得在有耗介质存在时的数值模拟结果与实测结果差别较大。

(a)实测 (b)数值模拟图4 Ф50 mm和Ф75 mm PVC管线成像比较

(a)实测 (b)数值模拟图5 Ф50 mm和Ф75 mm PVC管线成像比较

3.3 盐水中的Ф50 mm金属+Ф50 mm PVC 管线

GPR的探测深度与介质的导电率直接相关，3.2节的结果是在淡水条件下得到的，电导率很小。现在水中加入160 g盐，按照1 μS/cm=0.55～0.75 mg/l计算，电导率σ在0.2～0.3 S/m之间，取中间值σ=0.25 S/m作为电导率计算参数。

图6给出了实测所得的雷达剖面，金属和PVC管线直径均为50 mm，对应管线到前壁的距离(探测深度)分别从40 mm增大到150 mm.从图可以看出：金属管线的最大探测深度明显大于PVC管线，金属管线在120 mm处雷达剖面依然有较强的双曲线，而PVC管线在80 mm以上就无法成像。在较大的电导率的介质中，只有当管线在浅层(如图6(a))时才能在雷达剖面中观察到如前面所述的多个双曲线上下排列的现象；在当管线埋藏较深时，多次反射的信号基本上被介质吸收，这就是图6(b)～(f)一个管线在雷达剖面上对应一个双曲线的原因。

(a) h=40 mm (b) 60 mm

(c) h=80 mm (d) 100 mm

(e) h=120 mm (f) 150 mm图6 盐水中不同探测深度对应的实测雷达剖面图

图7是与图6一一对应的数值模拟雷达剖面图。根据前面的讨论，由于没有信噪比的限制，在管线与前壁任何间距下，数值模拟的雷达剖面都可以清晰辨别两个管线对应的两条双曲线。只是随着管线-前壁的距离增加，接收的电场强度减弱，雷达剖面上的双曲线变得暗淡。通过以下两种方法可以使数值模拟与实际相符：第一种方法是在数值模拟中对接收的电场强度设置某一个与实际噪音相关的临界值，当电场强度低于该临界值时表示电场强度太弱被探测物无法成像； 第二种方法是在激励Ricker子波中加入噪音。但受到FDTD数值色散的影响，第二种方法可能会受到限制。

(a) h=40 mm (b) 60 mm

(c) h=80 mm (d) 100 mm

(e) h=120 mm (f) 150 mm图7 盐水中不同探测深度对应的数值模拟图

4. 结 论

开展了中心频率为1.6 GHz的GPR水箱管线成像研究，并建立了二维的数值模拟程序和实验进行比较。针对有代表性的金属和PVC管线，在不同深度、不同管径组合以及介质损耗大小等情况下获得了实测的雷达剖面。相同模型下的数值模拟结果与实验符合得较好。讨论了一些以前较少关注的问题，如电磁波在被探物体和水箱前壁间多次反射形成的多条双曲线问题。对于在较强介质损耗情况下的实测与数值模拟不符现象，也给出了可能的解决办法。

[1] 蔚建斌, 陈自力, 江 涛. 探地雷达在地下介质中传播路径的计算[J]. 电波科学学报, 2009, 24(5): 925-933.

WEI Jianbin, CHEN Zili, JIANG Tao. Calculation of propagation path of GPR in media subsurface[J]. Chinese Journal of Radio Science, 2009, 24(5): 925-933. (in Chinese).

[2] BOURGEOIS J M, SIMITH G S. A fully three-dimensional simulation of a ground-penetrating radar: FDTD theory compared with experiment[J]. IEEE Trans Geosci Remote Sensing, 1996, 34(1): 36-44.

[3] CHEN Howwei, HUANG Taimin. Finite-difference time-domain simulationof GPR data[J]. Journal of Applied Geophysics, 1998, 40(1/3): 139-163.

[4] TEIXEIRA F L, WENG C C, STRAKA M, et al. Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils[J]. IEEE Trans Geosci Remote Sensing, 1998, 36(6): 1928-1937.

[5] 周 峰, 潘和平, 文国军, 等. 基于有限差分的井中激发极化三维正演[J]. 电波科学学报, 2010, 25(4): 785-792.

ZHOU Feng, PAN Heping, WEN Guojun, et al. Three-dimensional forward modeling of induced of polarization borehole using finite difference method [J]. Chinese Journal of Radio Science, 2010, 25(4): 785-792. (in Chinese)

[6] ROBERTS R L, DANIELS J. Finite difference time domain (FDTD) forward modeling of GPR data[C]// 5th Int Conf Ground Pentrating Radar. Waterloo, 1994: 185-204.

[7] 郑奎松, 葛德彪, 葛 宁. 三维电磁散射的网络并行FDTD计算和加速比分析[J]. 电波科学学报, 2004, 19(6): 767-771.

ZHENG Kuisong, GE Debiao, GE Ning. Parallel FDTD computing for 3D EM scattering problem and its speedup factor analysis[J]. Chinese Journal of Radio Science, 2004, 19(6): 767-771. (in Chinese).

[8] 冯德山, 戴前伟, 何继善. 探地雷达的正演模拟及有限差分波动方程偏移处理[J]. 中南大学学报:自然科学版, 2006, 37(2): 361-365.

FENG Deshan, DAI Qianwei, HE Jishan. Forward smulation of ground penetrating radar and its finite difference method wave equation migration processing[J]. Journal of Central South University: Science and Technology, 2006, 37(2): 361-365. (in Chinese).

[9] YEE K S. Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media [J]. IEEE Trans Antennas Propagat, 1966, 14(3): 302-307.

[10] 冯德山. 基于小波多分辨探地雷达正演及偏移处理研究[D]. 长沙: 中南大学信息物理工程学院, 2006.

FENG Deshan. Research Forward Smulation and Migration Processing of GPR Based on the Wavelet Multi-Resolution[D]. Changsha: Central South University, 2006. (in Chinese).

[11] BERENGER J P. A perfectly matched layer for the absorption of electromagnetic waves [J]. Journal of Computational Physics, 1994, 114 (2): 185-200.