各向异性介质中二阶弹性波方程模拟PML吸收边界条件 - 图文-

1、第27卷第5期2007年10月大地测量与地球动力学JOURNAL OF GEODESY AND GEODY NAM I CSVol.27No.5Oct.,2007文章编号:167125942(20070520054205各向异性介质中二阶弹性波方程模拟P ML吸收边界条件3朱兆林马在田(同济大学海洋与地球科学学院,上海200092摘要提出了基于弹性波动方程的二阶双曲系统的一种新的弹性波动方程模拟的最佳匹配层(P ML吸收边界条件。对于二阶系统方程,P ML吸收边界条件模型通常是以四分裂位移参量来构建的,这种方法需要求解时间的三阶导数,且占用较多的计算空间;作为另一种选择,非分裂的P ML算法可

2、以扩展到二阶系统中,但它需要求解二重时间积分。新方法可克服或简化上述问题。用交错网格有限差分方法加最佳匹配层边界条件新算法的模拟方法用数值模型作了试验,结果证实了此方法的有效性。关键词模拟最佳匹配层吸收边界二阶波动方程各向异性介质中图分类号:O343文献标识码:AS I M ULAT I O N O F P ML ABS O RB I NG BO UNDARY CO N D I T I O N W I TH SECO N D2O R D ER ELAST I C W AVE EQUAT I O N I N AN I SO TRO P I C M ED I AZhu Zhaolin and Ma

3、 Zaitian(School of O cean and Ea rth Science,Tongji U niversity,S hanghai200092AbstractA ne w alternative perfectly matched layer(P MLabs orbing boundary conditi on is devel oped t o at2 tenuate the artificial boundary reflecti ons generated in nu merical si m ulati on of the second2order elastic wa

4、ve equa2 ti on.The second2order equati on can be described by dis p lace ment,which is more app r op riate than the first2order one.Its P ML conditi on conventi onally needs t o s p lit the dis p lace ment int o f our parts,which occup ies a large a2 mount of me mory and requires s olving a third2or

5、der differential equati on in ti m e.A s f or the other choice,non2s p lit2 ting P ML method may be app lied t o the second2order equati on,but it requires s olving the dual integral in ti m e.The ne w method can s olve or si m p lity the above p r oble m s.Finally,a staggered2grid finite difference

6、 method with this P ML conditi on is used t o si m ulate an anis otr op ic media model and the results show that the method is efficient. Key words:si m ulati on,perfectly matched layer(P ML,abs orbing boundary,second2order elastic wave equati on, anis otr op ic media1引言由于计算区域的有限性,数值模拟需要有效限定吸收边界条件限制

7、人为边界的反射。在最近的30多年里,计算地球物理学家提出了许多地震波动方程的数值模拟吸收边界条件方法,主要分为两类:一类是傍轴近似的吸收边界条件方法,如Clayt on& Engquist1和Reynolds2等;另一类是阻尼层吸收边3收稿日期:2007205216基金项目:973项目“碳酸盐岩缝洞型油藏开发基础研究”(2006CB202402作者简介:朱兆林,男,1979年生,博士生,研究方向为地震波传播与反射地震学.E-mail:zhuzl2004第5期朱兆林等:各向异性介质中二阶弹性波方程模拟P ML吸收边界条件界条件,如Cerjan3等。当前研究比较广泛的吸收边界条件是最佳匹配

8、层(P ML吸收边界条件,属于第二类。P ML吸收边界条件最初被Berenger4应用于麦克斯韦电磁场方程的数值模拟,后来被许多研究者扩展到声波、弹性波等其他波场数值模拟过程中,且同样取得了较好的效果,例如Hastings5和Collino6把P ML吸收边界条件应用到了速度2应力的一阶弹性波方程中;Zeng7在粘弹性介质的波动模拟中做了进一步扩展。在弹性波的P ML吸收条件能够比较方便地应用到一阶双曲系统方程中的基础上,Komatitsch8提出了二阶弹性波动方程P ML吸收边界条件的分裂方法。另外,相对分裂方法而言,非分裂算法也是一种选择,W ang9就提出了一阶弹性波动方程的非分裂算法的

9、P ML模型。在各向异性介质中,P ML吸收边界条件可能出现不稳定现象,Becache10分析了不同各向异性介质的稳定性条件,Appel o11给出了各向异性介质中一种新的P ML吸收边界条件的实现方法。P ML吸收边界条件是在计算区域外加一定厚度的吸收层,并用有衰减因子的波动方程计算吸收层内的波场来实现的。波动方程一阶系统的P ML条件分裂算法是把位移场分为两项,而二阶系统分裂为4项,后者需要解一个三阶时间导数式且需较多的内存空间。如果把非分裂算法9扩展到二阶系统,会产生新的困难,既不能节省空间还需要解二重积分。本文提出了一种新的P ML分裂算法并做了相应的数值模拟实验。2原理和算法2.1二

10、阶弹性波动方程文中公式推导都基于下列二阶弹性波方程5U5t2=55x A5U5x+B5U5y+C5U5z+5 5y D 5U5x+E5U5y+F5U5z+55z G5U5x+H5U5y+I5U5z(1A=C11C16C15C16C66C56C15C56C55,B=C16C12C14C66C26C46C46C25C45C=C15C14C13C56C46C36C55C45C35,D=C16C66C56C12C26C25C14C46C45E=C66C26C46C26C22C24C46C24C44,F=C56C46C36C25C24C23C45C44C34G=C15C56C55C14C46C45C13

11、C36C35,H=C56C25C45C46C24C44C36C23C34I=C55C45C35C45C44C34C35C34C33其中,(x,y,z是空间坐标系统,U=(u,v,w是位移波场向量,u,v和w分别是沿x,y和z轴的波场分量;Cik(i,k=1,2,3,4,5,6是弹性参数,是介质密度,t表示时间。把方程(1变换到频率域,则-2U=55x A5U5x+B5U5y+C5U5z+55y D5U5x+E5U5y+F5U5z+55z G5U5x+H5U5y+I5U5z(2其中,是圆周频率。把方程(2的对应边界方向坐标变换到一个复数空间,然后就可以构建带有衰减因子的P ML公式。下面讨论P

12、ML公式的推导过程。2.2传统的P ML算法如图1所示,模型分为两部分:内部计算区域;P ML区。假设只考虑模型的左右边界(x方向P ML吸收边界条件构建过程。图1中,左右边界分别在x0=0和x0=15处。过程的第一步是坐标变换,新的x方向坐标为x=x-jxx0(xd x(3其中, x是新的复数空间坐标;是衰减因子,它是一个关于x变量的实函数,可以描述为12:(x=20f0x-x0L P ML2P ML域0计算域(4这里为指数常数,通常取值为1.79,f表示子波的主频,LP ML是P ML吸收层的厚度。根据式(3,得到算子5/5 x=(1/(5/5x,=1+/(j,然后用(1/(5/5x替代式

13、(2中的5/5x,得到-2U=155x1A5U5x+1255x A5U5x+155x B5U5y+C5U5z+155y D5U5x+155z G5U5x+55y E5U5y+F5U5z+55z H5U5y+I5U5z(5分裂波场U为4部分,表示为U=U1+U2+U3+U48。根据1/=(j/(j+和5(1/5x=-jx/(j+2,x表示5/5x。方程(5分裂为:55大地测量与地球动力学27 卷图1P ML模型(a及展开(bFig.1Model P ML(aand its devel opment(b(j+2U1=55x A5U5x1x (j+3U2=-A5U5x(j+jU3=5 5x B 5U

14、5y+C5U5z+55y D5U5x+55z G5U5x-2U4=5U 5y E 5U5y+F5U5z+5U5z H5U5y+I5U5x(6式(6反变换到时间域可以得到传统的二阶弹性波动方程的P ML吸收边界条件55t+2U1=55x A5U5x1x 55t+3U2=-A5U5x55t55t+U3= 55x B 5U5y+C5U5z+55y D5U5x+55z G5U5x52U4 5t2=5 5y E 5U5y+F5U5z+55z H5U5y+I5U5x(7从方程(7可知4项分裂P ML算法的形式,其需要求解U23次时间导数且需消耗额外的计算空间8。2.3新的P ML算法W ang9在2D柱坐

15、标下提出了速度-应力波动方程非分裂最佳匹配层(NMP ML方法。这种方法可以直接扩展到二阶系统方程,但是会带来新的问题,比如说二重积分问题,且因为增加了许多补偿量,所以也不能节省计算空间。本文提出了新的方法,能够合并上面分裂算法方程(6中的U1、U2和U3,这种方法能够避免时间方向的三阶导数也能节省部分计算空间。重写方程(6为:(j+2U1=55x A5U5x(j+2U2=-xj+A5U5x(j+2U3=1+j55x B5U5y+C5U5z+55y D5U5x+55z G5U5x-2U4=55y E5U5y+F5U5z+55z H5U5y+I5U5x(8因为式(8中第三个公式没有包括=0,所以

16、校正如式(9(可以直接从时间域反推导得到(j+2U3=1+j+(55x B5U5y+C5U5z+55y D5U5x+55z G5U5x(9如此,我们能够把原始波场分裂为两项U5=U1+U2+U3和U4,也可以分裂为3项U5=U1+U2,U4和U3。转换到时空域,式(8和式(9有两项傅立叶变换关系如下:F-1-xj+=-x e-tF-11+j+(=(t+s(t(10其中,(t是冲激函数,s(t是阶跃函数,即t>0,s(t=1;t<0,s(t=0。因此时间域二项分裂式为:55t+2U5=55x A5U5x-x e-t3A5U5x+55x B5U5y+C5U5z+55y D5U5x+55

17、z G5U5x+s(t355x B5U5y+C5U5z+55y D5U5x+55z G5U5x52U45t2=55y E5U5y+F5U5z+55z H5U5y+I5U5x(1165第5期朱兆林等:各向异性介质中二阶弹性波方程模拟P ML 吸收边界条件其中,3是褶积算子。下面求解式(11中的积分形式9。假设P =-x e-t3A 5U 5x (12和Q =2s (t 355x B 5U 5y +C 5U 5z +55y D 5U 5x +55zG 5U 5x (13由积分式(12和(13得到P n =e -t P n -1-x t 2e -t A 5U n -15x +A 5U n5x (14

18、Q n=Qn -1+t 455xB 5U 5y +C 5U5z +55y D 55x +55zG 55x ( U n + U n -1(15理论上,新的算法(两项分裂或者三项分裂P ML 和传统的四项分裂P ML 有同样的精度。下面用模型测试上述方法的有效性。3模型试验各向异性介质中弹性波的三维交错网格有限差分方法13用于求解式(11,(14和(15。震源函数为:f (x,y,z,t =g (x,y,z 1-2f 0t -0.5f 02e-f 0t -0.5f2g (x,y,z =e-7(r/r 02r20r =(x -x s 2+(y -y s 2+(z -z s 2其中,g (x,y,z

19、是高斯函数10,(x s ,y s ,z s 为震源位置。模型是一个3D 均匀介质立方体,长、宽、高都为1000m 。介质是TTI 介质(VTI 沿y 轴逆时针旋转20°,VTI 介质参数用Thom sen 形式14表示为V p 0=5800m /s ,V s 0=3400m /s ,=2.4g/c m 3,=0.2,=0.1,和=0.035。模型的空间采样间隔都为10m 。子波的主频取f 0=30Hz,震源点位置取为(500,500,500,P ML 宽度是10个采样点。图2(a 和2(b 分别是模型I 的70m s 和120m s 的波场快照,自上而下分别是x 、y 、z 方向的

20、分量。可以看到当波传播到70m s 的时候,波场能量还没有到达边界,外层为纵波,内层为横波;传播到120m s,纵波的波场完全透过边界,人为边界反射吸收相当成功。模型是个2.5D 多层介质,沿y 方向介质参数不变,x 、y 、z 方向的长度分别为2500m 、1500m 和1500m (图3。第一层、三层和五层都为各向同性介质,纵波参数分别为V p 0=5300m /s,6300m /s,6800m /s,横波分别为V s 0=3250m /s 、3800m /s 、4350m /s,密度大小分别为=2.4g/c m 3、2.5g/c m 3、2.7g/c m 3。第二层、四层都为TTI 介质

21、(VTI 介质沿y 轴逆时针旋转20°,VTI 介质参数用Tho m sen 形式表示分别为V p 0=5800m /s、6400m /s;V s 0=3400m /s 、4000m /s,=2.45g/c m 3、2.55g/c m 3;=0.16、0.16;=0.08、0.08;=0.035、0.035。模型离散为251×151×151,采样间隔为10m 。用纵波震源,主频为f 0=30Hz,震源位置 图2模型I 的波场快照Fig .2Snap shots of wave field of model 75大地测量与地球动力学27卷(1250、750、50。图

22、4是模型二中波场传播220m s处的快照,依次分别为x、y、z方向分量的波场。可以看到波场通过顶层和x方向两边时,P ML条件较好地吸收了人为反射能量。4结语对P ML条件新算法,用两个各向异性介质模型进行了测试,结果证明新算法不仅有效,对比传统的四阶分裂算法能够避免求解时间方向的三阶导数,且能够节省部分计算空间。图3模型(x×y×z=2500×1500×1500m3 Fig.3Model(x×y×z=2500×1500×1500m3图4模型在t=220m s的波场快照Fig.4Snap shot of wave

23、field of model(t=220m sReferences1Clayt on R and Enquist B.Abs orbing boundary conditi ons f or acoustic and elastic wave equati onsJ.Bulletin of the Seis2 mol ogical Society of America,1977,67:1529-1540.2Reynolds A C.Boundary conditi ons for the nu merical s olu2 ti on of wave p r opagati on p r ob

24、le m sJ.Geophysics,1978, 43(6:1099-1110.3Cerjan C,Kosl off D,Kosl off R and ReshefM.A nonreflect2 ing boundary conditi on for discrete acoustic and elastic wave equati onJ.Geophysics,1985,50(4:705-708.4Berenger J P.A perfectly matched layer f or the abs or p ti on of electr omagnetic wavesJ.Journal

25、of Computati onal Phys2 ics,1994,114:185-200.5Hastings F,Schneider J B and B r oschat S L.App licati on of the perfectly matched layer(P MLabs orbing boundary con2 diti on t o elastic wave p r opagati onJ.Journal of the Acoustic Society of America,1996,100(5:3061-3069.6Collino F and Ts ogka C.App li

26、cati on of the P ML abs orbing layer model t o the linear elast odyna m ic p r oblem in anis otr op2 ic heter ogeneous mediaJ.Geophysics,2001,66(1: 294-307.7Zeng Y Q,He J Q and L iu Q H.The app licati on of the per2 fectly matched layer in nu merical modeling of wave p r opaga2 ti on in por oelastic mediaJ.Geophysics,2001,66(4:1258-1266.8Komatitsch D and Tr omp J.A perfectly matched layer ab2 s orbing boundary conditi on f or the second-order seis m ic wave equati onJ.Geophysical Journal I nternati onal,2003, 154:146-153.9TsiliW a