2010, Vol. 53
探地雷达（Ground Penetrating Radar，GPR）数值模拟一直是该领域理论研究的热点，时域有限差分法（FDTD）以算法简单灵活而被广泛采用.何兵寿等[1]运用FDTD算法结合Mur超吸收边界条件进行了理想频散条件下的雷达正演合成；刘四新等[2]利用FDTD算法结合基于分裂场的Berenger PML边界条件实现了雷达波在有损耗介质中的模拟；肖明顺[3]研究了适合GPR正演模拟FDTD算法中的单轴各向异性理想匹配层（UPML）边界条件；李静等[4]为了提高模拟精度，采用高阶FDTD算法并结合UPML边界条件进行了探地雷达数值模拟；赵延文和聂在平[5]对比了UPML媒质中ADI-FDTD方法与常规FDTD算法精度及计算效率，刘四新[6]、刘磊[7]等考虑介电常数虚部对电磁波传播的影响，用含传导电流项的Maxwell方程来进行GPR正演，得到了能更精确地刻画地下介质特征的波形.
由此可见，运用FDTD进行探地雷达数值模拟中，提高计算效率与精度，合理选取吸收效果好、计算效率高、易于编程的边界条件是今后研究的重点.考虑到ADI-FDTD方法[8~11]将传统的一个时间步分成了两个分时间步，这两个分时间步分别采用前向和后向差分，其误差互相弥补，兼具隐式差分格式的无条件稳定性和显式差分格式计算相对简单的优点，可选取较大的时间步长，极大地提高了计算效率.
本文利用无条件稳定的ADI-FDTD算法结合高精度UPML吸收边界条件进行GPR二维全波场数值模拟，大大降低了边界上的反射，提高了GPR模拟效率与精度，除GPR模拟得到的剖视图之外，还得到了Ez、Hx、Hy分量不同时刻的波场快照，通过对这些全波场快照的分析，能更加直观全面地了解雷达波的空间传播特性，提高资料解释精度.
2 ADI-FDTD法差分公式由电磁波传播理论可知[12]，GPR遵循的Maxwell方程的6个耦合偏微分方程在二维情况下变为3个偏微分方程，其电磁分量仅为Ez、Hx、Hy三个，称之为TM波：
|
(1) |
|
(2) |
|
(3) |
其中Ez为电场强度（V/m）；Hx、Hy分别为x方向和y方向磁场强度（A/m）；ε为介电常数（F/m）；μ为磁导率（H/m）；σ为电导率（S/m）；σm为导磁率（Ω/m）.下面分两个过程分别对3个偏微分方程进行差分离散.
过程一 在n→n+1/2步中，Maxwell旋度方程左边的时间偏微分项仍旧采用中心差分格式，对x方向导数的差分离散取隐式，即在子时间步的末时刻（n+1/2）Δt作差分，对y方向导数的差分离散取显式，即在初时刻nΔt作差分离散，Hx分量在n和n+1/2时刻均取值，则式（1）各分量离散为
|
并定义：
|
|
同时，定义标号（m）=（i，j+1/2），则式（1）的差分格式为
|
(4) |
同理，对x方向导数差分离散取隐式，对y方向导数的差分离散取显式，Hy分量在n和n+1/2时刻均取值，则式（2）各分量离散为
|
如果定义标号（m）=（i+1/2，j），则式（2）的差分格式为
|
(5) |
再分析（3）式，对于该等式右端两项作差分离散时，对∂Hy/∂x取隐式差分格式，对∂Hx/∂y取显式差分格式，Ez分量在n和n+1/2时刻均取值，则式（3）各分量离散为
|
如果定义标号（m）=（i，j），则式（3）的差分格式为
|
(6) |
以上式（4）~（6）可用于子时间步n→n+1/2电磁场的时域推时计算，其中式（4）在计算n+1/2时刻场值时其右端只涉及n时刻场值，称为显式格式，但（5）式和（6）式都包含同时刻Ezn+1/2和Hyn+1/2场分量，不能构成显式时间推进计算，逻辑上讲时间步进是无法完成的，称为隐式格式.为克服此困难，重写（5）式，并把它代入到（6）式消去Hyn+1/2，得
|
(7) |
如果令：
|
|
则（7）式可改写为
|
(8) |
显然（8）式可写为矩阵形式
|
(9) |
其中
|
矩阵A为三对角条带矩阵，其首、末元素b1、c1及aimax、bimax需要用吸收边界条件得到.
过程二 在n+1/2→n+1步中，对x方向导数的差分离散取显式，对y方向导数的差分离散取隐式，则式（1）、（2）、（3）离散后的差分格式如下：
|
(10) |
|
(11) |
|
(12) |
以上式（10）~（12）可用于子时间步n+1/2→n+1电磁场的时域推时计算，其中式（12）在计算n+1时刻场值时其右端只涉及n+1/2时刻场值，称为显式格式，但（10）式和（11）式都包含同时刻Ezn+1和Hxn+1场分量，不能构成显式时间推进计算，逻辑上讲时间步进是无法完成的，称为隐式格式.为克服此困难，重写（11）式，并把它代入到（10）式消去Hxn+1，得
|
(13) |
如果令：
|
|
则（13）式可改写为
|
(14) |
显然（14）式可写为矩阵形式
|
(15) |
其中
|
矩阵A为三对角条带矩阵，其首、末元素b1、c1及ajmax、bjmax需要用吸收边界条件得到.
综上所述，二维GPR波ADI-FDTD的计算步骤如下：
在子时间步n→n+1/2
（1）利用式（4）求Hxn+1/2；
（2）计算Ezn+1/2，求解矩阵方程式（8）；
（3）利用式（5）求Hyn+1/2.至此，完成子时间步n →n+1/2的推进计算.
在子时间步n+1/2→n+1
（1）利用式（12）求Hyn+1；
（2）计算Ezn+1，求解矩阵方程式（14）；
（3）利用式（11）求Hxn+1.至此，完成子时间步n +1/2→n+1的推进计算.
3 PML边界条件的选取要达到真正有效将ADI-FDTD方法应用到精细模型的GPR正演模拟中，缩短计算时间，提高计算效率，有必要为ADI-FDTD方法设置正确有效的吸收边界.当然，现有的边界条件有很多，如：Mur吸收边界条件、完全匹配层（Perfect Matched Layer，PML）边界条件，但选择与ADI-FDTD差分格式相配合的吸收边界条件的原则是[13]：吸收边界条件的引入不破坏时间步进的无条件稳定性.考虑到PML边界条件的宽频带吸收特性，受到了人们极大关注，主要有以下几种形式：基于分裂场的Berenger PML[14]、非分裂场的CPML[15]、单轴各向异性吸收层（UPML）[16]等.CPML吸收层则对低频分析有一定的优越性；BerengerPML的理论体系是非Maxwell方程的，物理机制模糊，同时，其电磁场分量分裂技术增加了数值实现的难度、计算机内存的占用.而UPML不需要对电场和磁场进行分裂，波方程仍为Maxwell方程，节约计算机存储量，提高了计算精度和计算效率，迭代公式简单，便于编程.基于它的诸多优点，本文中选用UPML边界条件结合ADI-FDTD方法进行探地雷达数值模拟.
3.1 UPML边界条件的推导由电磁场与电磁波的理论可知，在单轴各向异性PML媒质中Maxwell两个旋度方程可以表示成：
|
(16) |
其中S具有单轴各向异性媒质的特征，它可以表示成
|
二维GPR波的UPML有4个平面区和4个棱边区，仅有Ez，Hx，Hy分量，这时UPML参数sx=κx+σx/jωε0，sy=κy+σy/jωε0，sz=1，式中sx和sy是有损耗单轴各向异性介质分别在x和y方向的相对介电常数张量和相对导磁率张量；参数κx和κy是用来吸收到达UPML层的凋落波；参数σx和σy是UPML区域的衰减因子.为了让波从模拟区域传播到UPML区域，没有反射，而且波全部被吸收，参数σ和κ在UPML区域应该有恰当的空间分布，例如，一个恰当的选择：
|
其中di为PML媒质厚度.根据Gedney的经验，一般取m=4，di=（8~10）Δ.
则二维GPR波情况下，UPML边界条件下的ADI-FDTD公式汇集如下：
|
(17a) |
|
(17b) |
|
(17c) |
|
(17d) |
|
(17e) |
|
(17f) |
式中Bx=μ1Hx/sx，By=μ1Hy/sy，Dz=εsyEz.其中式（17a），（17b）为Hx，Hy→Dz→Ez的时间推进计算.式（17c），（17d）为Ez→Bx，By的时间推进计算，式（17e），（17f）为Bx，By→Hx→Hy的时间推进计算，这样就完成了1个时间步的向前推进.
下面给出（17）式的ADI-FDTD形式.为了公式形式简明，设
|
(18) |
过程一 在子时间步n→n+1/2中，对x方向导数的差分离散取隐式，即在末时刻（n+1/2）Δt作差分离散；对y方向导数的差分离散取显式，即在初时刻nΔt作差分离散.于是，（17a）式的ADI-FDTD离散式为
|
(19) |
（17b）式的ADI-FDTD离散式为
|
(20a) |
即
|
(20b) |
（17c）式的ADI-FDTD离散式为
|
(21) |
（17d）式的ADI-FDTD离散式为
|
(22) |
（17e）式的ADI-FDTD离散式为
|
(23) |
（17f）式的ADI-FDTD离散式为
|
(24) |
注意（19）~（24）式中除（21）式以外，等式左右端都包含有n+1/2同时刻量.例如在（20）式中，Ezn+1/2和Dzn+1/2均为同时刻量，因此不能构成显式的时间推进计算.由（19）~（24）式消去Ezn+1/2以外的同时刻量，得到
|
(25) |
其中
|
(26) |
过程二 在子时间步n+1/2→n+1中，对y方向导数的差分离散取隐式，即在末时刻（n+1）Δt作差分离散；对x方向导数的差分离散取显式，即在初时刻（n+1/2）Δt作差分离散.于是，（17a）式的离散式为
|
(27) |
（17b）式的离散式为
|
(28a) |
即
|
(28b) |
（17c）式的离散式为
|
(29) |
（17d）式的离散式为
|
(30) |
（17e）式的离散式为
|
(31) |
（17f）式的离散式为
|
(32) |
为了得到n+1/2→n+1的时域推进计算，和n→n+1/2子时间步一样，由（27）~（32）式消去同时刻量，可得
|
(33) |
其中
|
(34) |
综上所述，二维GPR波UPML的ADI-FDTD的计算步骤如下：
在子时间步n→n+1/2
（1）计算Ezn+1/2，并求解矩阵方程式（25）；
（2）利用式（20b）求Dzn+1/2；
（3）利用式（21）、（22）求Bxn+1/2和Byn+1/2；
（4）利用式（23）、（24）求Hxn+1/2和Hyn+1/2.至此，
完成子时间步n→n+1/2的推进计算.
在子时间步n→n+1/2
（1）计算Ezn+1，并求解矩阵方程式（33）；
（2）利用式（28b）求Dzn+1；
（3）利用式（29）、（30）求Bxn+1和Byn+1；
（4）利用式（31）、（32）求Hxn+1和Hyn+1.至此，完成子时间步n→n+1/2的推进计算.
3.2 UPML边界条件应用实例为验证UPML边界条件的吸收效果，在二维空间中加入Ricker子波，对比它在有、无UPML边界条件下雷达波传播的空间特征.模拟区域为1.0 m ×1.0 m的正方形，空间为均匀素填土介质，其相对介电常数为10.0，电导率为0.002S/m，在模拟区域的正中心（x=0.5 m，y=0.5 m）位置加入一个Ricker子波.图 1为未加入边界条件1.0×10-8 s时的全波场快照，分析图 1中的（a）~（c）三个电磁场分量波场快照，无边界条件下波场分量在截断边界的四条边上产生了强的反射波，而四个角点附近尤其强烈；磁场分量Hx与Hy图形基本一致，把一个分量旋转90°，能很好地吻合，这与其实际物理意义是一致的.由图 2（a~c）波场快照可知，雷达波等值线是非常规则的同心圆，四个边角没有明显的反射，表明UPML层的吸收效果很好，雷达波无反射地进入完全匹配层后被迅速衰减掉.对比图 1，2表明：UPML是比较适合ADI-FDTD算法的边界条件，吸收效果良好，同时两图中也说明了均匀介质中雷达波前以球面波辐射传播的特点.
|
图 1 无边界条件的全波场快照 Fig. 1 The snapshot of full wave field without boundary condition |
|
图 2 加入UPML边界条件的全波场快照 Fig. 2 The snapshot of full wave field adding the UPML boundary condition |
目前探地雷达多采用在地表发射脉冲子波，同时由接收天线在地面接收的工作方式，为更真实地模拟雷达波在地下的传播特征，只考虑地下半空间情况.设定模拟区域为3.0 m×2.0 m，Ricker子波主频为900 MHz，位于模型深度0.1 m位置.模型中最上层为空气介质，厚度为0.1 m，空气层下面为素填土，其介电常数为10.0，电导率为0.002S/m，在素填土中存在一个混凝土矩状异常体，其两个对角点坐标为（1.425 m，0.6 m）、（1.575 m，0.9 m），其介电常数为6.0，电导率为0.0005S/m，详见图 3模型示意图.取网格空间步长为dx=0.005 m，dy=0.005 m，UPML边界设为10层网格，总的网格数为600×400，接收雷达波场快照的点位为（x=1.5 m，y=0.0 m）.受CFL条件的影响，应用常规的FDTD正演软件如（GPRmax，Author：Antonis Giannopoulos）取时间步长为1.17933×10-11 s，在Intel（R）CoreTM 2 Duo CPU，P8600@2.4GHz，2.96GB的内存物理地址扩展，Window XP操作系统的IBM X200笔记本上计算该雷达剖面50道数据，共计花费时间16 min46 s，而本文应用结合UPML边界条件ADI-FDTD算法采用的时间步长为2.0×10-11 s，采用更大时间步长的同时保证了算法的稳定性，尽管还需求解大型线性方程组，但本文方法把计算时间缩短了4 min，效率提高了近1/4（可能会因不同的编程风格存在出入）.分析wiggle图 4与扫描图 5可知，矩状混凝土的上界面在雷达正演图中仍然为平界面，矩状金属体的两个棱角位置出现绕射波，由于矩形的顶边不是很长，导致矩状体的上边在雷达剖面图中近似为双曲线形弧形，其实质是一个平界面与两个双曲线的半支组成，而矩状体的下界面与上界面类似，只是反射波能量较弱.分析图 6（a~i）三个电磁场波场快照可知，Ez、Hx、Hy分量的雷达波前以近似半圆形向下扩散传播，不同时刻雷达波前到达位置得到清楚的反映，当波前到达矩状异常体位置时，异常体的4条边产生了强的反射.通过分析这些不同时刻的雷达全波场快照，可以了解雷达波在地下半空间介质的传播详细过程，以便更好地认识雷达波的传播特性，进一步对工程资料作出精确解释.
|
图 3 模型一示意图 Fig. 3 Sketch map of model 1 |
|
图 4 模型一正演wiggle图 Fig. 4 Simulate wiggle map of model 1 |
|
图 5 模型一正演扫描图 Fig. 5 Simulate scan map of model 1 |
|
图 6 模型一的雷达波传播全波场快照 Fig. 6 GPR full wave field snapshot of model 1 |
模型二的模拟区域参数、介质参数、天线主频与模型一完全相同.混凝土与素填土的分界面存在两个凹的“V”形与一个三角形突起，在混凝土介质中存在两个球状异常体，半径均为0.02 m，其中位于（0.30 m，0.4 m）位置为金属球体，位于（2.50 m，0.4 m）位置为球状空洞，详见示意图 7.应用本文ADI-FDTD算法并取上例的模拟参数进行计算，模拟所得的正演wiggle图与扫描图如图 8、9所示，不同时刻全波场快照如图 10所示.分析wiggle图与扫描图可知，两个球状体都产生了强的反射波与绕射波，但两者产生的绕射波是有区别的，金属球体的绕射波能量更强，波幅更宽.在两个凹的“V”形与一个三角形突起位置也产生了强的绕射波，由于两个凹的“V”形角度大小的区别，在wiggle图与扫描图上体现为绕射波形态不一致，右边角度较平缓的凹“V”形下顶点位置产生了明显的回转波.
|
图 7 模型二示意图 Fig. 7 Sketch map of model 2 |
|
图 8 模型二正演wiggle图 Fig. 8 Simulate wiggle map of model 2 |
|
图 9 模型二正演扫描图 Fig. 9 GPR full wave field snapshot of model 2 |
|
图 10 模型二的雷达波传播全波场快照 Fig. 10 GPR full wave field snapshot of model 2 |
分析5.0×10-9 s时Ez波场快照（图 10a），雷达波前以半圆形向下传播，达到深度0.5 m位置；8.0×10-9 s时刻（图 10b）波前到达凸出的三角形附近，并在其上顶点产生绕射波，使雷达波一部分以该凸出的三角形顶点为新的波源向上传播，其波前主体部分继续向下传播；在1.0×10-8 s时刻（图 10c），原本以半圆形向下传播的波前，受两个凹的“V”形的影响，变成了带钜齿状半圆向下传播，而向上传播的波前呈三角形的形态传至顶边界处；在1.2×10-8 s时刻（图 10d），向下传播的波前已至深度1.2 m，受到两个凹的“V”形影响，图中产生了两个直立的菱形，而向上传播的波已超出边界.图 10（e~l）为磁场Hx、Hy的雷达波场快照，其传播特征与Ez基本类似，仅在磁场快照左上角区域中能更清楚地分辨出金属球体的位置.由于模型二相对复杂，各种异常体产生的绕射波、反射波、多次波混叠在一起，导致后期的波场快照相对复杂，波场细节特征难以分辨.通过分析雷达wiggle图、扫描图与不同时刻的雷达全波场快照，能更全面准确地推断模型特征，深刻认识雷达波的传播规律，指导GPR资料解释.
5 结论及建议（1）本文结合交替方向隐式差分的特点，推导了二维雷达（TM）波ADI-FDTD方法的迭代计算公式，包括第一个子时间步n→n+1/2对方程中的x方向导数取隐式差分格式，y方向导数取显式；第二个子时间步n+1/2→n+1的迭代过程中对x方向导数取显式差分格式，y方向导数取隐式，并分别给出了其计算步骤.
（2）成功地将UPML吸收边界引入到GPR数值模拟的ADI-FDTD方法中，而不破坏ADI-FDTD方法的无条件稳定性，消除了截断边界处的强反射，为探地雷达的正演模拟提供了一种新的模拟手段. ADI-FDTD方法是一种无条件稳定的交替方向隐式时域有限差分法，它克服了常规FDTD法受CFL约束条件限制，可选取较大的时间步长，尽管需要迭代求解大型线性方程组，仍然提高了计算效率.
（3）设计了1个复杂模型和1个简单模型，应用基于UPML边界条件的ADI-FDTD算法并编制了相应的计算程序，对这两个模型进行数值模拟，得到了正演模拟wiggle图、扫描图及雷达全波场快照.通过分析这些雷达剖面图与全波场快照，对于雷达波在空间的传播特征有了系统的了解和更清楚的认识，有助于雷达资料的解释.
| [1] | 何兵寿, 张会星. 地质雷达正演中的频散压制和吸收边界条件改进方法. 地质与勘探 , 2000, 36(3): 59–63. He B S, Zhang H X. The suppression of numerical dispersion and improvement of absorbing boundary conditions in forward modeling of GPR. Geology and Prospecting (in Chinese) , 2000, 36(3): 59-63. |
| [2] | 刘四新, 曾昭发, 徐波. 三维频散介质中地质雷达信号的FDTD数值模拟. 吉林大学学报(地球科学版) , 2006, 36(1): 123–127. Liu S X, Zeng Z F, Xu B. FDTD simulation for Ground Penetrating Radar signal in 3-Dimensional dispersive medium. Journal of Jilin University (Earth Science Edition) (in Chinese) , 2006, 36(1): 123-127. |
| [3] | 肖明顺, 昌彦君, 曹中林, 等. 探地雷达数值模拟的吸收边界条件研究. 工程地球物理学报(地球科学版) , 2008, 5(3): 315–320. Xiao M S, Chang Y J, Cao Z L, et al. Research on boundary conditions in forward modeling of Ground Penetrating Radar. Chinese Journal of Engineering Geophysics (in Chinese) , 2008, 5(3): 315-320. |
| [4] | 李静, 曾昭发, 吴丰收, 等. 探地雷达三维高阶时域有限差分模拟研究. 地球物理学报 , 2010, 53(4): 974–981. Li J, Zeng Z F, Wu F S, et al. Study of three dimension high-order FDTD simulation for GPR. Chinese J. Geophys (in Chinese) , 2010, 53(4): 974-981. |
| [5] | 赵延文, 聂在平. UPML媒质中无条件稳定的二维ADI-FDTD方法. 电波科学学报 , 2002, 17(6): 586–589, 603. Zhao Y W, Nie Z P. Unconditionally stable 2D ADI-FDTD method in perfectly matched uniaxial medium. Chinese Journal of Radio Science (in Chinese) , 2002, 17(6): 586-589, 603. |
| [6] | 刘四新, 曾昭发. 频散介质中地质雷达波传播的数值模拟. 地球物理学报 , 2007, 50(1): 320–326. Liu S X, Zeng Z F. Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium. Chinese J. Geophys (in Chinese) , 2007, 50(1): 320-326. |
| [7] | 刘磊, 昌彦君, 曹中林. 色散介质的探地雷达正演模拟. 资源环境与工程 , 2009, 23(2): 171–174. Liu L, Chang Y J, Cao Z L. Forward simulation of Ground Penetrating Radar in dispersive medium. Resources Environment & Engineering (in Chinese) , 2009, 23(2): 171-174. |
| [8] | 郑奎松. FDTD网络并行计算及ADI-FDTD方法研究.西安:西安电子科技大学理学院, 2005 Zheng K S. Study of FDTD parallel computing and ADI-FDTD method (in Chinese). Xi'an:Institute of Science, XiDian University, 2005. http://cdmd.cnki.com.cn/Article/CDMD-10701-2006055896.htm |
| [9] | Fu W M, Tan E L. Stability and dispersion analysis for higher order 3-D ADI-FDTD Method. IEEE Transactions on Antennas and Propagation , 2005, 53(11): 3691-3696. DOI:10.1109/TAP.2005.858588 |
| [10] | Namiki T. 3-D ADI-FDTD method-unconditionally stable time-domain algorithm for solving full vector Maxwell's equations. IEEE Trans. Microwave Theory and Techniques , 2000, 48(10): 1743-1748. DOI:10.1109/22.873904 |
| [11] | Namiki T. A new FDTD algorithm based on alternating-direction implicit method. IEEE Trans. Microwave Theory and Techniques , 1999, 47(10): 2003-2007. DOI:10.1109/22.795075 |
| [12] | 葛德彪, 闫玉波. 电磁波时域有限差分方法. (第二版). 西安: 西安电子科技大学出版社, 2005 . Ge D B, Yan Y B. Finite-Difference Time-Domain Method for Electromagnetic Waves (in Chinese). (2th Edition). Xi'an: XiDian University Press, 2005 . |
| [13] | Gedney S D, Liu G, Roden J A. Perfectly matched layer media with CFS for an unconditionally stable ADI-FDTD method. IEEE Transactions on Antennas and Propagation , 2001, 49(11): 1554-1559. DOI:10.1109/8.964091 |
| [14] | Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics , 1994, 114(2): 185-200. DOI:10.1006/jcph.1994.1159 |
| [15] | Roden J A, Gedney S D. Convolutional PML (CPML):an efficient FDTD implementation of the CFS-PML for arbitrary media. Microwave and Optical Technology Letters , 2000, 27(5): 334-339. DOI:10.1002/(ISSN)1098-2760 |
| [16] | Liu G, Gedney S D. Perfectly matched layer medium for an unconditionally stable three-dimensional ADI-FDTD method. IEEE Microwave and Guided Wave Letters , 2000, 10(7): 261-263. DOI:10.1109/75.856982 |