各向异性介质GPR非结构化网格有限元正演

石明，冯德山王洪华，戴前伟

（1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；2. 贵州省有色金属和核工业地质勘查局物化探总队，贵州 都匀，558004；3. 桂林理工大学 地球科学学院，广西 桂林，541004）

各向异性介质GPR非结构化网格有限元正演

石明1，2，冯德山1，王洪华3，戴前伟1

（1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；2. 贵州省有色金属和核工业地质勘查局物化探总队，贵州 都匀，558004；3. 桂林理工大学 地球科学学院，广西 桂林，541004）

为了更准确地认识电磁波在各向异性介质中的传播特征，从麦克斯韦方程组出发，推导电磁场在单斜各向异性介质中的波动方程。采用 Delaunay 非结构化网格有限元法进行空间域离散，采用中心差分公式进行时间域离散，导出单斜各向异性介质中GPR时空域有限元方程。在此基础上，建立各向异性介质GPR非结构化网格有限元正演算法。利用以上算法编制的程序对均匀各向异性介质模型进行正演计算，并与解析解对比，验证该算法的正确性和有效性；对均匀各向异性介质模型和圆形异常体模型进行计算，并与背景介质为各向同性介质的计算结果进行对比。研究结果表明：与FDTD计算结果对比，非结构化网格有限元法的拟合误差更小，精度更高；电磁波传播受介质各向异性的影响，波前面呈椭圆状向外扩散传播，其能量衰减较慢；与各向同性介质中圆形异常体反射波相比，各向异性介质的圆形异常体反射波的曲率、能量和双程走时，受不同方向上电磁波速度影响会发生显著变化。

各向异性介质；非结构化网格有限单元法；探地雷达；正演模拟

探地雷达（ground penetrating radar，GPR）作为一种重要的浅部地球物理勘探手段被广泛应用于工程勘察［1-2］、无损检测［3］等众多领域，具有广阔的应用前景［4］。目前，GPR实测数据的解释主要依赖各向同性介质的电磁波传播理论，然而，随着GPR应用领域的扩大，各向异性介质如冰川［5］、结晶体［6］、含水岩石、构造裂隙等已成为探测对象，与电磁波在各向同性介质中的传播规律相比，电磁波在各向异性介质中传播会发生畸变和衰减，因此，研究电磁波在各向异性介质中的传播和异常体回波特征对提高 GPR实测数据的解释精度有重要意义。GPR正演模拟是研究电磁波在各向异性介质中传播规律的有效手段。在众多的GPR正演算法中，时域有限差分法（finite difference time domain， FDTD）因其发展成熟和理论简单而被广泛应用。刘四新等［7-10］对FDTD算法在GPR正演中的应用进行了深入研究。但FDTD要求地电模型被剖分成如正方形、立方体等规则的单元，这严重制约着其在复杂介质模型的应用［11］。与FDTD相比，有限元法（finite element method，FEM）可以利用非结构化网格更精确的离散复杂介质模型［12］，在电磁波场梯度大的区域采用较密的网格，而在电磁波场变化平稳的区域采用较疏的网格。因此，采用较少的网格就可以构建复杂地质结构，具有更高的计算效率和模拟精度［13］。在目前的GPR正演模拟研究中，地下介质大都被假定为各向同性，人们对各向异性介质的GPR正演模拟研究较少，且大都采用FDTD算法［14-16］。为了更好地模拟电磁波在各向异性介质中的传播规律，本文作者应用非结构化网格有限元法对各向异性介质中的电磁波传播规律进行研究。通过对均匀各向同性介质模型、均匀各向异性模型、圆形异常体各向异性介质模型进行计算，并与各向同性介质的正演模拟结果进行对比，总结电磁波和异常体回波在各向异性介质中的传播特征。

1 TM模式和TE模式的电磁波理论

根据电磁波理论，频率域各向异性介质中的麦克斯韦方程可表示为

式中： εjj和σjj（j=xx，yy和zz）分别为j方向上的介电常数和电导率 。考虑二维情况，各向异性介质的主轴方向为z方向，将式（2）代入式（1），按照TE模式和TM模式分解为2组独立的方程组。

1） TM模式：

2） TE模式：

从式（3）可以看出：各向异性介质中TM模式满足的方程与介电常数为 εzz和电导率为σzz的各向同性满足的方程完全相同，而各向异性介质的 TE模式方程组与各向同性介质有较大的差别。将 TE模式方程中的前2式的Ex和Ey表达式代入第3式，经整理可得

将式（4）两边同除以εμεyyxx，代入磁激励源 Sh，并转化到时间域可得

其中：Hz和Hz分别为Hz对时间的一次和二次导数。式（5）即为 TE模式下各向异性介质中的磁场波动方程。

2 各向异性介质的 GPR有限元正演模拟算法

求解 GPR波动方程的实质是应用伽略金有限单元法，结合初始、边界条件，推导相应的有限元方程［17］。这里采用非结构化 Delaunay 网格［18］对计算区域进行剖分，如图1所示。非结构化网格可在激励源点附近与接收点附近采用密网格，而在大部分均匀介质区域采用疏网格剖分，因此，较少的网格数便可精确剖分复杂介质模型。根据有限单元法理论，三角单元内的磁场可用该单元上3个节点上的磁场表示为

式中：R为方程（7）的残差；Ω为计算区域面积；N为形函数向量。将式（5）和（6）代入式（7），可得

图1 网格剖分及节点编号示意图Fig.1 Sketch map of mesh division and point number

式（8）左边第1项：

式（8）左边第2项：

式（8）左边第3项：

根据高斯定理，式（8）左边第4项可写为

由电磁场理论可知，GPR波在模型边界上向外传播应该满足

其中：k为电磁波在介质中传播的波数；r为天线位置到边界上任一点的距离。对式（14）两边对 r求导，可得

将式（13）两边分别乘以k并与式（14）相加，可得

其中：n为模型边界的外法线方向。将式（16）代入式（12）右边第1式可得

则式（12）可写为

式中：

根据式（8）第 4项的推导方式，同理可推导式（8）中左边第5项为

其中：刚度矩阵K3的表达式为

式（8）左边第6项为

将式（9），（10），（11），（18），（20）和（22）代入式（8）可得

对时间域有限元方程式（23），采用中心差分法即可进行求解。将时间域［0 T］分为n等份，则采样间隔为nTt/=Δ，在t时刻，

其中：tI的时间迭代格式为

式（25）和式（26）即为GPR有限元正演模拟时间递推公式。由于零时刻的Hz，0，Hz，0和Hz，0以及tΔ-时刻的Hz，-Δt，Hz，-Δt和Hz，-Δt均为0，故根据式（25）和（26）可计算tΔ时刻的Hz，Δt，然后依次递推可得到各个时刻的电场强度Hz。在求解式（25）时，选用不完全LU分解预处理的双共轭梯度法BICGSTAB（biconjugate gradient stabilized）线性方程组求解算法［19］。该算法是一种高效迭代法，特别适用于求解方程组条件数很大的病态线性方程组。

3 数值算例

3.1 非结构化网格有限元正演算法的正确性验证

依上述方法原理，编制了各向异性介质的GPR非结构化网格有限元正演源程序。利用该程序分别对均匀各向异性和均匀各向同性介质模型进行计算，并与均匀各向同性介质的时域有限差分法的结果与解析解进行对比，以验证程序的正确性。

图2所示为均匀各向异性介质解析解、FDTD 和FEM正演单道波形对比结果。从图2可见： FEM 和FDTD的正演结果与解析解都非常吻合，但在波峰和波谷处，还存在一定差异，FEM的结果更靠近解析解，而FETD的结果在波峰和波谷处拟合较差。经过计算，在6.0～10.0 ns之间，FDTD的相对误差为1.2%，而FEM 的相对误差为 0.63%，由此可见基于非结构Delaunay网格的有限元法与FDTD算法相比具有更高的精度。这主要是由于非结构化Delaunay 相比FDTD中正方形网格剖分计算区域具有更小的离散误差，局部加密技术能精确逼近电磁场的剧烈变化。

图2 均匀各向异性介质解析解、FDTD和FEM正演单道波形对比Fig.2 Waveforms comparison among FETD results， FDTD results and analytical results for homogenous anisotropic media

3.2 各向异性介质圆形异常体模型

图3 4种不同相对介电常数分布条件下各向异性介质模型在9 ns时刻的波场快照图Fig.3 Snapshots of anisotropic models at 9 ns with four different permittivity distributions

图4 各向异性介质Ⅰ圆形异常体非结构化网格剖分示意图Fig.4 Schematic diagram of cylinder model discretized by unstructured meshes of Delaunay

图5 各向异性介质Ⅰ圆形异常体有限元正演剖面图Fig.5 Simulated profile of cylinder model with isotropic background media

图6 背景介质分别为各向同性介质Ⅰ和各向异性介质Ⅱ的GPR正演剖面图Fig.6 GPR simulated profiles of circle model with anisotropic background mediaⅠand isotropic background mediaⅡ

4 结论

1） 从Maxwell方程组出发，在介质满足单斜各向异性条件下，推导了时间域GPR波在各向异性介质中满足的波动方程；并阐述了非结构化Delaunay网格有限单元法和中心差分法求解该波动方程的详细解法。

2） 在同等条件下，非结构化网格有限元正演精度较传统有限元精度有显著提高。

3） 电磁波在各向异性介质中呈椭圆状向外扩散传播，其能量衰减较慢，与各向同性介质中圆形异常体反射波相比，各向异性介质中的反射波的曲率、能量和双程走时受不同方向电磁波速度的影响会发生显著变化。

［1］ 曾昭发， 刘四新. 探地雷达原理与应用［M］. 北京： 电子工业出版社， 2011： 168-169. ZENG Zhaofa， LIU Sixin. Ground penetrating radar theory and applications［M］. Beijing： Electronics Industry Press， 2010：168-169.

［2］ BATTISTA B M， ADDISON A D， KNAPP C C. Empirical model decomposition operator for Dewowing GPR data［J］. Journal of Environmental and Engineering Geophysics， 2009，14（4）： 163-169.

［3］ GOODMAN D. Ground-penetrating radar simulation in engineering and archaeology［J］. Geophysics， 1994， 59（2）：224-232.

［4］ SLOB E， SATO M， OLHOEFT G. Surface and borehole groundpenetrating-radar developments［J］. Geophysics， 2010， 75（5）：75A103-75A120.

［5］ GIFFORD C M， FINYOM G， JEFFERSON M， et al. Automated polar ice thickness estimation from radar imagery［J］. IEEE Transactions on Geoscience and Remote Sensing， 2010， 19（9）：2456-2469.

［6］ 魏兵. 各向异性介质电磁散射及参数反演研究［D］. 西安： 西安电子科技大学理学院， 2003： 1-2. WEI Bing. Study of electromagnetic scattering and reconstruction for anisotropic media［D］. Xi’an： Xiandian University. College of Science， 2003： 1-2.

［7］ 刘四新， 曾昭发， 徐波. 三维频散介质中地质雷达信号的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.

［8］ 刘四新， 曾昭发. 频散介质中地质雷达波传播的数值模拟［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.

［9］ GIANNOPOULOS A. Modeling ground penetrating radar by GprMax［J］. Construction and Building Materials， 2005， 19（10）：755-762.

［10］ JAMES I， ROSEMARY K. Numerical modeling of ground-penetrating radar in 2-D using MATLAB［J］. Computers & Geosciences， 2006， 32（9）： 1247-1258.

［11］ 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.

［12］ 戴前伟， 王洪华， 冯德山， 等. 基于三角形剖分的复杂 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.

［13］ 王洪华， 戴前伟. 基于UPML吸收边界条件的GPR有限元数值模拟［J］. 中国有色金属学报， 2013， 23（7）： 2003-2011. WANG Honghua， DAI Qianwei. Finite element numerical simulation for GPR based on UPML boundary condition［J］. The Chinese Journal of Nonferrous Metals， 2013， 23（7）： 2003-2011.

［14］ SCHNEIDER J， HUDSON S. The finite-difference time-domain method applied to anisotropic material［J］. IEEE Transactions on Antennas and Propagation， 1993， 41（7）： 994-999.

［15］ CARCIONE J M， SCHOENBERG M A. 3D ground-penetrating radar simulation and plane-wave theory in anisotropic media［J］. Geophysics， 2000， 65（5）： 1527-1541.

［16］ 王帮兵， 田钢， 孙波， 等. 南极冰盖内部结构特性研究： 基于三维各向异性电磁波时域有限差分方法［J］. 地球物理学报，2009， 52（4）： 966-975. WANG Bangbing， TIAN Gang， SUN Bo， et al. The study of the COF feature in the Antarctic ice sheet based on 3-D anisotropy FDTD method［J］. Chinese Journal of Geophysics， 2009， 52（4）：966-975.

［17］ 徐世浙. 地球物理中的有限单元法［M］. 北京： 科学出版社，1994： 266-275. XU Shize. Finite element method for geophysics［M］. Beijing：Science Press， 1994： 266-275.

［18］ KEY K， WEISS C. Adaptive finite element modeling using unstructured grids： the 2D magnetotelluric example［J］. Geophysics，2006， 71（6）： 291-299.

［19］ 柳建新， 蒋鹏飞， 童孝忠， 等. 不完全 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.

（编辑 陈灿华）

GPR numerical simulation for anisotropic medium by using finite element method based on unstructured meshes

SHI Ming1，2， FENG Deshan1， WANG Honghua3， DAI Qianwei1

（1. School of Geosciences and Info-Physics， Central South University， Changsha 410083， China；2. Geophysical and Geochemical Prospecting Team，Non-ferrous Metals and Nuclear Industry Geological Exploration Bureau of Guizhou， Duyun 558004， China；3. College of Earth Sciences， Guilin University of Technology， Guilin 541004， China）

In order to more accurately understand the propagation law of GPR （ground penetrating radar） wave in anisotropic medium， the wave equations in monoclinic anisotropic medium were deduced based on the Maxwell equation，and the finite element method of unstructured meshes and the central difference method were applied to discretize in space domain and time domain respectively. Then， the GPR finite element numerical simulation algorithm of monoclinicanisotropic medium was built and programmed， and its feasibility and effectiveness were tested by comparing the results of simulation and analytical solution of homogenous anisotropic model. The finite element method based on unstructured meshes was applied to simulate the circle anisotropic model and homogeneous anisotropic model and compared with calculated results of its isotropic medium model. The results show that the simulated results by finite element method based on unstructured meshes have lower fitting error and higher precision compared with the simulated results by FDTD，the electromagnetic wave propagates outward with elliptical circle and its energy attenuates slower duo to the different velocities of wave in different directions. And curvature， energy and two way travel time of reflected wave change significantly compared with those of isotropic medium due to different wave velocitis in different directions.

anisotropic media； finite element method based on unstructured meshes； GPR （ground penetrating radar）；numerical simulation

P631

A

1672-7207（2016）05-1660-08

10.11817/j.issn.1672-7207.2016.05.027

2015-10-10；

2015-12-22

国家自然科学基金资助项目（41574116，41004053）；中南大学创新驱动项目（2015CX008）；中南大学教师研究基金资助项目（2014JSJJ001）；教育部新世纪优秀人才支持计划项目（NCET-12-0551）；中南大学升华育英人才计划项目（2012）；湖湘青年创新创业平台培养对象项目（2013） （Projects（41574116， 41004053） supported by the National Natural Science Foundation of China； Project（2015CX008）supported by the Innovation Driven Program of Central South University； Project（2014JSJJ001） supported by the Teacher Research Fund of Central South University； Project（NCET-12-0551） supported by the New Century Excellent Talent Support Program of the Ministry of Education；Project（2012） supported by the Sublimation Yuying Talent Program of Central South University； Project（2013） supported by the Hunan Youth Innovation and Entrepreneurship Training Platform）

冯德山，教授，博士生导师，从事地球物理正反演与数据处理；E-mail： fengdeshan@126.com