DOI: 10.11817/j.issn.1672-7207.2016.05.027

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

石明^{1, 2}，冯德山¹，王洪华³，戴前伟¹

 (1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；

2. 贵州省有色金属和核工业地质勘查局物化探总队，贵州 都匀，558004；

3. 桂林理工大学 地球科学学院，广西 桂林，541004)

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

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

中图分类号：P631          文献标志码：A         文章编号：1672-7207(2016)05-1660-08

GPR numerical simulation for anisotropic medium by

using finite element method based on unstructured meshes

SHI Ming^{1, 2}, FENG Deshan¹, WANG Honghua³, DAI Qianwei¹

 (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)

Abstract: 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 monoclinic anisotropic 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.

Key words: anisotropic media; finite element method based on unstructured meshes; GPR (ground penetrating radar); numerical simulation

探地雷达(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模式的电磁波理论

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

            (1)

式中：E和H分别为电场强度和磁场强度；和分别为介电常数和电导率张量；为磁导率。根据各向异性介质理论，单斜各向异性介质的介电常数和电导率可分别表示为

，    (2)

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

1) TM模式：

2) TE模式：

          (3)

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

 (4)

将式(4)两边同除以，代入磁激励源，并转化到时间域可得

           (5)

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

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

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

         (6)

式中：i，j和k为三角单元上3个节点编号；和分别表示相应节点上的磁场和形函数。根据伽辽金有限单元法，式(5)的弱解形式可表示为

              (7)

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

图1  网格剖分及节点编号示意图

Fig. 1  Sketch map of mesh division and point number

               (8)

式(8)左边第1项：

          (9)

式(8)左边第2项：

     (10)

式(8)左边第3项：

       (11)

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

  (12)

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

                (13)

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

              (14)

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

               (15)

因为的物理意义是表示天线到边界上任一点的矢径与边界法线方向的夹角，所以，有。将其代入式(15)，可得

            (16)

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

      (17)

则式(12)可写为

 (18)

式中：

        (19)

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

                 (20)

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

    (21)

式(8)左边第6项为

   (22)

将式(9)，(10)，(11)，(18)，(20)和(22)代入式(8)可得

  (23)

对时间域有限元方程式(23)，采用中心差分法即可进行求解。将时间域[0 T]分为n等份，则采样间隔为，在t时刻，

 (24)

用中心差分代替微分，并令，式(24)可变为

        (25)

其中：的时间迭代格式为

              (26)

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

3  数值算例

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

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

图1所示为均匀各向异性介质非结构化Delaunay网格示意图。模型为1个长×宽为2.0 m×2.0 m的矩形区域，非结构化网格Delaunay网格用于离散计算区域。从图1可以看到：在激励源点和接收点附近采用密网格剖分以适应电磁波场的快速变化，而其他区域采用稀疏网格剖分。这样可用较少的网格去精确逼近电磁波的传播，克服激励源点和接收点产生的奇异性。介质的相对介电常数为，，电导率为0.001 S/m。GPR波脉冲激励源的中心频率为500 MHz，时间步长为0.02 ns，时窗为12 ns。由于各向异性介质的解析解求解非常困难，因此，将Y方向上的相对介电常数设置为5，这样可以方便求解解析解以及FDTD和FETD的正演结果对比。

图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

为了更详细地分析GPR波在各向异性介质中的传播规律，改变模型1中的相对介电常数，分别设置成均匀各向同性介质()模型、均匀各向异性介质(，)模型、均匀各向异性介质(，)模型、右边均匀各向异性介质(，)和左边各向同性介质()组合模型，得到如图3所示的4种不同相对介电常数分布条件下各向异性介质模型在9 ns时的波场快照。由图3可知：GPR波在各向同性介质中呈圆状向外传播，各个方向传播速度相同；而在各向异性介质中呈椭圆状向外扩散传播，各向方向传播速度不同；对比图3(a)，(b)，(c) 可知GPR波在各向异性介质中的传播能量比在各向同性介质中能量衰减慢。GPR波在各向异性介质中传播能量衰减慢，由图3(d)可知当GPR波传播到9 ns时，左、右两边的波前面已发生断开。这是由于GPR波在左边各向同性介质中传播速度快，而在右边各向异性介质中传播速度慢，导致相同时刻GPR波传播的距离不同。

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

图4所示为1个各向异性介质圆形异常体非结构网格剖分模型示意图。模型为长×宽为10 m×5 m的矩形区域，背景介质为各向异性介质Ⅰ，其相对介电常数为，，电导率为0.001 S/m。模型中间(5.0 m，2.2 m)处埋有1个半径为0.4 m的圆形异常体，其相对介电常数，电导率0.01 S/m。背景介质分别为均匀各向同性介质、各向异性介质，其电导率都为0.001 S/m。电磁波脉冲激励源的中心频率为100 MHz，时间步长为0.1 ns，时窗长度为80 ns。非结构化Delaunay网格用于剖分计算区域，可以看到在圆形异常体和地表激励源放置处，采用较密的非结构网格进行剖分以适应电磁场的剧烈变化，而在其他电磁波场变化平缓的区域，采用较疏的网格进行剖分，这样可用较少的网格精确剖分地电模型。图5所示为各向异性介质的圆形异常体正演剖面图。由图5可见：在40 ns和50 ns分别出现了2条非常明显的双曲线反射波，经过时深转换，刚好对应圆形异常体顶部和底部产生的反射，这也验证了各向异性介质GPR有限元算法的有效性。

为了更好地认识电磁波在各向异性介质中的传播特征，将图4 中模型中的背景介质参数分别设置为均匀各向同性介质()和各向异性介质Ⅱ(，), 将其计算结果与图5所示结果进行对比。图6所示分别为均匀各向同性介质和各向异性介质Ⅱ中圆形异常体正演剖面。由图6可见：当背景介质为各向异性介质Ⅰ时，圆形异常体反射波到达时间约为  40 ns，与均匀介质中圆形异常体的反射波到达时间相同；当背景介质为各向异性介质II时，其反射波到达时间增大，约为50 ns，这主要是电磁波在各向异性介质中不同方向上的传播速度不同所致。此外，图5中双曲线反射波的能量比图6中的能量高，而图6 (b)中反射波双曲线的能量要比图6(a)中的双曲线反射波能量低，这主要是电磁波在各向异性介质中传播速度不同，从而导致波前扩散损失不同。另外，图6(a)中双曲线反射波的曲率比图5和图6 (b)中的大，因此，各向异性介质主轴方向上相对介电常数会严重影响圆形异常体反射波的能量和曲率。

图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 ground- penetrating-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.

(编辑  陈灿华)

收稿日期：2015-10-10；修回日期：2015-12-22

基金项目(Foundation item)：国家自然科学基金资助项目(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[1]