双相TI介质中弹性波交错网格高阶有限差分法数值模拟

尹学爱，邱光辉

（1.滨州学院，山东 滨州 256603；2.中化地质矿山总局 山东地质勘查院，山东 济南 250013）

双相TI介质中弹性波交错网格高阶有限差分法数值模拟

尹学爱1，邱光辉2

（1.滨州学院，山东 滨州 256603；2.中化地质矿山总局 山东地质勘查院，山东 济南 250013）

应用双相介质波动方程，推导了双相横向各向同性介质（TI）中波动方程的有限差分格式，对双相TI介质中弹性波有限差分数值进行了模拟。结果表明，弹性波在双相TI介质中传播时，除了存在常规的快纵波（qP1）和横波以外，还存在慢纵波（qP2）。并且慢纵波的速度明显小于快纵波，而且受耗散系数的影响衰减地很快，所以在实际中很难观测到慢纵波。快纵波在固相和流相中相位相同，而慢纵波在固相和流相中的相位相反。慢纵波在流相中振幅大，而在固相中的振幅较小。

双相横向各向同性介质；交错网格；有限差分法；弹性波

0 前言

经典的地震波理论只适合于研究固体或流体单相介质中地震波传播规律。然而，无论是砂岩储层还是碳酸盐储层，或者是海底沉积物，都是由固体和流体两种部份组成，即由固体和流体组成的双相介质或多相介质。由于地震波在双相或多相介质传播时，其传播规律与在单相介质中不同，因此，根据双相介质地震波理论，通过正、反演数值模拟，才能可靠准确地确定储层的厚度、空间展布，以及储层的孔隙度、渗透率和油气饱和度等储层参数，所以对油气勘探开发等实际工作具有十分重要的意义。

但是，该领域的研究还存在大量凾待解决的问题，其中最突出的是计算速度和稳定性，国内、外学者仍在进行对于双相介质模型和高性能算法的研究。Gassmann［1］提出了关于弹性波在多孔介质中的传播理论，并建立了著名的Gassmann方程；Biot［2～4］根据潮湿土壤的电位特性和声学中声波的吸收特性，发展了Gassmann的流体饱和多孔隙双相介质理论，奠定了双相介质波动理论的基础，但由于Biot双相介质波动方程在复杂地质环境下没有确定的解析解，只能通过数值方法求得，所以国内、外对Biot双相介质的数值模拟做了大量研究；Zhu和McMechan［5］用有限差分法模拟了双相介质；N.Dai［6］等用一阶应力～速度波动方程模拟了各向异性双相介质；孙卫涛和杨慧珠［7］采用交错网格技术，建立了各向异性孔隙介质波动方程的高精度差分格式，并对这类差分格式的频散特性和稳定性作了详细分析讨论，解决了计算稳定性和边界反射问题；裴正林［8、9］基于Biot理论给出了三维双相各向异性介质应力～速度弹性波方程交错网格任意偶阶精度有限差分解法，对三维双相各向异性介质中弹性波场进行了模拟；王秀明等［10、11］利用基于Biot理论的非均匀孔隙弹性介质的高阶交错网格有限差分算法，模拟了地震波在其中的传播。

作者在Biot理论的基础上，建立了Biot双相介质波动方程交错网格差分格式，通过数值模拟对双相横向各向同性介质中的弹性波传播特征进行了观测分析。

1 原理

由饱和流体孔隙介质的Biot理论可知，声波传播的Biot线性理论基于以下六种假设：

（1）波长远大于颗粒尺寸，颗粒粒径又大于孔隙尺寸。

（2）固体与流体的质点运动相对较小，在线性弹性理论范围内。

（3）流体相在整个介质中是连续的，而不连通的孔道可以视为固体骨架。

（4）固体骨架认为是均匀和统计各向同性的，并且忽略了重力的影响。

（5）流体流动是Poiseuille类型，岩石骨架与孔隙流体之间存在相对运动，且流体的流动满足广义达西定律（这一限制只适用于低频情况），并忽略热效应。

（6）介质是完全饱和的。

基于上述假设，可以得到下面相关方程。

1.1 应变～位移关系

若用u3×1＝ ［uxuyuz］T表示固相的平均位移矢量，U3×1＝ ［UxUyUz］T表示流相的平均位移矢量。则固相应变分量可以用位移分量表示，固相正应变为：

固相切应变为：

固体部份的体应变θ，可由位移矢量场的散度表示：

同样地，流相部份的体应变ε为：

1.2 应力～应变关系

根据广义Hooke定律，双相TI介质中的应力～应变关系（即本构方程）表示为

其中 σ6×1＝［σxxσyyσzzτyzτxzτxy］T是应力矩阵的列矩阵表达式；e6×1＝［exxeyyezzeyz

exzexy］T是应变矩阵的列矩阵表达式；S是流相应力；ε是流相体应变；D6×6是固相的物性矩阵，它包含与固相有关的弹性参数；R1×1是与流相有关的弹性参数；Q1×6是与固相和流相之间耦合有关的弹性参数，且Q6×1＝Q1×6T。

对于横向各向同性双相介质，

其中 d66＝0.5（d11－d12），在D6×6中共有五个独立的物性参数，Q1×6为：

其中 Q1＝Q2，共有二个独立的物性参数，R1×1＝R中只有一个独立物性参数。

1.3 应力～位移关系

根据分析力学中的Lagrange方程，可导出双相TI介质中当有耗散存在时的平衡运动方程式：

式中

对于各向同性双相介质b11＝b22＝b33，对于垂向横向各向同性双相介质b11＝b22。在式（8）中ρ11、ρ12、ρ22为质量系数，ρ11表示单元体中固相相对于流相运动时固相部份总的等效质量；ρ22为流相相对于固体相运动时流体部份总的等效质量；ρ12表示流相和固相之间的质量耦合系数，是粘滞、摩擦等效应的综合反映，又称为视质量。

2 有限差分算法的实现

目前有限差分法在数值模拟中具有其它方法不可替代的优势，差分离散的途径主要有两种：

（1）对单变量的二阶波动方程直接利用时空二阶中心差分离散求解（规格网格差分）。

（2）将以位移表示的二阶波动方程化为以应力和质点速度表示的一阶双曲线方程组，用应力和速度的交错网格求解，交错网格是一种较为先进的差分格式。

作者在本文利用交错网格有限差分法原理，实现了双相介质相关微分方程向差分方程的转化，为下一步数值模拟打下了基础。

双相介质的运动平衡方程式为：

为了得到一阶速度～应力方程，首先要将固相位移分量和流相位移分量各自满足的方程分离开来。

令式 （11）×ρ22－式（12）×ρ12，并整理得到固相位移分量公式为式（13）。

令式（11）×ρ12－式（12）×ρ11，并整理得到流相位移分量公式为式（14）。

设固相速度矢量为v3×1＝ ［vxvyvz］，流相速度矢量为V3×1＝ ［VxVyVz］，将速度对时间的一阶导数替换成公式（13）和式（14）中位移对时间的二阶导数，我们可以进一步得到固相速度和流相速度的一阶时间导数方程形式。

下面推导应力分量的一阶方程形式，双相介质的本构方程为：

将方程式（16）两边同时对时间求导，我们就可以得到固相应力和流相应力的一阶时间导数方程形式。对于固相有式（17）。

对于流相有式（18）。

其中

我们建立了双相介质一阶速度～应力波动方程，接下来就可以建立双相介质交错网格差分离散格式［21］。首先看交错网格离散格式中双相介质相应的波场分量和弹性参数的空间位置，见下页表1和下页图1。

表1 弹性波场分量和弹性参数位置Tab.1 Elastic wave field components and elastic parameter

图1 三维双相介质交错网格示意图Fig.1 3Dtwo－phase medium staggered－grid

（1）固相位移各分量的交错网格差分离散格式为式（20）。

（2）流相位移各分量的交错网格差分离散格式为式（21）。

其中

（3）固相应力各分量的交错网格差分格式为式（22）。

（4）流相应力分量的交错网格差分格式为式（23）。

其中

上面离散格式中，DΔx、DΔy、DΔz是空间求导算子，分别代表了对x、y、z的导数；n表示时间上的离散点；i表示x方向上的离散点；j表示y方向上的离散点；k表达z方向上的离散点。

若令y方向的空间导数为零，则差分格式退化为x～z平面上的二维交错网格差分格式；若令y分量上的波场值都为零，则退化为只有x、z两个分量的差分格式。作者在本文的模拟都是在二维二分量上求取的。

3 数值模拟结果及分析

二维均匀双相横向各向同性介质模型弹性参数见表2。模型大小为400m×400m，震源为爆炸点源，高斯子波主频为30Hz，震源位于模型中心，时间步长为10－4s，网格大小为2.5m，交错网格差分精度为o（Δt2，Δx6）。

图2和图3（见下页）是利用表2中模型2所得到的二维均匀双相横向各向同性（VTI）介质弹性波波场模拟结果。

下页图4和后面的图5是利用表2中模型1所得到的二维均匀双相横向各向同性（HTI）介质弹性波波场模拟结果。

表2 均匀双相介质物性参数Tab.2 Parameters of medium physical

图2 均匀双相介质中弹性波快照Fig.2 Snapshots of TI media

图5 弹性波场快照与弹性波场记录（t＝0.09s）Fig.5 Snapshots and time series for elastic wave field

弹性波在双相各向异性介质中传播时存在四种波，即快纵波、快横波、慢横波和慢纵波。由于模拟为二维，且倾角和方位角为零度或者90°时，只能同时看到三种波，所以从图2中可以看到快纵波qP1、快横波qS1（SH波）和慢纵波qP2。从上页图4中可以看到快纵波qP1、快横波qS2（SV波）和慢纵波qP2。这三种波具有各向异性特性，只有慢纵波属于第二类波，其余为第一类波。

图2是利用表2中模型2所得到的，二维均匀双相横向各向同性（VTI）介质弹性波波场模拟结果。从图2中可以看到，弹性波在双相横向各向同性（VTI）介质中传播时存在三种波，即快纵波qP1、快横波qS1（SH）、慢纵波qP2。从图2中可知，当倾角和方位角都为零度时，快横波和慢横波在固流相中形成横波分裂盲点，速度一样，无法区分。在固相x、z分量上含有快纵波qP1、快横波qS1和慢纵波qP2；在流相x、z分量上也含有快纵波qP1、快横波qS1和慢纵波qP2，其中快纵波qP1和快横波qS1在z分量的能量比x分量要强一些。从图2中还可以看出，快纵波qP1和快横波qS1在方位上各向异性，并且在一个方位上是耦合的，同时可以看出横波波前产生尖角和三分叉现象。这一特殊现象使波场复杂化，并且容易引起误解。

根据物性参数与速度的关系，得到速度，用速度验证波前面，得知本模拟结果是正确的。

图3是z方向上双相横向各向同性（VTI）介质中弹性波波场快照与弹性波场记录的对比。其中检波器道间距为5m和最小偏移距为零，检波器排列在震源上方37.5m处水平接收。从图3中可以看出，各种波对应良好，只是在流相中的快横波有频散现象，这应该是快纵波与快横波耦合造成的。

图4是利用表2中模型1所得到的二维均匀双相横向各向同性（HTI）介质弹性波波场模拟结果。由图4可知，在固相x、z分量上含有快纵波qP1、快横波qS2（SV波）和慢纵波qP2，但是慢纵波qP2在x分量上能量很弱，在z分量上基本看不到。在流相x、z分量上也含有快纵波qP1、快横波qS2和慢纵波qP2，快横波qS2和慢纵波qP2是耦合的。从图4整体上看，快纵波qP1能量在z方向上要比x方向上要小一些。

图5是z方向上双相横向各向同性（HTI）介质中弹性波波场快照与弹性波波场记录的对比。其中检波器道间距为5m和最小偏移距为零，在震源深度水平接收。从图5中可以看出，波场快照与波场纪录对应良好。从图5中可知，当倾角和方位角都为零度时，快横波和慢横波在固流相中形成横波分裂盲点，速度一样，无法区分。慢纵波qP2在x、z方向上能量均较固相中要强很多，这进一步证明了在流相中更容易观测到慢纵波qP2。慢纵波qP2的振幅大小，主要受耗散系数b大小的影响。耗散系数越小，慢纵波越明显，反之，则越不明显。

根据物性参数与速度的关系，可得到速度，用速度验证波前面，得知本模拟结果是正确的。

4 结论

（1）根据Biot理论，作者对双相介质理论进行了讨论，推导了双相横向各向同性介质波动方程和双相横向各向同性速度～应力弹性波方程，通过交错网格高阶有限差分数值模拟，研究了弹性波在双相横向各向同性介质中的传播情况。

（2）在双相介质VTI和HTI中，除了存在常规的纵波（快纵波）和横波以外，还存在第二类纵波（慢纵波），慢纵波的速度明显小于快纵波，而且受耗散系数的影响衰减地很快，所以在实际中很难观测到第二类纵波。

（3）在双相介质中，快纵波在固相和流相中相位相同，而慢纵波在固相和流相中的相位相反。慢纵波在流相中振幅大，而在固相中的振幅较小，即在相同条件下，从流相中更容易观测到慢纵波。

从数值模拟的结果可以看出：交错网格高阶有限差分法是一种模拟精度高、计算效率好，实用性强的弹性波场数值模拟方法。

［1］ GASSMANN F.Elastic waves through a packing of spheres［J］.Geophysics，1951，16（4）：673.

［2］ BIOT M A.Theory of propagation of elastic waves in fluid filled porous solid I.higher requency range［J］.J.Acou.Soc.Am.，1956，28（2）：179.

［3］ BIOT M A.Mechanics of deformation and acoustic propagation in porous media［J］.J.appl.Phys.，1962，33（4）：1482.

［4］ BIOT M A.Generalized theory of acoustic propagation in porous dissipative media［J］.J.Acoust.Soc.Am.，1962，34（5）：1254.

［5］ ZHU X，MCMECHAN G A.McMechan.Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory［J］.Geophysics，1991，56（3）：328.

［6］ DAI N，VAFIDIS A，KANASEWICH E R.Wave propagation in heterogeneous，porous media：avelocity－stress，finite－difference method［J］.Geophysics，1995，60（2）：327.

［7］ 孙卫涛，杨慧珠.双相各向异性介质弹性波场有限差分正演模拟［J］.固体力学学报，2004，25（1）：21.

［8］ 裴正林.三维双相各向异性介质弹性波方程交错网格高阶有限差分法模拟［J］.中国石油大学学报，2006，30（2）：16.

［9］ 裴正林.双相各向异性介质弹性波传播交错网格高阶有限差分法模拟［J］.石油地球物理勘探，2006，41（2）：137.

［10］王秀明，张海澜，王东.利用高价交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播［J］.地球物理学报，2003，46（6）：842.

［11］王东，张海澜，王秀明.部分饱和孔隙岩石中声波传播数值研究［J］.地球物理学报2006，49（2）：524.

P 631.4

A

10.3969／j.issn.1001－1749.2012.05.06

1001—1749（2012）05—0533—08

滨州学院科研基金（BZXYL1107）

2011－10－21 改回日期：2012－05－30

尹学爱（1979－），女，汉，山东邹平人，教师，目前主要从事弹性波数值模拟方面的研究。