用伪谱法合成零偏移距地震记录

1、第28卷第4期长春科技大学学报V o l.28N o.41998年10月JOU RNAL O F CHAN GCHUN UN I V ER S IT Y O F SC IEN CE AND T ECHNOLO GYO ct.1998用伪谱法合成零偏移距地震记录张文生1何樵登2李树勇1(1中国科学院计算数学与科学工程计算研究所,北京100080;2长春科技大学,长春130026摘要零偏移距地震记录是叠后偏移的基础,常用射线追踪法来制作,但当介质复杂时,射线会出现交叉现象,对三维情况计算量更大。从二维标量波动方程出发,用伪谱法研究了零偏移距地震记录的合成,给出了不同情况下的数值计算公式,并对边界吸

2、收问题作了有效处理。数值计算结果表明了方法的正确性和有效性。该方法可用于偏移中正演模拟的快速计算,且易推广到三维情形。关键词伪谱法正演零偏移距记录反周期扩展法中图分类号P631.40引言在波动方程偏移中常需制作零偏移距地震记录,零偏移距记录通常用射线追踪法和波动方程法来制作。射线追踪法依据的是几何地震学,不易反映波场的动力学特点,另外对复杂介质会出现射线交叉现象,对三维情况,因需作走时计算,计算量更大。波动方程法基于波在介质中传播的波动性,通过数值求解波动方程来实现,能较全面地反映波场特征。波动方程数值求解通常有有限差分法、有限元法和伪谱法三种。有限差分法是最常用的一种有效方法,但精度低,网格

3、步长必须明显小于最短波长的一半,频散现象也较严重,如采用高阶差分法,能提高精度,但势必会增加计算量,目前该方法向高阶交错网格法、不规则网格法及大网格快速差分法等方向发展。有限元法也是正演模拟的有效手段,能有效地求解不规则区域,但计算量相对较大,目前该方法向多重网格的区域分解法和混合元方向发展。伪谱法是一种较新和发展较快的方法,不同于传统的有限差分法,其空间差分用快速傅氏变换来实现,精度可认为是无限阶的1,且理论上已证明,伪谱法在一个空间波长上最少仅需两个网格点,因此可采用较大的空间步长,减少计算量。有限差分法,有限元法和伪谱法在地震波动方程数值模拟中都得到了较广泛应用,但通常都用来合成共炮点地

4、震记录及作波场传播瞬时切片图,如要得到零偏移距即自激自收记录,固然可再作抽道集等处理得到,但显然是费时和不必要的。本文用伪谱法研究了零偏移距地震记录的合成,因其计算速度快,故可应用于C laerbou t偏移特别是三维偏移处理中。中国博士后科学基金资助第一作者简介张文生男31岁博士后地球物理专业已发表“最小二乘法与反问题误差最优模的选择”等论文收稿日期199711141方法原理1.1计算公式零偏移距地震记录的制作要用到爆炸反射界面原理2。该原理把地下反射界面假设成具有爆炸性的物质,计算时将震源函数放于界面上,并且假定在t =0时所有的爆炸反射界面同时起爆。当介质的速度用半速度代替后,由界面上各

5、散射源发出的初至波被地面检波器接收,其结果就是近似的CD P 叠加剖面即零偏移距记录。假定地震波满足如下二维标量波动方程52u 5t 2=v 2(52u 5x 2+52u 5z 2(1其中u (x ,z ,t 是位移,t 是时间,v (x ,z 是介质速度,x ,z 分别是水平和垂直距离。将(1式左端用二阶时间差商代替,右端用伪谱法计算,可得u (x ,z ,t +t 2u (x ,z ,t -u (x ,z ,t -t -t 242(k 2x +k 2z v 2(x ,z u (k x ,k z ,t e i (k x x +k z z d k x d k z (2其中t 为时间采样,k x

6、 、k z 分别为x 、z 方向的波数,u (k x ,k z ,t 为u (k x ,k z ,t =u (x ,z ,t e -i (k x x +k z z d x d z (3特别地,如速度v (x ,z 为常数c ,则由(1式可导出u (x ,z ,t 的解析解3。先对(1式作关于x 、t 的二维傅氏变换,整理得52u (k x ,z ,z 2=-(2c 2-k 2x u (k x ,z ,=-k 2z u (k x ,z ,(4其中u (k x ,z ,=u (x ,z ,t e -i (k x x +t d x d t (5因速度为常数,故在z 方向上的波数也为常数,由(4式并利

7、用爆炸反射界面条件,经推导可得u (k x ,z ,=u (k x ,k z ,t =0e ik z z c 1-(ck x 2(6再对(6式作关于k x 、的二维反傅氏变换,得u (x ,z ,t =142u (k x ,k z ,01c 1-(ck x 2e i (k x x +k z z +t d k x d (7由(7式即得零偏移距剖面u (x ,0,t =142u (k x ,k z ,0c 1-(ck x2e i (k x x +t d k x d (8其中k 2z =2c2-k 2x 。824长春科技大学学报1998年上面(2式是变速时的零偏移距剖面的数值计算公式,而(8式是常速

8、时的零偏移距剖面的解析表达式,这两式均使用了全波波动方程。对某些模型或复杂模型,记录中会出现多次波,为避免多次波的干扰,可采用无层间反射的单程波波动方程4。5u (k x ,k z ,t 5t=iv (x ,z k z 1+k 2x k 2z u (k x ,k z ,t (9该式也易从(1式导出(从略,其中u (k x ,k z ,t 如(3式所示。同推导(2式相类似的方法,采取时间中心差商格式,最后可得u (x ,z ,t +t u (x ,z ,t -t +ti 22u (k x ,k z ,t v (x ,z k z 1+k 2x k 2z e i (k x x +k z z d k

9、x d k z (10由(10式可算出波场u (x ,z ,t ,令z =0即得零偏移距剖面u (x ,z =0,t 。1.2震源、稳定性条件及吸收边界条件用方程(2、(8、(10作数值计算时,要用到初始条件,这里采用时间为雷克子波的源函数形式s (x ,z ,t ,即s (x ,z ,t =A f (x ,z co s2f 0(t -t 0e -(t -t 02f 20(11其中A 为振幅,为衰减系数,f 0为主频,f (x ,z 为界面函数。伪谱法计算的稳定性条件可从矩阵稳定性分析得到5,不难推得,(2式、(8式和(10式计算的稳定性条件为t 2d v m ax (12其中空间网格取d =

10、x =y ,v m ax 取成最大速度值。数值模拟中另一必须解决而又难以解决的问题是吸收边界条件。通常有限差分法采用粘滞型(或称海绵型吸收条件,它在扩充区域内逐渐将波场衰减吸收掉,从而不产生边界反射;该方法适用性广,但难以选择合适的吸收系数(一般用试验方法得到,而不合适的吸收系数仍会导致很强的反射。另一种要推导边界吸收方程,有限差分时将其结合在差分方程中计算,但它仅对一定角度的入射波有理想吸收效果。这些方法都难以应用于伪谱法中。在伪谱法中,根据傅氏变换的周期性条件,人们提出了一种反周期扩展的方法6,该方法已在各向异性介质共炮点地震记录正演中得到了成功应用5,本文采用了这种方法。它主要通过扩展计

11、算区域增加计算量来消除记录中的边界反射,效果良好。2数值计算前面导出了不同情况的用于零偏移距地震记录的公式(2式、(8式和(9式,数值计算时可根据情况利用。根据这些公式,给出如下三个模型的计算结果,加以说明和验证。模型1模型由一个点绕射体和一个水平界面组成,其空间采样大小为N x =64,N z =64。绕射点位于(32,7处,水平界面位于N z =20处,介质速度为3000m s 。图1是所得的零偏移距记录剖面,它由一个绕射双曲线和一个水平界面构成,图中的采样x =z =20m ,t =2m s ,时间样点数为250,所示单位均以空间或时间方向的采样点数来表示(下同,其长度可由给定的采样步长

12、换算得到。根据射线几何路径分析,不难验证,图1波场924第4期张文生,何樵登,李树勇:用伪谱法合成零偏移距地震记录 图1模型1的零偏移距记录F ig .1Zero -off set syn thetic record of m odel 1是正确的。模型2在模型1的基础上讨论一个稍复杂的模型(图2。它模拟了一系列向斜构造(45,速度仍假定为常数3000m s 。图3是用(2式算得的未加边界吸收的记录,图4是对应的加边界吸收后的结果。比较可知,吸收效果明显。图5是用射线追踪法所得的记录,比较图4与图5可知,从另一侧面说明了本方法及算法的可行性和正确性,图中的采样x =z =20m ,t =2m

13、s ,时间样点数为400,空间样点数为64。模型3如图6,该模型由三个界面四层构成,四层的层速度分别为2000m s 、2500m s 、3000m s 、2500m s 。采用无层间反射的波动方程计算,所得的零偏移距地震记录如图7所示。由于受第二界面中嵌型构造的影响。第三界面中央记录中出现弯曲;第二个界面所导致的绕射波形清晰可见。图中x =z =25m ,t =2m s ,时间样点为600,空间样点为64,所用CPU 时间约为0.42s ,所用计算时间约为2.86m in 。图2类似向斜构造的常速模型F ig .2A con stan t veloc ity m odel ana logue

14、d a se -r ies of syncli nes 图3模型2(图2的未加边界吸收的零偏移距记录F ig .3Zero -off set syn thetic record of m odel 2(f ig .2with no boundary absorbi ng con -d ition s034长春科技大学学报1998年 图4模型2(图2加边界吸收的零偏移距记录F ig .4Zero -off set syn thetic record of m odel 2(f ig .2af ter boundary absorbi ng 图5模型2(图2用射线追踪法合成的零偏移距记录F ig

15、.5Zero -off set syn thetic record of m odel 2(f ig .2usi ng the ray trac i ng m ethod图6层变速模型F ig .6A m odel with layer var ied veloc ity 图7模型3(图6的零偏移距记录F ig .7Zero -off set syn thetic record of m odel 3(f ig .63结论快速有效地合成零偏移距记录是进行叠后偏移的基础。本文给出了用全波波动方程和无层间反射的波动方程制作零偏移距地震记录的方法,对边界吸收用反周期扩展法作134第4期张文生,何樵登

16、,李树勇:用伪谱法合成零偏移距地震记录432 长春科技大学学报 1998 年 了有效处理 ( 当然不排除使用衰减型吸收边界条件 。 模型计算表明了所用方法及算法的 正确性和有效性, 效果良好。 由于傅氏积分采用了快速离散傅氏变换来实现, 提高了计算 速度, 尽管在边界吸收中采用了反周期扩展法扩展了计算域增加了计算量, 但计算速度还 是快的, 这些为推广到三维情形奠定了良好基础。 本文得到了张关泉老师的大力帮助和支持, 同韩立国副教授在傅氏变换方面也曾作 过有益讨论, 在此均表示衷心感谢。 参 考 文 献 1 Fo rnberg B. T he p seudo sp ectra l m etho

17、d: Com p a rison s w ith fin ite d ifferences fo r the ela stic w ave equa tion. Geop hysics, 1987, 52 ( 4 : 483 501 2 L oew en tha l D , L u L , Roberson R , et a l T he w ave equa tion app lied to m ig ra tion. Geop hysica l . P ro sp ecting, 1976, 24 ( 2 : 380 390 3贺振华主编 1 反射地震资料偏移处理与反演方法 1 重庆: 重

18、庆大学出版社, 1989 ( 6 : 854 859 1997 1 2 Zhang W en sheng , H e Q iaodeng , L i Shuyong 4 Gazdag J. M odeling of the acou stic w ave equa tion w ith tran sfo rm m ethod s Geop hysics, 1981, 46 . 5张文生 1 各向异性介质弹性参数反演理论方法与数值方法研究: 学位论文 1 长春: 长春科技大学, w ith the p seudo sp ect ra l m ethod. N o t on ly w e g iv

19、e the num erica l com p u ta t ion fo rm u lae . ca l ca lcu la t ion s show the co rrectness and effect iveness of m ethod s So w e supp ly a rap id m odelling m ethod fo r m ig ra t ion, w e a lso can ex tend it to th ree d i en sion s ca se ea sily. m p eriod ic ex ten sion m ethod ca ted m ed ia

20、, and com p u ta t ion is t i e 2con sum ing fo r a 32 ca se. In the p ap er, ba sing on m D 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. fo r d ifferen t ca ses, bu t a lso so lve the bounda ry ab so rb ing p rob lem effect ively. N um eri2 6 Fu rum u ra T , T akenaka H. A nw rap a round eli ina tion techn ique fo r the p seudo sp ectra l w ave syn 2 m ly, it is aecom p lished by the ray t racing m ethod, how ever, ray s m ay cro ss fo r a com p li2 the 22 sca la r