传统和旋转交错网格有限差分在双相介质中模拟对比

收稿日期：XXXX－XX－XX基金项目：国家重大科学仪器设备开发专项（2012YQ030126）；国家自然基金煤炭联合项目（U1261203）；中国地质调查局项目(12120115102101)传统和旋转交错网格有限差分在双相介质中模拟对比林朋，卢勇旭（中国矿业大学（北京）煤炭资源与安全开采国家重点实验室，北京100083）摘要：以Biot双相介质模型为背景，推导了双相各向同性介质二维三分量一阶速度—应力弹性波方程方程，建立了各向同性双相介质波动方程的二维三分量有限差分格式。分别采用传统交错网格有限差分技术和旋转交错网格有限差分技术对均匀和非均匀双相各向同性介质进行了波场模拟。结果表明：旋转交错网格有限差分技术能够有效模拟双相各向同性介质中弹性波的传播情况；通过传统和旋转交错网格有限差分技术的对比，说明了旋转交错网格有限差分算法在稳定性更强，避免了插值带来的误差，是一种有效的地震波场模拟方法。关键词：双相介质；传统交错网格；旋转交错网格；波场模拟随着人们对地球内部认识的逐步加深和地球物理学的快速发展，含饱和流体的孔隙介质被认为是与含油气储层最为接近的介质模型，其是由固体骨架颗粒和孔隙中的流体（如气、水等）组成双相或多相介质。双相介质中弹性波传播规律的正确认识，对于油气的开采具有重要意义。地震波数值模拟是人们用来描述和认识地震波传播规律的有效途径。交错网格有限差分技术作为地震波场模拟的常用方法之一，在地球物理中一直被广范应用。Biot理论描述了饱和流体空隙介质中地震波的传播，是双相介质波动理论的基础[1-3]；Schmitt讨论了柱坐标系下地震波在横向各向异性介质中的传播问题[4]；Crampin通过一系列实验研究发现，双相各向异性介质中存在横波分裂现象[5]；王尚旭研究了双相介质地震波传播规律，并利用有限元法实现了双相介质地震波场模拟[6]；牟永光应用有限差分技术对孔隙各向同性介质进行了波场分析[7]；刘洋等通过虚谱法对双相各向异性介质中弹性波的传播特征[8]；王秀明等使用高阶交错网格有限差分技术实现了非均匀孔隙介质的正演模拟[9]；裴正林通过交错网格有限差分法实现了双相各向异性介质和三维横向各向同性介质弹性波的高阶波场模拟[10-11]。本文使用传统和旋转交错网格两种有限差分技术对双相介质进行了波场模拟，通过对比，说明了旋转交错网格有限差分技术的有效性和优越性。1方法原理1.1双相各向异性介质二维三分量弹性波方程根据双相介质地震波传播理论[2]，可得到饱和流体孔隙介质运动方程[11]21111221122112212222112222111211121112()()()()()()iijiijiijiijiijiijvbvbVstVbvbVst（1）（2）式中，i、j表示x,y,z三个不同方向分量，bij表示耗散系数，ij为固相应力分量，s为作用在流体上的有效应力，,,xyzν为介质固体骨架速度矢量,,,xyzVVVV为介质流体速度矢量。令为孔隙度，11和22分别表示单位体积内固体骨架和流体部分的有效质量，12为流体相对固体骨架运动时的视质量[12]。三者与固相密度s和流相密度f之间满足[2]11121222=1=sf（3）（4）由饱和流体孔隙介质传播理论，易得双相各向同性介质二维三分量一阶速度—应力表达式。令Qi表示固体骨架和流体空隙之间体积变化的耦合参数，R表示描述流体的弹性参数，即将一定体积流体注入孔隙介质体积元，为保持总体积不变时在流体上施加的一种力的度量[12]。对于双相各向同性介质，固相部分速度与应力之间的关系为111211331112++00++xxxzxxxxxzzzzzzzVVbxzttbVVttxz（5）由广义达西定律可得Biot介质流相部分速度和应力关系式122211331222+00+xxxxzzzzsVVbxttbsVVztt（6）1.2交错网格有限差分技术交错网格有限差分技术是一种有效的地震波正演模拟方法。传统交错网格是将网格剖分成整网格点和半网格点，在相邻的两个时间层上的网格点处分别定义速度和应力分量，并且在空间分布上相邻的两个时间层上的物理量恰好交错半个网格，导数值在半网格点处计算，以实现时间和空间的交错[13]。由于速度和应力的相对关系，对于不同的介质模型，在计算过程中需要对部分场量和模型参数进行插值，增加了计算误差，降低了计算精度。网格定义、波场分量及弹性参数位置如图1（a）所示。旋转交错网格由传统交错网格发展而来，通过旋转对网格进行了重新划分和定义。与传统交错网格相比，不同之处在于，在同一网格点处仅定义同一物理量（速度、应力），通过计算沿网格对角线物理量的差分来计算微分，由于速度分量和应力分量各定义整网格点和半网格点上，避免了部分场量和模型参数的插值[14]，降低了计算误差，提高了计算精度[15-16]。网格定义、波场分量及弹性参数位置如图1（b）所示[17]。（a）传统交错网格（b）旋转交错网格图1交错网格示意图1.3稳定性分析稳定性问题是地球物理学科研究中所必须解决的难题，直接关系到数值模拟方法的成败。对于传统交错网格有限差分技术，在步长相等的情况下，时间域二阶、空间2M阶的稳定性条件如下：max11/nkktVDCh（7）式中，t是时间步长，D是空间维数，h为空间步长，Vmax为最大相速度，Ck时是空间差分系数。对于旋转交错网格有限差分技术，Saenger[14]在Neumann稳定性条件下，给出了在空间步长相等时，时间域二阶、空间2M阶的稳定性条件为：max11/nkktVCh（8）式中参数意义同（7）式。通过对比两种交错网格的稳定性条件可知，相较于传统交错网格，旋转交错网格有明显的优点，在理论上有更宽松的稳定性条件。2数值模拟及分析2.1均匀双相各向同性介质波场模拟为观察旋转交错网格有限差分技术在双相各向同性介质中的模拟效果，采用均匀各向同性介质为背景，模型大小为256256mm，网格间距1m，采样间隔0.1ms，震源采用Ricker子波，位于模型中间位置，加载于固相正应力处。模型如图2所示。取t=40ms时的波场快照，使用传统和旋转交错网格有限差分技术所得波场快照如图3、4所示。图2均匀双相各向同性介质模型（a）固相X分量（b）流相X分量（c）固相Z分量（d）流相Z分量图3均匀双相各向同性介质传统交错网格数值模拟波场快照（a）固相X分量（b）流相X分量（c）固相Z分量（d）流相Z分量图4均匀双相各向同性介质旋转交错网格波场快照结合图3和图4可知，使用传统和旋转交错网格有限差分技术对均匀双相各向同性介质模拟时，均可以得到清晰的波场快照，均匀双相介质中存在的快纵波、慢纵波和快横波都清晰可见，波场特征十分明显，波场分布符合地震波传播规律。说明了旋转交错网格有限差分技术可以较好的模拟地震波在均匀双相介质中的传播规律，是一种有效的地震波场模拟方法。2.2非均匀双相各向同性介质波场模拟为进一步验证传统和旋转交错网格有限差分技术在双相各向同性介质模拟中的差异，以非均匀双相各向同性介质为背景，模型分为A、B两个区域，B是一密度和速度极小的地质异常区块，其参数与均匀情况下相同。模型如图5所示。对非均匀双相各向同性介质分别使用传统交错网格技术和旋转交错网格技术进行波场模拟，取t=40ms时的波场快照进行对比，所得波场快照如图6、7所示。图5非均匀双相各向同性介质模型（a）固相X分量（b）流相X分量（c）固相Z分量（d）流相Z分量图6非均匀双相各向同性介质传统交错网格数值模拟波场快照（a）固相X分量（b）流相X分量（c）固相Z分量（d）流相Z分量图7非均匀双相各向同性介质旋转交错网格波场快照从图6可知，当使用传统交错网格有限差分技术对双相介质进行波场模拟时，除地质异常体外，波场均可以正常传播，符合双相介质波动理论，而在地质异常体附近，由于波场分量和模型参数的插值，导致波场出现了振幅异常现象，不能正确表示波场传播情况；从图7可以观察到，使用旋转交错网格有限差分技术模拟时，波场传播均符合双相介质地震波传播规律，在地质异常体附近有反射波和绕射波产生，可以正确表示波场在双相介质中的传播。3结论（1）相对传统交错网格，旋转交错网格有限差分技术在理论上具有更宽松的稳定性条件。（2）在均匀双相各向同性介质的波场模拟中，旋转交错网格有限差分技术能够得到清晰的波场快照，且波场特征十分明显，可以有效模拟双相各向同性介质中弹性波的传播情况。（3）对存在地质异常体的非均匀双相各向同性介质，传统交错网格有限差分技术稳定性不足，具有明显的局限性；而旋转交错网格有限差分技术适应性较强，稳定性更好，是一种有效的地震波场数值模拟方法。参考文献[1]BiotMA.Theoryofpropagationofelasticwavesinafluid-saturatedporoussolid,I:Low-frequencyrange[J].TheJournaloftheAcousticSocietyofAmerica,1956a,28:168-178.[2]BiotMA.Theoryofpropagationofelasticwavesinafluid-saturatedporoussolid,II:High-frequencyrange[J].TheJournaloftheAcousticSocietyofAmerica,1956b,28:179-191.[3]BiotMA.Mechanicsofdeformationsandacousticpropagationinporousmedia[J].J.Appl.Phys,1962a,33:1482-1498.[4]SchmittPD.Acousticmultipolelogginginthetransverselyisotropicporoelasticformations.TheJournaloftheAcousticSocietyofAmerica,1989,86(6):2397-2421.[5]CrampinS,YedlinM.Shear-wavesingularitiesofwavepropagationinanisotropicmedia.JGeophys,1981,49:43-46.[6]王尚旭.双相介质中弹性波问题有限元数值解和AVO问题[D].北京：中国石油大学（北京）博士学位论文，1990.[7]牟永光，裴正林.三维复杂地震数值模拟[M].石油工业出版社，2005-1.[8]刘洋，李承楚.双相各向异性介质中弹性波传播特征研究[J].地震学报，1999，21(4):367-373.[9]杨顶辉.双相各向异性介质中弹性波方程的有限元解法及波场模拟[J].地球物理学报，2002，45(4):575-583.[10]裴正林.三维各向同性介质弹性波方程交错网格高阶有限差分法模拟[J].石油物探，2005，44(4):308-316.[11]裴正林.双相各向异性介质弹性波传播交错网格高阶有限差分法模拟[J].石油地球物理勘探，2006，41(2):137-143.[12]孙卫涛，杨慧珠.双相各向异性介质弹性波场有限差分正演模拟[J].固体力学学报，2004，25(1):21-28.[13]孙瑞艳.TTI介质旋转交错网格有限差分及其组合边界条件[D].东营：中国石油大学（华东）硕士学位论文，2010.[14]SaengerE.H.,GoldN.,ShapiroS.A