频散介质中探地雷达有限元法正演模拟
王洪华1，戴前伟1, 2
(1. 中南大学 有色金属成矿预测教育部重点实验室，湖南 长沙，410083；
2. 中南大学 地球科学与信息物理学院，湖南 长沙，410083)
摘要：为了更准确地认识探地雷达(GPR)波在频散介质中的传播规律，从Maxwell方程组出发，在频散介质相对介电常数与频率满足Debye关系的条件下，经傅里叶变换，导出时间域GPR波电场和磁场的波动方程。对GPR波的电场波动方程采用伽辽金有限元法，推导频散介质中探地雷达二维有限元方程，采用透射吸收边界有效减弱截断边界的超强反射波，使得截断边界处的反射波充分吸收。利用依上述方法原理编制的程序分别对均匀模型和两圆状模型进行正演计算，验证该算法的正确性和可行性。研究结果表明：均匀频散介质中GPR子波的衰减比在非频散介质中的衰减要快得多，且GPR子波的宽度随着接收器与发射源距离的增加而变大；频散介质中异常体的反射波宽度比非频散介质中的反射波宽度大。
关键词：频散介质；有限单元法；探地雷达；正演模拟
中图分类号：P631 文献标志码：A 文章编号：1672-7207(2014)03-0790-08
Finite element method GPR forward simulation in dispersive medium
WANG Honghua1, DAI Qianwei1, 2
(1. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Changsha 410083, China;
2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)
Abstract: In order to more accurately understand the propagation law of GPR (ground penetrating radar) wave in dispersive medium, the Fourier transform was adopted under the condition that relative dielectric constant in dispersive medium and frequency satisfied the Debye condition, and the wave equation of electromagnetic field was obtained in time domain based on the Maxwell equations. The wave equation of GPR electric field was taken as an example, and the 2-D finite element equation of GPR in dispersive medium was gotten by means of Galerkin finite element method. By combining transmitting absorbing boundary, the strong reflected wave on truncation boundary was effectively weakened and the reflected wave on truncation boundary was absorbed adequately. On the basis of the theory mentioned above, the program of its feasibility and effectiveness was tested by homogeneous and two-rounded models. The results show that compared with the modeling results of non-dispersive medium, the GPR wave in dispersive medium attenuates faster than that in non-dispersive medium, and the width of GPR wave in dispersive medium gets larger with the increase of the distance.Width of the reflected wave with the abnormal body in dispersive medium is larger than that of the reflected wave in non-dispersive medium.
Key words: dispersive medium; finite element method; GPR; forward modeling
探地雷达(ground penetrating radar, GPR)是一种根据高频脉冲电磁波在地下不同物性介质之间的反射及绕射等波动规律来探测地下结构和特性的地球物理方法[1]。经过多年发展，探地雷达已成为工程和环境地球物理领域中一种重要的无损检测方法，在质量监测[2]、水文地质调查[3]等检测方面得到了广泛应用。由于一般地下介质比空气具有更强的电磁能量衰减特性，加之地下介质的分布复杂，电磁波在地下介质中的传播要比空气中的传播复杂得多，因此，利用探地雷达正演模拟研究GPR波在地下介质中的传播，对于提高探测的精度和解释的准确性具有重要意义[4]。在众多的探地雷达正演模拟方法中，时域有限差分法(finite difference time domain, FDTD)因发展成熟和使用简单而被广泛采用。国内外众多学者[5-8]对FDTD探地雷达正演模拟进行了系统而又深入的研究，但它不适合复杂的物性分界面。有限单元法(finite element method，FEM)不需要计算内部边界条件，适用于物性复杂问题，求解过程规范化，具有广泛的适应性[9-10]。在以往的探地雷达有限元正演模拟研究中，地下介质的物性参数被认为是与频率无关的，然而，当介质具有频率依赖时，介质的频散特性将显著地改变GPR波的传播特征[4]。Qing等[11]研究了计算频散介质的伪谱时间域有限差分法(pseudospectral time domain, PSFD)。Teixeira等[12]采用三维FDTD模拟了GPR波在频散、各向异性土壤中的传播，并与实验的结果进行了对比分析。刘四新等[13]在采用Sullivan理论的基础上，采用FDTD模拟了GPR在频散介质中的传播特征，并提出了一种简便的吸收边界条件，避免了Berenger完全匹配层中场分裂过程。为了进一步认识探地雷达(GPR)波在频散介质中的传播规律，本文作者在频散介质的相对介电常数与频率满足Debye关系的条件下，采用FEM法对GPR波在频散介质中的传播规律进行研究。对均匀频散模型和双圆状模型进行了正演模拟，并与非频散介质的正演模拟结果进行对比，总结GPR波在频散介质中的传播规律。
1 频散介质中GPR波动方程推导
根据电磁波场理论，在频率域Maxwell方程组的微分形式可以表示为[14]：
(1)
(2)
(3)
(4)
式中：
为角频率；
为相对介电常数(F/m)，它是频率的函数；
为真空中的介电常数(F/m)；
为磁导率(H/m)；
为电导率(S/m)；E为电场强度(V/m)；H为磁场强度(A/m)；
为电荷密度(C/m3)。
根据Maxwell方程组和本构关系，对式(1)两边求旋度，得
(5)
将式(2)代入式(5)得
(6)
式(3)为电磁场的Helmholtz方程，描述了GPR波在地下介质中的传播规律。由矢量乘法法则：
(7)
假设
，和为坐标缓变函数，其空间导数可以忽略不计。将式(7)代入式(6)可得
(8)
考虑电阻率的影响，相对介电常数随频率变化的Debye公式[15]可表示为
(9)
其中：
和
分别为频率为无穷和0时的相对介电常数值；
称为弛张时间；j为虚数单位。
将式(9)代入式(8)，可得
(10)
对式(10)两边同乘
，经整理可得
(11)
根据傅里叶时频变换原理，式(8)可转换到时间域，可得GPR电场满足的波动方程为
(12)
同理，磁场H满足的波动方程为
(13)
在式(12)和 (13)中代入激励源，如电场源SE或磁场源SH，则有
(14)
(15)
式(14)和(15)即为时空域中GPR在频散介质中传播满足的波动方程。
2 频散介质中GPR波动方程有限元求解
求解GPR波动方程其实质是应用有限单元法，结合雷达波所要满足的初始、边界条件求解偏微分方程。下面以电场满足的波动方程式(11)为例，采用伽辽金法推导频散介质中GPR有限元方程[16]。
2.1 网格剖分
首先将求解的二维区域剖分成矩形单元，如图1所示。
图1 网格剖分及节点编号示意图
Fig. 1 Sketch map of mesh division and node number
2.2 线性插值
图2(a)所示为母单元示意图，图2(b)所示为子单元示意图，这2个单元间的坐标变换关系为
；
(16)
其中：x0和y0为子单元中点的坐标；a和b为子单元的2个边长。微分关系为
；
； 
(17)
单元形函数为双线性插值形函数，可写为
(18)
其中：
和
为点i(i=1, 2, 3, 4)的坐标。形函数的分量可写为
；
；
；
(19)
2.3 单元积分
根据伽辽金法，将(14)两边同时乘以
，并求积分，有
(20)
其中：
；
；
；N1，N2，N3和N4为单元上各节点的形函数；E1，E2，E3和E4为各节点的电场值；
为单元面积。
式(20)左边第1项为
(21)
其中：M1e矩阵的计算公式为
(22)
式(20)左边第2项为
(23)
其中：
矩阵的计算公式为
(24)
式(20)左边第3项为
(25)
其中：M2e矩阵的计算公式为
(26)
式(20)左边第4项为
(27)
其中：K1e矩阵的计算公式为
(28)
式(20)左边第5项为
(29)
其中：K2e矩阵的计算公式为
(30)
式(17)右边项为
(31)
2.4 总体合成
根据式(21)，(23)，(25)，(27)，(29)和(31)得到单元积分：
(32)
将单元列向量
，
和
分别扩展成全体节点的列向量E，
和
，将4×4的系数矩阵M1e，，M2e，K1e及K2e分别扩展成ND×ND的矩阵M1，
，M2，K1及K2 (其中ND为节点总数)，将列向量Se扩展成ND维列向量SE，即
(33)
由于
不恒等于0，所以有
(34)
式(34)即为频散介质中GPR波电场的二维有限元方程。同理，也可以推导磁场的二维有限元方程为
(35)
式中：M1和M2为质量矩阵；为阻尼矩阵；K1和K2为刚度矩阵；和
为时间的一次导数；和
为时间的二次导数。
2.5 解GPR有限元方程
解常微分方程组(34)或(35)前，要将吸收边界条件进行加载(下面小节将具体阐述)采用中心差分法解方程(34)式。把时间域[0, T]分为几个相等步长
(其中，n为计算的步长数)，将t时刻的微分方程记为
(36)
用中心差分法代替微分，并令
，式(36)可变为
(37)
其中：
的时间迭代格式为
(38)
式(37)和(38)即为GPR有限元正演模拟时间递推公式。由于零时刻的E0，
和
以及
时刻的
，
和
均为0，所以，根据式(37)可计算出
时刻的
，然后，依次递推，可得到各个时刻的电场E。令
，
则式(37)可简化为
形式的线性方程组。在求解波动方程(37)时，本文中选用不完全LU分解预处理的稳定双共轭梯度BICGSTAB(biconjugate gradient stabilized)线性方程组求解算法[17]。该算法是一种高效迭代法，特别适用于求解方程组条件数很大的病态线性方程组。此外，系数矩阵A是质量矩阵M1和阻尼矩阵的线性组合，可以采用单元集中质量矩阵求解，以提高求解速度。
图2 母单元和子单元示意图
Fig. 2 Sketch maps of parent element and sub-element
3 透射吸收边界条件
透射吸收边界采用沿截断边界外法向衰减的单程波方程作为边界上得吸收边界。以左边界为例，用右行波作用于左行入射的波，两者相互抵消后实现了左边界上无反射的吸收效果。对于右边界使用左行波方程，底边界使用上行波方程，即可以实现消除截断边界反射波。透射吸收边界条件公式为[18-19]
(39)
式中：L为微分算子。则在整个计算区域中，有
(40)
其中：Ei为第i个节点的电场强度。残余量R为
(41)
利用Galerkin余量法，二维方程可写为
(42)
其中：
为计算区域的面积；r为残余量R的向量表示，由式(41)代入式(42)得到。由于内边界的积分相互抵消，所以，只需求解区域外边界积分。在均匀各向同性介质中，根据Clear.bout推导的傍轴近似[20]，其下行波、左行波、右行波方程分别为
；
；
(43)
由法向导数的定义知：
(44)
其中：v为电磁波在媒质中的传播速度；
和
分别为x和y边界外法线的方向余弦。将式(44)代入式(42)，并利用高斯定理可得
(45)
式中：
，
和
分别表示左、右及底边界；上边界为自由边界。将式(43)给出的边界条件代入式(45)，可以得到吸收阻尼矩阵F：
(46)
将边界阻尼加入到GPR有限元方程中，得到如下方程：
(47)
式(47)即为结合透射吸收边界条件的频散介质中GPR有限元方程。
4 计算实例
依方法原理编制了频散介质中GPR波二维有限元法源程序，对均匀模型和双圆模型进行正演计算，并与非频散介质中的结果进行对比分析。
4.1 均匀频散介质模型
图3所示为1个长×宽为10.0 m×7.0 m的均匀频散介质模型，模型的上部为空气层，厚度为1.0 m，下部为均匀频散介质，其相对介电常数为
，
，电导率为=0.001 S/m。
弛张时间为2 ns。发射源T为Riker子波，放置在地表3.0 m处，接收器分别放置在3.0，4.0，5.0，6.0，7.0和8.0 m处，如图3所示。图4所示为均匀介质中GPR子波振幅值衰减曲线。由图4可见：从各道的振幅看，不管是频散介质还是非频散介质，其振幅随着接收器与发射源距离的增大而减小；GPR波的衰减在频散介质中比在非频散介质中要快，并且随着接收器与发射源距离的变大，GPR波衰减地更快。图5所示为GPR子波均匀介质中模拟所得的信号。为了显示方便，图5(a)和图5(b)分别利用图中信号最大振幅进行标准化。事实上，各道信号的最大振幅是不相同的，频散介质中各道最大振幅分别为3.510 0×10-34，0.248 0×10-34，0.084 3×10-34， 0.036 2×10-34，0.018 2×10-34，0.010 7×10-34；非频散介质中各道最大振幅分别为4.330×10-34 V/m，0.571×10-34，0.328× 10-34，0.224×10-34，0.165×10-34和0.128×10-34 V/m。从相位上看，随着距离的增大，频散介质的作用越来越明显，GPR子波变得越来越宽；而在非频散介质中，GPR子波的宽度则基本上没有发生大的变化，说明在频散介质中GPR波的分辨率已有所降低。
4.2 双圆模型的正演模拟
两圆模型示意图如图6所示。该模型长×宽为10.0 m×6.0 m，上部为空气层，厚度为1.0 m；背景介质的相对介电常数，；电导率
=0.001 S/m；在背景介质中埋有2个相同的圆状异常体，其相对介电常数
，
，电导率
=0.001 S/m；整个区域被边长为0.1 m的矩形单元剖分为100×60的网格空间；GPR波脉冲激励源的中心频率为100 MHz，采样时窗长度为80 ns，采样时间间隔为0.1 ns。
图3 均匀模型示意图
Fig. 3 Sketch map of Homogeneous model
图4 GPR子波振幅值衰减曲线图
Fig. 4 Sub-wave of GPR amplitude attenuation graph
图5 频散介质和非频散介质中GPR波的比较
Fig. 5 Comparison of GPR wave in dispersive medium and non-dispersive medium
图 6 两圆模型示意图
Fig. 6 Sketch map of two circles
应用依本文建立的模拟频散介质的FEM算法编制的程序对两圆模型进行正演模拟，并得到非频散介质的正演模拟结果。GPR在非频散介质和频散介质中的正演模拟Wiggle图分别如图7和图8所示。从图7和图8可见：2个圆状异常体所在位置出现了双曲线绕射波，只是双曲线弧形的宽度比圆实际的直径要大得多，但双曲线的弧顶能准确地对应于2个圆状异常体的顶点；边界处的反射波能量很小，说明采用透射吸收边界条件，有效地减弱了边界处反射波的影响。对比图7和图8可知：频散介质中GPR波的衰减比在非频散介质衰减要快；在频散介质中，GPR子波和异常体双曲线反射波的宽度比在频散介质中要宽，增大了近1倍。通过分析这些GPR图件，能更全面认识GPR在频散介质和非频散介质中的传播规律，有益于指导GPR资料的解释。
图7 非频散介质中正演模拟Wiggle剖面图
Fig. 7 Wiggle section for forward simulation in dispersive medium
图8 频散介质正演模拟Wiggle剖面图
Fig. 8 Wiggle section for forward simulation in non-dispersive medium
5 结论
(1) 从Maxwell方程组出发，在频散介质的相对介电常数与频率满足Debye关系的条件下，采用傅里叶变换，推导了时间域GPR波在频散介质中满足的波动方程；并给出了伽辽金有限单元法求解该泛函变分问题的详细解法。
(2) 透射边界条件能很好地压制截断边界处的反射波能量，有效地减弱了截断边界处反射波的影响。
(3) 均匀频散介质中GPR子波的衰减比在非频散介质中要快得多，且GPR子波的宽度随着接收器与发射源距离的增大而变大；频散介质中异常体的反射波宽度比非频散介质中的反射波宽度要大。
参考文献：
[1] 曾昭发, 刘四新. 探地雷达原理与应用[M]. 北京: 电子工业出版社, 2011: 168-169.
ZENG Zhaofa, LIU Sixin. Ground penetrating radar theory and applications[M]. Beijing: Publishing House of Electronics Industry, 2010: 168-169.
[2] 魏超, 肖国强, 王法刚. 地质雷达在混泥土质量检测中的应用[J]. 工程地球物理学报, 2004, 1(5): 447-451.
WEI Chao, XIAO Guoqiang, WANG Fafang. The application of geo-radar in concrete quality detection[J]. Chinese Journal of Engineering Geophysics, 2004, 1(5): 447-451.
[3] 郭高轩, 吴吉春. 应用GPR 获取多孔介质水力参数研究进展[J]. 河海大学学报(自然科学版), 2005, 33(1): 18-22.
GUO Gaoxuan, WU Jichun. Advances in estimation of hydrogeological parameters of porous media with ground penetrating radar[J]. Journal of Hohai University (Natural Sciences), 2005, 33(1): 18-22.
[4] 刘四新, 曾昭发, 徐波. 三维频散介质中地质雷达信号的FDTD数值模拟[J]. 吉林大学学报(地球科学版), 2006, 36(1): 123-127.
LIU Sixin, ZENG Zhaofa, XU Bo. FDTD simulation for ground penetrating radar signal in 3-Dimensional dispersive medium[J]. Journal of Jilin University (Earth Science Edition), 2006, 36(1): 123-127.
[5] Goodman D. Ground-penetrating radar simulation in engineering and archaeology[J]. Geophysics, 1994, 59(2): 224-232.
[6] FENG Deshan, CHEN Chenshen, DAI Qianwei. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD[J]. Chinese J Geophys, 2010, 53(10): 2484-2496.
[7] Irving J, Knight R. Numerical modeling of ground-penetrating radar in 2-D using MATLAB[J]. Computers & Geosciences, 2006, 32(9): 1247-1258.
[8] Giannopoulos A. Modelling ground penetrating radar by GprMax[J].Construction and Building Materials, 2005, 19(10): 755-762.
[9] 底青云, 王妙月. 雷达波有限元仿真模拟[J]. 地球物理学报, 1999, 42(6): 818-825.
DI Qingyun, WANG Miaoyue. 2D finite element modeling for radar wave[J]. Chinese Journal of Geophysics, 1999, 42(6): 818-825.
[10] DI Qingyun, WANG Miaoyue. Migration of ground-penetrating radar data with a finite-element method that considers attenuation and dispersion[J]. Geophysics, 2004, 69(2): 472-477.
[11] QING Huoliu, GUO Xinfan. Simulations of GPR in dispersive media using a frequency-dependent PSTD algorithm[J]. IEEE transactions on Geoscience and Remote Sensing, 1999, 37(5): 2317-2324.
[12] Teixeira F L, Chen W C, Straka M, et al. Finite-difference time_domain simulation of Ground Penetrating Radar on dispersive in homogeneous and conductive soils[J]. IEEE Transaction on Geoscience and Remote Sensing, 1998, 36(6): 1928-1937.
[13] 刘四新, 曾昭发. 频散介质中地质雷达波传播的数值模拟[J]. 地球物理学报, 2007, 50(1): 320-326.
LIU Sixin, ZENG Zhaofa. Numerical simulation for ground penetrating radar wave propagation in the dispersive medium[J]. Chinese Journal of Geophysics, 2007, 50(1): 320-326.
[14] Sullivan D M. Electromagnetic simulation using the FDTD method[M]. New York: IEEE Press, 2000: 67-68.
[15] Tong X, George A M. GPR attenuation and its numerical simulation in 2.5 dimensions[J]. Geophysics, 1997, 62(1): 403-414.
[16] 戴前伟, 王洪华, 冯德山, 等. 基于三角形剖分的复杂GPR模型有限元法正演模拟[J]. 中南大学学报(自然科学版), 2012, 43(7): 2668-2673.
DAI Qianwei, WANG Honghua, FENG Deshan, et al. Finite element method forward simulation for complex geoelectricity GPR model based on triangle mesh dissection[J]. Journal of Central South University (Science and Technology), 2012, 43(7): 2668-2673.
[17] 柳建新, 蒋鹏飞, 童孝忠, 等. 不完全LU分解预处理的BICGSTAB算法在大地电磁二维正演模拟中的应用[J].中南大学学报(自然科学版), 2009, 40(2): 484-491.
LIU Jianxin, JIANG Pengfei, TONG Xiaozhong, et al. Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling[J]. Journal of Central South University (Science and Technology), 2009, 40(2): 484-491.
[18] 薛东川, 王尚旭, 焦淑静. 起伏地表复杂介质波动方程有限元数值模拟方法[J]. 地球物理学进展, 2007, 22(2): 522-529.
XUE Dongchuan, WANG Shangxu, JIAO Shujing. Wave equation finite-element modeling including rugged topography and complicated medium[J]. Progress in Geophysics, 2007, 22(2): 522-529.
[19] 冯德山, 陈承申, 王洪华. 基于混合边界条件的有限单元法GPR正演模拟[J]. 地球物理学报, 2012, 55(11): 3774-3785.
FENG Deshan, CHEN Chengshen, WANG Honghua. Finite element method GPR forward simulation based on mixed boundary condition[J]. Chinese Journal of Geophysics, 2012, 55(11): 3774-3785.
[20] Clear.bout J F. Imaging the earth’s interior[M]. London: Blackwell Scientific Publications, 1985: 267-271.
(编辑 陈灿华)
收稿日期：2013-03-21；修回日期：2013-05-27
基金项目：国家自然科学基金资助项目(41374118)；中南大学研究生自主探索创新基金资助项目(2013bjjxj011，2013zzts053)
通信作者：戴前伟(1968-)，男，湖南涟源人，博士，教授，从事电磁法方法及理论、工程地球物理勘探等研究；电话：13574816998；E-mail: qwdai@csu.edu.cn
