基于三角形剖分的复杂GPR模型有限元法正演模拟
戴前伟，王洪华，冯德山，陈德鹏
(中南大学 地球科学与信息物理学院，湖南 长沙，410083)
摘 要：
格剖分的时域有限差分法(FDTD)和有限单元法(FEM)，对于物性参数分布复杂或几何特征不规则的模型适应性差的问题，从雷达波所满足的Maxwell方程出发，推导探地雷达(GPR)有限元波动方程，通过采用三角形网格剖分和线性插值基函数，在满足时间步长与空间网格差分稳定性前提下，应用Galerkin有限单元法求解GPR波波动方程；同时为消除FEM进行GPR正演模拟时来自截断边界处的超强反射，采用透射边界条件把GPR波在截断边界处的反射波透射出去，进而压制来自截断边界处的反射波。然后，编制GPR有限元正演模拟的Matlab程序。应用该程序分别对起伏分界面、“V”字形2个复杂地电模型进行FEM正演模拟，得到基于三角形网格剖分的FEM正演模拟GPR剖面图，并把该正演模拟剖面图与常规的基于矩形剖分的FEM正演模拟剖面图进行对比，结果表明：基于三角形剖分的FEM对于复杂GPR模型的物性参数分界面拟合更好，其模拟所得的正演剖面与实际模型更相符，具有更高的模拟精度，更有利于指导雷达剖面的数据解译。
关键词：
中图分类号：P319 文献标志码：A 文章编号：1672-7207(2012)07-2668-06
Finite element method forward simulation for complex geoelectricity GPR model based on triangle mesh dissection
DAI Qian-wei, WANG Hong-hua, FENG De-shan, CHEN De-peng
(School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)
Abstract: For the problem that finite difference time domain method (FDTD) and finite element method (FEM) based on the rectangular mesh dissection have poor adaptability to the patterns that physical parameter distributions are complex or geometric features are irregular, ground penetrating radar (GPR) finite element wave equation was deduced; from the Maxwell equation which met with the radar wave, under the premise of meeting the time step and space grid difference stability, triangular mesh dissection and linear interpolation function and Galerkin finite element method were adopted to solve GPR wave undulation equation. Meanwhile, in order to eliminate the super strong reflection on the interceptive boundary when the FEM was conducting GPR forward simulations, the transmission absorbed boundary condition was adapted to transmit reflected waves that in the interceptive boundary of GPR wave, and then suppressed the reflected wave from the interceptive boundary. Afterwards, the Matlab program of GPR finite element forward simulations was established, and this program conducted FEM forward simulations for two complex geoelectric model of curl interface under-ground and “V” shape respectively. The GPR profiles of FEM forward simulations was obtained based on triangular mesh dissection. The profiles based on this forward simulations and the FEM forward simulations profiles based on conventional rectangular dissection were compared. The results show that the FEM based on triangular dissection fits better to the physical parameter boundary surface of the complex GPR model , the forward simulations profiles are more consistent with actual models, and it has higher simulations accuracy and is more advantageous to the guidance of the radar profile data interpretation.
Key words: ground penetrating radar (GPR); finite element; triangular dissection; transmission absorbed boundary condition
探地雷达(GPR)是采用高频电磁波在地下电性界面(目标体)上的反射与透射信息来探测地质目标体的一种地球物理方法[1]。具有高效、快速、无损、抗干扰能力强的优点，广泛应用于工程与环境地球物理探测的各个领域[2-3]，成为浅部勘探的重要方法技术。GPR正演模拟对数据解释具有重要的指导作用。通过对GPR模型的正演模拟，可以加深对GPR探测剖面的认识，提高解释精度[4]。时域有限差分法(FDTD)以其算法简单、灵活，程序易实现，而得到了广泛的应用[5-8]，但是FDTD法对于物性参数分布复杂或场域几何特征不规则的模型，适应性差。为此，引入有限单元法(FEM)，该算法能适用于求解各种复杂问题、求解过程规范，计算程序通用性强等优点具有广泛的应用性。它在GPR正演模拟中应用起步较晚，主要是借鉴声波、弹性波的有限元求解理论。如：沈飚[9]将声波方程代替电磁波的传播方程进行了二维介质的正演模拟；底青云等[10]进行了带衰减项GPR剖面的FEM仿真模拟及偏移处理；谢辉等[11]利用二十节点等参单元对Pulse探地雷达进行了有限单元正演模拟计算，建立了电磁波在Pulse探地雷达天线中的全反射模型；Arias等[12]将GPR有限单元模拟技术应用到考古成像中。以往利用FEM进行GPR正演都是基于矩形网格剖分，对复杂GPR模型拟合效果较差，难以实现复杂地电模型的正演模拟。为此，本文作者采用三角形剖分方式对复杂GPR地电模型进行FEM正演模拟，并把该结果与基于矩形剖分的FEM进行了对比。
1 雷达波波动方程
雷达波满足Maxwell方程组：
(1a)
(1b)
(1c)
(1d)
其中：E为电场强度，V/m；H为磁场强度，A/m；B为磁感应强度，T；D为电位移，C/m2；J为电流密度，A/m2；ρ为电荷密度，C/m3。
本构关系方程：
，
，
(2)
式中：ε为介电常数，F/m；μ为磁导率，H/m；σ为电导率，S/m。
把方程式(2)代入方程式(1a)，并求旋度，得：
(3)
由矢量恒等发展可知：
(4)
由于ε，μ和σ为坐标的函数，其空间导数可以忽略，则式(4)中右边第2项
把式(4)带入式(3)，并整理得：
(5)
将雷达波激励源电场源Se或磁场源Sh代入式(5)中，有：
(6)
同理对式(1b)两边求旋度，可得下式：
(7)
式(6)和(7)表明，磁场H和电场E及其分量满足相同的偏微分方程。
2 GPR波动方程有限元求解
求解GPR波动方程其实质是应用FEM结合雷达波所要满足的初、边值条件求解偏微分方程。根据文献[10]和[13]，GPR所满足的偏微分方程(6)和(7)写成FEM矩阵形式为：
(8)
(9)
式中：M为质量矩阵；
为阻尼矩阵；F为边界阻尼矩阵；K为刚度矩阵；
和
为时间的一次导数；
和
为时间的二次导数。以电场为例，采用只对空间域进行离散，并令
(10)
则式(8)对应的方程为：
(11)
式中：Nj(x, y)为关于节点j的插值函数；U(t)为节点j上的时间函数；nd为节点数；
，
，
表示沿外法向求导；
代表求解的区域；
为的外边界。
2.1 三角单元线性插值
在三角形单元e上，记三角形3个顶点按逆时针方向分别为i，j和m，每个顶点的坐标分别为(xi, yi)，(xj, yj)和(xm, ym)，顶点的场值为Ui，Uj和Um，则三角形单元中的场函数为：
(12)
式中：
Ni(x, y)，Nj(x, y)和Nm(x, y)为三角单元的形函数。ai, …, am为只与单元顶点坐标有关的常数；
为三角形单元的面积；令：
则三角形单元上的质量矩阵Me，阻尼矩阵为
，刚度矩阵Ke分别为：
，
(13)
(14)
2.2 透射边界条件及源函数
方程式(11)包含场函数外法向导数沿边界的积分，其内边界积分相互抵消，如果仅采用简单的自由边界条件，即令外边界上场函数法向导数为0，会在截断边界处产生全反射，严重影响模拟的精度。为消除反射波，在截断边界上用入射波方程代替全波动方程就是需要的吸收边界条件[14]。
在均匀各向同性介质中，根据Clear bout的推导的旁轴近似下行波波动方程[15]为：
(15)
左行波方程为：
(16)
右行波方程为：
(17)
由法向导数的定义可知：
(18)
式中：nx和ny为边界处法线的方向余弦。
如图1所示，上边界
为自由边界条件，左边界
、下边界
和右边界
采用透射边界条件，将(15)~(18)代入方程(11)的边界积分项得：
(19)
其中：v为电磁波在媒质中的传播速度，
，
，
N/A2，
，
C2/(N·m2)。
以左边界为例，左边界单元上的阻尼矩阵为：
(20)
采用如下子波作为GPR正演模拟的脉冲激励源
(21)
(22)
其中：(xs, ys)为发射天线位置的坐标；w0为发射天线中心频率；a为衰减系数。则方程(8)的源项为：
(23)
图1 区域剖分及透射边界条件示意图
Fig.1 Schematic diagram for finite element mesh
2.3 中心差分法显式时间外推
经空间域离散后，组合求解域内各单元的系数矩阵，得到线性方程组：
(24)
其中：
，
，
，
，
。
采用中心差商代替微分，把时间域[0, T]用相等步长进行剖分，得?t=T/n，则时刻t的微分方程记为：
(25)
在中心差分中
(26)
(27)
将式(26)和(27)代入式(24)，得中心差分法显式时间外推公式为：
(28)
式中：
(29)
(30)
(31)
并结合零时刻及零时刻以前，区域内没有初始波场；零时刻后，区域内的波场才有激励源产生的初始条件：
(32)
根据式(28)并结合初始条件式(32)，则可根据t时刻和t-?t时刻点上的波场值Ut，Ut-?t求得区域内t+?t时刻上的波场值Ut+?t。在实际正演模拟计算中，由于系数矩阵A往往是大型的病态稀疏矩阵，其条件数很大，如果采取对其直接求逆的方法求解，计算量巨大，其求解结果不可信，只有当系数矩阵A是对角阵时，其逆矩阵是对角元素求倒数，为此采取集中质量矩阵方法[19]来处理M矩阵与矩阵，即将每一行(或列)的元素都加到对角线元素上去，计算效率显著提高。此外，中心差分法显式时间外推公式是条件稳定的，时间步长与空间步长必须满足一定的条件才能保证计算的稳定，即：
?t ≤L/C (33)
式中：L为网格剖分后最小尺寸单元的最小边长；C为单元上的波数；?t为雷达波通过该单元的时间，因此，网格剖分要尽量均匀，避免出现个别尺寸很小的单元而使?t很小，计算量随之不合理的增加。
3 数值模拟实例
设计了起伏界面、“V”字形界面2种GPR地电模型，应用基于三角形网格剖分的FEM算法结合透射边界条件对设定的模型进行了正演模拟，并把该计算结果与基于矩形网格剖分的FEM算法进行了模拟效果对比。
3.1 起伏界面模型
图2所示为起伏界面模型示意图。模型为一个10 m×5 m的矩形区域，分为上下两层介质，中间是起伏分界面，上层介质的相对介电常数ε1为5.0，电导率σ1为0.001 s/m，下层介质的相对介电常数ε2为50.0，电导率σ2为0.01 s/m。空间网格步长为0.1 m，网格总数为100×50。GPR波脉冲激励源的中心频率为100 MHz，时间步长为0.25 ns，时窗长度为100 ns。分别应用基于三角形和矩形网格剖分的FEM对该模型进行正演模拟，其模拟所得的GPR正演模拟剖面如图3所示。由图3可见：在剖面中40 ns左右时刻，存在一能量较强反射界面，通过时深转换计算可知，该反射界面与图2模型中的上、下两层分界面位置相吻合，并且图中变化剧烈的角点、界面起伏波动大、欠光滑的断棱处，产生了角点绕射波形。但通过对比图3(a)与图3(b)剖面可知，由于矩形剖分拟合起伏界面主要是以阶梯状模型拟合，三角形剖分能以三角单元的某条边来更近似的拟合起伏的界面，因此三角形剖分拟合起伏分界面的能力强于矩形剖分，导致图3(a)中产生的类似于角点绕射的波形能量明显比图3(b)的要弱。
3.2 “V”字形模型
图4所示为“V”字形模型示意图，模型为一个10 m×5 m的矩形区域，分为上下两层介质，中间为一“V”字形界面，“V”处在下层的最中间位置，宽度为3 m，深度为1.5 m，上层介质的相对介电常数ε1为5.0，电导率σ1为0.001 s/m，下层位于1.5 m深的位置，相对介电常数设置为50，电导率设置为0.01 s/m。空间网格步长为0.1 m，网格总数为100×50。GPR波脉冲激励源的中心频率为100 MHz，时间步长为0.25 ns，时窗长度为100 ns。分别应用基于三角形剖分和矩形剖分的FEM对“V”字形模型进行正演模拟，其模拟所得的GPR正演模拟剖面如图5所示。由图5可见：在36 ns附近有一条能量较强的反射界面，通过计算可以得出它与图5模型中的上、下两层分界面位置相吻合；分界面中间，可以看到一个“V”字形，但比模型中的“V”字形要浅得多。对比图5(a)和(b)可知：图5(b)在“V”字形的三个角点绕射的现象比图5(a)的要强，特别是下角点处图5(b)产生了较强的回转波或绕射波，而图5(a)中则未见明显的绕射波(或能量较弱)。由此可见，基于三角形剖分的有限元法对于“V”字形模型拟合更好，模拟所得的正演剖面具有更高的分辨率、更利于指导雷达剖面的数据解译。
图2 起伏界面模型示意图
Fig.2 Sketch map of fluctuation interface model
图3 起伏界面GPR模型不同剖分方式的FEM正演模拟剖面图
Fig.3 GPR forward compose section map of fluctuation interface model based on different element meshing
图4 “V”字形模型示意图
Fig.4 Sketch map of “V” shape model
图5 “V”字形GPR模型不同剖分方式的FEM正演模拟剖面图
Fig.5 GPR forward compose section map of “V” shape model based on different element meshing
4 结论
(1) 推导GPR有限元波动方程，给出有限单元法求解该泛函变分问题的详细解法，并结合差分离散时，时间步长与空间步长须满足的数值稳定性条件，探讨集中质量矩阵算法在GPR有限元波动方程求解中的应用，解决大型病态稀疏矩阵求逆的问题，提高了有限元线性方程组求解的计算效率。
(2) 利用透射边界条件对GPR模拟区域的截断边界处强反射波进行处理，模型计算结果表明：透射边界条件能很好的压制截断边界处的反射波能量，有效地减弱了边界处反射波的影响。
(3) 2个复杂GPR模型算例结果表明：相比FEM的矩形网格剖分方式，本文采用的FEM三角形单元网格剖分方法，能更好地拟合复杂GPR模型的物性参数分界面，其数值模拟结果能更有效地消除类似角点处的绕射波，对界面反映更符合模型实际，具有更高的模拟精度。
参考文献：
[1] 曾昭发, 刘四新. 探地雷达原理与应用[M]. 北京: 电子工业出版社, 2010: 168-169.
ZENG Zhao-fa, LIU Si-xin. Ground penetrating radar theory and applications[M]. Beijing: Press of Electronics Industry, 2010: 168-169.
[2] 张伟, 李姝昱, 张诗悦, 等. 探地雷达在水利工程隐患探测中的应用[J]. 水利与建筑工程学报, 2011, 9(1): 34-38.
ZHANG Wei, LI Shu-yu, ZHANG Shi-yue, et al. Application of GPR in detecting potential hazards of water conservancy projects[J]. Journal of Water Resources and Architectural Engineering, 2011, 9(1): 34-38.
[3] 杨天春, 周勇, 李好. 超前探测中探地雷达应用与结果的处理分析[J]. 工程地质学报, 2010, 18(6): 971-1004.
YANG Tian-chun, ZHOU Yong, LI Hao. Use and analysis of ground penetration radar for casting geological conditions in highway tunneling[J]. Journal of Engineering Geology, 2010, 18(6): 971-1004.
[4] 郭成超, 王复明. 探地雷达电磁波正演模拟研究[J]. 公路交通科技, 2008, 25(8): 37-41.
GUO Cheng-chao, WANG Fu-ming. Study on forward modeling of GPR electromagnetic wave propagation[J]. Journal of Highway and Transportation Research and Development, 2008, 25(8): 37-41.
[5] Bergmann T, Robertsson J O A, Holliger K. Numerical properties of staggered finite difference solutions of Maxwell’s equations for ground-penetrating radar modeling[J]. Geophysical Research Letters, 1996, 23(1): 45-48.
[6] Carcione J M. Radiation patterns for 2-D GPR forward modeling[J]. Geophysics, 1998, 63(8): 424-430.
[7] Kowalsky M B, Dietrich, P Teutsch G, et al. Forward modeling of ground-penetrating radar data using digitized outcrop images and multiple scenarios of water saturation[J]. Water Resources Research, 2001, 37(6): 1615-1625.
[8] Good Man D. Ground-penetrating radar simulation in engineering and archaeology[J]. Geophysics, 1994, 59(2): 224-232.
[9] 沈飚. 探地雷达波波动方程研究及其正演模拟[J]. 物探化探计算技术, 1994, 16(1): 29-33.
SHEN Biao. Study of radar wave equation and its forward modeling[J]. Geophysical and Geochemical Computation Technology, 1994, 16(1): 29-33.
[10] 底青云, 王妙月. 雷达波有限元仿真模拟[J]. 地球物理学报, 1999, 42(6): 818-825.
DI Qing-yun, WANG Miao-yue. 2D finite element modeling for radar wave[J]. Chinese Journal of Geophysics, 1999, 42(6): 818-825.
[11] 谢辉, 钟燕辉, 蔡迎春. 电磁场有限元法在GPR正演模拟中的应用[J]. 河南科学, 2003, 21(3): 295-298.
XUE Hui, ZHONG Yan-hui, CAI Ying-chun. Application of finite element methods for electromagnetic fields in forward model of GPR[J]. Henan Science, 2003, 21(3): 295-298.
[12] Arias P, Armesto J, Capua D, et al. Digital photogrammetry GPR and computational analysis of structural damages in a mediaeval bridge[J]. Engineering Failure Analysis, 2007, 14(8): 1444-1457.
[13] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994: 266-275.
XU Shi-zhe. Finite element method for geophysics[M]. Beijing: Science Press, 1994: 266-275.
[14] 薛东川, 王尚旭, 焦淑静. 起伏地表复杂介质波动方程有限元数值模拟方法[J]. 地球物理学展, 2007, 22(2): 522-529.
XUE Dong-chuan, WANG Shang-xu, JIAO Shu-jing. Wave equation finite-element modeling including rugged topography and complicated medium[J]. Progress in Geophysics, 2007, 22(2): 522-529.
[15] Clear bout J F. Imaging the earth’s interior[M]. Blackwell Scientific Publications, 1985: 267-271.
(编辑 陈爱华)
收稿日期：2011-09-22；修回日期：2011-11-15
基金项目：国家自然科学基金资助项目(40804027，41074085)；湖南省自然科学基金重点资助项目(09JJ3084)；中南大学硕士生学位论文创新资助项目(2011ssxt055)
通信作者：冯德山(1978-)，男，湖南祁阳人，博士，副教授，从事探地雷达正反演及小波分析研究；电话：0731-88836145；E-mail: fengdeshan@126.com
摘要：针对基于矩形网格剖分的时域有限差分法(FDTD)和有限单元法(FEM)，对于物性参数分布复杂或几何特征不规则的模型适应性差的问题，从雷达波所满足的Maxwell方程出发，推导探地雷达(GPR)有限元波动方程，通过采用三角形网格剖分和线性插值基函数，在满足时间步长与空间网格差分稳定性前提下，应用Galerkin有限单元法求解GPR波波动方程；同时为消除FEM进行GPR正演模拟时来自截断边界处的超强反射，采用透射边界条件把GPR波在截断边界处的反射波透射出去，进而压制来自截断边界处的反射波。然后，编制GPR有限元正演模拟的Matlab程序。应用该程序分别对起伏分界面、“V”字形2个复杂地电模型进行FEM正演模拟，得到基于三角形网格剖分的FEM正演模拟GPR剖面图，并把该正演模拟剖面图与常规的基于矩形剖分的FEM正演模拟剖面图进行对比，结果表明：基于三角形剖分的FEM对于复杂GPR模型的物性参数分界面拟合更好，其模拟所得的正演剖面与实际模型更相符，具有更高的模拟精度，更有利于指导雷达剖面的数据解译。
