基于Biot机制的孔隙介质地震波传播特征研究

李倩宇，李红星，黄 宇，章晨望，徐文斌，樊嘉伟，梅振繁

（东华理工大学地球物理与测控技术学院，江西南昌330013）

为了满足日益增长的勘探需求，以及研究地震波传播过程中的衰减及频散等问题，国外早在20世纪中叶开展有关双相孔隙介质的相关科研工作。Gassmann于1951年提出了Gassmann理论[1]，将孔隙度这个概念引入；1956年Biot建立了Biot理论，并首次预言慢纵波的存在[2-3]；Panneton等人对其进行改进，实现了Biot模型的三维波场正演模拟[4]；Stoll等人对Biot模型进行改进，建立了Biot-Stoll模型[5]；Dvorkin等人将Biot流动机制和喷射流动机制综合考虑，即结合固流相相互作用的力学机制，提出了BISQ（Biot-squirt）模型[6]，其能够较好反映孔隙介质中弹性波的传播规律；Nicholas Chotiros在Biot和BISQ模型的基础上建立了BICSQS模型[7]。与国外相比，国内对双相孔隙介质的研究开展较晚，门福录于1965年最早开始相关研究[8]；随后唐应吾、刘银斌、牟永光等人也相继展开研究，这一系列孔隙介质理论为油气等资源的勘探提供了极其重要的理论支持[9]。本文基于Biot机制对地震波的传播进行数值模拟，取得较好效果，之后对其波场特征进行分析。并且基于Biot模型探讨波速、衰减与频率的相关性。

1 Biot模型地震波传播数值模拟

1.1 交错网格差分方法

平衡方程（应力部分）：

几何关系（应变部分）：

应力与应变的关系：

二维下，弹性介质中质点在x、z轴方向上的速度分量如果用Vx、Vz表示的话，那么，在外力为零时的介质部分：

将式（2）和Vx、Vz代入式（3）并对时间t求导可得一阶微分方程组：

将方程组（4）和（5）联合，便可得到质点速度和应力表示的一阶弹性波波方程：

在时间上，使用二阶中心差分对式（6）进行差分离散；在空间上，采用任意高阶差分：

其中：

其余各式同理，其中为交错网格差分系数。

从交错网格差分格式可看出，Vx、Vz在 (k+1/2)Δt时取得式(k-1)Δt时刻 σxx、σzz、σzx的值，而应力张量Vx、Vz则在 kΔt时取得式在 (k+1/2)Δt的值。在程序设计时，这点显得尤为重要，理解交错网格的意义不仅仅是空间上的交错，同时还有时间上的交错，这才是交错网格最特别的地方。

1.2 完全匹配边界条件

在利用计算机进行波动方程数值模拟时，有个不可避免的问题，那就是边界问题，因为需要在空间域选取一个有限范围。学者们通常会在有限区域的外围增加边界吸收条件，其目的是尽可能利用有限计算区域，实现更为真实的无限空间的波动特征，那么如何有效设置人为边界对数值模拟的可靠程度及精确度至关重要。

在理论意义上，可以完全吸收任意角度、任意频率入射的波动则称为完全匹配层。由JP.Berenger在计算电磁场有限域波动问题中首先引入。迄今为止，比较系统的基于二阶线性双曲型偏微分弹性动力方程是仍然缺少的。在这给出Collino的二维各向同性介质PML吸收边界条件与交错网格差分格式，这很重要，它能够在各向异性介质弹性波场模拟中运用，得到的效果比较好。

下列方程组为波动控制方程的一种形式，可以根据它得到不同形式的波动模型：

式中：A、B——m×m维矩阵；

ν——m维向量。

引入一个新变量，能够将待求解问题表达式改写成如下格式：

其中，‖为只保留与交界面方向平行的偏导数分量。如对y轴求偏导（而⊥则为只保留x轴的偏导分量）。

设下轴方向上的阻尼因子为d(x)，另外有个新变量，它对与交界面成法向的波动分量上施加阻尼，使左半空间的解答满足系统在右半空间新波动的解答，这样一个波动变量u同时被定义[10]。

于是，使u满足下列控制方程是非常重要的：

和：

在频域内，由（11）式可以得到：

由（12）式可以得到：

对自变量作复变换如下：

现在利用平面谐波解来研究该模型特性，设：

其中，iν0w-iAν0kx-iBν0ky=0 ，用相同方法，能够找到式（12）的一组特解：

式（12）的解，u‖和u⊥需满足条件如下：

取a‖+a⊥=ν0，即：

平面波动系统（11）的解可以写成如下形式：

并且满足下列条件：

（1）左半空间内，u≡v x≤0，可以保证在界面处无反射；

（2）右半空间内，u因为受到阻尼影响，发生衰减；

（3）完全匹配吸收层内，阻尼系数如下所示：

1.3 交错网格算法流程图

对有限差分法的技术路线（见图1）和实现步骤进行研究后，能发现有限差分法的实现相对简单并且模拟出的效果也比较好，因此本文在对其深入研究的基础上，得到了令人满意的效果。

图1 交错网格差分算法流程图[10]

1.4 Biot地震波传播数值模拟

根据Biot模型建立相关地质模型，数值模拟结果见图2。

从图2中的模拟结果可以清晰看出快纵波、慢纵波及横波。对固流两相的波场快照进行分析，可以得到两者快纵波的波形相位具有一致性，与单相介质中的纵波类似。而固流两相中慢纵波的波形呈现出相位相反的特征。

2 Biot模型波速和衰减特征分析

2.1 孔隙介质相速度和衰减

Biot机制基于孔隙介质的双相特性，提出了流体的粘滞性对孔隙流体相对运动产生影响，并造成弹性波在孔隙介质内传播过程中逐渐衰减。Biot理论模型利用变分原理得到的波动方程来自流体孔隙介质中波的传播，之后对地震波在流体孔隙介质中的传播进行理论预测，发现存在3种体波：快纵波、慢纵波及横波，三者中快纵波频散相对较小，而慢纵波在传播过程中衰减与耗散严重。

Stoll用势Φs、Φf、Ψs、Ψf定义固体骨架的位移向量u及液体的位移向量U：

图2 孔隙介质地震波场快照

ϕ为沉积物孔隙度。Biot关于标量势的方程（用平面波解exp[i(k·x-ωt)]）：

式中：k——波数

α——构造常数或曲度因子；

κ——渗透率；

η——孔隙液体粘滞系数；

ρs——沉积物颗粒密度；

ρf——孔隙液体密度。

并有：

式中：f——频率；

c——波速；

ρ——沉积物总体密度；

F——随频率增高的泊肃叶流偏移；

Kf——孔隙液体体积模量；

Ks——沉积物颗粒体积模量；

μ——框架剪切模量。

假定孔隙是圆筒状，Biot可得F的表达式：

式中：a——孔隙尺寸参数。

由上述方程可得：

其中，O=MH-C2，P=2BC-EM-AH ，Q=

由方程（23）可得：

其中：

则纵波相速度和衰减可表示为：

在Kb=0，μ=0的特殊情况时（海底沉积物的这两项值比 Ks小得多，大多情况下可以忽略），C=M=H ，形成EDF模型：

2.2 波速和衰减数值分析

通过对孔隙度、粘滞系数、孔隙气体含量等参数依次进行改变，分析其纵波衰减和相速度随频率发生的变化，可以就上述参数对地震波衰减特征的影响进行探讨，并且在理论上对地震波衰减特征进行分析。

起初对Biot模型中孔隙度、孔隙气体含量、粘滞系数对速度的影响进行探讨。数值分析中所用的频率范围为1～108Hz，孔隙度依次选择20%、40%及60%，分别选取0.1Pa·s、0.01Pa·s和0.001Pa·s作为粘滞系数，气体含量则取0、0.1及0.2，为了研究它们对地震波衰减的影响，因此对这些参数值进行改变。

从图3～图5中分析可知，Biot模型速度在低频段上基本保持稳定，不过也存在着不足之处，那就是频散相对来说较明显；同时，Biot模型速度在频率增大的同时也在增大，这能从速度与频率的关系中分析得到。频段相同的条件下，孔隙度、孔隙气体和粘滞系数都与速度成反比。不同于速度的是，孔隙度、孔隙气体含量和粘滞系数对衰减几乎没有影响，并且衰减与频率成正比，103Hz作为频率的一个分界点，衰减在不超过这个频率的范围内稳定增大，超过这个频率后，衰减增长的速度变缓慢。

3 结论

（1）本文利用较易实现的有限差分法数值模拟Bi⁃ot模型地震波场，快纵波、慢纵波及横波都可以从模拟结果中清晰看出。通过比较流相与固相的波场快照，可以看出两者快纵波的波形相位一致，而慢纵波在固

图3 不同粘滞系数下频率与衰减和速度关系

图4 不同孔隙度下频率与速度和衰减关系

图5 不同孔隙气体含量下频率与速度和衰减关系

相与流相中的波形则表现出相位相反的特征。在固相波场中，横波波形最明显，能量最强；在流相波场中，慢纵波比快纵波及横波更清晰，这是不言而喻的。

（2）分析Biot模型速度与衰减的关系，显而易见，Biot模型速度在低频段上基本保持稳定，不过也存在着不足之处，那就是频散相对来说较明显；同时，Biot模型速度在频率增大的同时也在增大，这能从速度与频率的关系中分析得到。频段相同的条件下，孔隙度、孔隙气体和粘滞系数都与速度成反比。不同于速度的是，孔隙度、孔隙气体含量和粘滞系数对衰减几乎没有影响，并且衰减与频率成正比，103Hz作为频率的一个分界点，衰减在不超过这个频率的范围内稳定增大，超过这个频率后，衰减增长的速度变缓慢。

[1]Gassmann F. Über die Elastizität Poröser Medien[J].Vier⁃teljschr.naturforsch.ges.zürich,1951,96(96):1-23.

[2]Biot M A.Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid.I.Low-Frequency Range[J].Journal of the Acoustical Society of America,1956,28(2):168-178.

[3]Biot M A.Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid.II.Higher Frequency Range[J].Journal of the Acoustical Society of America,1956,28(2):179-191.

[4]吴从辉,李红星,李勤武,等.基于波场分离的双相各向同性介质正演模拟[J].东华理工大学学报：自然科学版,2017(1):58-65.

[5]Stoll R D,Bryan G M.Wave Attenuation in Saturated Sedi⁃ments[J].Journal of the Acoustical Society of America,1970,47:1440-1447.

[6]Dvorkin J,Nur A.Dynamic poroelasticity:a Unified Model with the Squirt and the Biot Mechanisms[J].Geophysics,1993,58(4):524-533.

[7]Chotiros N P.An Inversion for Biot Parameters in Water-satu⁃rated Sand[J].Journal of the Acoustical Society of America,2002,112(5):1853-1868.

[8]门福录.波在饱水孔隙弹性介质中的传播[J].地球物理学报,1965:107-114.

[9]李红星.杭州湾沉积物原位声学特性分析及浅表低速层研究[D].吉林大学,2008.

[10]梅振繁.双相介质地震波场模拟与特征分析[D].东华理工大学,2014.