基于旋转交错网格的探地雷达正演数值模拟-有色金属在线

以旋转交错网格(RSG)差分为基础，实现探地雷达(GPR)在非均匀突变介质中的正演数值模拟。通过利用旋转交错差分算子，将雷达波场各场分量设置在非均匀突变介质的基本单元中，单元内同一物性参数定义在同一位置，将对角线方向上物性参数差分值的线性组合来计算坐标轴方向上的物性参数差分值，针对非均匀特征变化剧烈的网格区域，算法中的单元网格形式不再需要进行插值，扩宽了数值稳定性条件的限制。推导了适合于探地雷达的旋转交错网格差分格式及TM_y极化模式下的更新方程组，实现了均匀介质和非均匀突变介质中的TM_y雷达波场的数值模拟，并分别从解析解、模拟剖面和波场快照的角度，与标准交错网格差分算法进行了对比。结果表明，在模拟非均匀特征明显的介质时，旋转交错网格算法可以选取更大的时间步长，数值色散程度控制更好，提高了模拟效率和成像精度，可更有效地指导非均匀突变介质中探地雷达数据的解译。

关键词：

旋转交错网格；标准交错网格；探地雷达；正演模拟；数值色散；

中图分类号：P631　　文献标志码：A

GPR simulation based on rotated staggered grid

ZHANG Bin^{1, 2, 3}, DAI Qian-wei^{1, 2}, YIN Xiao-bo^{1, 3}

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,

Central South University, Changsha 410083, China;

3. Institute of Resources exploration and Environmental Geology, Central South University, Changsha 410083, China)

Abstract: Based on the rotated staggered grid finite-difference, forward simulation in inhomogeneous media with strong discontinuities for ground penetrating radar (GPR) was implemented. By using the rotating staggered difference operator, the GPR wave field components and other physical parameters were distributed in the elementary cells of RSG, in which all field components of one physical property were located at one elementary unit in computational domain, then the difference of field components and physical parameters along the coordinate axes were calculated by using the linear combination value of them across the diagonal coordinate axes, no averaging of elementary cell was needed even in grids-domain with strong heterogeneities, which relax the limitation of the numerical stability condition. The rotated staggered grid finite difference scheme for GPR and the corresponding update equations was deduced in TM_y polarization mode, then the numerical simulation of GPR wave field in TM_y polarization mode was implemented, in addition, the comparation between the standard staggered grid and RSG was presented from three respects of the analytic solutions, the simulated sections and the field snapshots. The results show that with the more relaxed limitation of the numerical stability condition and the better controlling of numerical dispersion degree, the RSG difference algorithm improves the efficiency and accuracy of simulation, which can effectively guide the GPR data inversion and interpretation in the inhomogeneous media with strong heterogeneities.

Key words: rotated staggered grid; standard staggered grid; ground penetrating radar; forward modeling; numerical dispersion

探地雷达正演数值模拟是研究雷达波波场传播规律的有效途径，诸多数值计算方法已成功应用于探地雷达的数值模拟中：射线追踪法正演成像^[1-2]、频率域法正演模拟^[3-4]、积分法数值模拟^[5]、伪谱法^[6-9]、时域有限差分法^[10-19]等。高效的有限差分方法已经发展成为了一种较成熟的数值模拟方法，然而在应对物性参数差异较大或非均匀特征明显的异常目标体时，无论是隐式还是显式的标准交错网格差分格式，均会造成不同程度的数值不稳定性^[20-21]。为解决该问题，国内外学者从不同方面着重研究非均匀介质中的正演模拟。FELLINGER等^[22]采用弹性有限积分法模拟了弹性波的传播和散射，改善了各向异性非均匀介质的模拟效果；底青云等^[23]则采用有限单元法，推导了含衰减项的有限元方程，实现了复杂介质的探地雷达有限元数值模拟；陈承申等^[24]提出了基于混合边界条件的GPR有限单元正演模拟；PITARKA等^[25]提出了不规则差分网格的思想，在粗细网格的过渡区采用插值算法来处理各项异性特征明显的网格区域，相应地提高了模拟精度也降低了计算机内存；张慧等^[26]在空间变网格差分法的基础上，提出了空间和时间的双变网格差分算法，提高了复杂介质数值模拟的精度；冯德山等^[27]结合单轴各向异性完全匹配层边界条件和交替方向隐式差分法，突破了探地雷达稳定性条件的限制；朱自强等^[28]采用卷积完全匹配层边界条件实现了钻孔雷达的全波场数值模拟；戴前伟等^[29]采用三角形剖分实现了对复杂探地雷达模型的有限元模拟，提高了模拟精度；冯德山、王洪华等^[30]则采用了只需节点而无需单元的无网格Galerkin算法，得到了高次连续的解；方宏远^[31]等则推导了一阶显式辛分块龙格库塔的辛算法，在与传统差分网格算法精度同等的情况下节约了迭代计算的时间。

旋转交错网格(Rotated staggered grid，RSG)差分法于2000年首次由波动反演技术项目组(Wave inversion technology consortium, WIT)的SAENGER 等^[20]提出，该差分格式在标准交错差分网格(Standard Staggered Grid，SSG)的基础上进行了改进。该算法将所有的介质参数定义在同一单元的合适位置，单元内同一物理量定义在同一位置，通过对角线方向物理量差分值的线性组合来计算坐标轴方向上的物理量差分，彻底解决了标准交错差分网格或传统同位网格在模拟各向异性介质时需插值的缺点。目前，未见相关文献和研究成果将该差分网格应用于探地雷达电磁波数值模拟领域，而在地震波和声波数值模拟领域有一些应用，相应的成果主要有：SAENGER等^[32]介绍了当模型中存在单一断裂结构异常体时所产生的弹性波波场散射和绕射情况，并分析了RSG差分法的数值稳定性和精确性；SAENGER等^[33]采用RSG差分法模拟了地震波在黏弹性介质和各向异性介质中的传播；陈浩等^[34]系统地介绍了RSG的差分格式，并给出了相应的完全匹配层吸收边界条件及实现步骤；张鲁新等^[35]则将不分裂卷积完全匹配层的RSG差分法应用于孔隙弹性介质的数值模拟中，解决了RSG中普通PML边界条件吸收不足的缺点；严红勇等^[36]推导了RSG任意偶数阶精度有限差分格式，实现了各向异性介质和黏弹性介质的地震波数值模拟，较好地反映地下介质的真实情况。

考虑到复杂地质体具有非均匀的广泛性，将弹性波和声波领域的旋转交错网格差分法应用于探地雷达正演模拟中。该方法继承了标准交错网格的高效迭代特点，利用旋转交错差分算子，将各雷达波场分量及物质参量(介电常数、磁导率、电导率)设置在非均匀介质的基本单元中，利用对角线方向上差分值的线性组合来计算非均匀特征变化剧烈的网格区域坐标轴上的差分值，扩宽了标准交错网格中数值稳定性条件的限制，解决了标准交错网格在选取极限时间步长后引起的数值不稳定或不精确等缺点，可提高探地雷达在非均匀突变介质中的模拟效率和成像精度。

1 旋转交错网格差分格式

与标准交错网格或传统同位网格差分算法相比，旋转交错网格差分法采用不同的空间交错策略，先计算单元格对角线方向上的相关物理量的差分，然后将这些对角线方向上的差分值的线性组合来计算沿坐标轴方向上的差分，即是将相同类型的物理量定义在同一位置^[20]。

在二维(x，z)介质情况下，旋转交错网格有限差分格式的算子可以写成如下形式：

(1)

(2)

(3)

(4)

式中：z′和x′分别代表对角线方向上的向量；z和x分别代表坐标轴方向上的向量；Δz和Δx则是沿坐标轴方向上的空间步长；Δr是对角线的长度，且有，和为离散差分算子，其不同的空间交错策略及单元网格场量分布如图1和图2所示。

图1 标准交错网格单元

Fig. 1 Schematic diagram of elementary cells of standard staggered grid

图2 旋转交错网格单元

Fig. 2 Elementary cells of rotated staggered grid

2 TM_y模式下RSG差分法的离散

探地雷达正演模拟是在特定的介质模型区域中求解Maxwell方程组的数值计算过程，在二维情况下只含E_y、H_x、H_z这3个分量的TM_y极化模式独立的方程组可以写为

(5)

(6)

(7)

式中：μ为磁导率，H/m；σ为电导率，S/m；上标m和e分别表示磁场和电场作用；ε为相对介电常数；E为电场强度，V/m；H为磁场强度，A/m。在标准交错网格的中心差分中，以t时刻(x，z)坐标中的电场分量E(x，z，t)为例(磁场分量类似)：

(8)

(9)

在旋转交错网格中，则有

(10)

(11)

结合式(3)、(4)、(10)和(11)，在旋转交错网格中式(8)和(9)可近似地写为

(12)

(13)

则式(5)、(6)、(7)可近似地写为

(14)

(15)

(16)

旋转交错网格中FDTD各场分量的位置如图3所示，其他物性参数(相对介电常数ε、磁导率μ、电导率σ)的标记方式与场分量相同，该图为场量网格单元向y(y′)轴方向投影至x′z′平面得到。

在TM_y极化模式下，从图3所示的场分量分布，以标准网格坐标标号来表示旋转网格坐标，可得到旋转交错网格中的FDTD更新方程：

(17)

(18)

(19)

图3 TM_y极化模式的旋转交错网格场分量

Fig. 3 Field component of RSG in TM_y polarization mode

取，，整理得：

(20)

(21)

(22)

在程序中令所有的系数项分别为

(23)

(24)

(25)

(26)

(27)

(28)

(29)

显然，RSG有限差分算法中系数项总数与SSG算法的相同，只是各更新方程所对应的系数项不同。

3 RSG算法的实现

从RSG差分算法的更新方程中可以看出，各场分量分别取前一时刻场值的相应半步长进行迭代，在更新方程(20)~(21)中，磁场分量H_x和H_z均在(n+0.5)Δt时刻获取更新，而电场分量E_y则利用更新式(22)在(n+1)Δt时刻获取更新。在各个边界处理上，本文作者采用普通的PML吸收条件，其实现方式与SSG中相同，最后将所有更新的当前场量进行存储，RSG差分算法的流程如图4所示。

图4 RSG模拟程序流程示意图

Fig. 4 Flow chart of RSG for simulation

4 数值算例与分析

4.1 RSG算法正确性和稳定性验证

为了验证RSG算法的正确性，建立一个有解析解的含薄层均匀介质模型。如图5所示，模型为2.0 m× 2.0 m的均匀混凝土介质，相对介电常数为9.0，电导率为0.0001 S/m，薄层介质为湿砂，厚度为0.20 m，相对介电常数为20.0，电导率为0.001 S/m。为更细致地分析RSG算法的稳定性，选取两种不同的时间步长，分别采用SSG和RSG算法进行模拟，两种条件下均选取空间步长dx=dz=0.005 m，即分别选取较大的时间步长和接近极限值的时间步长：dt₁=0.015 ns，dt₂=0.03 ns，模拟时天线的中心频率为800 MHz，边界均采用常规PML吸收边界条件，各占据8个网格。

图5 含薄层的均匀介质模型

Fig. 5 Homogeneous medium model with thin layer embedded

图6所示为时间步长取dt₁=0.015 ns的单道波计算结果，实细线为解析解，实粗线为RSG模拟解，虚线为SSG模拟解。与解析解相比，模拟解对薄层的上下界面响应一致，初步证明了RSG算法的正确性，由于选取了较大的时间步长，在25~30 ns时间段，SSG模拟结果开始产生不稳定现象，见图6中的放大图、图7和图8中箭头所示，说明在较大时间步长条件下，RSG模拟结果与解析解拟合得更好。

为考察取更大时间步长时的稳定性情况，针对该均匀介质模型，选取了接近SSG稳定性条件限制的极限时间步长dt₂=0.03 ns，分别采用SSG和RSG算法进行模拟，结果如图9所示。两种方法的模拟结果均产生了不同程度的数值色散现象，但是RSG计算的结果仍是以分辨薄层的上下界面，而SSG计算结果的色散现象更严重，已无法分辨薄层的界面。由可进行粗略估算(Δi为空间步长)，SSG算法在模拟均匀介质的时间步长限制条件约为dt≤0.03536 ns，而RSG算法约为dt≤0.05ns。显然，在该时间步长条件下，RSG计算的过程更稳定。

4.2 非均匀突变介质的模型

图6 SSG、RSG数值模拟解与解析解的比较

Fig. 6 Comparison of SSG, RSG modeling solution with its analytic solution

为验证算法对非均匀突变点模拟的有效性，设置了如图10所示的复杂模型，并分别采用SSG和RSG算法进行计算，模型为2.0 m×2.0 m的方形区域，区域中存在诸多非均匀的多边形异常突变介质，突变介质为湿砂，其他参数与上述均匀介质模型的一致，同样选取接近SSG稳定性条件限制的极限时间步长dt₂=0.03 ns，分别采用SSG和RSG算法进行模拟。

图7 SSG模拟的雷达剖面图

Fig. 7 GPR simulation scans with SSG method

图8 RSG模拟的雷达剖面图

Fig. 8 GPR simulation scans with RSG method

图9 时间步长为0.03 ns时分别采用RSG(a)和SSG(b)算法模拟计算的单道波信号

Fig. 9 Single scanning signal simulated by RSG (a) and SSG (b) with time-step of 0.03 ns

图11所示为SSG算法模拟的雷达剖面图，受突变介质相对介电常数的主要影响。由于在雷达波传播至突变介质内部时，0.03 ns的极限时间步长已超出了SSG稳定性条件的限制，在5个绕射反射波簇的尾部均产生了大量的杂波和多次波，严重影响了突变介质异常体的数目和分界面的分辨，数值色散现象严重，且该现象的严重程度随传播时间的增加而增强。图12所示为RSG算法模拟的雷达剖面，在5个绕射反射板簇的尾部也产生了诸多的杂波和多次波，使得突变介质异常体的分界面难以分辨，但并不影响其数目的分辨，模拟结果存在一定程度的数值色散现象，但该现象并未随传播时间而叠加，深部的信号仍能分辨突变介质异常体的绕射数目，数值色散现象得到明显改善。

图10 非均匀突变介质地电模型

Fig. 10 Geoelectric model with strong heterogeneities

图11 SSG模拟的雷达剖面图

Fig. 11 GPR simulation scans with SSG method

图12 RSG模拟的雷达剖面图

Fig. 12 GPR simulation scans with RSG method

4.3 波场快照分析

为更细致地比较SSG和RSG的数值色散程度和成像精度，分别给出了两种算法的波场快照结果，计算的结果为源位于x=1.0 m处，E_y分量在4个不同传播时刻的瞬时波场快照，图示能量已均一化，如图13所示。其中图13(a)、(b)、(c)和(d)所示为SSG在5 ns、10 ns、15 ns和20 ns的波场快照，图13(e)、(f)、(g)和(h)所示为RSG在5 ns、10 ns、15 ns和20 ns的波场快照。在t=5 ns时刻，波前面还未传播至最浅部的第1个突变点时，两种算法均满足对应的稳定性条件，数值色散程度一致，色散程度较小，如图13(a)和(e)所示；在t=10 ns时刻，随着波前面的扩散，当波前面传播至第3个突变介质内部时，SSG模拟的结果叠加了前一时刻的计算误差，干扰了该时刻3个绕射波的分辨，而在RSG模拟结果中仍能清晰地分辨3个绕射波，如图13(f)所示；随着波前面继续往深部扩散，计算误差的叠加，SSG算法的色散程度愈加严重，如图13(c)和(d)所示。由于雷达波在突变介质异常体的内部已并不满足数值稳定性条件，使得计算误差和色散程度持续叠加，在主计算区域内产生了大量的杂波和多次波，严重影响突变介质异常体的数目和界面的分辨，而该色散现象在RSG算法中得到有效改善，主计算区域内也产生了诸多杂波和多次波，但该现象并未愈加严重，数值色散程度较小，计算误差仍能分辨突变介质异常体的数目，如图13(g)和(h)所示，RSG算法较好地控制了数值色散现象和计算误差，提高了正演数值模拟的成像精度。

5 结论

1) 将RSG应用于探地雷达正演数值模拟中，介绍了RSG的差分格式，TM_y激化模式下的RSG差分离散，详细推导了探地雷达的RSG差分格式及TM_y极化模式下的更新方程组，并给出了探地雷达RSG算法正演模拟的实现过程。

2) 分析了SSG算法和RSG算法的不同交错策略和稳定性条件，在相同的空间步长前提下，RSG可以选取更大的时间步长，来获取与SSG相同的成像精度，拓宽了数值稳定性条件的限制，提高了计算效率。

图13 SSG和RSG算法不同时刻E_y分量的波场快照(空间步长dx=dz=0.005 m，时间步长dt=0.03 ns)

Fig. 13 Snapshots of E_y component at different time-steps with SSG and RSG scheme (with space interval dx=dz=0.005 m, time-step is dt=0.03 ns)

3) 模拟了探地雷达在均匀介质中的传播，比较了SSG和RSG的数值模拟解与解析解的拟合程度，验证了RSG算法的正确性，结果表明在较大的时间步长情况下，RSG算法的数值色散程度更低。

4) 模拟了探地雷达在非均匀突变介质模型中的传播，比较了SSG和RSG算法的成像精度。在接近SSG稳定性条件限制的极限时间步长条件下，SSG计算的结果数值色散现象严重，而RSG算法仍具备较高的分辨率，提高了成像精度，并且RSG的计算过程更稳定。

REFERENCES

[1] GOODMAN D. Ground-penetrating radar simulation in engineering and archeology[J]. Geophysics, 1994, 59: 224-232.

[2] CAI J. Ray-based synthesis of bistatic ground-penetrating radar profiles[J]. Geophysics, 1995, 60: 87-96.

[3] POWERS M H, OLHOEFT G R. Modeling dispersive ground penetrating radar data[C]// Proceedings of the 5th International Conference on Ground-Penetrating Radar, Waterloo, Ontario, 1994: 173-183.

[4] ZENG X, MCMECHAN, G A, CAI J, CHEN H W. Comparison of ray and Fourier methods for modeling monostatic ground-penetrating radar profiles[J]. Geophysics, 1995, 60: 1727-1734.

[5] ELLEFSEN K J. Effects of layered sediments on the guided wave in crosswell radar data[J]. Geophysics, 1999, 64: 35-38, 1698-1707.

[6] CARCIONE J M. Ground-penetrating radar: Wave theory and numerical simulation in lossy anisotropic media[J]. Geophysics, 1996, 61: 1664-1677.

[7] CASPER D A, Kung K S. Simulation of ground-penetrating radar waves in a 2-D soil model[J]. Geophysics, 1996, 61: 1034-1049.

[8] LUI Q H, FAN G. Simulations of GPR in dispersive media using a frequency-dependent PSTD algorithm[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37: 2317-2324.

[9] 李展辉, 黄清华, 王彦宾. 三维错格时域伪谱法在频散介质井中雷达模拟中的应用[J]. 地球物理学报, 2009, 52(7): 1915-1922.

LI Zhan-hui, HUANG Qing-hua, WANG Yan-bin. A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media[J]. Chinese Journal of Geophysics, 2009, 52(7): 1915-1922.

[10] YEE K S. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media[J]. IEEE Trans Antennas Propagate, 1966, 14: 302-307.

[11] UMRAN S I, ROBERT A. MARSHALL. Numerical Electromagnetics[M]. Cambridge: Cambridge University Press, 2011.

[12] ROBERTS R L, DANIELS J. Finite-difference time domain forward modeling of GPR data[C]//Proceedings of the 5th International Conference on Ground-Penetrating Radar, Waterloo, Ontario, 1994: 185-204.

[13] TEIXEIRA F L, CHEW W C. Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, conductive soils[J]. IEEE Transactions on Geoscience and Remote Sensing, 1998, 36: 1928-1936.

[14] HOLLIGER K, BERGMAN T. Numerical modeling of borehole georadar data[J]. Geophysics, 2002, 67: 1249-1257.

[15] GEORGAKOPOULOS S V, BIRTCHER C R. Higher-order finite-difference schemes for electromagnetic radiation scattering penetration[J]. IEEE Antennas and Propagation Magazine, 2002, 44: 134-142.

[16] 刘四新, 曾昭发, 徐波. 三维频散介质中地质雷达信号的FDTD数值模拟[J]. 吉林大学学报(地球科学版), 2006, 36(1): 123-127.

LIU Si-xin, ZENG Zhao-fa, XU Bo. FDTD simulation for Ground Penetrating Radar signal in 3-Dimensional dispersive medium[J]. Journal of Jilin University (Earth Science Edition), 2006, 36(1): 123-127.

[17] 李静, 曾昭发, 吴丰收, 黄玲. 探地雷达三维高阶时域有限差分模拟研究[J]. 地球物理学报, 2010, 53(4): 974-981.

LI Jing, ZENG Zhao-fa, WU Feng-shou, HUANG Ling. Study of three dimension high-order FDTD simulation for GPR[J]. Chinese Journal of Geophysics, 2010, 53(4): 974-981.

[18] 冯德山, 张彬, 戴前伟, 陈德鹏. 基于速度估计的改进型线性变换有限差分偏移在探地雷达中的应用[J]. 地球物理学报, 2011, 54(5): 1340-1347.

FENG De-shan, ZHANG Bin, DAI Qian-wei, CHEN De-peng. The application of the improved linear transformation of finite difference migration based on the velocith estimation in the GPR data processing[J]. Chinese Journal of Geophysics, 2011, 54(5): 1340-1347.

[19] FENG D S, DAI Q W. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD[J]. NDT&E International, 2011, 44(6): 495-504.

[20] SAENGER E H, GOLD N, SHAPIRO S A. Modeling the propagation of elastic waves using a modified finite-difference grid[J]. Wave Motion, 2000, 31: 77-92.

[21] SAENGER E H, SHAPIRO S A. Effective velocities in fractured media: A numerical study using the rotated staggered finite difference grid[J]. Geophysical Prospecting, 2002, 50: 183-194.

[22] FELLINGER P, MARKLEIN R, LANGENBERG K J, KLAHOLZ S. Numerical modeling of elastic wave propagation and scattering with EFIT elasto-dynamic finite integration technique[J]. Wave Motion, 1995, 21: 47-66.

[23] 底青云, 王妙月. 雷达波有限元仿真模拟[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.

[24] 冯德山, 陈承申, 王洪华. 基于混合边界条件的有限单元法正演模拟[J]. 地球物理学报, 2012, 55(11): 3774-3785.

FENG De-shan, CHEN Cheng-shen, WANG Hong-hua. Finite element method GPR forward simulation based on mixed boundary condition[J]. Chinese Journal of Geophysics, 2012, 55(11): 3774-3785.

[25] KAMAE K, IRIKURA K, PITARKA A. A technique for simulating strong ground motion using hybrid Green’s function[J]. Bull Seism Soc Am, 1998, 88: 357-367.

[26] 张慧, 李振春. 基于双变网格算法的地震波正演模拟[J]. 地球物理学报, 2011, 54(1): 77-86.

ZHANG Hui, LI Zhen-Chun. Seismic wave simulation method based on dual-variable grid[J]. Chinese Journal of Geophysics, 2011, 54(1): 77-86.

[27] 冯德山, 谢源. 基于单轴各向异性完全匹配层边界条件的ADIFDTD三维GPR全波场正演[J]. 中南大学学报(自然科学版), 2011, 42(8): 2365-2371.

FENG De-shan, XIE Yuan. Three dimensional GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD[J]. Journal of Central South University (Science and Technology), 2011, 42(8): 2365-2371.

[28] ZHU Z Q, PENG L X, LU G Y, MI S W. Borehole GPR numerical simulation of full wave field based on convolutional perfect matched layer boundary[J]. Journal of Central South University, 2013, 20(2): 764-769.

[29] 戴前伟, 王洪华, 冯德山, 陈德鹏. 基于三角形剖分的复杂GPR模型有限元法正演模拟[J]. 中南大学学报(自然科学版), 2012, 43(7): 2668-2673.

DAI Qian-wei , WANG Hong-hua, FENG De-shan, CHEN De-peng. Finite element method forward simulation for complex geoelectricity GPR model based on triangle mesh dissection[J]. Journal of Central South University (Science and Technology), 2012, 43(7): 2668-2673.

[30] 冯德山, 王洪华, 戴前伟. 基于无单元Galerkin法探地雷达正演模拟[J]. 地球物理学报, 2013, 56(1): 298-308.

FENG De-shan, WANG Hong-hua, DAI Qian-wei. Forward simulation of ground penetrating radar based on the element-free Galerkin method[J]. Chinese Journal of Geophysics, 2013, 56(1): 298-308.

[31] 方宏远, 林皋. 基于辛算法模拟探地雷达在复杂地电模型中的传播[J]. 地球物理学报, 2013, 2(56): 653-659.

FANG Hong-yuan, LIN Gao. Simulation of GPR wave propagation in complicated geoelectric model using symplectic method[J]. Chinese Journal of Geophysics, 2013, 56(2): 653-659.

[32] SAENGER E H, SHAPIRO S A. Scattering and diffraction by a single crack: an accuracy analysis of the rotated staggered grid[J]. Geophysical Prospecting, 2002, 50: 183-194.

[33] SAENGER E H, BOHLEN T. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[J]. Geophysics, 2004, 69: 583-591.

[34] 陈浩, 王秀明, 赵海波. 旋转交错网格有限差分及其完全匹配层吸收边界条件[J]. 科学通报, 2006, 51(17): 1985-1994.

CHEN Hao, WANG Xiu-ming, ZHAO Hai-bo. A rotated staggered grid finite-difference with the absorbing boundary condition of a perfectly matched layer[J]. Chinese Science Bulletin, 2006, 51(17): 1985-1994.

[35] 张鲁新, 符力耘, 裴正林. 不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用[J]. 地球物理学报, 2010, 53(10): 2470-2483.

ZHANG Lu-xin, FU Li-yun, PEI Zheng-lin. Finite difference modeling of Biot’s poroelastic equations with unsplit convolutional PML and rotated staggered grid[J]. Chinese Journal of Geophysics, 2010, 53(10): 2470-2483.

[36] 严红勇, 刘洋. 黏弹TTI介质中旋转交错网格高阶有限差分数值模拟[J]. 地球物理学报, 2012, 55(4): 1354-1365.

YAN Hong-yong, LIU Yang. Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media[J]. Chinese Journal of Geophysics, 2012, 55(4): 1354-1365.

(编辑龙怀中)

基金项目：国家自然科学基金资助项目(41374118)；中国铁路总公司科技研究开发计划项目(2014G005-B)；湖南省交通科技计划项目(201423)；湖南省住房和城乡建设厅科技计划项目(BZ201408，BZ201411)；中南大学研究生自主探索创新项目(2014zzts048)

收稿日期：2014-08-25；修订日期：2015-05-16

通信作者：戴前伟，教授，博士；电话：13574816998；E-mail：qwdai@mail.csu.edu.cn

摘要：以旋转交错网格(RSG)差分为基础，实现探地雷达(GPR)在非均匀突变介质中的正演数值模拟。通过利用旋转交错差分算子，将雷达波场各场分量设置在非均匀突变介质的基本单元中，单元内同一物性参数定义在同一位置，将对角线方向上物性参数差分值的线性组合来计算坐标轴方向上的物性参数差分值，针对非均匀特征变化剧烈的网格区域，算法中的单元网格形式不再需要进行插值，扩宽了数值稳定性条件的限制。推导了适合于探地雷达的旋转交错网格差分格式及TM_y极化模式下的更新方程组，实现了均匀介质和非均匀突变介质中的TM_y雷达波场的数值模拟，并分别从解析解、模拟剖面和波场快照的角度，与标准交错网格差分算法进行了对比。结果表明，在模拟非均匀特征明显的介质时，旋转交错网格算法可以选取更大的时间步长，数值色散程度控制更好，提高了模拟效率和成像精度，可更有效地指导非均匀突变介质中探地雷达数据的解译。