时域多分辨率法在探地雷达三维正演模拟中的应用
冯德山，戴前伟，瓮晶波
(中南大学 信息物理工程学院，湖南 长沙，410083)
摘 要：
摘 要：应用小波伽略金方法，对Maxwell的2个旋度方程进行离散化，导出DB2-MRTD算法的探地雷达3D差分公式、数值稳定性条件。在此基础上，开发探地雷达MRTD法正演模拟程序。所开发的对三维立方体小空洞体模型进行正演模拟。研究结果表明：采用该程序极大地提高了运算速度，改善了三维探地雷达正演方法；相应的正演合成三维剖视图及切片图指示了立方体小空洞体的雷达波空间的三维响应特征；模拟结果表明，采用时域多分辨率法可提高探地雷达探测的可靠性、准确度。
关键词：
中图分类号：P631 文献标识码：A 文章编号：1672-7207(2007)05-0975-06
Application of multi-resolution time domain method in three dimensions forward simulation of ground penetrating radar
FENG De-shan, DAI Qian-wei, WENG Jing-bo
(School of Info-Physics and Geomatics Engineering, Central South University, Changsha 410083, China)
Abstract: The wavelet-Galerkin method was used to disperse two Maxwell rotation equations. The GPR three dimension (3D) finite difference formula and the numerical stability of DB2-MRTD algorithm were deduced. MRTD forward simulation program of GPR was developed which can raise the operational speed greatly and improve 3D GPR forward simulation method. This self-made program processes forward simulation to the 3D small cubical pocket model and gets its 3D fence diagram and slice map of forward simulation synthesize. Multi-resolution time domain method can be used to improve the credulity and accuracy of the GPR detection.
Key words: ground penetrating radar; multi-resolution time domain method; three dimensions forward simulation; DB2 wavelets basis
探地雷达是根据高频电磁波在地下介质中的传播特性来探测地下结构和特性的地球物理方法。由于它具有高效、快速、无损、抗干扰能力强等优点，已广泛应用于工程中各个领域，成为浅层勘探的有力工 具[1]。目前，有关探地雷达探测与解释技术大都建立在二维的基础上，二维雷达探测基于地质体在垂直探测剖面的方向是均匀不变的假设，即地质体是无限延伸的，但在工程实际中，常见的地质体如岩溶洞穴是不规则的三维体，因此，二维雷达探测通常能探测到地质体是否存在及其定位，但难于对探测对象提供更为准确的信息，如探测对象的形状、空间信息。因此，为提高探地雷达探测的准确度，开展三维探测与解释技术就显得特别重要，而探地雷达的三维正演是三维探测与解释技术的基础，故研究探地雷达的三维正演显得格外重要。探地雷达的正演技术一直备受国内外很多学者的关注，在进行雷达正演模拟时常采用时域有限差分法(FDTD)进行数值计算[2]。传统的FDTD法由于受数值色散及稳定性条件的限制，空间网格步长选取λ/10~λ/20(其中，λ为波长)[3]，同时，为保证计算的准确性，用FDTD法进行三维探地雷达正演模拟时需进行更细的网格剖分，导致运行时间长，内存需求大。在此，本文作者改用时域多分辨率(multi-resolution time-domain, MRTD)法进行探地雷达三维正演计算，计算时允许选取较大的空间网格[4]。
1 MRTD法的基本原理
Maxwell方程组概括了宏观电磁场的基本规律，探地雷达为高频电磁波，故用MRTD法正演模拟探地雷达[5]，也应从Maxwell的2个旋度方程出发：
MRTD法求解电磁场计算的基本思想是：从矩量法的原理出发，将未知的场函数按选定的基函数展开，代入Maxwell方程，并将时间微分近似为中心差分，空间微分近似为空间点场值加权平均，从而建立差分迭代格式。M. Krumpholz等[6]使用MRTD法采用的Battle-Lemarie小波基，它最突出的特点是对空间进行网格划分以后，在每个波长采2个点可达到FDTD每个波长采10个点的精度，而在每个波长采2个点是Nyquist采样定理的极限，其对应的网格个数只有三维FDTD的1/125，但是，由于Battle-Lemarie小波是无限支撑的，计算中将涉及更大范围内更多的场量，计算非常复杂，编程也非常困难，为此，在实际电磁场计算中，常采用紧支撑的Daubechies小波来展开Maxwell方程。与DB2-MRTD算法相比，基于Daubechies的高阶小波基计算精度并没有明显提 高[4]，但计算时间成倍地增大，故本文作者采用DB2小波为基函数的电磁场多分辨分析(MRTD)算法，进行雷达三维正演计算。DB2-MRTD与FDTD算法相比需要更小的内存、更快的计算速度和更好的数值色散特 性[7-9]。这里着重推导基于DB2小波基的MRTD算法的探地雷达差分公式。为此，把电场和磁场利用图1所示的K.S.Yee[10]空间网格，基于DB2小波的尺度函数对电场和磁场在空间展开[11-12]，用Haar尺度函数在时间上展开，有：
图1 1个Yee氏网格单元及电磁场各分量在网格空间离散点的相互关系
Fig.1 A K.S.Yee grid cell and interrelation of electro- magnetic field components in discrete point of space grid
DB2小波尺度滤波器的系数分别为：
hk=0=0.341 506 350 946 22，
hk=1=0.591 506 350 945 87，
hk=2= 0.158 493 649 053 78，
hk=3=-0.091 506 350 945 87；
x=l?x；y=m?y；z=n?z；t=k?t。
?x，?y和?z为空间离散间隔；?t为时间离散间隔。下脚标采用了传统的Yee-FDTD结构。按照“电场与磁场在空间和时间上都是正交的”这一性质安排下脚标的顺序。
将方程(3)~(8)代入方程(1)~(2)中，应用小波伽略金方法[11, 13]，可以得到以下6个迭代方程：
2 MRTD法的稳定性条件
上面导出的探地雷达三维差分方程(9)~(14)是一种显式差分格式，与FDTD法一样，MRTD的执行也是通过按时间步推进计算电磁场在计算空间内的变化规律，这种差分格式存在稳定性问题，即时间变量步长?t与空间变量步长?x，?y和?z之间必须满足一定条件，否则将出现数值不稳定性。这种不稳定性表现为：随着计算步数的增加，被计算的场量的数值无限制地增大。因此，若用所导出的差分方程进行稳定的计算，就需要合理地选取时间步长与空间步长之间的关系，可用类似于FDTD方法[14]导出三维探地雷达DB2- MRTD算法的稳定性条件[15]，其形式为：
?t≤
。 (15)
当?x=?y=?z=?s时，稳定性条件变为：
c?x≤0.433 013?s。
3 三维探地雷达MRTD法正演实例
采用所推导的DB2-MRTD差分方程、稳定性条件再结合边界条件，就可以编程进行三维探地雷达正演计算。三维探地雷达正演实例选用立方体小空洞地电模型。立方体小空洞三维地电模型示意图如图2所示，脉冲波形为雷克子波，脉冲位置位于地表下 0.05 m处，其频率为600 MHz；模拟区域为长方体，其x，y和z 3个方向的长度分别为0.6 m，1.3 m和 0.65 m，长方体中充满了混凝土介质，其介电常数为6.0，电导率为0.001 S/m，该长方体中埋有1个立方体小空洞，立方体小空洞的长、宽、高都为0.2 m。计算过程中取网格空间步长为0.01 m，时窗长度为12 ns。测线平行于y轴，共布置了20条测线。设第1条雷达测线的位置为位于x=0.1 m处，其他每条线之间的间隔为0.025 m，最后1条测线的位置为x=0.575 m，在检测过程中采用发射天线与接收天线分离的方式，以便更好地耦合。
图2 立方体小空洞三维地电模型示意图
Fig.2 3D geoelectricity model sketch map of small cubical pocket
图3所示为立方体小空洞模型在x=0.48 m，y= 0.58 m，z=0.42 m的雷达3D立体剖视图。可见，立方体所在位置产生了绕射波，而最上面的高亮显示的分层为直达波出现的位置；同时，在立方体小空洞位置产生了2个交叉的绕射波，此交叉波为立方体小空洞x向和y向切片的局部显示；图4所示为立方体小空洞三维地电模型在x=0.20 m，y=0.38 m，z=0.39 m时x，y，z 3个方向的3D雷达切片图，可见，3个切片交叉处高亮显示的圆亦指示了立方体小空洞的位置及其空间的绕射特性；图5所示为立方体小空洞模型在x=0.21 m时的雷达切片图，此图为x向单个切片图，其中0.25~0.45 m之间的绕射波为立方体空洞体异常所致，它与二维矩形切片的响应特征一致；图6所示为立方体小空洞模型在z=0.31 m时的雷达切片及电场幅值的对比振幅图，其中“草帽形状的帽顶”形象地指示了空洞的位置及场值对比关系。图7所示为立方体小空洞三维地电模型在x=0.125，0.195，0.245，0.375和0.610 m时同时显示的雷达切片图，可见，不同位置的三维响应特征是波场从空洞位置向两端逐渐变弱。图8所示为立方体小空洞三维地电模型在z方向7个不同位置的雷达切片图，它明确地指示了图的z向的准确位置。
图3 立方体小空洞模型在x=0.48 m，y=0.58 m，z=0.42 m时的雷达3D立体剖视图
Fig.3 3D radar fence diagram of small cubical pocket model when x=0.48 m, y=0.58 m, z=0.42 m
图4 立方体小空洞三维地电模型在x=0.20 m，y=0.38 m，z=0.39 m时x，y和z 3个方向的3D雷达切片图
Fig.4 Radar slice map of xyz three directional small cubical pocket model when x=0.20 m, y=0.38 m, z=0.39 m
图5 立方体小空洞三维地电模型在x=0.21 m时的雷达切片图
Fig.5 Radar slice map of three directional small cubical pocket model when x=0.21 m
图6 立方体小空洞三维地电模型在z=0.31 m时的雷达切片及振幅图
Fig.6 Radar slice and amplitude map of three directional small cubical pocket model when z=0.31 m
图7 立方体小空洞三维地电模型不同位置时的x向的雷达切片图
Fig.7 Different position radar slice map of small cubical pocket model in x directional
图8 立方体小空洞三维地电模型不同位置时的z向的雷达切片图
Fig.8 Different position radar slice map of small cubical pocket model in z directional
通过这些雷达剖面、雷达切片的不同形式，能更加全面与细致地了解立方体小空洞的三维雷达反射特征，加深对三维雷达反射特征的认识，克服了二维雷达探测的局限性，提高了探地雷达探测的可靠性、准确度与精度，同时也说明应用MRTD法进行探地雷达三维正演模拟是可行的。
4 结 论
a. 提出了探地雷达时域多分辨(MRTD)正演模拟方法，并推导了三维DB2-MRTD法探地雷达差分方程及稳定性条件，编制了MRTD法雷达正演模拟程序。运用该程序对立方体小空洞三维模型进行正演模拟，清楚地显示了雷达波三维空间传播特性，突破了过去仅能模拟简单二维模型的局限。
b. 由于MRTD算法具有较好的数值色散特性及稳定性条件，故可以通过选取更大的空间网格步长，提高运算速度，改善三维探地雷达正演方法。
参考文献：
[1] 李大心. 探地雷达方法与应用[M]. 北京: 地质出版社, 1994.
LI Da-xin. Method and application of ground penetrating radar[M]. Beijing: Geology Press, 1994.
[2] DAI Qian-wei, FENG De-shan, HE Ji-shan. Finite difference time domain method forward simulation of complex geoelectricity ground penetrating model[J]. Journal of Central South University of Technology, 2005, 12(4): 478-482.
[3] 冯德山, 戴前伟, 何继善. 探地雷达的正演模拟及有限差分波动方程偏移处理[J]. 中南大学学报: 自然科学版, 2006, 37(2): 361-365.
FENG De-shan, DAI Qian-wei, HE Ji-shan. Forward simulation of ground penetrating radar and its finite difference method wave equation migration processing[J]. Journal of Central South University: Science and Technology, 2006, 37(2): 361-365.
[4] 李 康, 孔凡敏, 郭毅峰, 等. MRTD和高阶FDTD算法的数值色散特性的分析[J]. 系统仿真学报, 2005, 17(9): 2089-2091, 2095.
LI Kang, KONG Fan-min, GUO Yi-feng, et al. Analysis of numerical dispersion properties of MRTD and higher-order FDTD schemes[J]. Journal of System Simulation, 2005, 17(9): 2089-2091, 2095.
[5] 冯德山. 基于小波多分辨探地雷达正演及偏移处理研究[D]. 长沙: 中南大学信息物理工程学院, 2006.
FENG De-shan. GPR forward simulation and migration processing research with multi-resolution of wavelets[D]. Changsha: School of Info-Physics and Geometrics Engineering, Central South University, 2006.
[6] Krumpholz M, Winful H G, Katehi L P B. Nonlinear time-domain modeling by multiresolution time domain (MRTD)[J]. IEEE Trans on MTT, 1997, 45(3): 385-393.
[7] Saber K F, Katehi L P B. A study of dielectric resonators based on two-dimensional fast wavelet algorithm[J]. IEEE Trans Microwave and Guided Wave Letters, 1996, 6(1): 19-21.
[8] Hong I P, Yoon N, Park H K. Numerical dispersive characteristics and stability condition of the multiresolution time-domain (MRTD) method[C]//1997 Electromagnetic Compatibility Proceedings. Beijing, 1997: 455-458.
[9] 郭汉伟, 梁甸农, 刘培国, 等. 电磁场计算中FDTD与MRTD[J]. 国防科技大学学报, 2001, 23(5): 84-88.
GUO Han-wei, LIANG Dian-nong, LIU Pei-guo, et al. FDTD and MRTD in electromagnetic field numerical computation[J]. Journal of National University of Defense Technology, 2001, 23(5): 84-88.
[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] 成礼智, 郭汉伟. 小波离散变换理论及工程实践[M]. 北京: 清华大学出版社, 2005.
CHENG Li-zhi, GUO Han-wei. The engineering practice of wavelet and discrete transform theory theory[M]. Beijing: Tsinghua University Press, 2005.
[12] Cheong Y W, Lee Y M, Ra K H, et al. Wavelet-Galerkin scheme of time-dependent inhomogeneous electromagnetic problems[J]. IEEE Trans Microwave and Guided Wave Letters, 1999, 9(8): 297-299.
[13] Krumpholz M, Katehi L P B. MRTD: New time-domain schemes based on multiresolution analysis[J]. IEEE Trans Microwave and Guided Wave Letters, 1996, 44(4): 555-571.
[14] Bergmann T, Robertsson J O A, Holliger K. Numerical properties of staggered finite- difference solutions of Maxwell’s equations for ground-penetrating radar modeling[J]. Geophysical Research Letters, 1996, 23(1): 45-48.
[15] Fujiji M F, Hoefer W J R. A three-dimensional harr-wavelet—based multiresolution analysis similar to the FDTD method—derivation and application[J]. IEEE Trans on MTT, 1998, 46(12): 2463-2475.
收稿日期：2007-02-05；修回日期：2007-03-28
基金项目：“十一五”国家科技支撑计划重点资助项目(2006BAC07B00-2); 湖南省科技计划资助项目(2007FJ4055)
作者简介：冯德山(1978-)，男，湖南祁阳人，博士，讲师，从事探地雷达研究
通信作者：冯德山，男，博士，讲师；电话：0731-8836145；E-mail: fengdeshan@126.com
