地震波方程人工边界的两种处理方法

第42卷第4期2003年12月石　油　物　探GEOPHYSICAL PROSPECTIN G FOR PETROL EUMVol.42,No.4Dec.,2003文章编号:100021441(2003)0420452204地震波方程人工边界的两种处理方法崔兴福1,2,张关泉2(1.中国石油辽河油田分公司勘探开发研究院计算所,辽宁盘锦124010;2.中国科学院数学与系统科学研究院计算数学与科学工程计算研究所,北京100080)摘要:在分析以往人工边界处理优缺点的基础上,提出了利用波动方程的频散关系式,得到可以吸收任何方向入射波的自适应校正吸收边界条件和旋转校正吸收边界条件。

同时,采用波阵面法和能流密度法判别入射波方向,克服了Pade近似吸收边界只对正入射波具有较好吸收性,而对非正入射的波吸收性不好的缺点。

数值试验结果表明了本方法的有效性。

关键词:自适应校正;旋转校正;波阵面;能流密度中图分类号:P63114 文献标识码:ATwo processing methods for artif icial boundary of seismic w ave equationCui Xingfu1,2,Zhang Guanquan2(puting Center of Exploration&Development Research Institute of CNPC Liaohe Oilfield Com pany,Panjin124010,China;2.Institute of Com putational Mathematics and Scientific/Engineering Computing,Academy of Math2 ematics and System Sciences,Chinese Academy of Sciences,Beijing100080,China)Abstract:In this paper,two absorbing boundary conditions,adaptive correction condition and rotation correction condi2 tion,were derived to absorb incident waves coming from any directions by using dispersion relation,based on an analy2 sis of the advantages and disadvantages of existing boundary conditions.The determination of incident wave direction by wave front and energy flux density was also described.Numerical ex periments were performed and their results were presented,which indicated the efficiency of these methods.K ey w ords:adaptive correction;rotation correction;wave front;energy flux density 实际人工模拟地震勘探是在半无界空间中进行的,而在计算机上进行数值模拟地震波在介质中的传播,只能在有限的模型空间中实现。

地震波到达人工边界将产生虚假的反射波,干扰了实际地震波传播的机理,使仿真剖面变得模糊不清,不利于地下地层构造信息的解释。

自1969年Lysmer和Kuhlemeyer[1]首先提出人工边界处理问题,发展至今,已形成多种人工边界的处理方法,如阻尼边界[1～3]、吸收边界[4～6]、透射边界[7]、移动边界等,用Pade近似法得到的吸收边界条件[4,5]是目前边界处理效果比较好的一种。

我们正是在分析这种边界处理优缺点的基础上,从不同的侧面提出了这种边界的2种校正方法,即自适应校正法和旋转校正法。

1　声波方程自适应校正吸收边界条件1.1　二阶精度吸收边界的推导由二维声波方程u t t=c2(u x x+u zz)(1)进行三维Fourier变换得到频散关系式ω2=c2(k2x+k2z)=k2c2(2)式中,k2=k2x+k2z,k x=k cosθ,k z=k sinθ,θ是波前k和k z的夹角。

引进中间参量a(θ)=1-cosθ1-cosθ2(3)由图1表示的方向波数与波前关系可以推得cosθ2=(α-1)k+k xαk(4)利用k x=k cosθ,cosθ=2cos2θ2-1和2α2-4α+4=11+cosθ及方程(4)得到在频率波数域右行波的边界条件ωkx-1cω2+11+cosθk2z=0(5)收稿日期:20021113;改回日期:20030323。

作者简介:崔兴福(1965—),男,高级工程师,博士在读,现从事地震资料处理方法研究工作。

基金项目:本文得到国家973重点基础研究项目(G199903280)资助。

由Fourier 逆变换得到时间空间域中右行波的边界条件u x t +1c u t t -c1+cos θu zz =0(6)同理可以得到左行波、上行波、下行波的边界条件。

方程(6)中存在一个如何判别入射波方向的问题,下面给出判定波入射方向的第一种方法。

1.2　波入射方向的判定———波阵面法根据波阵面理论,在任何时刻,同一波阵面(波前)上各质点的运动状态相同,并且波阵面的外法线方向一定是波的传播方向。

具体地,为了确定在第(n +1)Δt 时刻边界点M 的位移u (n +1)(M )的值,用第n Δt 和(n -1)Δt 时刻所有各点以及第(n +1)Δt 时刻内部节点的位移值表示,搜索满足如下条件的点m 0。

min m ∈Ω(M)|u n (M )-u n (m )|=|u n (M )-u n (m 0)|(7)可以近似地认为m 0点和边界点M 在同一波阵面上,则二者连线垂直方向可以认为是入射于点M 的波传播方向。

实际应用表明,用这种方法确定入射方向,对平面波比较精确,而对于球面波相对差一点。

2　波动方程旋转校正吸收边界条件文献[8]中对各种近似的边界进行了详细的分析,用Pad e 近似得到的边界条件,具有良好的局部逼近性优点和不好的整体逼近性缺点。

鉴于此,我们通过旋转Pad e 近似边界得到旋转校正边界条件。

2.1　旋转校正吸收边界条件的建立在二维空间中旋转具有如下的形式x ′=x cos θ+y sin θy ′=-x sin θ+y cos θ(8)x =x ′cos θ-y ′sin θy =x ′sin θ+y ′cos θ(9)二阶Pad e 近似的右边界条件为u x t +1c u t t -c2u zz=0(10)方程(10)是方程(6)在θ=0的特殊情况,对正入射的地震波吸收较好,而对非正入射的地震波吸收不好。

将方程(10)变换到频率波数域中,有ωc k x -ωc 2+12k 2z=0(11)旋转方程(11),得ωc(k x cos θ+k z sin θ)-ωc2+12(k z cos θ-k x sin θ)2=0(12)对应到时间空间域得到旋转后的右边界条件-u xt -tan θu tz -1c cos θu t t +c2tan θsin θu x x -c sin θu xz +12c cos θu zz =0(13) 同理可以得到左边界、底边界条件。

用这种方法可以对弹性波的近似边界进行旋转得到弹性波旋转校正边界。

方程(13)中也存在一个如何判别入射波方向的问题,下面给出判定波入射方向的另一种方法。

2.2　入射波方向的确定———能流密度法能流密度矢量的方向和波传播方向相同,在数量上等于单位时间内流经垂直于波传播方向的单位截面的流量。

弹性波的能流密度为I =I 1i +I 2j +I 3k(14)式中,I 1=σ11u 1+σ12u 2+σ13u 3I 2=σ21u 1+σ22u 2+σ23u 3I 3=σ31u 1+σ32u 2+σ33u 3u 1,u 2和u 3为位移向量U 的各个分量。

在二维(i ・k )情况下得到的弹性波能流密度具体表达式如下I =(2μ+λ)5u 15x 1 u 1+λ5u 35x 3 u 3+μ5u 15x 3+　5u 35x 1 u 3i +μ5u 15x 3+5u 35x 1 u 1+(2μ+λ)5u 35x 3 u 3+λ5u 15x 1 u 3k (15)式中,λ和μ为拉梅系数。

声波的能流密度为ω=P U (16)式中,P 为声压。

只要确定出波在某一点能流密度矢量各个分量值的大小,就确定了传播到该点波的传播方向。

2.3　数值计算实现步骤在地震波正演的数值模拟中,不论是采用有限差分法还是有限元法,人工边界处理是关系到正演模拟效果好坏的关键。

在实际处理中,首先是解决如何判别入射波传播方向的问题,即方程(6)和(13)中的入射方向角θ,可以采用方程(7)或方程(15)来求取,方法(7)易于实现,但不十分准确,方法(15)是一种判别入射波方向的较好方法,但需要给出地质模型的拉梅系数λ和μ;采用自适应校正吸收边界方程(6)或旋转校正吸收边界方程(13)进・354・第4期崔兴福等1地震波方程人工边界的两种处理方法行边界处理,前者较为简明,后者烦琐,但效果较好。

下面分别给出方程(6)和方程(13)的差分离散格式。

将波场函数进行差分离散u (x ,z ,t )=u mi ,j ,i =1,2,…,I ;j =1,2,…,J ;m =1,2,…,M ,i ,j 和k 分别表示x ,z 和t 方向的离散节点号,则方程(6)和方程(13)的离散差分格式为u m +1i +1,j -u m +1i ,j -u mi +1,j +u mi ,jΔt Δx+ 1cu m +1i ,j -2u m i ,j +u m -1i ,jΔt 2-c1+cos θu mi ,j +1-2u mi ,j +u mi ,j -1Δz 2=0(17)-u m +1i +1,j -u m +1i ,j -u mi +1,j +u mi ,jΔt Δx -tan θu m +1i ,j +1-u m +1i ,j -u mi ,j +1+u mi ,jΔt Δz-1c cos θu m +1i ,j -2u m i ,j +u m -1i ,j Δt 2+c2tan θsin θu m i +1,j -2u m i ,j +u m i -1,jΔx 2-c sin θu mi +1,j +1-u mi +1,j -u mi ,j +1+u mi ,jΔx Δz +12c cos θu m i ,j +1-2u m i ,j +u m i ,j -1Δz 2=0(18)3　数值计算比较图1是方向波数与波前关系示意图。