基于复频移完全匹配层的辛龙格库塔算法探地雷达正演

冯德山，谭阳，杨军，张陆军，袁忠明，柳杰，王珣*

(1.中南大学地球科学与信息物理学院，长沙 410083；2.有色资源与地质灾害探查湖南省重点实验室，长沙 410083；3.广州市市政工程设计研究总院有限公司，广州 510060)

探地雷达(ground penetrating radar，GPR)是一种利用电磁波的穿透能力来测定地下结构中物质分布的一种无损检测方法，具有无破损、高效率、分辨率高、结果直观、工作方式灵活等特点[1-2]。探地雷达作为地球物理勘察和岩石工程行业中一种有价值的工具，目前被广泛应用于道路评估[3]、水利工程[4]、地雷检测[5]、土木工程[6]、考古学[7]、航空航天[8]、军事[9]等领域。GPR正演模拟是通过构建地下结构模型，模拟电磁波在地下目标体传播时的反射、透射、绕射现象，获取目标体回波信号的一种数值方法[10-11]。正演模拟是研究探地雷达电磁波在地下结构中传播规律的有效手段，也是解释探地雷达实测数据的重要依据。此外，正演模拟仍然是基于探地雷达实测信号反演地下目标结构参数的基础，准确高效的正演模拟对提高反演精度和速度具有重要意义。GPR电磁波数值模拟方法主要包括：射线追踪法(ray tracing method，RTM)、积分方程法(integral equation，IE)、时域有限差分法(finite difference time-domain，FDTD)、交替方向隐式(alternative direction implicit，ADI)法、有限单元法(finite element time-domain，FETD)、时域伪谱法(pseudo spectral time-domain，PSTD)、有限体积法(finite volume time-domain，FVTD)等。但这些方法在效率和精度等方面总存在一些问题。比如几何射线法[12]忽略了电磁波的运动特性，不能反映绕射现象；FDTD法[13-14]计算速度快，但计算效率受到Courant-Friedrichs-Lewy(CFL)稳定条件的约束，即时间步长的选取受到空间网格尺寸的制约；ADI-FDTD法[15]虽然摆脱了CFL稳定条件的约束，但较大的时间步长会增加数值频散性，受到数值误差或精度要求的限制；FETD法[16]虽然可以处理地表起伏等问题，但是计算量较大；PSTD法[17]虽然计算精度高，占用内存相对较小，但是不利于构建复杂地质体模型，不利于并行计算；FVTD法[18]虽然本身包含几何信息更容易处理复杂结构，但其算法较为复杂，不利于提高算法的精度。而辛算法是近年来发展起来的一种求解Hamilton系统的保结构算法，在其求解过程中不会引入人为耗散和虚假激励，在有关结构整体性、稳定性以及长期模拟计算方面，具有显著优势[19]。Fang等[20]利用Hamilton系统的辛龙格-库塔方法，扩展了复杂地下结构模型的雷达剖面和波传播研究。Tang等[21]研究了具有连续阶的辛分块龙格-库塔方法的构造。Liu等[22]基于辛几何理论和改进策略改进辛分块龙格-库塔格式来求解哈密顿系统。Ahmad等[23]推导了一类有效阶为三阶的辛分块龙格-库塔方法用于可分离哈密顿系统的数值积分。Yang等[24]导出了三维一阶辛分块龙格-库塔差分迭代格式，并将其应用于基于Higdon ABC边界条件的探地雷达路面检测。Sim等[25]提出了一种新的基于局部一维技术的近似解析辛分块龙格-库塔方法来数值求解二维声波方程。

完全匹配层(perfectly matched layer，PML)是由Berenger[26]提出的一种应用于电磁波截断计算的边界吸收条件。Kuzuoglu等[27]提出了复频移完全匹配层(complex frequency-shifted PML，CFS-PML)加强了对低频隐失波、大角度掠射波的吸收。Roden等[28]将CFS-PML应用到弹性波的FDTD计算中，其能够有效地应用于不均匀介质、有耗介质、频散介质、各向异性介质及非线性介质的截断计算中。Movahhedi等[29]将CFS-PML应用于Maxwell方程的时域矢量有限元的模拟中，并且开展了最优吸收效果的各项参数选取实验，有效地改善了截断边界反射。Movahhedi等[30]提出了一种新的公式来实现三维矢量时域有限单元法中用于矢量波动方程边界截断的复频移完美匹配层。李义丰等[31]将CPML吸收边界层应用于二维声波的截断计算中，并在软件COMSOL中完成频域和时域的有限元法(finite element method，FEM)计算。冯德山等[32]采用卷积方法将CFS-PML应用于时域有限元求解GPR波动方程的数值模拟中。Wang等[33]提出了基于波动方程的无网格法的递归卷积CFS-PML吸收边界条件。杨茜娜等[34]基于均匀速度模型参数选取方法，探究了复杂速度模型的CFS-PML参数设置。在求解地震波波动方程上，Ma等[35]将该CPML与辛分块龙格-库塔法和近似解析的辛方法进行结合，形成了适当同步的辛分块龙格-库塔(symplectic partitioned Runge-Kutta，SPRK)+CPML和近似解析的SPRK(NSPRK)+CPML。Zhao等[36]将CFS-PML边界条件扩展到频域有限差分地震模型中，相较于比时域模型，可以更方便地实现多源和直接扩展衰减因子。在求解电磁波波动方程上，雷建伟等[37]基于辛Euler算法和曲面共形技术，结合图形处理器(graphics processing unit，GPU)并行计算实现了非金属地下管道探地雷达电磁响应的高效精细化计算。综上，目前对于基于辛龙格库塔算法GPR正演模拟的研究较少，而且基于复频移完全匹配层的辛龙格库塔算法GPR正演模拟对吸收边界条件的关键参数的选取的研究也较少。

1 原理

1.1 哈密顿系统与SPRK方法

利用哈密顿函数H(p,q)来考虑以下正则化的哈密顿系统，公式为

(1)

式(1)中：H(p,q)为可分的哈密尔顿系统，p、q为实向量；f(q,p)、g(q,p)为向量值函数。

假设对式(1)中的第一式中p的分量用一种Runge-Kutta方法进行求解，而对第二式中q的分量用另外一种Runge-Kutta方法进行求解，这样的计算方法称为分块Runge-Kutta方法。每一个Runge-Kutta方法都具有独立的Butecher表达式。

(2)

对于辛分块Runge-Kutta方法，式(2)的系数满足以下关系：

(3)

将式(2)代入式(1)，结果表示为

(4)

式(4)中：ξ为时间增量；Pi和Qi为j级分块多项式。

(5)

1.2 控制方程

各向同性有耗介质中，Maxwell方程组表示为

(6)

式(6)中：E为电场向量；H为磁场向量；电流密度为J=σE，σ为电导率；ε为介电常数；μ为磁导率。

引入矢量磁位H=∇×A并令E=-U，则有耗介质中的广义Hamilton函数可以表示为

(7)

根据式(1)，Maxwell方程组可简化为

(8)

采用二阶中心差分近似拉普拉斯算子，则

(9)

对于二维TM波，式(8)可以写为

(10)

式(10)中：Az和Uz分别为场组分A和U在z方向的分量。

电场Ez、磁场Hx和Hy可表示为

(11)

二维频率域电磁波动为

(12)

式(12)中：ω为角频率；j为虚数单位。

在PML区域，复数伸展坐标下的波动方程表示为

(13)

其中x轴被拉伸为

(14)

式(14)中：σx为PML层内x方向电导率参数；κx的引入是为了改善PML对倏逝波的吸收特性；αx的引入则是为了改善PML对低频分量的吸收特性；σx和αx为大于零的正实数，κx≥1。在PML区域之外，αx=0，κx=1。对式(13)第二式右边的项应用复合函数求导公式，并代入式(13)可得

(15)

式(15)中：κ′x、σ′x、α′x分别为κx、σx、αx关于变量x的导数。式(15)转换至时间域比较复杂，为了方便推导，引入变量p，表达式为

p=jωε0+αx+σx/κx

(16)

将式(16)代入式(15)，可得

(17)

将式(17)展开，可得

(18)

引入辅助变量e1x，令

(19)

式(13)可以改写为

(20)

对式(19)两边同时乘以p可得

(21)

式(21)中：

(22)

式(22)两边同时乘以p，可以得到

(23)

式(23)中：

(24)

式(24)两边同时乘以p，可得

(25)

将式(16)代入式(21)、式(23)和式(25)，可得

(26)

把式(26)各项分别作拉普拉斯逆变换，即

(27)

式(27)中：

(28)

(29)

(30)

1.3 CFS-PML辛龙格库塔离散

为了方便对式(27)进行求解，对于CFS-PML所引起的辅助微分方程进行更新。

(31)

式(31)中：

(32)

式(32)中：v代表x或y。

2 CFS-PML吸收边界条件最佳参数选取实验

在计算细长结构或尖角的波相互作用时，CFS-PML吸收边界条件对倏逝波具有很强的吸收能力，并能显著提高边界截断的性能[29]。因此为了获得SPRK求解GPR波动方程的CFS-PML吸收边界条件的最佳参数，本文采用狭长模型测试不同参数对CFS-PML吸收边界条件吸收效果的影响。

图1为加载CFS-PML吸收边界条件的狭长模型的仿真区域示意图。仿真区域被离散成60×250的均匀网格区域，并且在仿真区域周围加载了30层CFS-PML环绕。网格间距为0.04 m，空气层厚度为0.4 m，介质层厚度为2.0 m。介电常数设为4.0，电导率设为0.001。采样时间间隔为0.08 ns。为了更加贴近探地雷达数值模拟实际情况，将激励源放在空气层进行数值模拟。S是激励点源；P1、P2、P3是参考点，其中P1和P3是大角度入射下设置的参考接收点，而P2是垂直入射下设置的参考接收点。假设左上角的网格坐标设为(1,1)，S、P1、P2、P3的网格坐标分别为(155,40)、(35,40)、(155,85)、(35,85)。仿真区域采用的激励源为100 MHz的雷克子波。

图1 狭长模型的仿真区域示意图

要衡量PML吸收边界条件的吸收效果，定义电场的总能量为

(33)

式(33)中：u(i,j)为网格点(i,j)内的电场值。定义反射误差为

Error=20lg(|ES-Eref|/|Erefmax|)

(34)

式(34)中：ES为观测点处的实测接收信号；Eref为参考模型信号；Erefmax为当仿真区域扩展的足够大(即在理想的仿真时间内，边界产生的误差不会被参考点接收到)时参考区域信号中振幅最大的值。

参考区域是在原有模型的基础上保证激励点源、参考点位置不变然后扩大计算区域得到的，采取12倍扩大区域作为参考区域。参考区域的作用是使参考点在设置好的理想仿真时间内接受到的波形可以看作是该时间范围内“无限大的理想模型”接收到的电磁波，达到测试模型的边界条件吸收效果作参考的目的。图2为加载CFS-PML吸收边界的参考区域示意图。参考区域被离散成720×3 000的均匀网格区域。S是激励点源；P1、P2、P3是参考点，左上角的网格坐标设为(1,1)；S、P1、P2、P3的网格坐标分别为(1 500,340)、(1 380,340)、(1 500,385)、(1 380,385)。

图2 参考区域示意图

图3为观测点不同参数m和反射系数R的最大反射误差分布图。通过分析观测点P1、P2、P3关于参数m和反射系数R的反射误差分布图可知，综合观测点P1、P2、P3，当m=2、R=10-4时，反射误差达到最小值。图4为观测点不同参数κmax和αmax的最大反射误差分布图，综合观测点P1、P2、P3，当κmax=4、αmax=0.008时，反射误差达到最小值。

图3 参数m和反射系数R的最大反射误差分布图

图4 κmax和αmax的最大反射误差分布图

图5为P1、P2、P3观测点加载CFS-PML吸收边界条件在不同时刻的反射误差对比图。可以看出，CFS-PML吸收边界条件的反射误差大部分在40 dB以内，说明CFS-PML吸收边界条件能更好地吸收反射波并减少了仿真过程中由于边界反射引起的误差。

图5 CFS-PML吸收边界反射误差图

为了更好地说明CFS-PML吸收边界条件在GPR正演截断边界处对反射波的吸收效果，采用图1所示模型，与透射边界条件进行对比实验进一步分析两种边界条件的吸收效果。图6为两种边界不同时刻的波场快照对比图。图6(c)、图6(d)、图6(e)、图6(f)中圈出的部分可以看出相同时刻加载CFS-PML的边界处几乎没有反射，但是加载透射边界条件的边界处的反射较为明显，波的形状也发生了畸变。图6(h)中圈出的部分的颜色明显比图6(g)中圈出的部分的颜色要深，说明透射边界对于角点所产生的反射的吸收效果明显不如CFS-PML。

图6 不同边界条件的波场快照对比图

图7为施加CFS-PML吸收边界条件和透射边界条件的仿真区域内电场总能量随时间的变化曲线。仿真的前期两类吸收边界的曲线基本重合，但是在仿真的后期，可以清楚地看到加载CFS-PML作为吸收边界条件的总能量低于加载透射边界作为吸收边界条件的总能量，这也表明CFS-PML吸收边界条件可以更好地消除边界反射波。

图7 仿真模型计算区域内总能量随时间的变化曲线图

图8为分别施加CFS-PML和透射边界的GPR辛龙格库塔方法计算得到的P1、P2、P3接收点的GPR信号与参考解的对比图。采用的参考解是在不施加边界条件下，同时保证在模拟时间内虚假反射不会到达接收点，将仿真区域扩大12倍后得到参考区域后采用相同的数值方法计算得到的，激励源和接收点的相对位置保持一致。由图8明显可见，垂直入射得到的参考接收点P2在加载不同吸收边界条件下的GPR信号都与参考解吻合很好，而大角度掠射下得到的参考接收点P1和P3在加载不同吸收边界条件下的GPR信号都与参考解吻合的不是很好，其中P1的吻合效果更差。而且P1点施加透射边界下的GPR信号中由于边界产生的虚假反射波清晰可见，而施加CFS-PML的GPR信号与参考解吻合很好，未见明显的虚假反射。

图8 不同边界条件下接收点获得的单道波形与参考解的对比

图9为不同边界条件下接收点获得的单道波形与参考解的相对误差对比。图9用具象化的相对误差来对比不同边界条件下的吸收效果之间的差异。以参考模型为例，参考波形的最大值即为归一化值，对参考波形进行归一化就是除以该值。将加载不同的边界条件的仿真模型归一化与参考模型归一化的差值，即为相对误差。加载不同的边界条件的误差不同，可以从相对误差值之间的差异来评估加载不同的边界条件的吸收效果。从图9明显可以看出入射角度不同对CFS-PML吸收边界的吸收效果影响效果较小，但是会很大程度上影响透射边界的吸收效果。

图9 不同边界条件下接收点获得的单道波形与参考解的相对误差对比

3 CFS-PML在SPRK探地雷达数值模拟中实例

为了更好地说明CFS-PML边界条件在GPR正演截断边界处对反射波的吸收效果，也为了验证该算法是否可以很好地拟合了起伏不平的地下界面，因此基于MATLAB平台编写了CFS-PML边界的辛龙格库塔GPR正演程序。建立如图10所示的起伏地形地电模型，采用规则网格剖分，剖分网格为160×390，网格间距为0.05 m，最上层的空气层厚度为0.5 m，第一起伏层介电常数为7.0，电导率为0.000 1，第二层起伏层介电常数为10.0，电导率为0.001，第二层起伏层中的存在两个圆形异常体，左侧异常体介电常数为6.0，电导率为0.000 1，右侧异常体介电常数为5.0，电导率为0.000 1，背景介质介电常数为15.0，电导率为0.001。考虑到二维GPR探测中经常采用剖面法与宽角法两种数据采集方式，因此，应用辛龙格库塔方法对该模型进行正演时也基于剖面法和宽角法进行展开。

图10 起伏地形模型

剖面法采用自激自收方式，自左向右同步移动。时间步长为0.1 ns，时窗长度为200 ns。激励源仍为100 MHz的零相位Ricker子波，为了更贴合实际，将激励源位置放在空气层中，距离地表5个网格。模拟区域的CFS-PML吸收边界层数为30层。图11(a)为基于辛龙格库塔算法加载CFS-PML边界条件正演接收到的200道波形组成的雷达剖面图。由图11(a)可见，CFS-PML对边界的截断效果比较理想，在正演剖面图中未出现人工截断边界的干扰。在30 ns的位置，从左至右有5处较为明显的绕射，分别对应起伏上界面隆起的4个脚点以及起伏上界面的最低点，上界面的起伏形态与位置在剖面图中亦清晰可辨。在70 ns的位置，更深部埋深为4 m的圆形异常体产生的绕射波都清晰可辨，左侧异常体上下界面的反射波较为清晰，且上界面能量相对较强，而右侧异常体的绕射波可能由于上界面起伏最低点的影响而产生了形变。还可以看到起伏下界面存在明显的多处饶射。

再考虑该地电模型宽角法的正演结果。将发射天线置于模拟区域的正上方，接收天线置于发射天线两侧，共记录200道波形。时间步长为0.1 ns，模拟时窗长度为200 ns，激励源仍为100 MHz的零相位Ricker子波，激励源的位置距离地表5个网格。采用相同的吸收边界，图11(b)即为宽角法接收到的正演剖面。

从图11(b)宽角法得到的正演剖面可以看出，由于激励源位于空气层中，而且激励源与地表之间存在一定的距离，因此不仅空气中的直达波界面反射清晰可见而且地表的直达波界面反射也清晰可见，但是该模型起伏地形的分界面无法像剖面法得到的分界面那样明显。

图11 起伏地形模型雷达正演剖面

图12为激励源位于中心点的不同时刻的波场快照。图12(a)中，电磁波在空气和地表的产生了分离，在空气中的波速较快，波长较长。图12(b)中，电磁波进入地下，但是由于激励源与地表之间存在一定的距离，因此电磁波波前的两侧比中部宽。图12(c)中，电磁波波前遇到起伏的上界面发生反射。图12(d)中，由于该地电模型上界面是起伏界面，电磁波传播至起伏界面的时间不同，因此透射波的波前是不规则的波形。图12(e)中，由左侧的圆形异常体引起的绕射波清晰可见。图12(f)和图12(g)中，由右侧的圆形异常体引起的绕射波清晰可见。图12(h)中，由于下层起伏界面的影响，波的形状呈不规则状。

图12 起伏地形模型不同时刻波场快照

4 结论

提出了一种基于CFS-PML吸收边界条件的SPRK算法，将其应用于GPR二维正演模拟中，详细对比分析了CFS-PML吸收边界条件与透射边界条件在入射角度不同的情况下的吸收效果，并将该算法成功应用于复杂地质模型。结论如下。

(2)以一个狭长模型为例，进行了CFS-PML吸收边界条件的最佳参数选取实验。从反射误差图中得出m=2，R=10-4，κmax=4，αmax=0.008时CFS-PML吸收边界条件吸收效果最好，而且不同入射角度下的接收参考点误差大都在40 dB以内。与透射边界条件的对比实验可以看出，入射角度不同不会影响CFS-PML吸收边界的吸收效果，但是会很大程度上影响透射边界的吸收效果。

(3)以1个起伏地形的GPR地电模型正演为例，得到了一个高精度的数值模拟结果，说明所提出的算法对于地下任意复杂介质有着较好的适用性。