基于波动方程有限差分算法的接收函数正演与偏移
傅里叶有限差分法三维波动方程正演模拟
( x , y , z + Δz ;ω) . P′
( 3)
最后 ,将瞬时波场 P ( x , y , z + Δz ;ω) 代入 ( 4 ) 式进 行有限差分校正 , 9 P ( x , y , z + Δz ;ω) 9z
1856
2 2
地 球 物 理 学 报 ( Chinese J . Geophys. )
[22 ] [19 ] [13 ,14 ]
可以进一步改进相应的算法 . 基于 ADIPI 格式的 FFD 法既消除了分裂误差 , 又不影响 x 方向和 y 方向的精度 , 而且保持了双向 分裂 FFD 法的高效性 , 因此它可以胜任高陡倾角 、 强横向变化介质的三维正演模拟 ,成为精度高 、 效率 高、 适应性强的三维单程波正演方法 .
2 基于 ADIPI 格式的 FFD 方法
211 三维 FFD 算子
横向非均匀介质中单程波方程的三维 FFD 算 子为
[16 ,17 ]
9 9 2 + 2 9 x 9 y 0 Δs + kz ≈ k z + ω 2 2 9 9 1 + a 2 + 2 9x 9y
b
2
2
,
( 1)
其中 ω 是角频率 , k z =
Zhang J H ,Wang W M , Zhao L F ,et al. Modeling 3 2D scalar waves using the Fourier finite2difference method. Chinese J . Geophys . ( in Chinese) , 2007 , 50 (6) :1854～1862
1 引 言
波动方程地震偏移成像方法的现状与进展-刘喜武
第17卷 第4期 地 球 物 理 学 进 展 Vol.17 No.4 2002年12月(582~591) PROGRESS IN GEOPHYSICS Dec.2002波动方程地震偏移成像方法的现状与进展刘喜武1,2 刘 洪1(1.中国科学院地质与地球物理研究所,北京100101; 2.中国科学院研究生院,北京100039)[摘 要] 综述了波动方程地震偏移成像方法研究的现状,指出了各种算法的问题和解决方案;给出了相空间小波分析步进算法偏移成像以及保幅保结构李群算法偏移成像的最新进展;对波动方程偏移成像发展趋势进行展望.[关键词] 波动方程地震偏移成像;算法;相空间;小波分析;李群;研究进展[中图分类号] P315 [文献标识码] A [文章编号] 1004 2903(2002)04 0582 100 引 言地震偏移是一种将地震信息进行重排的反演运算,以便使地震波能量归位到其空间的真实位置,获取地下真实构造图像[1].除了深度域构造成像外,地震偏移还为其它特殊处理提供振幅、相位等信息,用于速度估计和属性分析,建立在波动方程基础上的地震偏移成像技术代表了地震处理的极致[2].地震偏移最初是在水平迭加基础上进行的,目的是使倾斜界面共深度 映像聚焦,使绕射波归位,即将能量还原到它们正确位置上.早期人工偏移是按照偏移空间的时距关系作图;若将共深度点剖面看作一系列绕射点组成的源反射,可用计算机实现对这些绕射点的偏移,即建立在射线理论基础上的绕射扫描迭加方法以及后来的Kirchhoff偏移.20世纪70年代初美国斯坦福大学以J.F.Claerbout为首的SEP研究小组第一个对标量波动方程提出了有限差分近似解法[3,4],实现了地震偏移.此后建立在波动方程基础上的地震偏移成像方法如有限差分法、Kirchhoff积分法、F!K方法及其各种变形等方法广泛应用[5].爆炸反射面模型[6](Loewenthal et al.,1976)为波动方程偏移成像条件的建立奠定了理论基础.由于波动方程描述地震波地下传播规律,因而波动方程偏移一方面可以解决复杂介质条件下成像问题,另一方面保持了波场的动力学特征.地震偏移各种方法最初是作为时间偏移方法出现的,目的是满足二维时间叠加剖面成图需要,后来为满足横向变速情况下成像精度需要,发展了深度偏移方法[7,8](Hubral1977, Larner1981).近20年来,偏移方法又发展到了三维和叠前偏移[9],三维叠前深度偏移代表地震偏移的发展水平.当今各种各样的偏移技术方法极为丰富,如时间偏移、深度偏移、二维、三维、迭前、迭后;如使用共炮集、共方位角道集、面炮方法实现等等.偏移算法也多种多样.实际应用中根据具体情况和要求选取相适应的方法.考虑如此复杂而庞大的内容,本文只对波动方程地震偏移中核心算法!!!波场延拓和成像的现状与进展进行阐述,并对其发展趋势进行展望.[收稿日期] 2002 01 18; [修回日期] 2002 04 04.[基金项目] 中国科学院知识创新工程项目(KZCX L SW 18)资助.[作者简介] 刘喜武,1971年生,黑龙江讷河人,中国科学院地质与地球物理研究所博士研究生,主要从事地震信号处理和成像研究.1 波动方程地震偏移成像方法的研究现状波动方程地震偏移成像的各种方法都是建立在波场反向外推基础上,按照算法实现的原理可以分为两大类:基于射线理论的偏移方法和基于波场延拓的偏移方法[10].基于射线的Kirchhoff 积分类方法,依靠射线追踪获得成像所需的旅行时,不受反射界面顷角限制,计算效率高,灵活,但在复杂地质条件下,多值走时使射线追踪难于获得正确旅行时,导致成像效果较差.此外,基于射线的方法缺少动力学信息如振幅等.基于波场延拓的方法,如有限差分法、F !K 方法等,物理概念清晰,自然解决了多值走时问题,能够更为精确成像.这类方法包括由双程波动方程导出的逆时偏移,由单程波动方程导出的各种方法.实质上,建立在波动方程基础上的Kirchhoff 积分方法与波场外推的F !K 方法、有限差分方法数理基础相同[5].已证明在常速介质中Kirchhoff 积分方法与F !K 方法的波场外推公式完全等价,而有限差分方法使用波动方程的各种近似,其波场外推公式除相位精度外,形式与前两者基本相同,是它们的近似式.然而,由于波场反向外推的实现算法不同,导致它们各自不同的特点.1.1 Kirchhoff 积分类偏移方法Kirchhoff 积分法波动方程偏移建立在波动方程Kirchhoff 积分解[11]的基础上,把Kirch hoff 积分中的格林函数用它的高频渐进解(即射线理论解)来代替.其基本过程包括从震源和接收点同时向成像点进行射线追踪,然后按照相应的走时从地震记录中拾取子波并进行叠加,如果对所有的路径计算出的走时都正确,对所有记录数据的叠加结果会在某些部位产生相对较大的值,这些值给出地下界面(即反射体位置).Schneider(1978)建立了Kirchhoff 积分偏移的波动方程理论基础[11],Bleistein(1987)将Kirchhoff 积分方法拓展到求解反射系数,进一步推进到偏移后的参数估计[12].虽然Kirchhoff 积分公式是严格的波动方程分解,但它的实现是利用波动方程的零阶高频渐近近似(射线方程),这种近似只有在t 时刻的圆频率 较大时才合理,因而源点或接收点的几个波长以内的绕射点不能正确成像,存在焦散区.其次,在复杂介质中,由于速度不均匀和高频近似,绕射点与源点和接收点之间的传播距离要远大于几个波长,这种大距离反向外推波场,就存在多重路径问题,造成旅行时求取困难,近年来人们不断改进旅行时求取方法,Audebert e t al (1997)对这些方法进行了总结和对比[13].由于Kirchhoff 积分法成像时对绕射面扫描到的数据,没有考虑频率成份,绕射面较陡部分在它扫过未偏移数据平坦部分时,将对地震子波进行重采样,造成算子假频.Gra y (1992)和Lumley et al.(1994)提出减少绕射面段陡部分扫描到子波频率成份,克服了这个问题[15!17].归纳起来,Kirchhoff 积分法偏移成像局限[18]为:(1)分辨率随着深度的增加而变差,从而导致对深部结构分辨率降低,这一现象源于利用射线解近似格林函数时对菲涅尔带的影响;(2)成像信息中缺乏正确的振幅信息,这一现象源于射线近似在复杂介质中存在焦散、多重路径和干涉等问题.为了保持灵活性,同时提高Kirchhoff 积分偏移的精度,Hill(1990,2001)提出了高斯束偏移方法[19,20].该方法将源点和接收点的波场局部分解为 束 导回地下.几个束可能来自不同地面位置,且不同束指向不同的初始传播方向.每一束由各自射线管引导,独立于其他束∀583∀4期 刘喜武,等:波动方程地震偏移成像方法的现状与进展传播.射线管可以重叠,这样能量可以通过多个路径在像点位置与源点和接收点之间旅行,解决了多路径问题.Bevc(1997)也提出一种解决多路径方法[21].该方法先向下应用标准的非递归Kirchhoff 偏移方法将记录偏移到地表以下多个波长深度处.在这个深度范围内认为多值走时不很严重,Kirc hhoff 积分方法可以较精确进行.接着在该深度处采用Kirc hhoff 基准面方法[22](Berryhill,1984)计算一个向下延拓的波场.然后将这个波场用于下一个有限范围深度内Kirchhoff 偏移.经过这些Kirchhoff 偏移和向下延拓的结合,完成偏移.这种方法通过分级办法处理多路径问题,但只对二维实用[22].1.2 基于波场延拓的波动方程偏移方法Kirchhoff 积分偏移采用地表所有记录数据单一整体空间褶积计算每一个点的像,并且地表位置与像点位置之间通常只采用一个路径,基于波场外推的波动方程偏移方法递推地从前一个深度Z 的波场计算深度Z + Z 的波场,自然考虑到每个深度可能的绕射点与每个源点或接收点的多重路径,因而可以在较为复杂介质条件精确成像[23].随着偏移技术的发展,波场外推偏移成像从最早单一的T !X 域实现,发展为诸如T !K 、F !K 、F !X 、 !P 域、小波变换域等多种域实现偏移方法.偏移方法包括由双程波方程导出偏移方法,如T !X 域有限差分法,逆时偏移方法;以及基于单程波方程近似解的各种方法.1.2.1 T !X 域或波数域有限差分方法(F !D)有限差分法波动方程偏移是最早提出的一种波场延拓波动方程偏移方法[3!5,25].这种方法直接对T !X 波动方程进行坐标变换并略去二阶导数项,得到变换后的简化的波动方程,然后再利用有限差分方法求解波动方程,进行成像.当地面已知波场P(x ,0,t)向地面下延拓至反射点Q(x #,z #)时,波场函数P (x #,z #,t #)的旅行时t #为波从震源到反射点Q (x #,z #)的下行传播时间[23].F !X 域有限差分偏移方法存在一些固有的困难[18]:∃由于空间离散化造成的数值频散,导致不同频率的波以不同的速度传播,从而造成成像的误差和人为假象;%由于矩形网格划分导致的三维数值各向异性,造成沿不同方向波传播速度不同.用差分方程近似替代偏微分方程,可分为隐格式和显格式两大类.隐格式求解困难,但精度高,对倾角较大情况偏移效果好.隐格式有限差分偏移自Claerbout(1971)提出15&有限差分偏移以来得到充分发展,马在田(1982)和张关泉(1986)分别提出高阶方程降阶方法,有效地解决了隐式有限差分方法难以高角度成像问题,但是没有考虑到降低计算量的问题[35,37].李志明的双线性变换和三维隐式有限差分多方向分裂算法中各项异性的补偿是隐式有限差分偏移方法的出色工作[27].显格求解较容易,快速,但精度较低,对大倾角适应性差.显式有限差分法由Berkout 提出[26],Holberg 和Hale(1999)在二维情况下有效解决了该算法的稳定性[2],Blacquiere 进一步将其推广到三维[2].对近似方程中x 做Fourier 变换到波数域,可进一步得到近似方程,再利用有限差分求解.基于波场外推的波动方程偏移方法(尤其是F !D 方法)与Kirchhoff 积分偏移方法在偏移孔径处理上存在显著差异[2].Kirchhoff 偏移方法通过绕射曲线对输入道叠加或将输入采样散布到输出道.波场外推偏移从记录地表将整个波场向下延拓,如果需要较大孔径,这个∀584∀ 地 球 物 理 学 进 展 17卷孔径必须作为一部分包含在整个计算内,换句话说即使记录地表这些道位置处没有能量,在所有深度偏移运算也必须包含所有的输出道.这就导致额外巨大计算量,在三维叠前深度偏移中,难于经济适用.1.2.2 逆时偏移(Reverse time migration)逆时偏移也是一种应用有限差分求解波动方程实现波场延拓的方法,它不是深度域外推,而是进行时间外推,求解双程波声波或弹性波方程,并且允许波向各个方向传播.这种方法没有倾角限制,精度较高.它的计算方法正好与地震正演模拟的计算顺序相反,以最大时间开始向最小时间计算.Baysal et al (1983)和Mcmechan(1983)给出这种方法的详细内容并描述了高精度全倾角成像能力[24,25].之后许多学者提出了不同的差分格式,不同计算方法实现逆时偏移方法,如双线性变换逆时偏移[27](李志明,1991).逆时偏移法完全遵循全波波动方程,是最为精确的方法,但同有限差分模拟一样,存在稳定性和数值频散问题,同时计算量非常大.1.2.3 基于单程波方程的偏移方法从理论上讲,波场延拓应该用波动方程的边值问题来解决,但由于地面上的观测范围有限,这种边值问题是不稳定的,为此,地球物理学家提出用单程波方程来做波场延拓.通常的叠后偏移(也就是零炮检距剖面偏移)就要用到单程波方程波场延拓.叠前深度偏移的共炮集偏移,是通过炮点波场和接收点波场向下延拓实现的.在叠后偏移中的成像运算是取 零时刻的波场值 .炮集叠前深度偏移的成像运算是炮点波场与接收点波场做互相关.基于单程波方程的方法有 !K 域stolt 法、 !K 域相移法(phase shift)和相移加插值法(phase shift plus interpolation), !X 域有限差分算法(finite difference)、 !K 域与 !X 域交替的裂步Fourier 算法(split step Fourier), !K 域与 !X 域Fourier 有限差分算法(Fourier finite difference), !K 域与 !X 域广义相位屏(generalized screen propagator)方法.以二维声波方程为例,这些算法的共同基础[10]是k 2x p (k x ,z , )- 2p (k x ,z , ) z 2= 2v 2p (k x ,z , ).(1)其中P (k x ,z , )为压力波场,z 为深度, 是圆周率,v (x ,z )是速度函数,k x 是x 方向波数,进一步对z 做Fourier 变换得到频散关系k 2x +k 2z = 2v 2.(2)其中k z 是z 方向波数,相应的单程波方程为p z=∋ik z p .(3)其中∋号分别表示上行波和下行波方程.由解析解可以得出 K 域波场延拓公式为p (k x ,z + z , )=p (k x ,z , )e∋ik z z .(4)(1)Stolt 偏移算法[28]Stolt(1978)提出的 !K 域偏算法.该算法不采用波场延拓途径而是通过已知叠加数据p (x ,t,z =0)的谱p ((k x , ,z =0)一步求出偏移剖面p (x ,t =0,z )的谱p ((k x ,t =0,k z ),再经过二维Fourier 反变换形成偏移剖面.该算法运算极快,效率高,无倾角限制和频散.该算法∀585∀4期 刘喜武,等:波动方程地震偏移成像方法的现状与进展要求每次外推时全部使用一个平均速度,不能适应横向变速,是一个时间偏移方法.(2)相移法(PS)和相加插值算法(PSPI)[29!31]Gazdag(1978)提出相位移方法.(3)式是其基本公式.该方法无倾角限制,无频散,但由于在 !k x 域中延拓,k x 与V(x )关系无法确定,因而不能适应横向变速情况.Gazdag(1984)提出相移加插值算法克服这个问题,在每一个深度用多个不变的参考速度计算几个外推波场,再变换到 !x 域,真正波场根据V(x )与参数速度的关系通过线性插值得到.该方法纯粹是一种数学手段,无物理本质意义,计算量较大,且不能完全改善对横向变速的适应.(3)裂步Fourier 方法(SSF)[32,33]针对相移法不能适应横向变速情况,Stoffa et al(1990)提出一种分裂步Fourier 方法.该方法基于小扰的理论,将速度场分为背景速度和扰动项之和,推导出波场延拓公式,交替在 !k x 域和 !x 域进行,实现对横向速度变化的处理.该方法无倾角限制和频散现象.对强烈的横向变速情况,扰动理论不成立.为解决这一问题,类似于相移加插值思想,Kessinger(1992)将多参考速度对数(MRVL)引入SSF 方法,对两参考速度之间重叠区进行插值,但增加了较大计算量.(4) !X 域有限差分方法(FD)[34!37]也是针对(3)式不能处理横向变速,想办法将其变换到 !x 域通过求解差分方程,实现波场延拓.这种方法的最大好处是可以处理任意方向任意变化的速度.缺点是特别陡倾角的轴无法成像.但通过对平方根算子使用高阶近似(Glaerbout,1985)可以改善倾角限制.马在田(1982)提出的平方根算子连分式展开,将高阶近似方程分裂为多个低阶方程,在 !x 域用有限差分方法求解这些低阶方程实现波场延拓,从而进一步解决了大倾角有限差分偏移问题.(5)Fourier 有限差分方法(FFD)[38,39]Ristow,et al (1994)提出结合相移法和有限差分方法各自优点在双域操作的一种混合算法.通过加入有限差分校正项对Stoffa et al (1990)的SSP 方法进行一般化并进行改善,以适应剧烈横向变速.该方法将速度分裂为常速的背景场和横向变化速度场,将向下延拓算子分裂为适于常速度的相位移算子和适于变速度的有限差分算子.单程波方程双域(波数!空间)解法,首先将问题变换到频率域,然后利用快速Fourier 变换(FF T)在空间域和波数域之间往返转换,实现双域操作.波的自由传播是在波数域通过一个具有某种参考速度的均匀介质进行的,参考速度可以纵向变化;而波和介质非均匀部分的相互作用是在空间域中依靠修改波前来进行的.双域操作的方法没有离散化造成的格点频散,计算量比时间!空间域的有限差分小得多.(6)广义屏方法(GSP)[40!42]广义屏方法(GSP)同FFD 方法一样,是一个结合双域操作得混合算法,两者的思想和做法一致.目前GSP 传播算子主要基于两条途径导出:一是基于散射理论和Born 近似[43],另一途径是基于Hamilton 路径积分[44].对双域传播算子的改进提高了大角度波响应得精度,并相应给出这些算子的命名,如拟屏传播算子(Pesudo screen propagator),复屏传播算子(C omple x screen propagator ),加窗传播算子(Windowed screen propagator ),Pade 屏传播算子(Pade screen propagator),高阶广义屏传播算子(High order screen propagator).这些算子在实际∀586∀ 地 球 物 理 学 进 展 17卷应用中得到很大得发展,能够在非常复杂的地区成像.三维波动方程叠前深度偏移的基础是三维波场延拓.在将二维单程波方程波场延拓推广到三维时,分裂步方法和相移加插值方法都没有什么限制,FD 和FFD 则遇到矩阵求逆的困难.提高大角度成像效果,人们探索了两个方向,一个是绕开FD 和FFD,找出快速实现的显格式短算子;另一个是探索矩阵求逆的快速算法设法实现三维FD 和FFD.在短算子构造方面张剑峰提出一种优化方法[45],可以较快构造出二维短算子,然后引入Mcclellan 变换进行一个二维到三维的变换,求出三维短算子.Claer bout1997年提出了螺旋坐标系谱分解[46]的方法,把多维数据体看成一维数据体,引入螺旋坐标系使得拉普拉斯算子的矩阵表述具有良好性质,采用谱因式分解方法可以比较快地求出矩阵的LU(下三角和上三角)分解,从而找到一种至少在均匀介质中最准确的求逆方法.Ris to w 和Ruhl 则提出用多方向分裂近似求算子的逆[39].1.3 基于射线和基于波场延拓的波动方程偏移方法的结合有学者结合射线追踪和有限差分各自的优点提出射线+有限差分偏移方法[47].该方法利用射线追踪速度快的优点,替代炮记录上波动方程偏移中源子波的正向传播过程,确定地震子波在各个时刻的位置,又利用有限差分法保持振幅的优点对原始记录作反向传播成像.2 波动方程地震偏移成像方法研究的进展现代地震偏移技术的发展已从叠后偏移到叠前偏移,从二维偏移到三维偏移,从时间偏移到深度偏移,但就偏移成像的核心技术仍是波场延拓和速度估计.在围绕继续解决陡倾角和变速问题,构造稳定、快速、精确算法的同时,保结构和保幅波动方程偏移算法成为重点研究方向.在人们致力于对传统波场表示(傅立叶变换)的波场延拓算子进行分析同时,数学和物理领域新兴的理论和工具为波动方程偏移提供了新的思路和方法,如前文提到的高斯束方法,以及小波分析和现代群论逐渐被应用于地震波场研究中,取得重要进展.2.1 相空间小波分析地震波场偏移成像的步进算法[48!53]研究波传播与成像的常用方法是利用特定的基函数(如Fourier 谐波,Green 函数)等对波场作分解.这条途径是否有效其关键在于:变换后的基函数所满足的方程是否易解,边界条件和初始条件是否容易处理.若对波场的展开能满足这个条件,称这种展开为有效表示,例如求均匀介质中的Helmholz 方程解时,波场对u(x ,z )对x 做Fourier 展开,这时在变换域u )(k x ,z )满足的方程易于求解.但是在非均匀介质中,u )(k x ,z )所满足的方程与u(x ,z )所满足的方程一样难于求解.在解决非均匀介质中波传播问题时,Fourier 变换之所以失效,其根源是Fourier 变换中基函数没有局部化性质.应用具有良好时频局部化性质的基函数对波场进行展开,将对波场求解有益.Weyle Heisenberg 框架(W !H 框架)、短时Fourier 变换(WFT)、小波变换(WT )等近年来发展起来的数学分析方法,由于其基函数具有时空局部化特性,成为函数和算子表示的有力工具.Steinberg(1993)在频率!空间域从标量波动方程出发研究了非均匀缓变介质中波场的步进算法,使用短时Fourier 变换和扰动方法,通过对非均匀介质作部分均匀近似的途径,推导出相空间地震波传播算子,给出步进算法,但计算量较大.在理论上,Steinberg(1993)提出∀587∀4期 刘喜武,等:波动方程地震偏移成像方法的现状与进展的方法开辟了研究非均匀介质中波传播的新的重要途径.高静怀(2000)对Steinberg 步进算法进行改进,提出基于W H 框架的波场步进算法,并将这种算法用于地震成像,减少了计算量.吴如山等(2000)将非均匀介质分解为局部背景介质和局部扰动介质,推导出了离散表示的局部化传播算子,将波场分解为一个紧标架(Garbor Daubechies 标架)或局部余弦基的表示,并且均有快速算法.用于叠后SEG/EAE G 盐丘模型成像,取得好的效果.相空间小波分析类地震偏移试图用具有局部化特征的小波类的基函数取代传统的Fourier 变换基函数,从根源上解决非均匀介质波传播和成像问题,值得深入研究和探讨.2.2 哈密顿体系下地震波场延拓的辛群算法和李群算法[54!62]在地震波传播理论和偏移成像研究中,应用辛群和李群算法是中国科学院地质与地球物理研究所以刘洪为首的研究小组探索的一个方向,获得了比较准确、快速的二维和三维波场延拓算法.牛顿体系、拉格朗体系与哈密顿体系是力学的三个等价的体系.牛顿体系将运动表示为二阶微分方程;拉格朗体系将力学体系表示为一个变分极值问题;哈密顿体系则给出在相空间(位置和动量)的一阶微分方程并指出相空间物体运动满足辛几何.当将它们表示成数值算法时,计算过程和结果是不等价的,有限差分方法是牛顿力学体系表述下的算法,有限元是拉格朗日体系表述下的算法,哈密顿表述下的算法是辛几何算法.辛几何就是体积不变的几何,而欧氏几何是长度不变.冯康院士在研究拉格朗日体系后,独立地发现了有限元方法.冯康又进一步发现,当把连续的Hamilton 体系离散化时,存在一种离散化方法,可以得到离散的Hamilton 系统,这种系统的运动方程满足辛几何性质,这种离散化的方法,就是辛几何方法.地震波传播过程从本质上说是一个具有无穷维自由度的哈密顿体系随时间演化过程,系统在任意时刻状态由初始状态通过单参数辛变换群变换得到.为正确计算地震波传播和偏移成像过程必须采用保持系统基本性质的算法.对(3)式中平方根算子的不同处理方法是波动方程偏移不同算法的基本原理.对于横向均匀介质在 !K 域可对单程波方程算子e 指数进行准确计算,为处理横向变速介质,相继开展了对单程波方程平方根算子的不同近似的方法研究.单程波深度延拓算子是由辛算子不对称分解得到的李群算子,将其表示成一个辛算子与一个窗算子的乘积,用辛几何方法来改进辛算子数值计算精度,用短算子非平稳滤波来加快窗算子计算.已有方法(如FD 、FFD)中辛算子的原有计算方法与冯康辛几何方法中的一阶Pade 近似相同.将冯康辛几何方法中二阶Pade 近似用于该算子,就导致了一种新方法,称为二阶辛几何方法(Yang Hui,2001).辛算法对单程波方程算子中平方根算子和e 指数提出了高阶辛Pade 近似辛几何格式,同其它非辛方法相比提高了成像精度,而且比Claerbout 用以提高精度的两步法要精确,降低了频散,适应横向变速,具有保幅特点.用短算子非平稳滤波来加快窗算子计算,能够根据速度的变化最佳切除消逝区干扰波,比SSP 保留更多有用信息.Li Bing(2001)提出一种优化方法,可以较快构造出二维短算子,然后进行一个二维到三维的Mcclellan 变换,求出三维短算子.Luo Mingqiu(2001)将二阶辛几何方法用分裂法推广到三维,并用OVER THRUST 模型(即法国石油研究院三维MARMOUSI 模型)作了检验.进一步研究表明,单程波方程波场随深度延拓过程,实际上是波场在一个以Z 为单参数的单参数李群作用下的变换过程.保持波场演化的李群变换性质的离散化算法为李群算∀588∀ 地 球 物 理 学 进 展 17卷。
有限差分波动方程正演模拟中的吸收边界条件
有限差分波动方程正演模拟中的吸收边界条件王开燕1，周妍1，刘丹1，郝菲2【摘要】在地震波传播的数值模拟过程中，在有限的区域内建立吸收边界条件是一个很重要的问题。
主要运用有限差分的方法对二维声波方程进行正演模拟，介绍并分析了利用有限差分的方法进行波动方程正演模拟过程中的几种吸收边界条件。
先通过理论阐述，然后通过建立均质模型和层状介质模型来研究不同吸收边界条件下的边界吸收效果，得到对应的波场快照和单炮记录，并加以比较。
通过实际验证得知当运用完全匹配层(PML)吸收边界条件时吸收效果最好，基本上不产生虚假反射。
【期刊名称】当代化工【年(卷),期】2014(000)005【总页数】4【关键词】关键词：有限差分法；正演模拟；吸收边界条件；二维声波方程；虚假反射模拟与计算地震数值模拟是地震勘探和地震学的重要基础，并已经在地震勘探和天然地震勘探中得到广泛的应用。
地震勘探过程中，我们只能得到地表和地下很少部分的数据，不可能得到波场的全部信息，只能通过波场正演模拟来获得波场的全部信息，从而全面地反映地震波在地下介质中的分布与传播情况。
地震数值模拟[1]是在已知地下介质结构情况下，研究地震波在地下各种介质中传播规律的一种地震模拟方法，其理论基础就是表征地震波在地下各种介质中传播的地震波传播理论。
本文主要采用有限差分的方法进行正演模拟[2]，但实际地震波是在无限介质中传播的，由于受计算机内存和计算时间的限制，有限差分法只能得到有限数量网格点上的波场值，所有就必须截断计算空间并设置边界条件，得到有限的计算模型，所以边界吸收条件就非常重要，如果处理不好就会产生虚假反射，影响得到的结论。
近年来，国内外许多学者在吸收边界条件方面做了大量的工作，提出了各种边界条件[3-6]。
本文通过声波方程有限差分方法，验证不同吸收边界条件下的正演模拟效果，优选出效果好的吸收边界条件。
1 二维声波方程二阶精度有限差分算法二维声波波动方程的表达式为:其中：c—声波波速；u（x，z，t）—声波波场值；f（x，z，t）—震源项。
波动方程叠前偏移与波形反演研究.
波动方程叠前偏移与波形反演研究【中文摘要】本文就当前反射地震勘探的两个热门题目——波动方程叠前深度偏移与基于波动方程的全波形反演进行了多方面的探讨。
深度偏移的核心是波场延拓算子的求取,我们编程实现了目前使用较多的几种算子:相移、相移插值、分步傅立叶、傅立叶有限差分以及广义屏算子。
数值例子说明,这些算子均有各自的优点和局限,应该根据介质的特点选择成像精度与计算时间可接受的方法。
共炮记录的偏移是最直观、易实现、精度高的叠前偏移技术,然而其计算量很大。
近些年出现了旨在进步叠前偏移计算效率的多种非共炮记录的偏移方法。
其中双平方根偏移基于沉降观测概念,将震源和记录同时向下延拓,计算迅速,而平面波偏移则先对震源和记录做平面波分解,再偏移叠加,是一种有效降低多次覆盖地震记录冗余性的好方法。
这两种方法均可为速度分析等输出重要的共成像道集。
为给叠前偏移实验提供必要的输进,我们基于惠更斯原理天生了共炮记录,与其它正演方法相比,该方法计算效率高、各种绕射波完整清楚、信噪比高,可以灵活的选择波场延拓算子,并且对观测系统有很好的适应能力,能够方便处理地表起伏时的记录正演题目。
起伏地表的处理是困扰复杂地区成像质量的一个关键因素,本文鉴戒Reshef“逐步—累加”的波场外推概念,在起伏地表上充填常速度,变不平坦地形为平坦地形,实现了共炮记录的天生和偏移,结果表明,我们的方法是处理复杂地表的简单有效策略。
偏移是构造成像的基本手段,然而随着油气勘探的不断深进,人们越来越希看由地震资料获得更多的岩石物性信息,其中高精度的介质速度场亦是偏移成像必要的输进,所以,各种地震反演方法迅速发展起来。
由于基于波动方程的波形反演方法直接采用微分方程模型,能够充分利用地震信息,可以得到更加精确可靠的反演结果,近些年景为国内外研究的热门。
然而,传统的梯度法反演存在收敛速度慢、不易重建速度场低波数成分的严重缺点,而高斯—牛顿法的计算量又很大。
我们采用近些年国外学者提出的将虚震源法与互易定理相结合计算雅可比矩阵的方法实现了高斯—牛顿法波形反演,并将结果与梯度法做了比较,深进熟悉了波场反演的规律。
面向波场分析的波动方程有限差分正演方法
面向波场分析的波动方程有限差分正演方法随着勘探对象的日益复杂化,使得研究的模型介质也越来越复杂,为了更详细地了解和认识复杂介质波场,进行地震波场数值模拟是非常必要的。
波动方程有限差分法是近几年最流行和最受欢迎的数值模拟方法之一,其优点就是模拟波场信息丰富、实现简单和运算速度较快等。
波动方程有限差分法虽然有很多优点,但是也存在一些不足之处,其中稳定性条件就是一个很重要且急需解决的问题。
稳定性条件直接决定正演模拟结果的成败,在粘弹介质波动方程的有限差分解中,介质的粘弹特性会增加数值解的不稳定性。
本文从粘弹性理论入手,研究了Kelvin-Voigt和Maxwell在不同精度差分格式下的稳定性条件,而对于标准线性粘弹模型,只给出了稳定性条件与速度和品质因子的数值解。
当空间网格步长一定时正演模拟中的时间步长是由模型中的最大速度决定的,速度越大时间步长越小,所需的运算量也越大,为了在不影响结果的前提下提高波动方程有限差分法模拟的效率,本文提出了波动方程变时间步长有限差分数值模拟。
近地表的主要特点就是低速度、低密度和低品质因子。
随着山地、沙漠等地表复杂区地震勘探的发展,近地表地震数值模拟技术受到了地球物理学家的广泛关注。
结合上述得到的结果进行地震波场近地表效应的模拟及分析,对这一部分的模拟采用自由边界条件,对近地表的瑞雷面波进行了模拟和定量分析,确定模拟出的是瑞雷面波。
粘弹性介质是更符合实际介质,进行了含有低速层的粘弹性介质在自由边界条件下的数值模拟。
还从Q值大小、震源子波频率等方面定量讨论粘弹介质地震波衰减影响因素。
实际介质(和油气藏有关的介质)往往是非均匀的,非均匀介质中存在很多微小异常且分布极不规则(通过测井数据和岩心样品就可以发现)。
对于上述情况实现了随机介质的构造和模拟,得到非均匀介质中的波场。
最后对几个实际速度模型进行数值模拟,得到了复杂储层(含有孔洞、裂缝等)下的全波场和储层在不同含水饱和度下的波场,并且通过频谱和时频分析对结果进行了分析,还研究了低速层和近地表对波场模拟的影响。
波动方程叠前深度偏移原理DOC
3 波动方程法叠前深度偏移常用的叠前深度偏移方法包括射线法和波动方程法。
射线法主要指基于绕射旅行时计算的Kirchhoff 积分法，在绕射旅行时计算方法上可以采用基于函程方程的变速射线追踪法、基于费马原理的二维有限差分法和稳健高效的三维迎风有限差分法；而波动方程叠前深度偏移是复杂介质成像的有效手段，能够解决强横向变速条件下复杂地质体的地震波成像问题。
基于共炮集的波动方程叠前深度偏移的基本思路是，首先对每一炮进行单炮偏移成像，然后再把各炮成像结果在对应地下位置上叠加，从而得到整个剖面成像。
从计算角度而言，成像过程是很简单的步骤，波场外推算子决定了偏移方法的效率、成像精度及其适应范围。
一般要求偏移算子能够适应陡倾角反射的成像及剧烈的横向速度变化，同时具有较高的计算效率。
3.1 波动方程叠前深度偏移的基本思路基于共炮集的波动方程叠前深度偏移的基本思路是，首先对每一炮进行单炮偏移成像，然后再把各炮成像结果在对应地下位置上叠加，从而得到整个成像剖面。
对于每一炮，标准的波动方程叠前深度偏移可以分为三步：震源波场的正向延拓、炮集记录波场的反向延拓和应用成像条件求取成像值（Clearbout, 1971）。
为了方便叙述基于共炮集的波动方程叠前深度偏移的基本过程，我们引入基于单程波方程的波场传播算子（Berkhout, 1987），并以频率域二维波场为例加以说明。
对震源波场);,(ωz x u s 和炮集记录波场);,(ωz x v s 做如下定义：（1）);0,(ωx u s ：它是炮点s 处频谱为)(ωf 的点源激发产生的震源波场，有)()();0,(ωδωf s x x u s -= （4-1）（2）);0,(ωx v s ：它是点s 处激发，排列接收到的记录波场，该波场可以写成：dr x v x v r s s ⎰=);0,();0,(,ωω （4-2） 其中，);0,(,ωx v r s 含有一非零道，即在接收点r 处的记录道，它满足： );0,()();0,(,ωδωx v r x x v s r s -= （4-3）（3）);,(ωz x u s ：它表示在深度0>z 处的正向延拓波场，如果引入表征波场从地面传播到深度z 的传播算子)0(z W →，则有：);0,()0();,(ωωx u z W z x u s s →= （4-4）（4）);,(,ωz x v r s ：它表示记录波场);0,(,ωx v r s 在深度z 的反向延拓波场： []);0,()0();,(,1,ωωx v z W z x v r s r s -→= （4-5） 其中，[]1)0(-→z W 为记录波场的反向传播算子。
三维波动方程有限差分正演方法
三维波动方程有限差分正演方法摘要：本文主要介绍了三维波动方程有限差分正演方法的基本原理和实现过程，并对其优缺点进行了分析。
通过数值实验验证了该方法的可行性和准确性，为地震勘探、地下水文学等领域的研究提供了参考。
关键词：三维波动方程、有限差分、正演方法、数值实验一、引言三维波动方程是地震勘探、地下水文学等领域的基础理论，其求解方法对于相关领域的研究具有重要意义。
有限差分正演方法是求解三维波动方程常用的数值方法之一，其具有计算速度快、精度高等优点，在实际应用中得到了广泛的应用。
本文将介绍三维波动方程有限差分正演方法的基本原理和实现过程，并对其优缺点进行分析，同时通过数值实验验证该方法的可行性和准确性。
二、三维波动方程有限差分正演方法的基本原理三维波动方程可以表示为：其中，u(x,y,z,t)为波动场，c(x,y,z)为介质速度，ρ(x,y,z)为介质密度，t为时间。
有限差分正演方法通过将空间和时间离散化，将三维波动方程转化为差分方程，进而求解波动场在不同时刻的数值解。
具体而言，有限差分正演方法将空间和时间分别离散化，将空间网格点和时间网格点相结合，构成一个四维网格空间，其中每个网格点对应一个波动场的数值解。
有限差分正演方法的求解过程可以分为以下几个步骤：1. 空间离散化对于三维空间，可以将其分为x、y、z三个方向，分别进行离散化。
假设在x方向上，空间步长为Δx，在y方向上，空间步长为Δy，在z方向上，空间步长为Δz。
则可以将空间网格点表示为(xi,yj,zk)，其中i=1,2,...,nx，j=1,2,...,ny，k=1,2,...,nz，nx、ny、nz分别为x、y、z方向的网格数。
2. 时间离散化假设时间步长为Δt，则在t时刻波动场的数值解可以表示为ui,j,k^n，其中n表示时间步数，i、j、k表示空间网格点的索引。
3. 有限差分近似将三维波动方程中的导数项用有限差分近似表示，例如：其中，ui,j,k^n表示在t时刻，在(xi,yj,zk)处的波动场数值解，c(xi,yj,zk)表示在(xi,yj,zk)处的介质速度，ρ(xi,yj,zk)表示在(xi,yj,zk)处的介质密度，Δx、Δy、Δz、Δt分别为空间和时间步长。
波动方程有限差分正演模拟误差来源分析
１２ 波动 方程 的高 阶差分 近似 ．
在做 有限 差分波 动方程 正演 数值模 拟 时主要 是 要综 合考 虑差 分格式 、 震源 的选择 、 边界 条件 的处 理
一 △ ） 开有 ｆ 展
电
０ ，
利用 泰勒公 式将 函数 （ ） ｔ 在 （ ＋△ ） ｔ ｆ 和 （ 一 ｔ
６ 。 ，。 ｘ
，
，
ｙ
，
ｚ
＝：
．
。
｛ 。 ：
动 方程 时 ， 常会产 生 不 期望 的数 值 频散 或 称 网格 常 频散 ， 导致 数值 模 拟结 果 分 辨 率 的降 低 。所谓 数值
频散 实质 上是 一种 因离散化求 解 波动方 程而 产生 的 伪 波动 ， 它是有 限差 分 方 法求 解 波 动 方程 时 所 固 有
０ ｕ １ Ｚ （ ０ ｕ ０（１ ）＋ ０（１ ）＝ Ｏｕ＋６
ｐ ｐ
㈩ （） ２ ｆ（ Ａ
ｆ ）。 （） １
取 到展开 式 的二阶项 ， 令式 （ ）＋ （ ） 到二 阶有 ３ 式 ４得
限差 分计 算公式
，ｆ ： ，） （ （ｔ Ａ） 。（ ） ８
假设 物质 的密 度不 随空 间变化 ， 方程 可表示 为 则
（． １ 东华理 工学 院 ， 西 抚 州 ３４ ０ ； ． 江 ４ ００ ２ 中南大 学 ， 南 长 沙 湖
湖南 长沙 ４０ １） １０ １
４０ ８ ； ． 南 省核 工 业地 质局 ， １０ ３ ３ 湖
摘 要： 有限差分波动方 程正演数值模拟 主要是解决差分格 式 、 震源选 择 、 边界 条件处 理等 问题 。这些 问题 的合 理 解决对得到正确和准确的波动方程数值解 有着 十分 重要 的作 用。笔者通 过对波动 方程计算 机有 限差 分正演模拟 过程的叙 述 ， 分析了在做正演模 拟时的误差来源 以及 如何有效地减小误差 。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
王红落,常 旭,陈传仁.基于波动方程有限差分算法的接收函数正演与偏移.地球物理学报,2005,48
(2)
:415～422
WangHL,ChangX,ChenCR.Receiverfunctionforwardmodelingandmigrationbasedonwave2equationfinitedifferencemethod.Chinese
J.Geophys.(inChinese),2005,48(2):415～422
基于波动方程有限差分算法的接收函数正演与偏移王红落1,常 旭1,陈传仁21中国科学院地质与地球物理研究所,北京 1000292长江大学地球物理与石油资源学院,荆州 434023
摘 要 针对接收函数正演与偏移,本文采用波动方程有限差分算法.借鉴成熟的勘探地震学方法,引入等效速度概念,建立接收函数转换波与地震勘探反射波的等效走时方程,实现了基于波动方程有限差分算法的接收函数正演与偏移.数值计算表明,波动方程有限差分叠后偏移方法可以对点绕射和穹隆构造模型实现高精度成像.本文利用数值计算讨论了波动方程有限差分叠后偏移与Kirchhoff叠后偏移对于接收函数偏移的适用性,还对偏移过程中速度模型的误差进行了分析.
关键词 有限差分正演,接收函数阵列,共转换点(CCP)叠加,有限差分偏移文章编号 0001-5733(2005)02-0415-08 中图分类号 P631 收稿日期 2004-04-13,2004-12-08收修定稿
基金项目 国家自然科学基金项目(40174017,40235055)资助.
作者简介 王红落,男,1975年生,1998年毕业于大庆石油学院,现为中国科学院地质与地球物理研究所博士研究生,主要从事勘探地震数字处理方面的工作.E2mail:whl@mail.iggcas.ac.cn
Receiverfunctionforwardmodelingandmigrationbasedonwave2equationfinitedifferencemethod
WANGHong2Luo1,CHANGXu1,CHENChuan2Ren21InstituteofGeology&Geophysics,ChineseAcademyofSciences,Beijing100029,China
2TheCollegeofGeophysicsandPetroleumResources,YangtzeUniversity,Jingzhou434023,China
Abstract Thewave2equationfinitedifferencealgorithmisappliedhereasasolutionforforwardmodelingandmigrationofreceiverfunction.Withthehelpofexplorationseismologyconcepts,theequivalentvelocityandtraveltimeequationareestablished,thereforetheforwardcalculationandmigrationofreceiverfunctionareimplementedbywave2equationfinitedifferencemethod.Thenumericalcalculationindicatesthatimagingofpointdiffractoranddomecanbeaccomplishedwithhighaccuracyemployingpoststackwave2equationfinitedifferencemigration.ThentheapplicabilitiesforreceiverfunctionimagingbetweenpoststackfinitedifferencemigrationandKirchhoffmigrationarediscussedaccordingtonumericalresults.Finally,theeffectsonvelocityerrorforimagingareanalyzed.Keywords Finitedifferencemodeling,Receiverfunctionarray,Commonconvertedpointstack,Finitedifferencemigration
1 引 言地震台站记录是震源时间函数、地下介质结构以及仪器特性耦合作用的结果.接收函数是消除震源效应以及仪器特性后,台站下方介质结构的径向响应.早期的接收函数方法是用单个固定台站远震体波资料来反演台站下方地壳的一维速度结
第48卷第2期2005年3月地 球 物 理 学 报
CHINESE JOURNAL OF GEOPHYSICSVol.48,No.2Mar.,2005
© 1994-2007 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net构[1～3],后来逐渐用于调查和证明岩石圈深部间断面的存在[4～6].将不同点的一维速度结构在横向上连成剖面,就可以追踪解释地下速度界面分布[7].随着全球范围内数字地震台网的进一步加密,就有可能提取近于勘探地震22D甚至32D观测系统的接收函数阵列.尤其是最近流动地震台网的大范围应用,使得台站按需要布置成任意的观测系统成为可能.勘探地震在地面高密度数据采集的前提下,利用数字信号分析和波传播理论,得到22D或32D的地下精细结构[8].借鉴成熟的勘探地震学的方法,直接从接收函数阵列对地下界面成像,是最近地震学研究的新方向.对于高密度空间采样接收函数阵列数据,将勘探地震中更精确的波动方程偏移成像方法推广到接收函数阵列,就具有明显的现实意义.近年来,在研究接收函数阵列PS震相运动学规律的基础上,前人已经做了很多开创性的工作.发展了对接收函数抽取共反射面元阵列并进行叠加的方法,大大提高了接收函数的信噪比[9～11].按照偏移成像的概念,要对转换界面成像,就必须把观测数据中的转换震相按波传播的相反方向聚焦到转换点上.Jonesetal.[12]提出了把叠加后的接收函数进行射线路径逆向反传的技术,这种沿射线逆传的成像是一种高频近似方法,而非波动方程偏移.Rybergetal.[13]提出了接收函数转换波与地震勘探反射波的等效时距关系,使反射地震数字处理方法和流程在接收函数上的应用成为可能.本文主要研究波动方程有限差分偏移算法在接收函数阵列上的应用.本文考虑地震波动力学特征,引入相对传播时间的概念,采用波动方程有限差分算法正演接收函数,这样得到的接收函数保存了地震波的动力学信息,更接近于实际的接收函数波形记录.根据反射地震共中心点叠加原理,引入等效于反射地震的时距关系,在此基础上,本文对接收函数阵列实现了波动方程有限差分偏移,该方法与Kirchhoff积分偏移结果比较,成像精度明显改进.数值计算结果证明了本研究的可行性.2 基于波动方程有限差分的接收函数正演 一个地震台站接收到的三分量远震P波地震记录可表示为[3]DV(t)=I(t)3S(t)3EV(t),DR(t)=I(t)3S(t)3ER(t),
DT(t)=I(t)3S(t)3ET(t),(1)其中D
V(t),DR(t)和DT
(t)分别表示地震记录中
垂直分量、径向分量及切向分量,S(t)表示入射平面P波震源时间函数,I(t)表示仪器特征响应,
EV(t)、ER(t)及ET(t)分别表示介质结构脉冲响应
的垂直、径向及切向分量.由径向分量D
R
(t)、切向
分量D
T(t)分别与垂直分量DV
(t)反褶积,获得的
地球介质的径向、切向脉冲响应E
R(t)、ET
(t)即为
通常所说的接收函数.因此,要正演接收函数,首先要正演三分量地震记录.
本文着重研究二维情况下两分量地震记录接收函数的计算.给定二维弹性介质模型,用二维弹性波动方程模拟P2SV波场,并从中提取模型表面各点波场随时间的变化,构成水平、垂直两分量地震记录.二维P2SV波在非均匀各向同性介质中的速度-应力弹性动力学方程为
ρv
・
x=9σxx9x+9σxz9z,
ρv
・
z=9σxz9x+9σzz9z,
σ・xx=(λ+2μ)9vx9x+λ9vz9z,
σ・zz=λ9vx9x+(λ+2μ)9vz9z,
σ・xz=μ9vx9z+μ9vz9x,(2)
其中(v
x,vz)为质点两分量速度矢量,ρ(x,z)为密
度,
(
σ
xx,σxz,σzz)为应力张量,v和σ分别表示速度
及应力对时间的导数,λ(x,z)和μ(x,z)为拉梅系数[14].接收函数是从远震记录中得到的,在正演接
收函数的计算中,考虑到远震的传播距离,可以假设远震波前进入到模型空间时是平面波,并将这一时刻作为震源的相对初始时刻,在模型空间范围内计算波场随时间的变化.
为了保证接收函数一致性并提高接收函数信噪比,最终得到波场的径向分量,需要对弹性波正演获得的两分量地震记录进行坐标旋转[15].如图1所
示,把从XZ坐标系观测的水平、垂直两分量记录投影到LQ坐标系,其中L方向为直达P波的射线路径方向,Q方向与L方向正交,L与Q的分量分别对应(1)式中的D
V(t)和DR
(t).对各向同性介
质,L方向分量主要是P波能量,Q方向分量是SV
波能量,坐标旋转其实质是一个波场分离的过程.
614地球物理学报(ChineseJ.Geophys.)48卷
© 1994-2007 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net