基于小波多分辨探地雷达三维正演模拟及偏移处理

冯德山^{1, 2}，戴前伟^{1, 2}

 (1. 中南大学  信息物理工程学院，湖南 长沙，410083；

2. 有色资源与地质灾害探查湖南省重点实验室，湖南 长沙，410083)

摘  要：通过对Maxwell方程进行离散化，导出DB2-MRTD算法的探地雷达3D差分公式，并开发探地雷达MRTD法正演程序，利用该程序对三维雷达模型进行正演模拟，得到其相应的正演合成三维剖视图及切片图；通过对这些模拟结果进行分析，验证MRTD法在探地雷达三维正演中的有效性；利用爆炸反射原理和浮动坐标变换，推导出三维探地雷达波动方程差分格式，并把波场外推矩阵在小波域进行求解，进而得到探地雷达小波域三维偏移算法及偏移处理程序，把该程序应用于正演结果中。对比偏移处理前、后的雷达资料发现，该三维偏移算法能使3D正演剖面中的反射波归位、绕射波收敛，极大地提高了雷达剖面的分辨率。

关键词：探地雷达；时域多分辨法；正演；偏移处理

中图分类号：P631         文献标识码：A         文章编号：1672-7207(2007)04-0758-06

GPR three dimension forward simulation and migration processing with multi-resolution of wavelets

FENG De-shan^{1, 2}, DAI Qian-wei^{1, 2}

 (1. School of Info-physics and Geomatics Engineering, Central South University, Changsha 410083, China;

2. Hunan Key Laboratory of Non-ferrous Resources and Geological Hazard Detection, Changsha 410083, China)

Abstract: The GPR three-dimension finite difference formula of DB2-MRTD algorithm was deduced through dispersing the Maxwell equations. Based on this, MRTD forward simulation program of GPR was developed. Making use of this program, forward simulation to the 3D radar model was made and its 3D fence diagram and slice map of forward simulation synthesize were obtained. By analyzing these analog results, the validity of MRTD 3D forward simulation method in GPR was illuminated. Making use of the bursting reflection theory and floating coordinate transform, the finite difference format of GPR three-dimension wave equation was deduced. After the extrapolate matrix of wave field in the wavelets domain was solved, the migration algorithm of GPR three-dimension wave equation in wavelets domain and migration processing program was gotten. Comparing the radar data before and after the migration processing, the results show that this three-dimension migration algorithm can make the reflection wave return to the original position, and make the diffraction wave converge in the three-dimension sections. The lateral resolution of radar sections can be highly enhanced.

Key words: ground penetrating radar; multi-resolution time domain; forward simulation; migration processing

探地雷达的正演模拟是雷达理论研究的重要内容之一，而探地雷达的偏移处理是探地雷达数据处理(反演算法)的重要方法之一^[1]。通常利用正演结果进行反演，能验证反演算法的正确性，提高解释精度。但是，目前有关探地雷达正演模拟与偏移处理的数据解释技术大都建立在二维的基础上，二维雷达探测通常能探测到地质体的存在及定位，但难于对探测对象提供更准确的信息^[2]，如探测对象的形状、空间信息。而探地雷达三维探测技术可以克服二维探测技术的不足，有助于更加详细了解地下目标的分布和相对位置关系。在此，本文作者对作为雷达三维技术基础的三维雷达的正演模拟与偏移处理技术进行研究。

1  MRTD法的基本原理

考虑Daubechies的紧支撑性^[3]，引入以DB2小波基函数的时域多分辨分析(MRTD)算法，进行雷达三维正演。DB2-MRTD与FDTD算法相比具有更好的数值色散特性和稳定性条件，所以，可以选取更大的网格空间步长，从而需要更小的内存、更快的计算速度。M. Krumpholz等^[4-5]对MRTD与FDTD这2种算法进行了比较，在此，本文作者着重推导基于DB2小波基的MRTD算法的探地雷达差分公式。其方法是：从Maxwell的2个旋度方程出发，将时间微分近似为中心差分，空间微分近似为空间点场值加权的平均，基于DB2小波的尺度函数对电场和磁场在空间展开，用Haar尺度函数在时间上展开^[6]，并以E_x和H_x电磁场分量为例(其他分量类似可得到)：

下脚标采用传统的Yee-FDTD结构，按照电场与磁场在空间和时间上都是正交这一性质安排下脚标的顺序，应用小波伽略金方法^[7-8]，可以从Maxwell旋度方程中推得下列6个迭代方程：

利用上述探地雷达差分方程，结合边界条件，对图1所示的地电模型进行正演合成，其合成的切片图如图2和图3所示。

混凝土：介电常数为6.0，电导率为0.001 S/m

图1  雷达三维正演地电模型示意图

Fig.1  Radar geoelectricity model sketch map of 3D forward simulation

图2  雷达三维正演不同位置时的y向切片图

Fig.2  y directional slice map of 3D forward simulation

图3  雷达三维正演不同位置时的x向切片图

Fig.3  x directional slice map of 3D forward simulation

2  探地雷达小波域三维偏移算法基本原理

由于雷达波在地下传播时，在高频情况下，介质的位移电流远大于传导电流^[9]，故对于探地雷达在三维空间^[10]传播的Maxwell方程就可以写成：

偏移的目的就是根据式 (9)中地面波场值E(x, y, z)=E(x, y, z=0, t)外推成像。应用地震偏移^[11]处理理论中的爆炸反射成像原理，可取式(1)中的V=v/2，使波动方程按单程时间计算，得：

而在地层倾角较小时，可略去对τ的2阶偏导数项，作傍轴近似处理，从而得到三维的差分偏移方程：

采用维数分裂法^[12] 对式(15)求解，可得到交替方向上的2个分裂方程：

对于这2个分裂方程，相当于分别在2个不同方向做2步偏移。在不同方向，2步偏移的方法一样，只是第2步的偏移是在第1步偏移结果的基础上进行的。首先对式(16)求解(式(17)的求解类似)。先假定速度只在垂向上变化而在横向上没有变化，并且将时间方向和水平方向上微分算子利用对应的有限差分算子代替，则有

将每个深度上的雷达波场表示成二维矩阵的形式，行表示时间方向，列表示空间方向，即雷达波场矩阵E={E_{i, l}}, i=0, 1, …, nx-1, l=0, 1, …, nt-1。并且将时间方向上的一阶有限差分算子表示成算子矩阵A，空间方向上的二阶有限差分算子表示成算子矩阵B，那么，深度方向上波场外推过程的矩阵形式^[13]为：

如果仅仅在时间—空间域中利用式(21)进行波场外推，计算效率会非常低。其原因是：一方面，矩阵A是由有限差分算子得到的算子矩阵，因为它的正则值随着矩阵维数增加而很快增大，造成A^-1的求解效率低；另一方面，雷达波场外推过程^[14]涉及3个矩阵相乘，其算法的复杂度大致为O(N_xN_t)。可以将时间-空间域中的波场外推过程在小波域中进行求解，其效率大大提高。通过阈值运算，可以将雷达波场标准形式的多分辨分解压缩成O(N_xlogN_t)个元素的稀疏矩 阵^[15]，而且算子A^-1和B标准形式的多分辨表示也是稀疏的^[15]，这样，三者之间乘积的算法复杂度大约为O(N_xlogN_t)，其小波域中的具体求解过程为：

同理，对于应用式(23)偏移得到的雷达波波场，在1个外推步长?τ内，在y方向上进行小波域差分外推运算，其结果即为小波域中在深度方向外推1个步长?τ的值。在此外推结果的基础上，再进行下一个步长的外推，直到t=τ，即达到所需要的深度为止。

3  探地雷达三维正演及偏移处理实例

三维地电模型如图1所示，脉冲波形为雷克子波，其频率为600 MHz；模拟区域为：在x，y和z方向的长度分别为0.60，1.30 和0.65 m；其中充满了混凝土介质，其介电常数为6.0，电导率为0.001 S/m，混凝土介质中埋有3个相同的球体空洞，球心的坐标分别为(0.20, 0.30, 0.20)，(0.30, 0.60, 0.35)，(0.45, 1.00, 0.40 )，半径为0.05 m。计算时取网格空间步长为0.01 m，时窗长度为15 ns。测线平行于y轴，共布置20条测线。在检测过程中，采用发射天线与接收天线分离的方式，以便更好地耦合。利用MRTD法对图1中3个球体小空洞雷达模型进行正演模拟，然后，应用小波域偏移对正演结果进行处理，以检验偏移算法的效果。

图2所示为3个球状空洞体不同位置的y向雷达切片图，图3所示为3个球状空洞体不同位置的x向雷达切片图。显然，由于3个小空洞的位置与埋深的变化导致切片图中绕射波形发生变化。图4所示为3个球状空洞体三维正演x=0.24 m，y=0.55 m和z=0.37 m的雷达立体剖视图，图5所示为3个球状空洞体三维正演x=0.24 m，y=0.55 m和z=0.37 m的雷达切片图，由于模型中存在3个球体小空洞，使得雷达三维正演剖面中3个球体小空洞所在的位置产生双曲线绕射波。正是由于存在这些绕射波，雷达正演合成图不能真实地反映3个球体小空洞的空间形态。为了提高雷达资料的解释精度，必须对正演合成的三维雷达剖面进行偏移处理。图6和图7所示为3个球状空洞体正演剖面经小波域偏移处理后三维效果图，它们分别对应偏移处理前的三维效果图图4和图5。通过对比偏移处理前、后的三维雷达图片可知，该偏移算法对坐标分别为(0.20, 0.30, 0.20)，(0.30, 0.60, 0.35)，(0.45, 1.00, 0.40)的3个球状空洞体，对所产生的绕射波进行聚焦处理，其发散的能量收敛到球体空洞所在的位置，基本恢复了球体小空洞的形状。由此可见，小波域三维偏移算法可以提高三维雷达图像的分辨率。

图4  雷达三维正演立体剖视图

Fig.4  Radar fence diagram of 3D forward simulation

图5  雷达三维正演的xyz向切片图

Fig.5  Radar xyz direction slice map of 3D forward simulation

图6  经小波域偏移处理后的三维雷达立体剖视图

Fig.6  Radar fence diagram of 3D after wavelets domain migration

图7  经小波域偏移处理后的xyz向雷达切片图

Fig.7 Radar xyz direction slice map of 3D after wavelets domain migration

4  结  论

a. 提出了探地雷达时域多分辨(MRTD)正演模拟方法，并推导了三维DB2-MRTD法探地雷达差分方程及稳定性条件，编制了MRTD法雷达正演模拟程序，该程序能快速实现探地雷达三维地电模型的正演模拟，克服过去仅能模拟简单二维模型的局限，并且实现了雷达三维正演数据的可视化，改善了三维探地雷达正演方法。

b. 提出了探地雷达小波域探地雷达三维偏移处理算法，并编制了相应的偏移处理程序。该程序能用于快速实现探地雷达三维偏移处理，极大地提高运算速度，改善三维探地雷达偏移方法。通过对比偏移处理前、后的雷达资料，得知该小波域三维偏移算法能使三维正演剖面中反射波归位、绕射波收敛，提高了雷达探测的可靠性、准确度及雷达剖面的分辨率，更有利于雷达资料的地质解释。

参考文献：

[1] 李大心. 探地雷达方法与应用[M]. 北京: 地质出版社, 1994.
LI Da-xin. Method and application of ground penetrating radar[M]. Beijing: Geology Press, 1994.

[2] 冯德山, 戴前伟, 何继善. 探地雷达的正演模拟及有限差分波动方程偏移处理[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.

[3] 崔锦泰. 小波分析导论[M]. 西安: 西安交通大学出版社, 1995.
CUI Jin-tai. An introduction to wavelets[M]. Xi’an: Xi’an Jiaotong University Press, 1995.

[4] 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.

[5] 郭汉伟, 梁甸农, 刘培国, 等. 电磁场计算中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.

[6] 成礼智, 郭汉伟. 小波离散变换理论及工程实践[M]. 北京: 清华大学出版社, 2005.
CHENG Li-zhi, GUO Han-wei. The engineering practice of wavelet and discrete transform theory theory[M]. Beijing: Tsinghua University Press, 2005.

[7] 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.

[8] 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.

[9] 李成方. 偏移处理在探地雷达数据处理中的研究[D]. 成都: 成都理工大学信息工程学院, 2002.
LI Cheng-fang. The application of migration technique to processing GPR data [D]. Chengdu: College of Information Engineering, Chengdu University of Technology, 2002.

[10] 雷文太, 粟 毅, 黄仕家. 探地雷达近场三维距离偏移成像算法[J]. 电子与信息学报, 2003, 25(12): 1641-1646.
LEI Wen-tai, SU Yi, HUANG Shi-jia. Ground penetrating radar near field 3-D range migration imaging technique[J]. Journal of Electronics & Information Technology, 2003, 25(12): 1641-1646.

[11] 马在田. 地震成像技术有限差分法偏移[M]. 北京: 石油工业出版社, 1989.
MA Zai-tian. Finite-difference method migration of seismic image technique[M]. Beijing: Petroleum Industry Press, 1989.

[12] 金胜汶, 陈必远, 马在田. 三维波动方程有限差分正演方法[J]. 地球物理学报, 1994, 37(6): 805-810.
JIN Sheng-wen, CHEN Bi-yuan, MA Zai-tian. Three- dimensional wave equation forward modeling by the finite-difference method[J]. Chinese Journal of Geophysics, 1994, 37(6): 805-810.

[13] 冯德山. 基于小波多分辨探地雷达正演及偏移处理研究[D]. 长沙: 中南大学信息物理工程学院, 2006.
FENG De-shan. GPR forward simulation and migration processing research with multi-resolution of wavelets[D]. Changsha: School of Info-physics and Geomatics Engineering, Central South University, 2006.

[14] 张安学, 蒋延生, 汪文秉. 探地雷达扫频三维成像方法[J]. 电波科学学报, 2000, 15(3): 313-316.
ZHANG An-xue, JIANG Yan-sheng, WANG Wen-bing, et al. 3D image by sweeping frequency in ground penetrating radar system[J]. Chinese Journal of Radio Science, 2000, 15(3): 313-316.

[15] 张海江. 小波多分辨波动方程正演模拟与偏移成像[D]. 北京: 中国石油天然气总公司石油勘探开发科学研究院, 2000.
ZHANG Hai-jiang. Wave equation forward simulation and migration imaging with multi-resolution of wavelets[D]. Beijing: Research Institute of Petroleum Exploration and Development, China National Petroleum Corporation, 2000.

收稿日期：2006-11-24

基金项目：国家自然科学基金资助项目(60672042)

作者简介：冯德山(1978-), 男, 湖南祁阳人, 博士, 讲师, 从事探地雷达研究

通讯作者：冯德山, 男, 博士；电话: 0731-8836145；E-mail: fengdeshan@126.com