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

收稿日期：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() ()() () ()() i ij i ij i ij i ij i ij i ij v b v b V s t V b v b V s tρρρρρρρρρτρρρρρρρρρτ∂-=+-++∇-∇∂∂-=+-++∇-∇∂ （1）（2） 式中，i 、j 表示x,y,z 三个不同方向分量，b ij 表示耗散系数，ij τ为固相应力分量，s 为作用在流体上的有效应力，(),,x y z υυυT =ν为介质固体骨架速度矢量,(),,x y z V V V T =V 为介质流体速度矢量。

令φ为孔隙度，11ρ和22ρ分别表示单位体积内固体骨架和流体部分的有效质量，12ρ为流体相对固体骨架运动时的视质量[12]。

三者与固相密度s ρ和流相密度f ρ之间满足[2]()11121222 = 1 = sf ρρφρρρφρ+-+ （3）（4） 由饱和流体孔隙介质传播理论，易得双相各向同性介质二维三分量一阶速度—应力表达式。

令Q i 表示固体骨架和流体空隙之间体积变化的耦合参数，R 表示描述流体的弹性参数，即将一定体积流体注入孔隙介质体积元，为保持总体积不变时在流体上施加的一种力的度量[12]。

对于双相各向同性介质，固相部分速度与应力之间的关系为()()()()111211331112++00++xx xz x x x x xz zz z z z z V V b x z t t b V V t t x z ττρυρυττρυρυ∂∂∂∂⎡⎤⎡⎤⎡⎤-⎢⎥⎢⎥⎢⎥⎡⎤∂∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥∂∂∂∂⎢⎥⎢⎥⎢⎥⎣⎦-⎢⎥⎢⎥⎢⎥∂∂⎣⎦⎣⎦∂∂⎣⎦ （5） 由广义达西定律可得Biot 介质流相部分速度和应力关系式 ()()()()122211331222+00+x x x x z z z z s V V b x t t b s V V z t t ρυρυρυρυ∂∂∂⎡⎤⎡⎤⎡⎤-⎢⎥⎢⎥⎢⎥⎡⎤∂∂∂=-⎢⎥⎢⎥⎢⎥⎢⎥∂∂∂⎢⎥⎢⎥⎢⎥⎣⎦-⎢⎥⎢⎥⎢⎥∂∂∂⎣⎦⎣⎦⎣⎦ （6）1.2 交错网格有限差分技术交错网格有限差分技术是一种有效的地震波正演模拟方法。

传统交错网格是将网格剖分成整网格点和半网格点，在相邻的两个时间层上的网格点处分别定义速度和应力分量，并且在空间分布上相邻的两个时间层上的物理量恰好交错半个网格，导数值在半网格点处计算，以实现时间和空间的交错[13]。

由于速度和应力的相对关系，对于不同的介质模型，在计算过程中需要对部分场量和模型参数进行插值，增加了计算误差，降低了计算精度。

网格定义、波场分量及弹性参数位置如图1（a ）所示。

旋转交错网格由传统交错网格发展而来，通过旋转对网格进行了重新划分和定义。

与传统交错网格相比，不同之处在于，在同一网格点处仅定义同一物理量（速度、应力），通过计算沿网格对角线物理量的差分来计算微分，由于速度分量和应力分量各定义整网格点和半网格点上，避免了部分场量和模型参数的插值[14]，降低了计算误差，提高了计算精度[15-16]。

网格定义、波场分量及弹性参数位置如图1（b ）所示[17]。

（a ）传统交错网格（b ）旋转交错网格 图1 交错网格示意图 1.3 稳定性分析稳定性问题是地球物理学科研究中所必须解决的难题，直接关系到数值模拟方法的成败。

对于传统交错网格有限差分技术，在步长相等的情况下，时间域二阶、空间2M 阶的稳定性条件如下：max 11/n k k tV C h =∆⎫≤⎪∆⎭ （7）式中，t ∆是时间步长，D 是空间维数，h ∆为空间步长，V max 为最大相速度，C k 时是空间差分系数。

对于旋转交错网格有限差分技术，Saenger [14]在Neumann 稳定性条件下，给出了在空间步长相等时，时间域二阶、空间2M 阶的稳定性条件为：max 11/n k k tV C h =∆⎛⎫≤ ⎪∆⎝⎭∑ （8）式中参数意义同（7）式。

通过对比两种交错网格的稳定性条件可知，相较于传统交错网格，旋转交错网格有明显的优点，在理论上有更宽松的稳定性条件。

2数值模拟及分析2.1 均匀双相各向同性介质波场模拟为观察旋转交错网格有限差分技术在双相各向同性介质中的模拟效果，采用均匀各向同性介质为背景，模型大小为256256m m ，网格间距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]Biot M A. Theory of propagation of elastic waves in a fluid-saturated porous solid, I: Low-frequency range [J]. The Journal ofthe Acoustic Society of America, 1956a, 28:168-178.[2]Biot M A. Theory of propagation of elastic waves in a fluid-saturated porous solid, II: High-frequency range [J]. The Journal ofthe Acoustic Society of America, 1956b, 28:179-191.[3]Biot M A. Mechanics of deformations and acoustic propagation in porous media [J]. J. Appl. Phys, 1962a, 33: 1482-1498.[4]Schmitt P D. Acoustic multipole logging in the transversely isotropic poroelastic formations. The Journal of the AcousticSociety of America, 1989, 86(6): 2397-2421.[5]Crampin S, Yedlin M. Shear-wave singularities of wave propagation in anisotropic media. J Geophys, 1981, 49: 43-46.[6]王尚旭. 双相介质中弹性波问题有限元数值解和A VO问题[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]Saenger E. H., Gold N., Shapiro S. A.. Modeling the propagation of elastic waves using a modified finite-difference grid[J].Wave Motion, 2010, 31:77-92.[15]王亚妮，李长江，李庆春. 旋转交错网格VTI介质波场模拟与波场分解[J]. 物探化探计算技术，2015，37(2):198-202.[16]李长江，李庆春，王亚妮. 旋转交错网格TTI介质波场模拟与波场分解[J]. 物探与化探，2015，39(3):553-557.[17]Bohlen T, Saenger E H. Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves [J].Geophysical Prospecting, 1995, 43(6):805-829.The contrast of simulation in the two-phase media betweenthe traditionaland rotated staggered gridLIN Peng, Lu Yong-Xu(State Key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology(Beijing),Beijing 100083,China)Abstract:Starting from the model of the Biottwo-phase media theory, it is deduced that the 2D/3C velocity-stress elastic wave equation of the two-phase isotropic media and the finite difference time domain scheme is set up. The wave field simulation is based on the two-phase homogeneous and inhomogeneous isotropic media using the traditional staggered-grid method and the rotated staggered-grid method. The results show that the rotated staggered grid difference scheme performs perfect for the two-phase isotropic media. By comparing rotated stagger-grid and traditional stagger-grid, it can be proved that the finite difference algorithm of rotating staggered-grid is more stable, avoiding errors caused by interpolation, so the rotated staggered-grid difference scheme is a numerical simulation of seismic wave field method with a strong effectiveness.Key words: two-phase media; traditional staggered-grid; rotated staggered-grid; wave field simulation作者简介：林朋（1990—），男，在读硕士研究生，地球探测与信息技术专业，中国矿业大学（北京）煤炭资源与安全开采国家重点实验室。