当前位置：文档库 › 200509-Kirchhoff偏移法在探地雷达正演图像处理中的应用

200509-Kirchhoff偏移法在探地雷达正演图像处理中的应用

第20卷第3期地球物理学进展Vol .20 N o .32005年9月(页码:849～853)

PROGRE SS I N GEOPHYS ICS

Sep t . 2005

K irchhoff 偏移法在探地雷达正演图像处理中的应用

戴前伟, 冯德山, 何继善

(中南大学信息物理工程学院,长沙410083)

摘要本文首先通过对波动方程的分析,得出了声波波动方程和雷达波波动方程形式一致性,从而说明了把广泛应用在地震数据处理中的偏移技术引入到GPR 资料处理中的可行性;其后说明了时域有限差分法(FDTD )的原理,并用它合成了几种常见的雷达正演剖面;最后利用K irchhoff 积分偏移法对正演所得的雷达剖面进行偏移处理,通过对比偏移处理前后的雷达正演剖面,可知K irchhoff 积分偏移法能使雷达正演剖面中的反射波的归位,绕射波收敛,从而大大地提高了雷达正演剖面的分辨率,更好地指导GPR 剖面的地质解释和验证偏移方法的有效性.关键词克希霍夫(K irchhoff)积分法,探地雷达(GPR ),有限差分法,正演模拟,偏移

中图分类号 P631 文献标识码 A 文章编号 1004-2903(2005)03-0849-05

Th e app lica tion o f K irchh o ff ′s m igra t ion m e th od in th e i m age

pro c e s s in g o f th e grou nd pen e tra t ing radar fo rw ard s i m u la te

DA I Q ian -w e i , FENG D e-shan , HE J i-shan

(School of I nf o-Phys i cs and Geo m atics Engi neeri ng ,C entr al Sout h Un i vers it y ,C hangsha 410083,China)

Ab s tra c t F irstly ,this a rticle ana lysis the w ave equa tion ,and gets the for m a t cons istency of the sound wave equa tion and the radar w ave equa tion.Accord ing to th is ,it illustra tes the fea sib ility of inducting the m igra tion techn ique co m -m on ly used in the se is m ic da ta process ing into the GPR data processing .S econd ly ,th is a rticle exp la ins the princip les of FDTD,and syn thesizes so m e radar for wa rd si m u la te sections w ith it .A t last ,th is a rticle utilizes K irchhoff 's m ethod of integ ra tion m igra tion to do the da ta process ing to the rada r section co m ing fro mthe for w ard s i m ula te .B y co m pa ring the rada r for wa rd si m u la te sections before and a fter the m igra tion process ,it is sho wed tha t K irchhoff 's m ethod of integ ra tion m igra tion cou ld m ake the reflect ion wave return to orig ina l pos ition ,the d iffraction w ave con -vergen t .S o that ,the resolut ion of radar f or wa rd s i m u la te sections cou ld be enhanced h igh ly ,and it cou ld instruct the geology exp lana tion of GPR section and va lida te the va lid ity of m igra tion m ethod better .

K e yw ord s K irchhoff 's m ethod of integra tion m igra tion ,g round penetra ting rada r ,fin ite d ifference ti m e do m a in m ethod ,for wa rd si m u late ,m ig ra tion

收稿日期 2004-10-28; 修回日期 2004-11-10.

基金项目国家自然科学基金重大项目(50099620-03-02)资助.

作者简介戴前伟,男,汉族,1968年生,博士,副教授,就职于中南大学信息物理工程学院,从事工程地球物理勘探工作.

(E -m a il :q wda i @https://www.wendangku.net/doc/2812446728.html,)

0 引言

在地震勘探的反射波法中,各种数据处理理论与方法都已相当成熟,而探地雷达(GPR)剖面观测法与反射地震剖面的单发单道接收具有相似的特性,只是由于二者的观测系统的差异,在方法上略有不同;电磁波与弹性波两者的传播波动方程形式上也非常一致,只是波动方程中代表的物理意义略有不同,如地震中的声学近似方程

Δ2

p -1v 2５2

p ５t

2=0,

(1)

式中p 为压力,v 为纵波波速,t 为地震波在介质中传播的时间.而探地雷达的高频电磁波在二维介质传播的麦克斯韦方程可以写成

Δ2

E =μ0ε*

５2E ５t 2=1v 2５2E

５t 2

(2)

式中

v =

1μ0ε

地球物理学进展

20卷

ε*为复介电常数,它与介电常数、介质的导磁率以

及频率有关.E 为电场强度.μ0为空气中的导磁率.由式(1)、(2)可见,电磁波与弹性波两者的传播波动方程形式上也非常一致,正是由于雷达波与地震波的这种相似性,使得地震资料处理中的许多技术,可用于GPR 资料的处理,如偏移技术、振幅恢复、道间振幅均衡等技术.

文中首先利用有限差分法(FDTD )法正演模拟合成了几种地电模型,并得到了雷达正演剖面.剖面中因有岩性突变点(断层断棱、地层尖灭点、不整合面的突起点等),这些突变点就会成为新震源,再次向四周发出球面波,即绕射波,当然如果存在凹界面,则还会产生特殊的回转波.由于这些异常波的存在,使探地雷达正演剖面不能准确地反映地下介质的真实情况.例如,对于GPR 剖面,来自点状反射体上的波的同相轴为一双曲线形态;来自倾斜界面的反射波的同相轴的倾角比实际反射界面的倾角变得要小.要消除这些影响,必须对GPR 正演剖面作偏移处理.

1 雷达正演原理

M ax w ell 方程组概括了宏观电磁场的基本规律.时域有限差分法是在时域计算电磁场的一种数值方法,自然应该从含时间变量的两个M axw ell 旋度方程出发.于是在无源区域,我们可把M ax w e ll 方程的两个旋度方程表示为如下的形式

Δ×E =-μ５H ５t -σ

m H ,

(3)Δ×H =ε５E

５t

+σe E ,

(4)

其中E 为电场强度(V /m);H 为磁场强度(A /m);ε为介电常数(F/m);μ为磁导率(H/m );σe 为电导率(S /m);σm 为等效磁阻率(ω/m).

在二维T M 波情况下,运用K.S.Yee 氏网格模型,利用中心差商代替微商,把连续变量离散化,把电场进行规约化,得出二维空间的时域有限差分方程,即我们所要推导的探地雷达正演模拟方程

E ～

n+1

(i ,j)=CA(i ,j)?E ～

z (i ,j)

+CD ?C B(i ,j)?H n+1

i +12,j

-H n+1

i -

,j +H n+1

i ,j -12-H

12x

i ,j +

(5)

n+1

i ,j +

n+1

i ,j +

+CD ?[E ～n z

(i ,j)-E ～

n z (i ,j +1)],(6)

n+1

i ,j +

n+1

i ,j +

+CD ?[E ～n z

(i +1,j)-E ～

n z (i ,j)],(7)

其中

CA(i ,j )=1-σ

(i ,j)Δt

2ε(i ,j)

σ(i ,j)Δt 2ε(i ,j),

CB(i ,j)=

ε0

ε(i ,j)+

σ(i ,j)Δt

CD =

Δt Δs ?1ε0μ0

,而Δt 与Δs 为差分时间步长和空间步长.利用上述探地雷达正演方程,结合边界条件,我们得到了两个雷达正演剖面,如图1、图3所示.

2 GPR 偏移的主要分类和应用

偏移技术作为现代数据处理的三大主要技术(叠加、反褶积和偏移)之一,一直是雷达处理研究的一项重要内容,它主要分为两大类:基于射线理论的绕射偏移方法和基于波场延拓的偏移方法.基于射线理论的偏移方法有波前模糊法、绕射曲线叠加法和偏移叠加法等,但是射线理论只是近似描述的电磁波的传播规律,在速度变化平缓的情况下,能使绕

射波、弯曲界面上的回转波、陡倾界面或断层面上的反射波收敛到正确的位置.当在复杂地质条件下反射层归位效果差,同时缺少动力学信息,如振幅.基于波场延拓的方法,如有限差分法、文中用到的建立在波动方程基础上的K irchhoff 积分法、波场外推的F -K 方法、由双程波动方程导出的逆时偏移等偏移方法.这类方法建立在波动方程的基础上的,即使在复杂的地质条件下,偏移归位过程中也可以保持反射波的特征,因此偏移成像精度较高.

3 K irchhoff 积分偏移法

3.1 K irchho ff 积分偏移法的简介

克希霍夫积分偏移法是1883年著名德国学者克希霍夫(K irc hhoff )根据惠更斯—费涅尔原理提出,它建立在物理学的基础上,利用K irchhoff 绕射积分公式把分散在地表各道上来自于同一个绕射点能量收拢在一起,置于地下相应的物理绕射点上,

3期戴前伟,等:K irchhoff偏移法在探地雷达正演图像处理中的应用

K irc hhoff积分法偏移法能适应任意倾斜角度的反

射界面,对剖面网格要求比较灵活.缺点是难以处理

横向速度变化,偏移背景噪声大;确定偏移参数比较

困难,参数选取不合适会大大影响偏移剖面的质量.

3.2K irchho ff积分偏移法的原理

由惠更斯—费涅尔原理可知,如果在围绕场源

所在的某一闭合面Q上,已知位移φ(x,y,z,t)及其

导数,且这些值是连续的(没有奇点),就可以算出Q

面以外任意一点M(x,y,z,t)上场源引起的位移φ.

因此当场源产生的波动已传至Q面上M(x,y,z,t)

点,且不知道场源的情况时,可以使用前一时刻t1,

Q面上的位移[φ]及其导数[φ]′值来表示t时刻Q

内任意M(x,y,z)点的位移φ.

φ(x1,y1,z1,t)=-

4π

[φ]

５

５n

-1

５φ

５n

v r

５r

５n

５φ

５t

d Q,(8)

式中用方括号[]表示位移不是在时刻t,而是在

t1=t-r v

的,[φ]称为延迟位,r表示点M1(x1,y1,z1)至Q面上各点M(x,y,z)的距离,n表示Q面的外法线方向.

3.3K irchho ff积分偏移法的算法研究

在均匀各向同性介质中,电场E满足波动方程

５2E ５2x +

５2E

５2y

５2E

５2z

-k2

５E

５t

=0,(9)

k为传播常数,当不考虑介质对电磁波的吸收,则

k2=1

(v电磁波速度).

围绕观测点p(x p,y p,z p)取一个闭合曲面,该闭合曲面上半面取地表观测面s1,下半面为半径无限大半球面,则p点电场强度E(x p,y p,z p,t)可由地表观测面上的电场强度表示为

E(x p,y p,z p,t)=

2π

５

５n

-1

５r

５n

５

５t

×E x,y,0,t-

d s1,(10)

式中

r=[(x-x p)2+(y-y p)2+z2p]12, n为地表面的法线方向,

５r ５n =-５r

５z

(z垂直地面向下).

当发射天线和接收天线合二为一(自激自收),发射天线发射的电磁波以速度向地下传播,遇到反

射界面时,波被反射回地面,被接收天线接收.根据惠更斯原理,可以把反射界面上的每一点当做二次波源,由这些二次波源辐射出来的二次波在地面各记录点叠加的结果,训是地面所观测到的波场,其定量描述由K irchhoff积分公式给出.对于GPR资料处理为这一逆过程,即由地面上任意点的电磁场记录,要求确定未知的反射界面上的二次波源的位置.

波动方程

ΔE=

５2E

５t2

的解E(x p,y p,z p,t),另

W(x p,y p,z p,t)=E(x p,y p,z p,-t),

W(x p,y p,z p,t)仍然满足同样的波动方程.对于E(x p,y p,z p,t)是时间“向前”的问题,而对于E(x p, y p,z p,-t)则是时间“倒退”的问题.所以在实测偏移处理中,测量数据定义在z=0平面上,地下不均匀点的电场为待求,待求电场在时间上应当较z=0平面上的记录超前r/v,所以将

E x,y,0,t-

改写为

E x,y,0,t+

(10)式变为

E(x p,y p,z p,t)=1

2π

５

５n

５r

５n

５

５t

×E x,y,0,t+

d s1.(11)

3.4K irchho ff积分偏移法的实现步骤

克希霍夫积分法偏移简单过程:

(1)将非零源距的探地雷达记录,转化为零源距的探地雷达记录.

(2)根据反射界面成像原理,偏移速度取真实速度的一半并在t=0时成像.同时,根据GPR在实际中沿一条测线观测,接收天线与发射天线同步移动,每移动一次只能获得一个记录道.当GPR测线的走向与地下地质体的走向近似垂直时,可以将地下地质断面近似看作二维情形.有

E(x p,z p)=

2π

５

５n

-1

５r

５n

５

５t

×E x,

d s1,(12)

(3)令

T=2z p/v,

t0=2r/v,

地球物理学进展20卷

并考虑

５５n =-５５z =-z p

r =-vT/2r ,

５５n 1r

=vT

2r 3

,利用中心差分

５

５t E (x p ,t 0)=E (x ,t 0+Δt)-E (x ,t 0-Δt)2Δt ,其中Δt 为时间采样间隔.上式可以写成

E (x p ,z p )=12π∫

T Δt ?1v 2t 2[E

(x ,t 0+Δt)-E (x ,t 0-Δt)+

2Δt

t 0

E (x ,t 0)]d x ,(13)

(4)使用地面规则网格,若Δx 为点距,Δz 为深度采样间隔,有

x p =m Δx ,z p =l Δz ,x =αΔx ,

r =[(m -α)2Δx 2+l 2Δz 2]

1/2

式中m,α,l 为正整数.上式化为离散形式

E(m,l)=

2πΔtv 2t 20∑a

[E (α,t 0+Δt)-E (α,t 0-Δt)+2Δt t 0

E (α

,t 0)]Δx .

(14)

4 应用实例

4.1 圆状地电断面

应用时域有限差分法对圆状地电断面模拟进行了模拟,并利用了M ur 吸收边界条件对各边界进行了处理.模拟过程中取雷达正演的天线中心频率为60MH z ,周围的介质的介电常数ε1为10.0,电导率σ1为0.005S/m,待模拟的圆形异常体为空洞,空洞中富含水和其它填充物,其介电常数ε2为30.0,,电导率σ2为0.05S/m,取网格的空间步长为0.05m,时间步长为3.0e -10s ,模拟时间为512个时间步长,该圆的圆心距上边界为5.5m,位于水平刻度5.0m 处,半径长为0.5m.其合成的雷达剖面图如图1所示.图中可见该圆的在雷达正演模拟中呈双曲线的形状,深度约5.0m 处为绕射双曲线的弧顶.但是从图中我们可以看出,绕射波的双曲线弧线的幅度比实际模型圆的半径大了将近一倍,这与实际模型大小有很大的出入,这样直接用时域有限差分法模拟出来的雷达正演剖面图,并不能反映出模型的真

实形状和大小.图2是利用K irchhoff 积分法偏移后的剖面,图中可见它较好地消除了绕射波的影响,使反射波正确归位,归位后的圆的大小及位置都与

模拟模型十分吻合.

图1 圆梯状地电模拟模型的雷达正演剖面图F ig.1 The rada r for wa rd co m pose section

m ap of circina l geoe lectricity m

odel

图2 K irchhoff 积分法偏移处理后的

阶圆状模型雷达正演剖面图

F ig.2 The radar for w ard co m pose sect ion m ap of circina l geoe lectricity m odel a fter K irchhoff

in tegra tion m ig ra tion process ing

4.2 阶梯状模型

应用时域有限差分法对阶梯状地电断面模拟进行了模拟,并利用了M ur 超吸收边界条件对各边界进行了处理.模拟过程中取雷达正演的天线中心频率为100MHz ,阶梯状地电模型的上下层介质的介电常数、电导率分别为ε1,σ1,ε2,σ2,数值大小与上面模型中的一致,其雷达正演合成的剖面图如图3所示,图中深处为3.6m 和7.0m 处,清晰可见两层层位,图中水平坐标8.0m 至12.5m 位置处的斜面,为上层断棱右端点处产生绕射波所致.而原阶梯状地电模型中没有斜边,可见未经处理的正演模拟

3期戴前伟,等:K irchhoff 偏移法在探地雷达正演图像处理中的应用

剖面与实际模型的形状是有一定的区别的,它不能很好地反映模型的特性.图4是利用K irchhoff 积分法偏移后的剖面,图中水平坐标8.0m 处为阶梯状模型的垂直断边,断距约3.4m 左右,而原来水平坐标8.0m 至12.5m 位置处的斜面绕射波在此处得到收敛并正确归位.可见经K irchhoff 积分法

偏移后的正演剖面与模拟模型吻合更好.

图3 阶梯状地电模拟模型的

雷达正演剖面图

F ig.3 The rada r for wa rd co m pose section

m ap of ladder geoelectricity m ode

图4 K irchhoff 积分法偏移处理后的阶梯状模型雷达正演剖面图

F ig.4 The radar for w ard co m pose sect ion m ap of ladder geoe lectricity m odel a fter K irchhoff

in tegra tion m ig ra tion process ing

5 小结

(1)雷达偏移是一种将雷达信息进行重排的反演算法,以便将雷达波能量归位到其空间的真实位置,获取地下真实构造图像.作者首先采用时域有限

有差分法对两种地电模型进行了地电合成,然后采用K irchhoff 积分偏移法,对探地雷达正演剖面进行了偏移处理,较好地实现反射波的归位,绕射波收敛,从而更加精确地反映地下介质的真实位置、大小,这可以让我们更好地认识雷达剖面的特征和指导GPR 剖面的地质解释.

(2)对于用K irchhoff 积分法对雷达剖面的偏移处理,虽然能适应任意倾角的反射界面,但仍然存在一些缺点和问题,如干扰噪音背景较强,偏移参数确定比较困难,有待于我们引进一引些更好、更先进的保幅保结构偏移算法,如地震偏移处理中的哈密顿体系下的辛群、李群等算法.参考文献(R eferences):

[1] 冯德山,戴前伟.用V is ua l C ++开发地质雷达正演模拟软件

[J].物探化探计算技术,2004,26(3):231～234.

[2] 冯德山.地质雷达二维时域有限差分正演[D].长沙:中南大

学,2003.

[3] 李成方.偏移处理在探地雷达数据处理中的研究[D ].成都:

成都理工大学,2002.

[4] 薛桂霞,邓世坤,刘秀娟.吸收边界条件在探地雷达偏移处理

中的应用[J].物探化探计算技术,2004,26(3):235～237.

[5] 陈树文,刘洪,李幼铭.三维叠前深度偏移的准三维算法研究

[J].地球物理学进展,2001,16(4):23～28.

[6] 耿建华,黄海贵,马在田.K irchhoff 积分波场延拓基准面静校

正方法研究[J].地球物理学进展,1996,24(6):665～669.

[7] 刘喜武,刘洪.波动方程地震偏移成像方法的现状与进展[J].

地球物理学进展,2002,17(4):582～591.

[8] 黄德济,贺振华,包吉山.地震勘探资料数字处理[M ].北京:

地质出版社,1990.

[9] 李世华,杨有发.物探数据处理[M ].北京:地质出版社,1995.

[10] 陆基孟.地震勘探原理(第三版)[M ].北京:石油大学出版

社,2001.

[11] 李大心.探地雷达方法与应用[M ].北京:地质出版社,1994.[12] 赵永辉,吴健生,万明浩,谭春.有限差分法探地雷达波动方

程偏移成像[J ].物探化探计算技术,2001,23(1):47～51.

[13] 邓世坤.克希霍夫积分偏移法在探地雷达图像处理中的应用

[J].地球科学—中国地质大学学报,1993,18(3):303～308.

[14] 底青云,王妙月.雷达波有限元仿真模拟[J ].地球物理学报,

1999,42(6):818～825.

[15] Jeffrey E.Lucius ,M ichael H.Po wers .GPR da ta process i ng

co mputer soft ware f or the PC [R].USGS ,2002.