探地雷达小波域三维波动方程偏移

冯德山,戴前伟.探地雷达小波域三维波动方程偏移.地球物理学报,2008,51(2):566～574Feng D S ,Dai Q W.The migration of G PR three dimension wave equation in wavelets domain.Chinese J .G eophys .(in Chinese ),2008,51(2):566～574探地雷达小波域三维波动方程偏移冯德山,戴前伟中南大学信息物理工程学院,有色资源与地质灾害探查湖南省重点实验室,长沙　410083摘　要　阐述了矩阵多分辨分析理论中的标准形式与非标准形式,并以Hilbert 算子为例,说明了算子多分辨表示的压缩效果,为小波域偏移算法奠定了理论基础.从三维雷达波动方程出发,利用爆炸反射原理和浮动坐标变换,推导出三维探地雷达波动方程差分格式,并通过方程分裂算法及小波多分辨算法,在小波域求解波场外推矩阵,进而得到探地雷达小波域三维波动方程偏移算法,在此基础上,开发了探地雷达小波域偏移处理程序,并把该程序应用于三个球体空洞的32D 正演结果及实际的雷达数据中,通过对比偏移处理前后的雷达资料,得知该三维偏移算法能使32D 正演剖面中的反射波归位、绕射波收敛,极大地提高了雷达剖面的分辨率,有利于探地雷达资料的地质解释.关键词　探地雷达,小波域,波动方程,偏移文章编号　0001-5733(2008)02-0566-09中图分类号　P631收稿日期2007-02-25,2007-12-11收修定稿基金项目　国家自然科学基金项目(60672042)与湖南省科技计划项目(2007F J4055)联合资助.作者简介　冯德山,男,1978年生,博士,讲师,主要从事探地雷达与地震勘探的研究.E 2mail :fengdeshan @1261comThe migration of GPR three dimension w ave equation in w avelets dom ainFE NG De 2Shan ,DAI Qian 2WeiSchool o f In fo 2Physics and G eometrics Engineering ,Central South Univer sity ,Changsha 410083,China Non 2ferrous Resources and G eologic Disaster s Pro specting Laboratory o f Hunan ,Changsha 410083,ChinaAbstract This paper elaborated on the standard and nonstandard forms of matrix multi 2res olution analysis theory in detail ,and taking Hilbert operator as an exam ple ,explained the com pression result of operator with multi 2res olution form.It laid a g ood theoretical foundation for migration alg orithm of wavelets domain.Then ,setting out from the three 2dimension radar wave equation ,and making use of the bursting reflection theory and floating coordinate trans form ,we have deduced the finite difference format of G PR three 2dimension wave equation.Through equation splitting and multi 2res olution wavelets theory ,and s olving the extrapolate matrix of wave field in the wavelets domain ,we have als o g ot the migration alg orithm of G PR three 2dimension wave equation in wavelets domain.Based on this ,the authors developed the migration program of G PR three 2dimension wave equation ,and applied this program to three 2dimensional forward m odeling result of three spheric caverns and practical G PR data.Through com paring the radar data before and after the migration processing ,it is known that this three 2dimension migration alg orithm could make the reflection wave return to original position ,and make the diffraction wave converge in the three 2dimension sections.The lateral res olution of radar sections could be highly enhanced ,and the migration alg orithm could make the radar three 2dimension detection m ore reliable and precise ,which is propitious to the geology explanation of G PR data.K eyw ords G round penetrating radar ,Wavelets domain ,Wave equation ,Migration第51卷第2期2008年3月地　球　物　理　学　报CHI NESE JOURNA L OF GE OPHY SICSV ol.51,N o.2Mar.,20081　引　言 探地雷达波偏移是利用波动方程及其给定的初始、边界条件,通过反向外推波场重建地下地质构造[1].偏移的目的就是使雷达剖面中岩性突变点产生的绕射波,凹界面产生的特殊回转波,真实归位[2],使地下界面得到真实明确的反映,从而更好地指导G PR剖面的地质解释和验证偏移方法的有效性[3].目前,有关探地雷达探测与解释技术大都建立在二维的基础上[4],二维雷达探测通常能定位出地质体的存在,但难以给出探测对象的三维空间信息[5].如Armando等[6]在处理有耗介质中传播的雷达数据时,应用了分裂步傅里叶偏移算法,E lizabeth Fisher[7]采用逆时偏移算法对单通道的雷达实例数据进行了处理等.而探地雷达三维探测可以解决二维的不足[5],如S zerbiak[8]在研究水库模型时应用了三维的雷达探测技术,Streich和K ruk[9]等给出了Vector2Migration技术,它可对三维雷达数据进行偏移处理,这些技术的发展势必促进三维雷达偏移技术的更深入研究[10].而探地雷达的小波域三维偏移,可以通过研究二维矩阵的多分辨分析理论,把三维雷达偏移通过分裂算法,转化为两步在小波域中求解二维雷达外推矩阵问题,从而在小波域中提高偏移算法效率.2　矩阵的多分辨表示 雷达波场外推过程主要是利用波动方程,将地表接收到的雷达波场递归地外推到地下某一深度处,它在频率域中是一个广义的空间卷积运算,外推算子可以表示成矩阵的形式,通常情况下这个矩阵是稠密的,而小波多分辨分析可使雷达波场外推矩阵变成稀疏矩阵[11],从而极大地提高了雷达偏移效率.在二维雷达波场外推矩阵的小波多分辨分解[12,13]中,由小波理论可知[14,15],二维矩阵存在两种形式的多分辨分解:标准形式(长方形)和非标准形式(正方形).尽管这两种多分辨分解存在差别,但它们之间可相互转换[16].由于探地雷达的小波域中求解要用到多分辨分解,故这里对于这两种多分辨分解进行详细的阐述.211　非标准形式设最粗和最细尺度对应的空间分别为Vj 和V,最粗的小波空间为Wj,则V0=Jj=0W j+V j,又记P j和Qj=P j-1-P j分别是从L2(R)到V j和W j的投影算子,其中Pj:L2(R)→V j为光滑投影算子,Q j: L2(R)→W j为细节投影算子,即:P j f(x)=∑l〈f,φj,l〉φj,l(x),Q j f(x)=∑l〈f,ψj,l〉ψj,l(x).(1) 那么对于任何给定的算子T=K(x,y):L2(R)→L2(R),可用多分辨表示为T=∑j∈Z(QjTQ j+Q j TP j+P j TQ j),(2)则T在最细分辨尺度的近似或离散为[17]T≈T0=P0TP0=∑Jj=0(Pj-1TP j-1-P j TP j)+P J TP J=∑Jj=0[(P j-1-P j)T(P j-1-P j)+(P j-1-P j)TP j+P j T(P j-1-P j)]+P J TP J=∑Jj=0[(Q j TQ j+Q j TP j+P j TQ j]+P J TP J.(3) 记Aj=Q j TQ j,B j=Q j TP j,Γj=P j TQ j,T j= P j TP j,则算子T的非标准形式可表示为T={A j,B j,T j}j∈Z,(4)算子的多分辨率表示的非标准形式可以通过递归描述为T j=A j+1B j+1Γj+1T j+1,(5)如果j为最大尺度,则T={{Aj,B j,T j}j∈Z,j≤n,T n},其中Tn=P n TP n.如果j=n有限,则非标准形式可以写成如图1的分块矩阵形式[18].212　标准形式根据小波多分辨分析,Vj可以表示为V j=j′>jW j′,在尺度j上{B j,Γj}就可以表示为算子序列{B j′j,Γj′j}j′>j,其中B j′j:W j′→W j,Γj′j:W j→W j′,这样就得到与非标准形式T={Aj,B j,Γj}j∈Z相应的标准形式[19]:T={A j,{B j′j}j′>j,{Γj′j}j′>j}j∈Z,(6)设最粗和最细尺度分别为j=n,j=0,那么Πj∈[0,n],V j=V nnj′=j+1W j′,所以算子的标准形式中除了765　2期冯德山等:探地雷达小波域三维波动方程偏移图1　矩阵的多分辨率分解(a)非标准形式;(b)标准形式.Fig.1　Multi2res olution decom position of a matrix(a)N on2standard form;(b)S tandard form.{B j′j,Γj′j}j=j+1,…n之外还要增加算序列{B n+1j}和{Γn+1j},其中B n+1j:V n→W j,Γn+1j:W f→V n.特别地,B n+1n 和Γn+1n就是非标准形式中的Bn和Γn.综上所述,最细尺度上的算子T=P0TP0可以表示为如下的标准形式:T0={A j,{B j′j}n j′=j+1,{Γj′j}n j′=j+1,B n+1j,Γn+1j,T n}j=1,2,…,n.(7) 式(7)可以用如图1b所示的分块矩阵形式表示. 3　小波多分辨的压缩特性 二维雷达图像的一级多分辨分解可以分成4个部分:原始雷达图像的低频逼近(LL),原始雷达图像对角线上的细节(HH),原始雷达图像水平方向上的边缘细节(LH),原始雷达图像垂直方向上的边缘细节(H L).如果对二维雷达数据进行多级小波多分辨分解就可得到两种形式的多分辨分解:标准形式和非标准形式.其分解实质上是从一维数据多分辨分解的张量积得到的,其算法为:先对二维雷达数据的行(或列)做一维小波多分辨分解,后对得到的中间矩阵的列(或行)做一维小波多分辨分解.二维雷达数据小波域中多分辨表示形式是通过对二维雷达数据的低频逼近部分LL做进一步分解得到的.把雷达数据进行小波多分辨分解的目的是为了通过阈值运算使二维雷达波场外推矩阵变得稀疏,大大压缩矩阵的元素,从而提高偏移处理速度.小波在时空域和频率域都具有局部性.解可能在大部分区域是光滑的,只是在少量点发生突变,这样便为利用小波进行数据压缩提供可能[20～22].图2(a～c)分别为Hilbert算子在非标准形式表示下的无数据截断、截断小于10-12、10-1不同阈值数据情况下,DB6小波基的多分辨表示及其重构组图,从这些图中可以总结出,Hilbert变换算子的非标准矩阵在不同的阈值截断情况下,将呈现不同程度的稀疏形式,图2c中可见,当截断小于1×10-1数据时, Hilbert变换算子的非标准矩阵已变成了极度稀疏的矩阵,由此可见小波对算子的压缩效果.针对雷达波场外推矩阵,首先进行小波多分辨分解,分解后大部分元素近似为零,不为零的元素非常稀疏,而且大部分非零元素集中分布在各分辨空间的对角线附近,此时再对雷达波场函数的小波多分辨分解设置很小的阀值,将小于此阀值的雷达波场值充零,而将大于此阀值的雷达波场值保留,这样雷达波场外推稠密矩阵将变为非常稀疏的矩阵,估计小波域中的雷达波偏移算法应该具有很快的速度.4　探地雷达小波域偏移算法的基本原理 雷达波在地球内部的传播满足三维Maxwell方程,所以根据地面接收到的反射波场准确地重建地下界面的形态位置,必须把二维偏移问题扩展到三维偏移.由于探地雷达的高频电磁波在地下传播时,介质的位移电流远大于传导电流,所以对于探地雷达在三维空间传播的Max well方程[23]就可以写成:Δ2E=1V292E9t2=92E9x2+92E9y2+92E9z2=1V292E2.(8)865地球物理学报(Chinese J.G eophys.)51卷　图2　Hilbert算子多分辨表示及其重构组图(a)无数据截断;(b)截断小于1×10-12数据;(c)截断小于1×10-1数据.Fig.2　The group map of Hilbert operator multi2res olution decom position and reconstruction (a)W ithout data cut off;(b)Cut off the data less than1×10-12;(c)Cut off the data less than1×10-1.其中E=E(x,y,z,t)仍表示雷达波场值,V=V(x,y,z)为介质速度,x,y,z和t分别为空间和时间的坐标分量,其相对应的空间和时间采样间隔分别为Δx,Δy,Δz和Δt.偏移的目的就是根据式(8)中地面波场值E(x,y,z)=E(x,y,z=0,t)外推成像.首先应用地震偏移[24]处理理论中的爆炸反射成像原理,可取(1)式中的V=vΠ2,这样把波动方程改成按单程时间计算,可得:92E9x2+92E9y2+92E9z2=4v292E9t2.(9)然后,进行浮动坐标变换:x′=x,y′=y,τ=2zΠv,t′=t+τ,(10)式中τ为拟深度分量,由于坐标变换不改变实际波场,故原坐标系中的波场E(x,y,z,t)和新坐标系965　2期冯德山等:探地雷达小波域三维波动方程偏移中的波场^E(x′,y′,τ,t′)是相等的,即E(x,y,z,t) =^E(x′,y′,τ,t′),由复合函数微分法得[25]:9E 9x=9^E9x′,　9E9y=9^E9y′,　9E9x2=92^E9x′2,92E 9y2=92^E9y′2,　92E9t2=92^E9t′2,(11)因为9τ9z=2v,t′=t+τ=t+2zv,9t′9z=2v,所以有:9E 9z=9^E9x′9x′9z+9^E9y′9y′9z+9^E9τ9τ9z+9^E9t′9t′9z=9^E9τ9τ9z+9^E9t′9t′9z=2v9^E9τ+2v9^E9t′,(12)92E 9z2=2v99z2v9^E9τ+2v9^E9t′=2v2v9^E9τ9^E9τ+9^E9τ9^E9t′+2v9^E9t′9^E9τ+9^E9t′9^E9t′=4v292^E9τ2+9^E9τ9^E9t′+9^E9t′9^E9τ+92^E9t′2=4v292^E9τ2+292^E9τ9t′+92^E9t′2.(13) 把式(11)、(13)代入式(9)中可以得到:92^E 9x′2+92^E9y′2+4v292^E9τ2+292^E9τ9t′+92^E9t′2=4v292^E9t′2.(14)对(14)式化简,略去肩标,并把波场^E写成原来的E 形式,可得浮动坐标系下的偏移方程[26]:v2 892E9x2+92E9y2+1292E9τ2+92E9τ9t=0.(15) 而当地层倾角较小的情况下,可略去对τ的两阶偏导数项,作傍轴近似,从而得到三维的15°差分偏移方程:v2 892E9x2+92E9y2+92E9τ9t=0.(16) 使用维数分裂法对式(16)求解,可得到交替方向上的两个分裂方程:92E9τ9t=-v2892E9x2;(17)92E9τ9t=-v2892E9y2.(18) 对于这两个分裂方程,相当于分别在两个不同方向做两步偏移[27],两步偏移的方法是一样的,只是第二步的偏移是在第一步偏移结果的基础上进行的.首先讨论对(17)式的求解,(18)式的求解类似,先假定速度只在垂向上变化而在横向上没有变化,并且将时间方向和水平方向上微分算子利用对应的有限差分算子代替[28],则有1Δt99τ(δtE(τ))=-v(τ)8Δx2δxxE(τ).(19) 将每个深度上的雷达波场表示成二维矩阵的形式,行表示时间方向,列表示空间方向,即雷达波场矩阵E={Ei,l},i=0,1,…,nx-1,l=0,1,…,nt-1.并且将时间方向上的一阶有限差分算子表示成算子矩阵A,空间方向上的二阶有限差分算子表示成算子矩阵B,那么深度方向上波场外推过程的矩阵形式为1Δt99τ(EA)=-v(τ)8Δx2B E.(20)将(20)式中关于垂直方向上的一阶微分算子离散表示,得到(E(τ+Δτ)-E(τ))A=-Δtv(τ)8Δx2B・E(τ).(21)对矩阵A求逆,可得E(τ+Δτ)=-Δtv(τ)8Δx2B・E(τ)・A-1+E(τ).(22) 如果仅仅在时间-空间域中利用(22)式进行波场外推,则计算效率将会非常低[29,30].这一方面是,由于矩阵A是由有限差分算子得到的算子矩阵,它的正则值随着矩阵维数增加而很快增大,造成A-1的求解效率低下;另一方面是由于雷达波场外推过程[20]涉及三个矩阵相乘,其算法的复杂度大致为O(N x N t).我们可以将时间-空间域中的波场外推过程在小波域中进行求解,其效率有很大的提高,因为通过阈值运算,可以将雷达波场标准形式的多分辨分解压缩成O(Nxlog N t)个元素的稀疏矩阵[24],而且算子A-1和B标准形式的多分辨表示也是稀疏的,这样三者之间乘积的算法复杂度大致为O(N x log N t),其小波域中的具体求解过程为WE(τ+Δτ)W T=-Δt・v(τ)8Δx2WBW T・WE(τ)W T・WA-1W T+WE(τ)W T.(23)设W为小波变换矩阵, E=WEW T, B=WBW T, A-1=WA-1W T,则(23)式可以重写为E(τ+Δτ)=-Δt・v(τ)8Δx2B・ E(τ)・ A-1+ E(τ).(24)075地球物理学报(Chinese J.G eophys.)51卷　其中 E 表示雷达波场向量E 的小波多分辨分解, B 表示矩阵B 的标准形式的多分辨分解, A 表示矩阵A 的标准形式的多分辨分解,这些都可由前述的多分辨分解方式得到,这些矩阵在一定程度上表现为极度稀疏的形式,只在对角线区和其他条带区域分布着几个重要的非零元素,而在远离对角线的地方,在些小波系数的值变得非常小,通过阀值运算将其略去之后,对偏移处理的波场影响不大,但偏移处理速度可以大大地提高,因此可以对 E 、 B 、 A 设置一个阀值系数ε1,将其中小于此阀值的元素充零,而保留大于此阀值的元素.同理,对于应用(24)式偏移得到的雷达波波场值,在一个外推步长Δτ内、y 方向上进行小波域差分外推运算,其结果即为小波域中在深度方向外推一个步长Δτ的值.在此外推结果的基础上,再进行下一个步长的外推,一直到t =τ时,达到所需要的深度为止,这时可得到所要求的最终的三维雷达偏移成像结果.图4　经小波域偏移处理后的三维雷达正演的剖视图与切片图(a )剖视图;(b )xyz 向切片图;(c )xy 向切片图.Fig.4　The 32D radar fence diagram and slice map after wavelets domain migration processing(a )Fence diagram ;(b )xyz direction slice map ;(c )xy direction slicemap.图3　三维雷达正演的剖视图与切片图(a )剖视图;(b )xyz 向切片图;(c )xy 向切片图.Fig.3　The radar fence diagram and slice map of 3D forward simulation(a )Fence diagram ;(b )xyz direction slice map ;(c )xy direction slice map.5　探地雷达小波域三维偏移处理实例511　三维正演实例的偏移处理利用MRT D 法对三个球体小空洞雷达模型进行了正演模拟[16],然后应用本文介绍的小波域偏移算法对正演结果进行处理,以检验偏移算法的效果.图3a 为三个球状空洞体三维正演x =0124m ,y =0155m ,z =0137m 的雷达立体剖视图,图3b 为三个球状空洞体三维正演x =0124m ,y =0155m ,z =0137m 的雷达切片图,图3c 为三个球状空洞体三维正演x =0124m ,y =0155m xy 向的雷达切片图,由于模型中三个球体小空洞的存在,使得雷达三维正演剖面中三个球体小空洞所在的位置产生了双曲线绕射波,正是这些绕射波的存在,使得雷达正演合成图不能真实地反映三个球体小空洞的空间形态,为了提高雷达资料的解释精度,必须对正演合成的三维雷达剖面进行偏移处理.175　2期冯德山等:探地雷达小波域三维波动方程偏移 图4a 为三个球状空洞体三维正演剖面经小波域偏移处理后x =0124m ,y =0155m ,z =0137m 的雷达立体剖视图,图4b 为三个球状空洞体三维正演剖面经小波域偏移处理后x =0124m ,y =0155m ,z =0137m 的雷达切片图;图4c 为三个球状空洞体三维正演剖面经小波域偏移处理后x =0124m ,y =0155m xy 向的雷达切片图.图6　小波域偏移处理后的雷达实测数据三维立体剖视图与切片图(a )剖视图;(b )xyz 向切片图.Fig.6　The 32D radar fence diagram and slice map of practical data after wavelets domain migration processing(a )Fence diagram ;(b )xyz direction slicemap.图5　雷达实测数据三维立体剖视图与切片图(a )剖视图;(b )xyz 向切片图.Fig.5　The 32D radar fence diagram and slice map of practical data(a )Fence diagram ;(b )xyz direction slice map.通过对比偏移处理前后的三维雷达图片可知,该偏移算法对坐标分别为(012m ,013m ,012m )、(013m ,016m ,0135m )、(0145m ,110m ,014m )的三个球状空洞体所产生的绕射波进行了聚焦处理,其发散的能量收敛到了球体空洞所在的位置,基本恢复了球体小空洞的形状;由此可见,小波域三维偏移算法可以提高三维雷达图像的分辨率.512　雷达实测数据的三维探地雷达小波域偏移处理图5(a ～b )为雷达实测数据x =2164m ,y =11185m ,z =3128m 时的三维雷达立体剖视图及其切片图,由于该雷达数据中存在着许多的双曲线绕射波,故对资料解释带来了一定的难度.图6(a ～b )为小波域偏移算法处理后x =2164m ,y =11185m ,z =3128m 的三维雷达立体剖视图及切片图.通过对比偏移处理前后的雷达剖面得知,偏移处理后的雷达剖面中的双曲线绕射波基本收敛,其能量也归位良好,说明了该小波域偏移算法在处理三维实际雷达资料方面的有效性. 由此可见,可把此偏移算法应用到资料处理中去,它可使绕射波收敛,使断点更清晰,并可消除多次波产生的混淆,对提高雷达记录的横向分辨率也275地球物理学报(Chinese J.G eophys.)51卷　具有重要作用,为地质解释提供更可靠的雷达剖面成果图.6　结论与讨论 (1)详细地阐述了二维多分辨分析理论中的标准形式与非标准形式的分解方式,并针对小波在时空域和频率域都具有局部性的特点,以Hilbert算子为例,说明了应用小波进行矩阵多分辨表示的压缩效果,为小波域偏移算法奠定了理论基础.(2)针对雷达波场外推矩阵,首先进行小波多分辨分解,分解后大部分元素近似为零,不为零的元素非常稀疏,且大都集中分布在各分辨空间的对角线附近,此时再对雷达波场函数的小波多分辨分解设置很小的阀值,将小于此阀值的雷达波场值充零,保留大于此阀值的雷达波场值,这样雷达波场波场外推稠密矩阵将变为非常稀疏的矩阵,显然此算法思路适用于二维雷达偏移算法.(3)提出了探地雷达小波域探地雷达三维偏移处理算法,并编制了相应的偏移处理程序,该程序能快速实现探地雷达三维偏移处理,改善了三维探地雷达偏移方法,通过对比偏移处理前后的雷达资料,得知该小波域三维偏移算法能使三维正演剖面中反射波归位、绕射波收敛,提高了雷达探测的可靠性、准确度及雷达剖面的分辨率,更有利于雷达资料的地质解释,同时该算法在实际的三维雷达资料处理中也具有广阔的发展前途.参考文献(References)[1]　冯德山,戴前伟,何继善.探地雷达的正演模拟及有限差分波动方程偏移处理.中南大学学报(自然科学版),2006,37(2):361～365 Feng D S,Dai Q W,He J S.F orward simulation of ground penetrating radar and its finite difference method wave equation m igrationprocessing.Journal o f Central South Univer sity Technology:Scienceand Technology(in Chinese),2006,37(2):361～365[2]　底青云,许　琨,王妙月.衰减雷达波有限元偏移.地球物理学报,2000,43(2):257～263 Di Q Y,Xu K,W ang M Y.The attenuated radar wave m igration with finite element method.Chinese J.G eophys.(in Chinese),2000,43(2):257～263[3]　周　辉,韩　波,邱东玲等.地质雷达资料的偏移速度分析和叠前偏移.吉林大学学报(地球科学版),2005,35(2):248～252,256 Zhou H,Han B,Qiu D L,et al.M igration velocity analysis and prestack m igration of ground penetrating radar.Journal o f JilinUniver sity(Earth Science Edition).(in Chinese),2005,35(2):248～252,256[4]　黄　伟,邓世坤.偏移方法用于探地雷达图像处理的有效性研究.地质科技情报,2001,20(4):99～102 Huang W,Deng S K.S tudy on effect of m igration method in G PR data processing.G eological Science and Technology Information(inChinese),2001,20(4):99～102[5]　卢成明,秦　臻,朱海龙等.探地雷达检测公路结构层隐含裂缝实用方法研究.地球物理学报,2007,50(5):1558～1568 Lu C M,Qin Z,Zhu H L,et al.Practical methods for detection of concealed cracks in highway pavement using ground penetration radardata.Chinese J.G eophys.(in Chinese),2007,50(5):1558～1568 [6]　Armando R S,Paul L S,Mrinal K S,et al.S plit2step F ourierm igration of G PR data in lossy media.G eophysics,2006,71(4):77～91[7]　E lizabeth F,G eorge A M,Annant A P,et al.Exam ples of reverse2time m igration of single2channel,ground2penetrating radar profiles.G eophysics,1992,54(7):577～586[8]　S zerbiak R B,M cM echan G A,C orbeanu R,et al.32Dcharacterization of a clastic reserv oir analog:From32D G PR data to a32D fluid permeability m odel.G eophysics,2001,66(4):1026～1037 [9]　S treich R,K ruk J V D,G reen A G.Vector2m igration of standardcopolarized3D G PR data.G eophysics,2007,72(5):65～75 [10]　张安学,蒋延生,汪文秉.探地雷达扫频三维成像方法.电波科学学报,2000,15(3):313～316 Zhang A X,Jiang Y S,W ang W B,et al.3D image by sweeping frequency in ground penetrating radar system.Chinese Journal o fRadio Science(in Chinese),2000,15(3):313～316[11]　张海江,刘雯林.波场外推算子的多分辨分解和压缩.石油地球物理勘探,2000,35(3):290～297 Zhang H J,Liu W L.Multi2res olution decom position and com pression of wave2field extrapolation operator.OGP(in Chinese),2000,35(3):290～297[12]　Beylkin G,M ohlenkam p M J.Numerical operator calculus in higherdimensions.Proceedings o f the National Academy o f Sciences,2002,99(16):10246～10251[13]　Beylkin G,K eiser J M.On the adaptive numerical s olution ofnonlinear partial differential equations in wavelet bases.Journal o fComputational Physics,1997,132:233～259[14]　崔锦泰.小波分析导论.西安:西安交通大学出版社,1995　Cui J T.An Introduction to W avelets(in Chinese).X iπan:X iπanJiaotong University Press,1995[15]　Ingrid Daubechies.T en Lectures on W avelets.United S tates ofAmerica:S ociety for Industrial Applied M athematics Press,1992 [16]　袁修贵,侯木舟.系统工程中矩阵与向量积计算的快速小波方法.系统工程,2002,20(6):94～96 Y uan X G,H ou M Z.The rapid wavelet com puting method of matrix multiplied vector in system engineering.Systems Engineering(inChinese),2002,20(6):94～96[17]　成礼智,郭汉伟.小波离散变换理论及工程实践.北京:清华大学出版社,2005　Cheng L Z,G uo H W.The Engineering Practice of W avelet and375　2期冯德山等:探地雷达小波域三维波动方程偏移Discrete T rans form Theory(in Chinese).Beijing:Tsinghua UniversityPress,2005[18]　Jin S,Wu R S.C omm on2offset pseudo2screen depth m igration.In:SEG69th Annual M eeting.Expanded abstracts:SEG,199911516～1519[19]　Wu Y afei,G eorge A M cM echan.W ave extrapolation in the spatialwavelet domain with application to poststack reverse2time m igration.G eophysics,1998,63(2):589～600[20]　M iller D,Oristaglio M,Beylkin G.A new slant on seism ic imaging:M igration and integral geometry.G eophysics,1987,52(7):943～964[21]　Beylkin G.On the representation of operators in bases of com pactlysupported wavelets.SIAM Journal on Numerical Analysis,1992,29(6):1716～1740[22]　Wu R S,Y ang F,W ang Z,et al.M igration operator com pression bywavelet trans form:beam let m igration.In:SEG67th AnnualM eeting.Expanded abstracts:SEG,199711646～1649[23]　雷文太,粟　毅,黄仕家.探地雷达近场三维距离偏移成像算法.电子与信息学报,2003,25(12):1641～1646 Lei W T,Su Y,Huang S J.G round penetrating radar near field32D range m igration imaging technique.Journal o f Electronics&Information Technology(in Chinese),2003,25(12):1641～1646 [24]　冯德山.基于小波多分辨探地雷达正演及偏移处理研究[博士论文].长沙:中南大学信息物理工程学院,2006 Feng D S.G PR forward simulation and m igration processing research with multi2res olution of wavelets[Ph.D.thesis](in Chinese).Changsha:School of In fo2Physics and G eometrics Engineering,CentralS outh University,2006[15]　金胜汶,陈必远,马在田.三维波动方程有限差分正演方法.地球物理学报,1994,37(6):805～810 Jin S W,Chen B Y,M a Z T.Theree2dimensional wave equation forward m odeling by the finite2difference method.Chinese J.G eophys.(in Chinese),1994,37(6):805～810[26]　赵永辉,吴健生,万明浩等.有限差分法探地雷达波动方程偏移成像.物探化探计算技术,2001,23(1):47～51 Zhao Y H,Wu J S,W an M H,et al.W ave equation m igration image with finite2difference for ground2penetrating putingTechniques for G eophysical and G eochemical Exploration(inChinese),2001,23(1):47～51[27]　马在田.地震成像技术有限差分法偏移.北京:石油工业出版社,1989 M a Z T.Finite2Difference M ethod M igration of Seism ic Image T echnique(in Chinese).Beijing:Petroleum Industry Press,1989 [28]　张海江.小波多分辨波动方程正演模拟与偏移成像[博士论文].北京:中国石油天然气总公司石油勘探开发科学研究院,2000 Zhang H J.W ave equation forward simulation and m igration imaging with multi2res olution of wavelets[Ph.D.thesis](in Chinese).Beijing:Research Institute of Petroleum Exploration and DevelopmentChina National Petroleum C orporation,2000[29]　W illiam A,Schneider J.Analyzing operators of32D imaging s oftwarewith im pulse responses.G eophysics,1999,64(4):1079～1092 [30]　陈树文,刘　洪,李幼铭.三维叠前深度偏移的准三维算法研究.地球物理学进展,2001,16(4):23～28　Chen S W,Liu H,Li Y M.The pseudo232D method operator for32Dprestack depth m igration.Progress in G eophysics(in Chinese),2001,16(4):23～28(本文编辑　汪海英)475地球物理学报(Chinese J.G eophys.)51卷　。