一阶弹性波方程的变网格高阶有限差分数值模拟

1、一阶弹性波方程的变网格高阶有限差分数值模拟2019年12月第43卷第6期正演技术一阶弹性波方程的变网格高阶有限差分数值模拟李振春3张慧张华（中国石油大学地球资源与信息学院,山东东营257061）李振春,张慧,张华.一阶弹性波方程的变网格高阶有限差分数值模拟.石油地球物理勘探,2019,43（6）:711716摘要使用可变网格的有限差分法进行地震模拟有许多独特的优点,主要表现为对地质模型的离散化更为合理,在低速带和复杂构造区域,可将局部网格划分得相对精细些,不仅提高了模拟精度,消除了因采样不足导致的频散现象,而且可以减少计算机内存需求,保持模型计算的灵活性。本文提出一种新的基于高阶交错网格技术的

2、弹性波数值模拟方法,通过改变网格的空间步长实现了局部网格加密技术,弥补了常规网格的缺陷和不足。试算结果表明,本文提出的算法稳定性较好,且能够提高模拟精度,减少计算时间,提高计算效率。关键词变网格交错网格高阶有限差分弹性波1引言法技术（地震采集。为适应复杂勘探条件的需要,、高效率的地震数值模拟技术变得越来越重要。自Virieux1率先采用交错网格有限差分法模拟P2SV波在非均匀介质中的传播以来,交错网格成为了有限差分法正演模拟中普遍使用的网格剖分技术。同二阶位移方程相比,使用交错网格的一阶速度应力方程,无论是模拟精度还是模拟的稳定性都有了较大提高。但是如果局部地区构造较为复杂,或者存在低速层时,

3、为了提高复杂构造的正演精度,通常对整个模型使用较小的采样网格,以保证对频散现象的抑制和对稳定性的保障。这对保持模型计算的灵活性非常不利,还会对模型的其他部分造成过采样,增加惊人的计算量。所以,设计一种更合理的离散化方式进行正演模拟是非常必要的。从空间采样的角度考虑,既要提高模拟精度、又要降低计算机内存需求的最有效方法就是在模型的不同区域采用不同步长的网格，即变网格（也称为不规则网格）。为此Jastram等提出了基于二维声波方程针对某一深度变网格步长的算法2;Jastram等=vx3;4Pitarka实现了低阶精度的不5。本文在前人研究的基础上提出一种高阶变网格数值模拟方法,实现了高精度的局部网

4、格加密技术,另外还改进了前人对变化网格过渡带采用的波场插值算法,避免了在过渡带区域进行大量的插值计算,也不需要平滑函数进行平滑,同时还保持了算法的稳定性。2方法原理文中所提出的算法是基于Virieux1和Graves6提出的一阶速度应力方程,对于各向同性介质,有txx+5t5tvx=b(5xzxz)+fxtxz+5t5tvz=b(5xzzz)+fz入+2u)5xvx+入5t5zvztxx=(入+2u)5zvz+入5t5xvxtzz=(5zvx+5xvz)5Ttxz二口(1)tt式中:vx,vz分别表示速度；txx,zz,xz分别表示应力向量；5x、5z、5t分别表示一阶偏微分算子5/5x、5/

5、5z、5/5t;b表示密度的倒数；fx,fz为体力向量；入及卩为拉梅常数。3山东省东营市中国石油大学地球资源与信息学院,257061本文于2019年1月8日收到,修改稿于同年7月30日收到。本课题由国家973专题(2019CB209605)、国家863专题(2019AA06Z206)和CNPC物探重点实验室中国石油大学(华东)研究室资助。712石油地球物理勘探2019年通过对时间导数采用二阶近似,得到弹性波一阶速度应力方程的离散形式为vxvn+2i,jn-2xxi+2,jni,jtxx+Dt+Atbx(Dxzxz)|ni,jnA1=dxi/2,A2=dxi/2n-1n-1k=1n+2zi+2,

6、j+2=vn-2zi+2,j+2t+Atbz(Dxxz+Dtzzz)|i+,j+22n+2i+2,jn+2i+2,jAn-1=k=1dxi+k+dxi/2,An=Edxi-k+dxi/2Ettn+1xxi+2,jn+1zzi+2,j=tnzzi+2,jDxg(x,z)入+2口)Dxvx+入+At(Dzvz图1计算对称点为i+的差分系数所需的变网格节点2n-1=T)Dzvz+入+At(入+2uDxvx同理，以i为对称点的差分算子为An-1=n-1k=1n+1nTxz(Dzvx+Dxvz)xzi,j+2=txzi,j+2+At(n+2i,j+2dxEi+(k-1)+dxi+n-1/2,An=k=1

7、dxEi-k+dxi-n/2。(2)式中:Dx,Dz分别表示对x,z的一阶微分算子；bx=(bi,j+bi+1,j)/2,bz=(bi,j+bi,j+1)/2;uxz=(1/ui,j+1/ui+1,j+1/ui,j+1+1/ui+1,j+1)/4。用变量g来表示速度VX、vz或应力向量TXX、ttzz、xz,则在常步长的交错网格技术中,用下式计算式(2)中的g关于x的2n阶空间偏导数7Cg(x+AxAxEnii=1下面给出变网格算法的任意偶数阶精度差分近似式及差分系数的求取方法,即nDxg(x,z)=i-1cE2i-1g(x+A2i-1,z)+(4)+c2i-A2i,z)n,)-i；i是空间差

8、分算子，它dg关于x。令g(x,z)=gzexp(ikx),则式(4)可写成n(3)-g(x-Axi-lzx2n式中，系数Cin由x(x2iT,z)和g(x-Ax2i-1,z)。与常网格不同,变网格的差分算子是空间变化ik=i-lcE2i-lA2iT)+c2iexp(-ikA2i)exp(+ik(5)的。给定地质模型,在确定网格剖分后就可求出各网格节点的差分算子供递推计算时调用。在交错网格技术中差分算子有两个对称点，即i+或i,不2同的对称点对应的Ai也不同。下面以对称点i+2为例求取Ai,如图1所示。对式(5)中的指数项进行泰勒展开，它的2n阶泰勒展开式为2Ai)1+ikAi+(ik)2Ae

9、xp(iki+22n-12n(kAi)2n-1+o(A(6)+i)(2n-1)!把式(6)代入式(5)并整理,可以得到2222AA)()(c1AA2ik=(c1+c2+-+c2n)+ik(c1A-c+c+i1222n2n1+c2A2+c2n2n)+2+(i)336332n-1(c1AA31+C2A2+c2n2n)+(i)n-12n-12n-1(clA-c2A+c2nA2n)12(2n-l)!(7)2n-1写成矩阵的形式1iA11-iA231iA2n-11-iA2nc1c2二2n-22n-10i(i)2n-2A12n-22n-1(i)2n-1A12n-2(i)2n-2A22n-12n-lA-(i

10、)2(i)2n-2A2n2n-1A-(i)2n2n-2(i)2n-2A2n-12n-1-(i)2n-1A2n-1c2n-1c2n00(8)第43卷第6期李振春等:一阶弹性波方程的变网格高阶有限差分数值模拟713通过解方程组(8),可以得到有限差分算子Dx的系数ci。同样地，差分算子Dz的差分系数也是通过这种方法求得。对于交错的网格，由于变量所定义的网格点的位置不同，对于每一个差分算子Dx或Dz都有两组不同的差分系数。在使用变网格时，网格步长的变化可能会导致数值反射的出现。究其原因是因为波场离散后，相速度是网格步长的函数。当相速度梯度较大时，即使速度和密度都没有变化，入射波的能量也会部分反射回来

11、，导致了数值反射现象。因为频散程度与频率有关，所以当入射波的频率较低时，产生的数值反射较高频时弱。在网格变化的区域对波场进行插值计算虽可以在某种程度上压制数值反射，但是插值算法计算量大，且效果并不令人满意，所以本文没有采用插值计算。本文提出的算法是根据各点所对应的步长的变化计算其差分系数，在一个循环内完成所有点的波场计算，不会产生数值反射，并且减小了计算量。式中:h=max(dxi,dzi)表示最大的空间网格步长；vmin是地震波速的最小值。本文所提算法需要满足的稳定性条件为Wmin(dxi，dzi)ni=1(-1)i-1Ci式中:Vmax为最大速度值；P为密度；Ci为式(3)中做2n阶空间展

12、开时的系数。4数值模拟实例为了验证提出的变网格差分算法的稳定性及有效性,设计两个模型正演实例来说明文中所述方法的实际效果。4.1。模X3000m,纵波速度为3500m/s,横0125来计算，密度为常数。P波震源激发，震源主频为30Hz,位于(1500m,1500m)处。变网格步长为10m,5m。与其进行对比的常规交错网格步长为5m。数值模拟中，差分解法的精度为o(At2,Ax10)。弹性波波场快照和单炮记录分别记录于图2和图3中。从上面的弹性波波场快照和单炮记录中可以看出,变网格的过渡区没有产生数值边界反射,说明本文提出的算法成功地消除了空间网格变化造成的数值反射,进而说明了该算法的合理性和良

13、好的稳定性。3法较好地解决了边界反射问题。使用透射边界能够吸收绝大部分边界反射的能量,剩余能量通过衰减边界再次减小,这样计算量增加不大,吸收也较彻底,可以较好地消除人为边界反射。稳定性条件是有限差分数值模拟中十分重要的问题9。经过验证,利用本文所提出的算法能够精确模拟的最大的地震波频率为fmax图2a变网格波场(左)和常规网格波场(右)水平分量快照(t=360ms)714石油地球物理勘探2019年第43卷第6期李振春等:一阶弹性波方程的变网格高阶有限差分数值模拟715412背斜模型模型大小为3000mX2000m,模型中横波速度取泊松比为0.25来计算；P波震源激发，震源子波主频为20Hz。采

14、用精度为。(At2,Ax10)的常规交错网格有限差分法和变网格有限差分法进行模拟对比。模型参数及单炮记录如图4和图5所示。从上面的模拟效果对比分析得知:在用常规交错网格有限差分法模拟的炮记录中,由于背斜的影响,自背斜顶界面反射以下出现了严重的数值频散;而用变网格算法进行模拟时,由于在背斜低速区域采用小网格步长计算,故由背斜影响导致的频散得到了较好地压制,充分显示了此算法较高的模拟精度。图4背斜模型716石油地球物理勘探2019年5结论本文在网格过渡区域,没有使用插值计算,而是根据各点所对应的步长的变化计算其差分系数,在一个循环内完成对所有点的波场计算,成功地消除了数值反射,且同波场插值计算相比

15、减少了计算量。对于常规网格算法难以解决的地下复杂构造的问题,采用变网格算法可较好地解决,且能较灵活地适应地层模型的变化,对低速体的适应性较常规网格算法有了较大改善另外,相对常规网格算法而言,变网格算法节省了计算机内存,减少了计算时间,且模拟的区域越大,节省的时间就会越多,避免了对整个模型的精细采样造成的资源浪费。本文提出的算法精度可达任意偶数阶精度,高阶变网格算法的提出,尤其是与交错网格有限差分法相结合,可以有效地提高模拟精度和计算效率。参考文献1VirieuxJ.P2media:velocity2stressfinite2differencemethod.Geo2physics,1986,5

16、1,889012JastramCandBehleA.Acousticmodelingonaverti2callyveryinggrid.GeophysProsp,1992,40,1571693JastramCandTessmerE.Elasticmodelingonagridofverticallyvaryingspacing.GeophysProsp,1994,42,357370FalkJ,TessmerEandGajewskiD.Tubewavemodel2ingbythefinite2differencemethodwithvaryinggridspacing.Prosp,1996,14

17、8,7793PitarkaA.3Delasticfinite2differencemodelingofseis2micmotionusingstaggeredgridswithnonuniformspacing.BulletinofSeismologicalofAmerica,1999,89,5468GravesRW.Simulatingseismicwavepropagationin3Delasticusingstaggered2gridfinite2differences.BullSeismSocAm,1996,86,109111077董良国,马在田,曹景忠.一阶弹性波方程交错网格高阶差分

18、解法.地球物理学报,2000,43(3):4114198廖振鹏,周正华,张艳红.,2002,45(4):534545,地球物理学报,2000,43(6):(本文编辑:冯杏芝)(上接第688页)虚线框在横向上的分布与图8中虚线框在横向上的分布比较可以得知,阻抗低值处并非都富集天然气。通过时频分析技术和阻抗反演技术相结合,最终确定储层的分布范围为图8中白色实线框内。从对实际资料的应用和分析中发现,将时频分析技术和波阻抗反演技术相结合可以更加精细地刻画礁滩相储集层的空间分布、油气的富集范围,计算油气储层的厚度,估算油气的储量,从而有效地对礁滩相储层进行预测,降低钻探风险。只对垂直剖面进行了探讨,建议

19、在实际的生产研究中也应该考虑对沿层切片、水平时间切片及波阻抗剖面和单频剖面叠合图进行分析,从而更加准确地对礁滩相油气藏进行预测。参考文献1ThomasEnglehart;石天译.墨西哥湾中部气云四分量地震数据的解释.石油勘探开发情报,2001(6):33432李等华,唐跃,殷积峰等.川东黄龙场构造上二叠统长兴组生物礁特征与潜伏礁预测.中国地质,2019,33(2):4274343CastagnaJP,SunShengjie,SiegfriedRW.Instantane2ousspectralanalysis:Detectionoflow2frequencyshad2owsassociatedw

20、ithhydrocarbons.TheLeadingEdge,2019,1201274董臣强,王军,张金伟.时频分析技术在三角洲层序分析中的应用.断块油气,2002,9(2):18205张宏彬,何樵登等.宽带约束反演.石油物探,1995,34(1):1106张霖斌，姚振兴等快速模拟退火算法及应用.石油地球物理勘探,1997,32(5):6546597王赟，王妙月，彭苏萍.地球物理随机联合反演.地球物理学报,1999,42(Supp.1):1421518寇晓丽，刘三阳.基于模拟退火的粒子群算法求解约束优化问题.吉林大学学报(工学版),2019,37(1):1361405结束语通过对礁滩相油气储集

21、层的特征分析和研究,提出了将时频分析技术和波阻抗反演技术相结合对礁滩相油气藏进行预测的方法。文中分析了单独应用两种技术对礁滩相油气藏预测的优缺点。综合两种技术的优点,可以更加有效地对礁滩相油气藏进行预测。通过对实际资料的应用发现,采用文中提出的礁滩相预测研究方法和思路,可以有效识别储集层的位置,精细刻画储集层的范围,计算储集层的厚度,减少钻井的风险,更好地指导油气开发。文中(本文编辑:冯小球)Vol.43No.6Abstracts町IVariable2gridhigh2orderfinite2differencenumericsimulationoffirst2orderelasticwave

22、equation.LiZhen2chun1,ZhangHui1andZhangHua1.OGP,2019,43(6):711716Usingvariable2gridfinite2differencealgorithmtocarryouttheseismicsimulationhasmanyu2niqueadvantagesthatmainlyaremorereasonablefordiscretizationofgeologicmodel,andrelativelyfinedivisionoflocalgridcanbeimplementedinlow2velocityzoneandcomp

23、lexstructuralarea,whichnotonlyimprovesthesimulatedprecisionandeliminatesthedispersionresultedfrominsuffi2cientsampling,butalsocanreducetherequire2mentofcomputermemory,keepingtheflexibilityofmodelcomputation.Thepaperpresentedanewelasticwavenumericsimulationmethodbasedonhigh2orderstaggered2girdtechniq

24、ue,whichreal2izedthetechniquedensifyinglocalgridbychangingthespatialstepofgrid,remedyingtshortcomingofordinaryshowedthatthethepaperischaracterizedby,improvedpreci2sionofsimulation,reducingcomputationaltimeandimprovingcomputationalefficiency.Keywords:variablegrid,staggeredgrid,high2or2derfinitediffer

25、ence,elasticwave1.CollegeofEarthResourceandInformation,Chi2naUniversityofPetroleum,DongyingCity,Shan2dongProvince,257061,ChinaElasticwavesimulationofP2andS2wavesseparationby22Dstaggeredgridandapplication.ZhangJian2lei1,TianZhen2ping1andWangCheng2xiang1.OGP,2019,43(6):717722Throughsystematicallystudy

26、ingonthemeth2odseparatingP2wavefromS2wavebysecond2orderelasticwaveequationandfirst2orderstaggeredgridalgorithm,thepaperpresentedtheformularepre2sentationsimulatingP2andS2wavesseparationinfirst2orderstaggeredgrid.TheP2andS2wavesarenaturallyseparatedbutdonotbecoupledtogetherinasimulatingprocessoftheme

27、thod,producingfullP2waveandfullS2wavesimulatedrecords.Theresultsofnumericsimulationshowedthecorrect2nessofthemethod.ThereisimportantpracticalmeaningforknowingpropagationruleofelasticwaveandguidingnumericsimulationofreservoirbystudyoffullP2waveandfullS2wavesimulation.Keywords:elasticwavesimulation,P2

28、andS2waveseparation,staggeredgrid1,TechnicalResearchCenterofGeophysicalPros2pecting,BGP,ZhuozhouCity,HebeiProvince,072751,ChinaStaggered2gridpseudo2spectrumsimulationoffrac2tureandcavereservoir.LiuJiong,MaJian2wei1,YangHui2zhu1andBaJing1.OGP,2019,43(6):723727Thepapercarriedouttheanalysisandcorrela2t

29、ionoftraditionalpseudo2spectrumandimprovedstaggered2gridpseudo2spectrumfromtheoreticalmethodcases,thesuperiorityof2method.Themoreofdifferentiallocalizationofpseudo2spectrumisshowedbyfirst2orderdifferentialcomputationalcaseofpulsese2quence;thecorrectnessofstaggered2gridpseudo2spectrummethodisprovedby

30、22Dacousticwavecomputationalcase;theprecisionofthestaggered2gridpseudo2spectrummethodbeinghigherthanthatoftraditionalpseudo2spectrumisprovedbythecomputationalcaseofelasticwaveinnon2hom2ogeneousmedium.Finally,thepropagationissueofelasticwaveinnon2homogeneousfractureandcavereservoirissimulatedbystagge

31、red2gridpseu2do2spectrummethod.Keywords:numericsimulationofwaveequation,traditionalpseudo2spectrum,staggered2gridpseu2do2spectrum,fractureandcavereservoir1.InstituteofSeismicWaveExplorationandDe2velopment,TsinghuaUniversity,BeijingCity,100084,ChinaMethodforidentifyingfluidpropertyinvolcanitereservoi

32、r.ZhangLi2hua1,PanBao2zhi1,ShanGang2yilandHeLi2.OGP,2019,43（6）:728730Thevulcanitereservoirischaracterizedbycomplexlithologyandvarietyofporestructure,whichbroughtaboutdifficultyinidentifyingande2valuatingthelithologyofvolcanitereservoir.Tak2ingdeepvolcanitereservoirinDaqingOilfieldas2019年12月第43卷第6期作者

33、介绍吕公河教授级高工,1963年生;1984年毕业于长春地质学院物探专业,1996年获得青岛海洋大学海洋地质专业硕士学位,2019年毕业于西安交通大学工程力学专业,获得博士学位;现任胜利油田地球物理勘探高级专家,长期从事滩海以及高分辨率地震勘探技术研究工作。陈浩林高级工程师,1969年生;1991年毕业于中国地质大学（武汉）地球物理专业,获学士学位;2001年毕业于中国石油大学（北京）地球物理专业,获工程硕士学位;2019年毕业于中国矿业大学（北京）地球物理专业,获博士学位。长期从事地震资料采集、处理方法及理论研究。王振华硕士,1958年生;1982年毕业于华东石油学院勘探地球物理专业,198

34、5年毕业于华东石油学院北京研究生部应用地球物理专业。主要从事地球物理勘探技术的开发与应用。现在北京华阳捷泰技术有限公司工作。单联瑜高级工程师,1959年生;1982年毕业于华东石油学院地球物理勘探专业,究和地震资料处理研究所在读博士王玉贵高级工程师;年毕业于石油大学（华东）,2019年获中国石油大学应用地球物理专业硕士学位,现为中国地质大学地球物理与信息技术学院在读博士。一直从事VSP资料处理、解释及方法研究工作。裴正林博士,副研究员,1962年生;1984年毕业于西安地质学院物探专业,2000年毕业于中国地质大学（北京），获博士学位,20002019年在中国石油大学（北京）从事博士后研究。主要从事地震波传播理论及数值模拟,井间地震层析成像,小波域地震属性提取和油气识别,以及小波理论应用研究。获省部级科技进步一等奖1项,二等奖1项,三等奖2项。合著1部,发表论文50余篇。段云卿副教授,1958年生;1982年毕业于华东石油学院物探专业,现在中国地质大学（北京）从事教学和地震资料处理及解释方法研究。已发表论文多篇。何兵寿副教授,1973年生;1996年毕业于中国矿业大学（徐州）物探专业,2002年获石油大学（北京）应用地球物理专业