2016, Vol. 59
2. 中国科学院地质与地球物理研究所, 北京 100029;
3. 有色金属成矿预测教育部重点实验室, 长沙 410083
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Changsha 410083, China
探地雷达(Ground Penetrating Radar，GPR)是一种对地下或物体内部不可见目标体进行定位的电磁无损探测技术，它根据地下介质的电性差异，在雷达剖面上进行成像，可分析地下介质的电性结构及空间形态(曾昭发等，2010；冯德山等，2014).为了更好地刻画雷达波动现象、理解其传播规律与物理过程，有必要开展GPR波动方程正演.GPR正演的模拟方法有许多，包括李展辉等(2009)将基于Runge Kutta方法的时间步迭代求解错格时域伪谱法应用于三维井中GPR模拟，有效解决了时域伪谱法的Gibbs现象；方宏远等(2012, 2013)利用Hamilton系统的辛分块龙格库塔方法开展了道路结构层雷达波传播特征研究；底青云和王妙月(1999)推导了含衰减项的雷达波有限元方程，实现了GPR波的有限单元法正演；冯德山等(2012)提出了一种结合透射边界与Sarma边界的混合边界条件，改善了有限单元法GPR正演的边界反射.尽管GPR数值模拟方法种类许多，但时域有限差分法(Finite Difference Time Domain，FDTD)因具有直接时域计算、编程容易、节约存储空间和计算时间、适合并行计算等优点仍占据最重要的地位，得到了广泛的应用.Xu (1997)、Giannopoulos (2005)、Irving和Knight (2006)、Cassidy和Millington (2009)、李静等(2010)、冯德山等(2010)应用FDTD法进行GPR数值模拟时，由于计算机内存有限，模拟网格区域必须截断为有限区域，为了在截断边界处不引起虚假反射，通常在网格外围加上吸收边界条件(Absorbing Boundary Condition，ABC)来削弱边界的反射(Berenger，1999).早期提出的Taflove & Brodwin (1975)边界、Engquist-Majda边界(1977)、Liao边界(1984)有效地减弱了边界反射；Mur (1981)在单行波的基础上提出了Mur二阶吸收边界，得到了广泛的应用；Mei和Fang (1992)提出了超吸收边界条件，它是对单行波方程为基础的吸收边界条件的一种有效改进，但是在模拟区域的角点处未能达到很好的改善.这些边界条件通常在模拟区域的外边界仍有0.5%~5%的反射系数，他们逐渐被完全匹配层(Perfectly Matched Layer，PML)所取代.
PML边界条件最初是由Berenger (1994)提出的，它在FDTD网格外边界加载一种非物理的吸收媒质，可以使电磁波在匹配层中无反射地通过并按指数规律衰减，相比于Mur二阶吸收边界，完全匹配层的吸收效果可以提高几个数量级，将边界条件的研究向前推进了一大步.Katz等(1994)将Berenger PML的使用范围延伸到三维模拟区域；Moerloose和Stuchly (1995)通过分析和数值实验指出Berenger PML对倏逝波的吸收毫无效果，且在低频情况下，会存在很大的数值误差；Berenger (1997)提出由倏逝波产生的数值反射主要存在于与PML电导率有关的截止频率之下，并且指出了Moerloose和Stuchly论文中仅考虑了传输波平行于PML界面的情况，理论上，倏逝波越强，PML对其的吸收效果也越好，但是由于时间上和空间上的离散，效果不尽人意；Berenger (1999)对波导应用中的倏逝波的反射进行了系统的讨论，并且对电导率递增的PML层的倏逝波反射系数进行了研究.Berenger PML的理论体系是非Maxwell方程的，物理机制模糊，同时，其电磁场分量分裂技术增加了计算内存与数值实现的难度(葛德彪和闫玉波，2011；2005)，而且只对行波有吸收效果，对空域中衰减的隐失波、低频波以及入射角度较小的掠角波无吸收效果.为了改善边界条件的吸收效果.Sacks (1995)和Gedeny (1996)提出了单轴各向异性完全匹配层(Uniaxisal PML，UPML)，UPML边界无须对电磁场分裂，计算更简洁、易编程实现，因其在吸收系数中引入了线性的吸收因子，UPML边界对隐失波、凋落波等干扰波有一定的吸收效果，能适用于有耗、无耗、色散等介质中.詹应林(2008)将UPML应用于二维FDTD探地雷达正演中；冯德山等(2011)将UPML边界应用于ADI-FDTD差分中实现了GPR的三维数值模拟.尽管UPML边界对于凋落波、隐失波、低频波的干扰波有一定的吸收效果，但仍然不是特别完美.为了克服UPML低频区性能变差的缺点，Kuzuoglu和Mittra (1996)考虑到复频移拉伸函数的PML技术可以吸收低频信号在边界表面的反射，同时对隐失波等大角度掠射波也具有好的吸收效果，为此，提出了复频偏移技术(Complex Frequency Shifted，CFS)加强了对低频、隐失波的吸收；Roden和Gedney (2000)将复频移拉伸函数加入到电磁波模拟的PML方程，有效压制了信号的假反射；Berenger (2002)验证了CFS PML能够对低频区倏逝波取得非常优越的吸收效果；Wang和Liang (2006)将复频移PML边界(Complex Frequency Shift Perfectly Matched Layer, CFS-PML)应用到ADI-FDTD中；李展辉和黄清华(2014)将CFS-PML应用在瞬变电磁法正演中，以替代传统狄利克雷边界条件；Drossaert和Giannopoulos (2007)以及Dimitri和Roland (2007)提出了基于递归积分的波场非分裂的复频移PML边界(Recursive Integration CFS-PML, 简称CFS-RIPML)；Zhang和Shen (2010)及Martin等(2010)将基于辅助微分方程的CFS-PML应用到弹性波数值模拟中，该方法易于扩展到高阶差分形式；张显文等(2009)将递归复频移PML应用到交错网格高阶差分弹性波波动方程正演中；Giannopoulos (2008)及Li等(2012)将基于CFS-RIPML吸收边界技术应用在电磁波波动方程的有限差分正演中.Roden和Gedney (2000)在CFS-PML基础上应用卷积技术提出了卷积完全匹配层(Convolutional PML，CPML)，大大提高了运算效率.张鲁新等(2010)将不分裂的完全匹配层(CPML)结合旋转交错网格有限差分技术应用到孔隙弹性介质模拟中；相比常规PML、UPML吸收边界，CPML对于隐失波、低频波等具有更好的吸收效果，所以CPML具有非常广阔的发展前景.
本文分析了PML边界中各参数的设置，采用波场快照着重探讨了二、三维UPML、CPML对倏逝波的吸收效果对比.然后，采用GPR二维剖面法和宽角法正演剖面直观地展示了UPML和CPML两种边界条件下倏逝波的反射对正演剖面的影响，同时在三维情况下展示了复杂模型的波场快照，直观地显示了倏逝波的形态.结果显示，CPML可有效减少边界反射误差，并能取得满意的精度，为GPR的正演提供了合适的吸收边界.
2 卷积完全匹配层原理CPML吸收边界条件是建立在伸缩坐标PML基础上.在伸缩坐标中，z方向的安培环路定律为
|
(1) |
其中si为坐标伸缩因子：
|
(2) |
把公式(1)从频域转换到时域，由于坐标伸缩因子与频率有关，因此转换到时域后公式右边会存在卷积形式，即
|
(3) |
其中，sx和sy分别为1/sx和1/sy的逆拉普拉斯变换.公式(3)右侧的卷积采用Luebbers提出的递归卷积来加速运算(李建雄，2007).为了改善CPML对倏逝波的吸收特性，Kuzouglu和Mittra (1996)提出伸缩坐标因子的选择为
|
(4) |
式(4)中αi和σi为正实数，κi为大于1的实数.κi值的引入是为了改善PML对倏逝波的吸收特性，αi值的引入则是为了改善PML对低频分量的吸收特性.根据拉普拉斯变换理论，si的脉冲响应为
|
(5) |
其中，δ(t)和u(t)分别是单位冲击函数和单位阶跃函数，将(5)式代入到(3)式中，有
|
(6) |
直接在递推方程(6)中计算时域卷积的效率是很低的，为了提高卷积计算效率，将ξi(t)离散冲击响应定义为
|
(7) |
式(7)中
|
(8) |
根据公式(7)和(8)，按照Yee氏网格将(6)式进行时间和空间上的离散，可得
|
(9) |
公式(9)中的离散卷积计算起来十分复杂，但是，由于Z0i(m)是简单的指数形式，从而它们的和可以使用递归卷积来递归得到(葛德彪和闫玉波，2011；2005).因此，引入一组新的辅助表达式ψi.
|
(10) |
式(10)中
|
(11) |
|
(12) |
|
(13) |
其中ai已在式(8)中给出，进一步整理公式(10)，可得
|
(14) |
式(14)中标号m=(i, j, k+1/2)，CA、CB与常规FDTD中表示的内容相同.其余场量可类似推出.由式(14)可看出，除κ项以及后面的ψ项外，公式(14)前面部分与常规FDTD公式十分相似.故编程时可令模拟区域内κ=1，而PML层内令κ渐进变化，从而将模拟区域与PML边界层中采用统一的差分公式，然后，再在PML区域中加上后面的ψ项，这种编程方式可使程序更简洁, 有效减少计算量.为了减少反射误差，在PML层内σ、κ、α参数分布为非均匀的，各参数单调变化，以z方向为例，通常取
|
(15) |
其中d为PML层网格数，z0为PML与模拟区域分界处的网格编号.由于电场和磁场在Yee网格中空间位置相差半个网格，故σ、κ、α等参数也需在PML网格上根据电场和磁场的对应位置来分别设置，同时，z0代表的分界处对于电场和磁场也有所不同，需要相差半个空间网格.
3 完全匹配层对倏逝波的吸收效果在PML边界条件中，Berenger PML、UPML、CPML这3种完全匹配层应用较广泛，对于传输波都具有很好的吸收效果，但对倏逝波的吸收，尤其是低频倏逝波的吸收，却有着截然不同的效果.倏逝波(evanescent wave)又称凋落波、隐失波(Berenger，1998).以二维TM波为例对比这3种边界条件对倏逝波的吸收能力.假设真空中雷达传输波的传播方向为X，倏逝波的方向为Y，X与坐标轴x轴夹角为θ，其中θ的选取为负数，任一场量的场值都可写为(Berenger，2000)
|
(16) |
指数项e-(ωsinχ/c)Y为倏逝波的自然衰减项.在Berenger PML中，假设PML层垂直于x轴，即σy=0，任一分量的场值可写为
|
(17) |
式(17)中为引入的衡量倏逝波强度的参数，
|
(18) |
分析衰减项表达式，结合倏逝波coshχ >1可知，理论上，PML对倏逝波的吸收能力较传输波更强，倏逝波越强，衰减越快.同时衰减项与入射角θ存在联系，入射角越大，对倏逝波的吸收能力越弱.当倏逝波足够强时，即coshχ足够大，倏逝波可以在极小的距离x内充分吸收.但是在FDTD离散网格中，在cosh χ足够大时，倏逝波甚至可在小于一个网格距离内被吸收.显然，有限的FDTD网格无法完成变化巨大的倏逝波的采样工作，所以必然会造成剧烈的反射误差，尤其是在低频区.Berenger (2000)通过数值计算验证指出，在低频区，真空与均匀PML分界面上，传输波相对于倏逝波的反射系数更小，而倏逝波越强，其反射系数也越大.所以倏逝波强烈反射的关键原因在于PML层的离散化.在Berenger PML中，可通过增加PML层的厚度来压制倏逝波.
对于UPML，其关键参数为s=κ+σ/jωε0，当κmax=1时，也即κ=1，可以证明，与Berenger PML完全等价.在离散FDTD中计算时，UPML采用的是半隐式差分，而Berenger PML在PML中采用的是指数差分，而且编程平台的精度有限，故在数值计算时，难免会存在舍入误差.但两者的数值误差可以忽略不计，所以在数值上，κmax=1时的UPML与Berenger PML也等价.在UPML中，当κmax增大时，对倏逝波的吸收能力也加强，但随之而来的代价是，由于κ的离散化，导致它在PML网格间递增值也逐渐变大，UPML的全域误差也剧烈增大.
CPML中，坐标伸缩系数为si=κi+σi/(αi+jωε0)，很显然，fα=αi/2πε0为一关键的分界参数，当f>>fα时，si趋近于si=κi+σi/jωε0，此时CPML与UPML无异.特别地，当κi=1时，丨ψPML丨变为
|
(19) |
可以看到公式(19)与公式(17)相同，即在高频时，CPML与Berenger PML、UPML的作用相当.当f < < fα时，si趋近于si=κi+σi/αi，即si变为实坐标伸缩因子.CPML中任一场分量的值可写为
|
(20) |
特别地，当κi=1时，将fα的值代入，丨ψPML丨可变为
|
(21) |
倘若将ψvacuum中的自然衰减项合并进去，则式(21)中总的衰减相当于倏逝波在距离(1+σx/αx)x上的自然衰减.分析式(21)可知，当f < < fα时，CPML并不能对传输波(sinhχ=0)起到吸收作用，但是对倏逝波能起到出色的吸收作用(sinhχ>0)，这是其他PML所没有的优点.公式(21)中的指数衰减项要小于公式(19)中的指数衰减项.对于Berenger PML来说，正是因为倏逝波衰减太快而不能被网格合理采样，所以在低频区吸收性能较差，而CPML中的指数衰减项的衰减能力适中，虽然还是存在一定的反射，但是基本上能够被网格合理采样，所以CPML对倏逝波的吸收能力要强于Berenger PML.注意到式(21)中，指数系数与频率相关，即f和sinhχ存在一定的变动，但实际上，在频率在fα之下的区域，fsinhχ接近于常数，故CPML在低频区能取得稳定的倏逝波吸收效果.
4 倏逝波的吸收数值实验从Berenger PML、UPML、CPML在高频区时的波场幅值公式中可看出，3种PML对倏逝波的吸收性能都与cosθ相关.当θ趋近于π/2时，cosθ趋近于0，此时传输波的方向与PML表面接近于平行时，PML对倏逝波的吸收效果最弱.为了清楚地反应倏逝波的反射，分析和对比倏逝波吸收的相关参数，开展大角度入射下的TM波实验.将激励源置于模拟区域的左上角距左分界面和上分界面一个网格处，模拟区域统一设为真空.激励源为主频900 MHz的blackman-harris脉冲，取空间步长为0.012 m，时间步长为0.012 ns，UPML中参数κmax分别等于1、3、5、7、9、11，CPML中κmax=1、αmax=0.008.为了能够更好地呈现出低频倏逝波异常反射，图中右边的颜色棒(colorbar)统一采用lines显示方式.图 1为不同PML参数条件下400时间步的电场Ez波场快照，在图 1(a-g)图中，左分界面和上分界面处都出现比较强烈的倏逝波的反射.从图 1(b-g)中看出，随着κmax的增大，UPML对倏逝波的吸收效果逐渐改善，但并不是κmax越大越好，如果κmax过大，会因为在PML中分布时增加速率太大而造成强烈的后期虚假反射，通常取κmax=5~11为宜.注意到，在κmax=1时，UPML与Berenger PML在数学与物理是完全等价的，所以在图 1a和图 1b中Ez波形完全一致.
|
图 1 Berenger PML、UPML、CPML三种吸收边界在400时间步时Ez的波场快照 (a) Berenger PML；(b) κmax=1的UPML；(c) κmax=3的UPML；(d) κmax=5的UPML；(e) κmax=7的UPML；(f) κmax=9的UPML；(g) κmax=11的UPML；(h) κmax=1、αmax=0.008的CPML. Fig. 1 The Ez snapshot of Berenger PML、UPML、CPML in the timestep of 400 |
图 1h为边界条件采用CPML时同时刻的雷达波场快照，其中主要参数设置为κmax=1、αmax=0.008.可以看到CPML对倏逝波取得了非常出色的吸收效果.通过分析这个实验可知，Berenger PML是UPML参数κmax=1时的特例，通过选取合适的κmax值可以改善UPML边界对倏逝波吸收效果，而CPML中参数α对倏逝波吸收效果具有非常重要的作用，通过选取合适的κmax及αmax的参数组合，可以达到对倏逝波的充分吸收.
为了进一步对比Berenger PML、UPML、CPML对倏逝波的吸收效果，设置图 2所示的简单模型，模型网格数量为200×200，通过设置一个检测点来分析3种不同边界条件下的反射误差.脉冲源点A距上侧和左侧PML层仅一个网格，检测(接收)点B在A点右侧100个网格处，激励源为900 MHz的Blackman-harris脉冲，空间步长设为0.012 m，时间步长为0.025 ns.
|
图 2 三种PML对倏逝波吸收性能的检验模型及检测点分布示意图 Fig. 2 The testing model and distribution of detection point for absorption properties of evanescent wave in three kinds of PML |
在检测点B记录下600个时间步内的电场值，记为Ezn.接下来，设置一个扩大6倍网格数量为1200×1200的参考空间，其激励源点A与接收点B距离不变，由于模型的扩大，此时点A与B近似处于模拟区域中心，由于空间足够大，在观察时间内，边界的反射尚未到达，因此可将该模型中观测点B处接收到的信号设为参考模型信号Ez, refn，并求取参考模型信号中600个时间步内振幅最大的值为Ez, refn，通过以下计算公式获得反射误差(单位为dB)：
|
(22) |
应用式(22)计算3种PML所得B点反射误差如图 3所示，其中UPML中参数κmax=5，CPML中κmax=5，αmax=0.008.由图 3可看出，CPML对倏逝波的吸收效果有很大的改善，大约提高了一个数量级.至于后期反射误差剧烈抖动，推断是κ和α离散化所致.对比于Berenger PML和UPML可知，κmax=5时，UPML对倏逝波的吸收有一定的改善，但是仍达不到CPML边界条件的吸收效果.
|
图 3 B点不同边界条件的反射误差对比 Fig. 3 Relative error at point B of different boundary conditions |
全局反射误差是反映吸收边界吸收性能的一项重要参数，它能反映一种边界条件在整个区域模拟计算中的误差精度(Berenger，2002).虽然κ和α的引入有效地改善了对倏逝波的吸收效果，特别是对低频倏逝波的吸收，但是引入κ和α需要在PML中进行离散化，从而会引入全域误差.κmax和αmax越大，引入的全域误差越大.仍设模拟区域网格数为200×200，将每个时刻的不同位置的电场值记为Ezn(i, j)，参考模型网格数为1200×1200的，其每个时刻的不同位置的电场值记为Ebn(i, j)，则全域误差的计算方法为
|
(23) |
Berenger PML、κmax=5的UPML和κmax=5，αmax=0.008的CPML的全域误差如图 4所示.图中可见，Berenger PML的全域误差最小.在相同的κ情况下，CPML和UPML的全域误差在同一个数量级，CPML的全域误差要比UPML的大，但都在合理的范围之内.因此在CPML与UPML之间权衡时，考虑到对倏逝波的吸收效果，CPML作为优先考虑的边界.
|
图 4 不同边界条件的全域误差对比 Fig. 4 Contrast of global error under different boundary conditions |
为了更好地说明CPML边界条件在GPR正演截断边界处对反射波的吸收效果，基于Matlab平台编写了CPML边界的二维雷达FDTD正演程序.设置模拟区域网格数量为200×200空间步长，空间步长为0.006 m，时间步长为0.015 ns，CPML边界层厚度设为10个网格.为了验证程序的正确性，采用最简单的均匀介质模型，模拟区域内介质的相对介电常数为3.5，电导率为2.5 S·m-1.激励源为主频900 MHz的Blackman-harris脉冲，位于模拟区域正中心，CPML参数取κmax=5，αmax=0.008.
图 5为应用FDTD算法正演所得的不同时间步电场Ez传播快照.由图 5(a-b)可见，当激励源位于模拟区域中心时，电磁波以圆形方式向外传播，这与电磁波传播规律是一致的，从而证明了程序的正确性.从图 5c中可见，当雷达波传播至边界处时，由于CPML边界条件的加载，使截断边界得到了较好的处理，波形得以传播出模拟区域而没有产生强反射.图 5d说明雷达波形传播出去后，区域内部后期也未见截断边界处的反射波，说明了CPML边界条件在正演中取得了满意的吸收效果.
|
图 5 加载CPML边界的不同时刻FDTD法正演雷达波场快照 Fig. 5 The snapshot of radar wavefeild in difference moment by FDTD with CPML boundary |
由于二维GPR探测中常采用剖面方式显示雷达数据，而雷达数据采集方式主要分为两类：剖面法与宽角法.为了加深CPML对二维雷达正演剖面中倏逝波的吸收效果的理解，建立图 6所示的矩状地电模型，背景介质介电常数为11.0，电导率为0.1 mS·m-1.矩形模型上顶面埋深为0.52 m，宽度为0.16 m，厚度为0.16 m，矩状小异常体的介电常数6.0，电导率0.5 mS·m-1.应用FDTD算法对该模型进行剖面法和宽角法两种采样方式开展正演模拟.
|
图 6 矩状地电模型 Fig. 6 Rectangular geoelectric model |
为此，采用发射天线(激励源)和接收天线(记录点)距为0.24 m，自左向右同步移动.网格步长为0.004 m，时间步长为0.02 ns，时窗长度为24 ns.激励源仍为900 MHz的Blackman-harris脉冲，模拟区域的外边界分别采用κmax=1的UPML和κmax=5，αmax=0.008的CPML，激励源的位置均距上侧PML内表面1个网格.图 7分别呈现了UPML与CPML波场快照，图中白色线框即为异常体.从图 7a中可看出，在边界处存在着明显的倏逝波反射，而在图 7b中，倏逝波的反射几乎不可见.
|
图 7 矩状模型正演波场快照 (a) UPML边界Ez波场快照; (b) CPML边界Ez波场快照. Fig. 7 The snapshot of Ez with rectangular geoelectric model (a) Ez snapshot with UPML boundary; (b) Ez snapshot with CPMLboundary. |
二维情况下，倏逝波的反射开始出现的位置是在PML内表面附近且距离激励源约1个波长处，故在剖面法中，若自激自收或收发天线位置靠近时，则记录剖面上并没有倏逝波的反射出现.但倘若收发天线水平方向上存在着一定距离，则雷达剖面上会存在着明显的倏逝波反射.
图 8a为采用UPML的雷达正演剖面图，图中可见，在直达波之后存在着连续的倏逝波反射干扰，如图 8a中黑色虚线框所示.而CPML对倏逝波的压制效果比较理想，在正演剖面图 8b中未出现倏逝波反射的干扰.对比图 8(a-b)可知，由于倏逝波的反射的存在会一定程度上影响正演剖面的正确性，所以在大收发距时CPML边界要优于UPML边界.
|
图 8 矩状模型剖面法探地雷达正演剖面 (a) UPML边界; (b) CPML边界. Fig. 8 The forward modelling profile of GPR using profile method to rectangular geoelectric model (a) UPML boundary; (b) CPMLboundary. |
再考虑该矩状模型宽角法的正演结果，将发射天线置于矩状模型的正上方，接收天线置于发射天线两侧，共记录200道波形.设置网格步长为0.004 m，时间步长为0.02 ns，模拟时窗长度为24 ns，激励源仍为900 MHz的Blackman-harris脉冲，激励源的位置距离上侧PML内表面1个网格.图 9(a-c)分别为加载κmax=1与κmax=5的UPML和κmax=5，αmax=0.008的CPML边界条件雷达正演剖面.图 9a中UPML内表面处的倏逝波的反射依次到达记录点，在正演剖面中形成黑色虚线所圈的“扫帚形”的干扰波，即使采用图 9b中κmax=5的UPML倏逝波反射有一定的减弱，但在正演剖面中仍然存在；分析图 9c可见，加载了CPML边界的宽角法雷达正演剖面中未出现倏逝波反射的干扰，说明CPML边界条件对倏逝波具有良好的吸收特性.
|
图 9 矩状模型宽角法探地雷达正演剖面 (a) κmax=1的UPML边界; (b) κmax=5的UPML边界; (c) κmax=5，α=0.008的CPML边界. Fig. 9 The forward modelling profile of GPR using wide angle method to rectangular geoelectric model |
实例分析三维GPR正演中CPML对倏逝波的吸收效果，设置三维模拟区域大小为200×200×200个空间步长，介质为真空介质，空间步长为0.01 m，时间步长为0.0167 ns.激励源与上例一致，为了较好地体现倏逝波的反射情况，方便分析边界条件对倏逝波的吸收效果，将源置于x、y平面中心，且距上侧PML层1个网格，便于雷达波在大角度情况下入射到PML边界.PML吸收边界分别加载κmax=1的UPML和κmax=5，αmax=0.008的CPML.
图 10(a-c)为加载UPML边界条件下180、260、340三个时间步的Ez波场快照；图 10(d-f)为加载CPML边界条件下180、260、340三个时间步的Ez波场快照.分析图 10(a-c)可知，随着雷达波前向外扩散，加载UPML边界的三维雷达波形中出现了明显的倏逝波反射现象，并且倏逝波的反射沿着UPML的内表面向外扩散.而图 10(d-f)加载了CPML边界的三维波形中倏逝波能量非常弱、几乎不可见.两者对比说明：相比于UPML，CPML对倏逝波具有更好的吸收效果.
|
图 10 加载UPML、CPML边界在不同时间步时电场Ez的三维雷达波场快照 (a) 180个时间步的UPML波场快照; (b) 260个时间步的UPML波场快照; (c) 340个时间步的UPML波场快照; (d) 180个时间步的CPML波场快照; (e) 260个时间步的CPML波场快照; (f) 340个时间步的CPML波场快照. Fig. 10 The Ez wavefield snapshot in different moment with UPML, CPML boundary (a) The UPML wavefield snapshot in the moment of timestep 180; (b) The UPML wavefield snapshot in the moment of timestep 260; (c) The UPML wavefield snapshot in the moment of timestep 340; (d) The CPML wavefield snapshot in the moment of timestep 180; (e) The CPML wavefield snapshot in the moment of timestep 260; (f) The CPML wavefield snapshot in the moment of timestep 340. |
为突出在复杂介质下CPML吸收边界对倏逝波的吸收效果，建立如图 11所示的三维复杂模型.背景介质的介电常数为5，电导率为0，其大小为1 m×1 m×1 m.模型中存在两个异常体，分别为边长0.03 m的正方体和直径0.06 m的球体.异常正方体与异常球体的介电常数为3.0，电导率为3 mS·m-1.激励源与上例一致，同样将源置于x、y平面中心.PML吸收边界分别加载κmax=1的UPML和κmax=5，αmax=0.05的CPML.图 12a为加载UPML边界条件下的Ez波场快照；图 12b为加载CPML边界条件的Ez波场快照；图 12c为z方向上的两道切片，可以看到第一道切片中倏逝波的反射呈规则的同心圆状，在第二道切片中可以看到两个异常体的明显散射；图 12d中同样有两道切片，在第一道切片中倏逝波的反射明显减弱.由此可见，雷达波遇到正方体异常和球体异常时会产生散射.UPML边界的波场快照中，上部存在着明显的倏逝波反射.而CPML边界的波场快照中，倏逝波的反射明显要小.所以在GPR正演中，倏逝波的反射会严重影响波场快照和正演剖面的正确性，有必要采取更好的吸收边界来压制倏逝波的反射.
|
图 11 三维情形下正方体与球体异常的复杂模型 Fig. 11 Cube and sphere objection in three-dimension |
|
图 12 加载UPML、CPML边界下复杂模型的三维雷达波场快照 (a) UPML波场快照; (b) CPML波场快照; (c) UPML波场快照，z方向的2个切片; (d) CPML波场快照，z方向的2个切片. Fig. 12 Three-dimensional radar snapshots of complex models under UPML and CPML boundary (a) The UPML wavefield snapshot; (b) The CPML wavefield snapshot; (c) The UPML wavefield snapshot; (d) The CPML wavefield snapshot. |
(1)探讨了Berenger PML、UPML、CPML三种PML边界条件的内在关系，对比了它们对倏逝波的吸收效果.Berenger PML由于阻抗匹配条件的特殊性，主要应用于真空介质.理论上，Berenger PML对倏逝波的吸收能力要强于雷达传输波，但实际中由于FDTD的离散，特别是在低频区，Berenger PML对倏逝波的吸收能力要远差于传输波，但可以通过加大PML层的厚度来压制倏逝波的反射.相比于Berenger PML，UPML不需分裂场量，而且编程相对简单，其适用范围更广.CPML为基于复频移和卷积计算技术的完全匹配层，对倏逝波，特别是低频倏逝波，有着出色的吸收效果.
(2)详细探讨了PML边界条件中关键参数κ和α对倏逝波吸收效果的影响.实验表明：UPML中，参数κmax通常在1~11范围内取值较为合适，当参数κmax=1时的UPML与Berenger PML在数学上和物理上完全等价.随着κmax的增加，UPML对倏逝波的吸收能力加强，但是当κmax过大超过取值范围时，其离散所引起的全域误差会影响计算精度.而CPML中有两个关键参数κ和α，其中αmax通常在0~0.05范围内，当αmax过大时，也会增大模拟区域中的全域误差，与UPML相同的κ值情况下，CPML的全域误差要大于UPML，但两者全域误差在同一个数量级.通常二维情况下CPML在κmax=1，αmax=0.008左右时可取得不错的倏逝波吸收效果.
(3)分别将UPML与CPML加载于GPR正演的边界条件处理中，了解到倏逝波的反射常出现在PML内表面的特性，验证了倏逝波反射的存在，而且倏逝波的反射通常与雷达传输波入射到PML层的角度有关，入射角越大，倏逝波的反射越强.倏逝波的存在会对雷达正演剖面的正确性造成干扰，因此，有必要压制倏逝波的反射.实验证明：CPML不仅对雷达传输波有着良好的吸收能力，对倏逝波也能取得不错的吸收效果.故在GPR正演的截断处，推荐加载CPML边界条件.
| Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2): 185-200. DOI:10.1006/jcph.1994.1159 | |
| Berenger J P. 1996. Perfectly matched layer for the FDTD solution of wave-structure interaction problems. IEEE Transactions on Antennas and Propagation, 44(1): 110-117. DOI:10.1109/8.477535 | |
| Berenger J P. 1997. Numerical reflection of evanescent waves from perfectly matched layers. IEEE Antennas and Propagation Society International Symposium, Montreal, Canada, July, (2): 1888-1891. | |
| Berenger J P. 1998. An effective PML for the absorption of evanescent waves in waveguides. IEEE Microwave and Guided Wave Letters, 8(5): 188-190. DOI:10.1109/75.668706 | |
| Berenger J P. 1999. Evanescent waves in PML's: Origin of the numerical reflection in wave-structure interaction problems. IEEE Trans. Antennas Propag, 47(10): 1497-1503. DOI:10.1109/8.805891 | |
| Berenger J P. 2000. Numerical reflection of evanescent waves by PMLs: Origin and interpretation in the FDTD case. Expected consequences to other finite methods.Int.J.Numer. Model, 13(2): 2-3. | |
| Berenger J P. 2002. Numerical reflection from FDTD-PMLs: a comparison of the split PML with the unsplit and CFS PMLs. IEEE Transactions on Antennas and Propagation, 50(3): 258-265. DOI:10.1109/8.999615 | |
| Berenger J P. 2002. Application of the CFS PML to the absorption of evanescent waves in wave guides. IEEE Microwave and Wireless Components Letters, 12(6): 218-220. DOI:10.1109/LMWC.2002.1010000 | |
| Cassidy N J, Millington T M. 2009. The application of finite-difference time-domain modelling for the assessment of GPR in magnetically lossy materials. Journal of Applied Geophysics, 67(4): 296-308. DOI:10.1016/j.jappgeo.2008.09.009 | |
| Di Q Y, Wang M Y. 1999. 2D finite element modeling for radar wave. Chinese J. Geophys., 42(6): 818-825. | |
| Dimitri K, Roland M. 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation. Geophysics, 72(5): 255-167. | |
| Drossaert F H, Giannopoulos A. 2007. A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves. Geophysics, 72(2): 9-17. | |
| Drossaert F H, Giannopoulos A. 2007. Complex frequency shifted convolution PML for FDTD modeling of elastic waves. Wave Motion, 44: 593-604. DOI:10.1016/j.wavemoti.2007.03.003 | |
| Engquist B, Majda A. 1977. Absorbing boundary conditions for the numerical simulation of waves. Math.Comput., 31(139): 629-651. DOI:10.1090/S0025-5718-1977-0436612-4 | |
| Fang H Y, Lin G. 2013. Simulation of GPR wave propagation in complicated geoelectric model using symplectic method. Chinese J. Geophys., 56(2): 653-659. | |
| Fang H Y, Lin G. 2012. Symplectic partitioned Runge-Kutta methods for two-dimensional numerical model of ground penetrating radar. Computers & Geosciences, 49: 323-329. | |
| Feng D S, Chen C S, Dai Q W. 2010. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD. Chinese J. Geophys., 53(10): 2484-2496. | |
| Feng D S, Chen C S, Wang H H. 2012. Finite element method GPR forward simulation based on mixed boundary condition. Chinese J. Geophys., 55(11): 3774-3785. | |
| Feng D S, Xie Y. 2011. Three dimensional GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD radar. Journal of Central South University: Science and Technology, 42(8): 2363-2372. | |
| Feng D S, Chen J W, Wu Q. 2014. A hybrid ADI-FDTD subgridding scheme for efficient GPR simulation of dispersion media. Chinese J. Geophys., 57(4): 1322-1334. | |
| Giannopoulos A. 2005. Modelling ground penetrating radar by GprMax. Construction and Building Materials, 19(10): 755-762. DOI:10.1016/j.conbuildmat.2005.06.007 | |
| Gedney S D. 1996. An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices. IEEE Transactions on Antennas and Propagation, 44(12): 1630-1639. DOI:10.1109/8.546249 | |
| Gedney S D. 1996. An Anisotropic PML absorbing media for the FDTD simulation of fields in lossy and dispersive media. Electromagnetics, 16(4): 399-415. DOI:10.1080/02726349608908487 | |
| Ge D B, Yan Y B. Finitedifference Time Domain Method for Electromagnetic Waves.Xi'an: Xidian University Press, 2005. | |
| Giannopoulos A. 2008. An improved new implementation of complex frequency shifted PML for the FDTD method. IEEE Transactions on Antennas and Propagation, 56(9): 2995-3000. DOI:10.1109/TAP.2008.928789 | |
| Irving J, Knight R. 2006. Numerical modeling of ground-penetrating radar in 2-D using MATLAB. Computers & Geosciences, 32(9): 1247-1258. | |
| Katz D S, Thiele E T, Taflove A. 1994. Validation and extension to three dimensions of the Berenger PML absorbing boundary condition for FD-TD meshes. IEEE Microwave and Guided Wave Letters, 4(8): 268-270. DOI:10.1109/75.311494 | |
| Kuzuoglu M, Mittra R. 1996. Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers. IEEE Microwave and Guided Wave Letters, 6(12): 447-449. DOI:10.1109/75.544545 | |
| Liao Z P, Wong H L, Yang B P, et al. 1984. A transmitting boundary for transient wave analysis. Scientia Sinica (Series A), 27(10): 1063-1076. | |
| Li J, Zeng Z F, Huang L, Liu F S. 2012. GPR simulation based on complex frequency shifted recursive integration PML boundary of 3D high order FDTD. Computers & Geosciences, 49: 121-130. | |
| Li J, Zeng Z F, Wu F S, et al. 2010. Study of three dimension high order FDTD simulation of GPR. Chinese J.Geophys.(inChinese), 53(4): 974-981. | |
| Li J X. 2007. Research on algorithms for implementing perfectly matched layers in the finite difference time domain method [Doctor's thesis]. School of Electronic Information Engineering, Tianjin university. | |
| Li Z H, Huang Q H, Wang Y B. 2009. A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media. Chinese J. Geophys., 52(7): 1915-1922. | |
| Li Z H, Huang Q H. 2014. Application of the complex frequency shifted perfectly matched layer absorbing boundary conditions in transient electromagnetic method modeling. Chinese J. Geophys., 57(4): 1292-1299. | |
| Martin R, Komatitsch D, Gedney S D, et al. 2010. A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML). CMES, 56(1): 17-40. | |
| Mei K K, Fang J. 1992. Superabsorption-a method to improve absorbing boundary conditons. IEEE Trans. on Antennas and Propagation, 40: 1001-1010. DOI:10.1109/8.166524 | |
| Mur G. 1981.Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations. IEEE Trans. Electomagn Compat, EMC-23(4): 377-382. | |
| Moerloose J D, Stuchly M. 1995. Behavior of Berenger's ABC for evanescent waves. IEEE Microwave and Guided Wave Letters, 5(10): 344-346. DOI:10.1109/75.465042 | |
| Roden J, Gedney S D. 2000. Convolution PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media. Microwave and optical technology letters, 27: 334-339. DOI:10.1002/(ISSN)1098-2760 | |
| Sacks Z S, Kingsland D M, Lee R, et al. 1995. A perfectly matched anisotropic absorber for use as an absorbing boundary condition. IEEE Transactions onAntennas and Propagation, 43(8): 1460-1463. | |
| Taflove A, Brodwin M E. 1975. Numerieal solution of steady-state EM scattering Problems using the time-dependent Maxwell's equations. IEEE Trans. MicrowaveTheory Tech, 23: 623-630. DOI:10.1109/TMTT.1975.1128640 | |
| Wang L N, Liang C H. 2006. A new implementation of CFS-PML for ADI-FDTD method. Microwave and Optical Technology Letters, 48(10): 1924-1928. DOI:10.1002/(ISSN)1098-2760 | |
| Xu T, McMechan G A. 1997. GPR attenuation and its numerical simulation in 2.5 dimensions. Geophysics, 62: 403-414. | |
| Zeng Z F, Liu S X, Feng X. Theory and Application of Ground Penetrating Radar.Beijing: Electronic Industry Press, 2010. | |
| Zhang L X, Fu LY, Pei Z L. 2010. Finite difference modeling of Biot's poroelastic equations with unsplit convolutional PML and rotated staggered grid. Chinese J. Geophys., 53(10): 2470-2483. | |
| Zhang W, Shen Y. 2010. Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling. Geophysics, 75(4): 141-154. DOI:10.1190/1.3463431 | |
| Zhang X W, Han L G, Huang L, et al. 2009. A staggered-grid high-order difference method of complex frequency-shifted PML based on recursive integration for elastic wave equation. Chinese J. Geophys., 52(7): 1800-1807. | |
| Zhan Y L, Chang Y J, Cao Z L. 2008. Modeling of groundpenetrating radar based on UPML boundary conditions. Resources Environment and Engineering, 22(2): 235-238. | |
| 底青云, 王妙月. 1999. 雷达波有限元仿真模拟. 地球物理学报, 42(6): 818–825. | |
| 方宏远, 林皋. 2013. 基于辛算法模拟探地雷达在复杂地电模型中的传播. 地球物理学报, 56(2): 653–659. | |
| 冯德山, 陈佳维, 吴奇. 2014. 混合ADI-FDTD亚网格技术在探地雷达频散媒质中的高效正演. 地球物理学报, 57(4): 1322–1334. | |
| 冯德山, 陈承申, 戴前伟. 2010. 基于UPML边界条件的交替方向隐式有限差分法GPR全波场数值模拟. 地球物理学报, 53(10): 2484–2496. | |
| 冯德山, 陈承申, 王洪华. 2012. 基于混合边界条件的有限单元法正演模拟. 地球物理学报, 55(11): 3774–3785. | |
| 冯德山, 谢源. 2011. 基于单轴各向异性完全匹配层边界条件的ADI-FDTD三维GPR全波场正演. 中南大学学报(自然科学版), 42(8): 2363–2372. | |
| 葛德彪, 闫玉波. 电磁波时域有限差分法.西安: 西安电子科技大学出版社, 2005: 67-81. | |
| 李建雄. 2007.时域有限差分法中完全匹配层的实现算法研究[博士论文].天津:天津大学电子信息工程学院. http://cdmd.cnki.com.cn/Article/CDMD-10056-2008185254.htm | |
| 李静, 曾昭发, 吴丰收, 等. 2010. 探地雷达三维高阶时域有限差分法模拟研究. 地球物理学报, 53(4): 974–981. | |
| 李展辉, 黄清华, 王彦宾. 2009. 三维错格时域伪谱法在频散介质井中雷达模拟中的应用. 地球物理学报, 52(7): 1915–1922. | |
| 李展辉, 黄清华. 2014. 复频率参数完全匹配层吸收边界在瞬变电磁法正演中的应用. 地球物理学报, 57(4): 1292–1299. | |
| 曾昭发, 刘四新, 冯晅. 探地雷达原理与应用.北京: 电子工业出版社, 2010. | |
| 张鲁新, 符力耘, 裴正林. 2010. 不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用. 地球物理学报, 53(10): 2470–2483. | |
| 张显文, 韩立国, 黄玲, 等. 2009. 基于递归积分的复频移PML弹性波方程交错网格高阶差分法. 地球物理学报, 52(7): 1800–1807. | |
| 詹应林, 昌彦君, 曹中林. 2008. 基于UPML吸收边界的探地雷达数值模拟研究. 资源环境与工程, 22(2): 235–238. | |