基于CFS-PML 的高阶有限元在地质雷达数值模拟中的应用研究

韩晓冰 郭涛 李航 梁冰洋 周远国

（西安科技大学通信与信息工程学院,西安 710600）

引言

地质雷达(ground penetrating radar，GPR)数值模拟是用电磁计算来研究地下介质分布规律的一种方法，通过将数值计算获得的模拟结果与实际研究区域的检测结果进行对比，能够有效地分析GPR 技术的实际应用效果，可为实际应用中检测数据的地质解释奠定理论基础，该方法逐渐成为近年来地球物理勘探领域的研究热点[1-3].常用的数值模拟方法包括有限差分法[4]、有限元法(finite element method,FEM)[5-6]等.在GPR 建模过程中，为了减少计算量，通常会采用完美匹配层(perfectly matched layer，PML)[7-8]对模型进行截断，PML 的吸收效果直接影响整个模型的计算精度.

然而，PML 边界条件在数值模拟过程中，仍存在一些复杂问题，例如，对倏逝波都不能很好地吸收[9]，因此在文献[10]中提出了基于单轴完美匹配层(uniaxial perfectly matched layer，UPML).UPML 不需要对吸收介质中的波场进行分裂，因此，计算简洁、内存消耗降低.并且在衰减系数中引入了吸收因子，使得UPML 边界能够较好地衰减倏失波、低频波等干扰波.文献[11]利用高阶时域有限差分(finitedifference time-domain,FDTD)法结合UPML 边界条件实现GPR 正演数值模拟，优化迭代UMPL 吸收边界，完善了三维模拟在计算效率与吸收效果的双重要求.文献[12]采用了基于交替方向隐式FDTD(alternating direction implicit-FDTD,ADI-FDTD)法，引入UMPL 消除了截断边界处的强反射现象.但UMPL 媒质的角形区域仍需频时转换，迭代求解大型线性方程组，增加了内存消耗.

虽然UPML 对于倏失波、低频波等干扰波具有一定的吸收效果，但吸收效果依然不是特别完美.复频移PML(complex frequency shifted-PML，CFS-PML)技术主要利用坐标展宽来增强倏失波在场中的衰减，其可以吸收低频信号在边界表面的反射，同时对大角度掠射波有很好的吸收效果[13](主要通过频移参数 α降低掠射波，入射到PML 界面时产生的非均匀波).文献[14]利用双线性变换技术增加CFS 因子，以避免场分裂情况，并且实现了柱坐标系下的CFS-PML吸收边界，与传统UMPL 边界相比，CFS-PML 边界条件具有更好的吸收性能.文献[15]实现了CFSPML 在三维海洋可控源电磁建模中的应用，有效地抑制住了人工边界效应，与Dirichlet 边界相比，降低了内存消耗，提高了计算效率.文献[16]改进了CFSPML 边界条件，并成功将其应用于瞬变电磁法的三维有限差分虚拟波域(finite difference-fictitious wave domain,FD-FWD)建模中.采用均匀半空间模型，数值结果表明在FD-FWD 中，CFS-PML 边界条件具有比PML 更好的吸收效果.

GPR 数值模拟时，吸收边界产生的人为反射波会影响整个模型的计算精度.为进一步改善GPR 三维数值模拟边界条件的吸收效果，同时优化计算效率，本文引入高阶有限元法(higher order finite element method，HO-FEM).HO-FEM 是FEM 的一种改进形式，不同之处在于其采用高阶正交多项式作为基函数[17].一方面提高了算法的效率，另一方面提高了离散网格内的建模精度，降低建模时对网格的依赖，且兼顾了常规FEM 对复杂地质模型模拟的灵活性[18]，该方法可使计算精度较高且计算量大幅减少[19].相较于HO-FEM 利用PML 作为吸收边界[20]，CFSPML 能更有效地处理倏失波衰减的问题.HO-FEM算法虽然已经应用于多种电磁场数值模拟中[21-22]，但在GPR 三维数值模拟中的研究尚未涉及.因此，本文做了以下两部分工作：1)将复频移坐标拉伸系统引入到PML 域内，使得频移参数 α对倏失波进行有效衰减；2)将CFS-PML 算法与HO-FEM 法进行结合，以此来求解三维GPR 数值模拟问题，提高计算精度与效率.

1 理论方法

1.1 HO-FEM 方法

从微分形式麦克斯韦方程组出发，利用时变因子 eiwt，得到如下控制方程：

式中：E为 电场强度，V/m；ω为角频率；B为磁通量密度，Wb/m2；H为磁场强度，A/m ；D为电通量密度，C/m2；J为电流密度，A/m2；μ 为磁导率；ε为介电常数.对式(1)中第二个方程取旋度，得到关于磁场矢量赫姆霍兹方程

式中：εr为相对介电常数；为真空中的波数；μr为相对介磁常数.与FEM 相同，对整个计算域进行区域离散，得到很多子区域，称之为单元.采用伽辽金加权残差法处理式(2)，令e单元内的加权残差积分为零，即

式中：Φi为加权基函数；Ωe是单元计算域.对式(3)展开如下：

式中：Nb表示单元外边界棱边数.使用伽辽金加权残差法时，规定采用的加权基函数和近似解展开函数相同，这里采用高斯-洛巴托-勒让德(Gauss-Lobatto-Legendre,GLL)高阶插值多项式作为加权基函数，一维参考空间上的N(N≥2)阶GLL 基函数

式中：PN(ξj)为Legendre 多项式；ξj为GLL 插值节点；j=0,1,···,N.

将所有单元进行集成，离散系统矩阵如下：

由式(1)控制方程推导得到磁场矢量赫姆霍兹方程，结合式(3)伽辽金加权残差法，引入式(5)中的GLL 高阶插值多项式，得到式(6) HO-FEM 的离散系统矩阵形式，最终我们将HO-FEM 算法应用于整个计算域内.

1.2 HO-FEM-CFS-PML 计算方法

CFS-PML[23]的主要思想是在人工边界内水平和垂直方向引入一个复杂的坐标拉伸系统.例如，x方向拉伸因子e(x)为

在电磁建模中主 要考虑到 参数k(x)和dx(x)，且k(x)和dx(x)均取正数，在本文中k(x)取值为25，与 文献[16]一致.

α(x)参数的本质是吸收倏逝波并改善掠射角处的吸收性能[24]，表达式为

式中：LPML为 每侧PML 区域厚度；pd取值一般为2 或3[16]；f0为电磁波主频率.

dx(x)作为x方向的衰减阻尼因子，表达式为

式中：cp为 电磁波传播速度；R是反射系数.式(10)中的1 为计算域x方向长度缩放为1.从式(10)看出，上述公式中dx(x)须在PML 域内取值.对于三维模型而言，PML 域内拥有三种媒质分区，分别是平面区、棱边区和角顶区，每个分区都应有d(x)作为衰减阻尼因子进行电磁波的衰减.

在HO-FEM 中引入CFS-PML 后，赫姆霍兹方程可以表示为

式(12)只针对PML 域内赫姆霍兹方程式.为使PML内边界与计算域外边界形成阻抗匹配，引入CFSPML 方法后三维空间中的相对介电张量以及相对介磁张量为：

结合HO-FEM 方法，采用伽辽金加权残差法处理，将式(6)中的质量矩阵和刚度矩阵修正为：

通过以上公式，我们将CFS-PML 复频移坐标拉伸系统引入到HO-FEM 中的PML 域内，对坐标空间中的电性参数进行复频移拉伸，从而改善介质本构张量的矩阵特性，增强了整个PML 域对倏逝波的吸收性能.本文并未描述参考空间与物理空间的转换，具体可参考文献[20].

2 数值算例

在这一节中，我们设计了两个GPR 模型，第一个模型采用三维分层介质结构，第二个模型采用三维复杂异常体介质结构.将本文的CFS-PML 方法与传统PML 方法在数值精度上做对比，为方便起见，两种算例都采用线源作为发射源，极化方向都为y方向，并与商业软件WCT 进行比较，来验证该方法的准确性.同时，对比特定频率下反射误差、相对误差，来验证本文算法的有效性和精确性.所有算例的仿真都在Intel(R) Xeon(R) Gold 6248R 的CPU 工作平台进行.

2.1 分层介质GPR 模型

算例1 三维分层介质GPR 模型如图1 所示，计算区域为4 m×4 m×2 m.以o为坐标原点，z=1 m 为分界面，规定三个方向均为正方向.上层空间为空气，空气层厚度1 m.下层空间为土壤，并且在分界面下方0.5 m 处设置水平层厚0.1 m 的矿体.在分界面上方0.125 m 处，沿x方向两端设置50 个接收点，每点相距0.08 m，接收点中心位置坐标为(2，2，1.125).发射频率设置为40 MHz，电流强度为25 A，长度为1 m的发射线源沿y轴放置于分界面上方0.01 m 处，其中心位置坐标为(2，2，1.01).表1 给出了算例1 具体的电性参数[25].

图1 分层介质GPR 模型Fig.1 Layered media GPR model

表1 算例1 电性参数Tab.1 Electrical parameters of case 1

由图2 可知，两种边界条件的数值模拟结果均表明HO-FEM 算法应用在GPR 数值模拟时具有较高精度.其中LMGF 求解器为WCT 中分层介质格林函数，主要应用于分层介质的电磁仿真.

图2 算例1 各接收点处磁场强度实部和虚部模拟结果Fig.2 Simulation results of the real and imaginary parts of the magnetic field strength at each reception point(case 1)

本文算法反射误差理论计算公式采用文献[26]方法，如式(18)所示：

式中：Hi为 观测点对应的磁场强度；为解析解中每个观测点对应的磁场强度；为解析解中磁场强度的最大值.

观察图3 的反射误差可知，由于频移参数 α的影响，CFS-PML 相比于传统PML 计算的各分量能够降低10～30 dB.由图4 可以看到，Hx、Hy、Hz相对误差分别小于0.72%、0.63%、0.46%，表明CFS-PML 吸收边界有效地抑制了人工边界效应，使得吸收性能更加优异.因此，本文所使用的CFS-PML 吸收边界在整体的接收点处所得的效果要优于PML 吸收边界.

图3 算例1 各接收点处磁场强度实部和虚部反射误差对比Fig.3 Comparison of real and imaginary reflection errors of magnetic field intensity at each reception point(case 1)

图4 算例1 各接收点处磁场强度实部和虚部相对误差对比Fig.4 Comparison of the relative errors of the real and imaginary parts of the magnetic field strength at each reception point(case 1)

与此同时，在本次算例中，CFS-PML 方法与传统PML 方法的计算耗时分别为18 044 s 和21 734 s，时间节省约为17%；内存消耗分别约为9.8 GB 和11.0 GB，内存节省约为11%.

2.2 起伏地岩GPR 模型

算例2 三维起伏地岩GPR 模型如图5 所示，计算区域为4 m×4 m×4 m.以o为坐标原点，z=2.5 m为分界面，规定三个方向均为正方向.上层空间为空气，空气层厚度1.5 m，下层空间土壤厚度为1.1 m，地岩隆起部分上方宽度为0.5 m，隆起高度为0.12 m，地岩未隆起的厚度为1.4 m.在模型结构上采用的是不规则网格剖分，假设在岩层下方存在一个棱柱形异常体，用以模拟地质空洞.发射中心频率设置为62 MHz，电流强度为25 A，长度为1 m 的发射线源沿y轴放置在分界面上方0.01 m，发射线源中心位置坐标为(2，2，2.51).接收点中心位置坐标为(2，2，2.625)，采用50 个接收点，间隔为0.08 m，沿着x方向展开长度为4 m.表2 给出了算例2 具体的电性参数[25].采用商业软件WCT 中CN-FDTD 算法与本文提出的算法进行对比.

表2 算例2 电性参数Tab.2 Electrical parameters of case 2

图5 起伏地岩GPR 模型Fig.5 GPR model of uphill terrane

算例2 虽然受地形改变的影响，但通过对比CNFDTD 算法结果，HO-FEM 算法仍能保持高精度.在测试过程中可以看到，PML 吸收效果并不能达到100%的吸收效率，但从实验结果来看，图6 观察到磁场强度变化较快的区域，CFS-PML 仍能很好地吻合CN-FDTD 结果，表现出较高的精度.图7 展示了传统PML 与CFS-PML 在观测点的反射误差，从中可以看到CFS-PML 的反射误差相比PML 较小.三维起伏地岩模型含有复杂有耗的结构，并且采用的不规则网格剖分也会使得传统的PML 边界产生一定的误差，然而CFS-PML 中的频移参数 α使截断区域产生频率相关的衰减，提高掠射角处吸收性能，使得CFS-PML 在传统的PML 基础上反射误差能够降低10 dB 左右，并且相对误差能够控制在合理的范围.从图8 可以看出CFS-PML 边界条件下Hx、Hy、Hz相对误差分别小于0.99 %、0.52 %、0.74 %，也体现了较高的精度.与此同时，算例2 中，CFS-PML 方法与传统PML 方法的计算耗时分别为54 243 s 和77 211 s，时间节省约为30%；内存消耗分别约为14.6 GB 和16.8 GB，内存节省约为14%.

图6 算例2 各接收点处磁场强度实部和虚部模拟结果Fig.6 Simulation results of the real and imaginary parts of the magnetic field strength at each reception point(case 2)

图7 算例2 各接收点处磁场强度实部和虚部反射误差对比Fig.7 Comparison of real and imaginary reflection errors of magnetic field intensity at each reception point(case 2)

图8 算例2 各接收点处磁场强度实部和虚部相对误差对比Fig.8 Comparison of the relative errors of the real and imaginary parts of the magnetic field strength at each reception point(case 2)

3 结论

本文在频域HO-FEM 算法的电磁场模拟中成功实现CFS-PML 人工边界条件，极大地改善了PML 的吸收性能，提高了电磁场GPR 数值模拟的计算效率.与商业软件WCT 数值模拟结果对比表明，CFSPML 边界条件的应用具有更高的计算精度和效率，验证了本文所提出算法的准确性和可靠性.两个算例的数值结果表明至少能够降低10 dB 的反射误差，计算效率上分别提升17%、30%.

值得注意的是，本文公式推导仅考虑各向同性介质，并未考虑各向异性介质材料.此外，本文基于CFS-PML 的HO-FEM 并未采用加速并行等技术，以后这将成为我们研究的方向.