基于三维孔隙网络模型的纵波频散衰减特征分析

魏乐乐， 甘利灯*， 熊繁升， 孙卫涛， 丁骞， 杨昊

1 中国石油勘探开发研究院， 北京 100083 2 清华大学周培源应用数学研究中心， 北京 100084

0 引言

地震波在通过孔隙介质时会发生频散和衰减现象.研究地震波衰减和速度频散特征对岩性分类、油气识别，以及油气田开发过程中的油气藏管理都具有重要的价值(何润发等，2020；凌云，2021；赵平起等，2021).由于地下介质条件复杂，地震波在地下含流体孔隙介质中传播时的衰减和频散特征一直是研究的热点问题.

近几十年来，国内外学者开展了大量岩石物理性质的研究，特别是孔隙介质流体的流动引起波的频散衰减过程，发展出许多理论模型.

Biot(1955)建立了完全饱和孔隙介质中波传播的动力学方程.在Biot理论中，黏性的孔隙流体和固体骨架之间的相互作用导致波的速度频散和衰减.Geertsma和Smit(1961)应用Biot理论推导出纵波的速度频散和衰减的近似解.然而，越来越多的研究发现，Biot理论的这种宏观的流体流动机制很难解释实际生产中地震波的高频散和高衰减现象(Carcione et al.，2010).

为了突破Biot理论的局限性，学者们着眼于研究微观尺度的流体流动机制.Mavko和Nur(1979)建立了喷射流模型来描述弹性波在流体部分饱和的裂缝/孔隙介质中的传播，而且认为喷射流动会造成更强的衰减.欧阳芳等(2021)对经典喷射流模型进行了扩展，结合等效介质理论和孔隙结构模型，通过数值模拟研究了微观孔隙结构下的速度频散和衰减特性.

Dvorkin和Nur(1993)认为，在含有流体的岩石中，波的速度频散和衰减除了受Biot机制的影响，还受到喷射流动机制的影响.他们将上述两种机制统一起来，提出了BISQ模型，并推导了纵波相速度和品质因子的表达式，研究了渗透率等对纵波速度频散和衰减的影响.

在Biot理论的框架下，White(1975a，b)首先引入中观尺度的耗散机制，提出非均匀部分饱和的孔隙介质模型，也称为斑状饱和模型.Dutta和Odé(1979a，b)通过严格求解Biot方程组，研究了在盐水充填的岩石中，衰减系数随频率、气体饱和度等的变化.Carcione等(2003)深入研究了White模型的物性及流体变化对衰减和频散的影响.王峣钧等(2014)应用斑块饱和岩石物理模型，从理论上分析了不同固结程度岩石中含气饱和度对速度和衰减的影响，结果表明局部含气储层在地震频带内会发生频散和振幅频变效应.在此基础上，他们建立了岩石物理量板，用于估计储层的含气饱和度.

Sun(2021)提出了一种BIPS(Biot-Patchy-Squirt)模型，用于描述非混相流体饱和裂缝孔隙弹性介质中的波的频散衰减特性.研究结果表明，BIPS模型与实验室数据吻合较好，这一发现将推动上述三种频散衰减机制在波速预测中的潜在应用.

随着研究的深入，学者们开始建立和研究多重孔隙介质模型.巴晶等(2012)建立了双重孔隙结构模型，利用描述非饱和双重孔隙介质地震波传播的Biot-Rayleigh方程求解平面波解，得到纵波、横波的相速度及逆品质因子的表达式，研究了非饱和岩石中的纵波频散与衰减特征.通过对三个地区砂岩储层的分析，发现在地震频带内纵波对储层中含气较为敏感，但对含气饱和度的定量表征效率不高，低孔隙度砂岩的纵波频散和衰减在地震频带更为显著.此后，郭梦秋等(2018)进一步研究了基于此模型的含流体致密砂岩的纵波频散和衰减特征.Ba等(2017)通过分析流体非均质性与岩石组构非均质性的叠加效应，提出了一种双重双孔隙理论.该模型描述了波在不同尺度下的具有组构非均质性的斑状饱和岩石中的传播.建模和实验数据表明，岩石组构的非均质性尺度会造成波的多频频散和衰减效应，特征频率的大小也与岩石组构的非均质性尺度有关.Zhang等(2020，2021)考虑到岩石组构具有多尺度的非均质性和自相似的分形特征.在此基础上，他们建立了多孔介质模型，通过数值模拟研究了波的频散和衰减特性与分形维数的关系.

在含油气复杂储层勘探时，不可避免的要遇到含裂缝和孔隙的储层岩石，这些岩石孔隙结构复杂、孔隙流体性质多变，导致波场信号出现频散衰减等特征.建立合理的裂缝-孔隙介质模型，是地震正演等技术的重要基础.

Chapman等(2002)提出了在微观尺度上建立孔隙弹性介质模型的方法，将模型设置成规则的立方体网格，每个网络节点代表孔隙/裂缝空间.但是此模型并未建立孔隙介质渗透率与裂缝参数等之间的几何关系.Tang(2011)、Tang等(2012)提出了含裂缝孔隙介质的波动方程模型，在推导过程中直接使用了Johnson等(1987)提出的动态渗透率模型，并没有在渗透率与裂缝参数等影响因素之间建立关系，也没有给出波频散衰减与渗透率之间的关系.Song等(2020)假定裂缝的形状可以是矩形的，建立了一个介质模型来估计饱和多孔岩石中包含一个随机分布的、无限小厚度的矩形裂缝时的纵波衰减和速度频散特征.

在天然岩石中，流体流动发生的空间是全局性的裂缝-孔隙网络，而不是局部的单一孔隙或孔管束.Xiong等(2020)、熊繁升等(2021)基于岩石内部的裂缝连通网络，用椭圆截面纵横比的变化模拟从扁裂缝、软孔隙到硬孔隙的多种情况，提出了具有椭圆形截面的三维裂缝/软孔隙网络模型，建立了宏观可测观量与裂缝参数之间的关系，并给出了渗透率的计算方法.该模型可以实现孔隙介质中跨尺度的整体建模，能更好地反映孔隙系统的整体效应.

考虑到熊繁升等(2021)建立的流体全饱和的三维裂缝/软孔隙网络模型没有考虑孔隙的存在，因此本文在此模型的基础上建立同时包含裂缝和孔隙的三维裂缝-孔隙网络模型，并推导出渗透率的计算方法.基于三维裂缝/软孔隙网络模型和三维裂缝-孔隙网络模型，运用体积平均的方法推导出基于渗透率的波动方程，利用平面波分析方法得到纵波频散和衰减曲线的表达式，详细研究了总孔隙度、裂缝孔隙度、裂缝纵横比、裂缝数密度以及孔隙流体黏度对快纵波频散衰减特征的影响，以及分析了慢纵波的频散衰减特征.

1 三维孔隙网络模型

1.1 三维裂缝/软孔隙网络模型

熊繁升等(2021)提出三维裂缝/软孔隙网络模型(图1).模型中单个管道横截面是纵横比可变的椭圆形，可以模拟从圆形孔隙到狭窄裂缝的不同类型的裂缝/软孔隙网络空间.

图1 三维裂缝/软孔隙网络模型Fig.1 3D fracture/soft pore network model

设单元立方体边长为l，沿X、Y、Z方向上分别排布有M、N、L个椭圆截面管道，用R1和R2表示裂缝椭圆截面的长短轴半径.定义如下参数：

(1)微裂缝的纵横比：a=R2/R1.

(2)裂缝的孔隙度：

(1)

(3)单元体内的裂缝数密度：

(2)

熊繁升等(2021)计算三维裂缝/软孔隙介质渗透率的步骤如下：

第一步：采用椭圆柱坐标系，基于质量守恒和NS方程推导出含不可压缩牛顿流体在椭圆截面微管中的流量表达式.

第二步：利用微管两端每个节点的流量守恒条件，得到全部网络节点满足的流量-压力线性方程组.该方程组中未知量为各节点的压力，入口和出口端的压力边界条件是方程的非齐次项.

第三步：求解上述流量-压力线性方程组，得到各节点压力，进而计算得到整个三维裂缝/软孔隙网络的流量.

第四步：结合达西定律计算得到岩石样本的渗透率.

1.2 三维裂缝-孔隙网络模型

在三维裂缝/软孔隙网络模型的基础上，本文建立了同时包含孔隙和裂缝的三维裂缝-孔隙网络模型(图2).

图2 三维裂缝-孔隙网络模型Fig.2 3D fracture-pore network model

同样地，设单元立方体边长为l，沿X、Y、Z方向上分别排布有M、N、L个椭圆截面管道.节点上的球形表示孔隙，连接节点的椭圆柱体通道表示裂缝.用rpore表示孔隙半径，Rasp表示平均孔隙半径与平均裂缝截面长轴半径之比.定义连接孔隙i、j的裂缝截面长轴半径为

(3)

定义如下参数：

(1)微裂缝的纵横比a：裂缝椭圆截面的长短轴半径之比的倒数；

(2)总孔隙度：φ=φf+φpore；

其中，裂缝所占孔隙度为

(4)

孔隙空间所占孔隙度为

(5)

(3)单元体内的裂缝数密度：表达式同公式(2).

对于三维裂缝-孔隙网络模型的渗透率计算，仍采用1.1节的推导方法.

对于同时含有裂缝和孔隙网络的介质，考虑任意一段微管，微管两端的网络节点为i和j，节点处压力分别是pi和pj，微管中压力与流量的关系式为

Qij=cijpi+dijpj，

(6)

其中：

(7)

(8)

这里：

(9)

(10)

(11)

接下来，对各节点列流量平衡方程式，得到全部网络节点满足的流量-压力方程组.求解该方程组得到各节点压力，进而计算得到整个三维裂缝-孔隙网络的流量.最后，结合达西定律计算得到岩石样本的渗透率.

2 基于体积平均法的波动方程推导

2.1 体积平均法

体积平均法是推导波动方程的重要方法.Whitaker(1999)详细介绍了采用体积平均法推导波动方程的原理和步骤.由于该方法通过对微观物理量在单元体上进行积分平均得到宏观物理量，因此可以考虑不同流体类型和孔隙空间的具体几何特征，可以认为是一种更一般化的建模方法.

记单元体Ω的体积为V，孔隙空间由固体骨架的体积和孔隙介质中流体的体积构成，即：

V=Vs+Vf，

(12)

记ψf为孔隙介质中与流体有关的物理量，规定在流体的外部，ψf的值为0.定义在整个区域Ω上对ψf作体积平均的方法为

(13)

定义本征体积平均量：

(14)

根据以上定义式，可知：

(15)

2.2 波动方程推导

在下述的推导过程中：假设流体为牛顿流体，不考虑流体发生正应变过程中黏性的作用，仅考虑剪切黏性；考虑孔隙介质含一种固体骨架、一种流体(或两种流体等效成一种).

用u、U分别表示固体、流体质点的位移向量；各物理量上方的一点表示对时间求导；以上下标中的s、f分别表示固体骨架、孔隙介质中所含的流体.

下面运用体积平均法来推导基于渗透率的三维孔隙网络模型(三维裂缝/软孔隙网络模型与三维裂缝-孔隙网络模型)的波动方程.

第一步，建立固体和流体的微观运动方程(动量守恒方程).

固体骨架部分的微观运动方程为

(16)

流体的微观运动方程为

(17)

式中，ρs为固体骨架的质量密度；σ为固体应力；ρf，pf为等效后流体的密度、压力；⊗表示张量积；s为流体应力.

第二步，给出固体部分和流体的本构方程.

固体部分的微观本构方程为

(18)

式中，Ks为固体颗粒的体积模量；G为固体骨架的剪切模量；e为固体的体应变.

流体本构方程为：

(19)

式中，η为等效后流体的黏度.

第三步，对微观方程做体积平均得到宏观大尺度方程.

对流体微观运动方程取体积平均后得到：

(20)

对固体微观运动方程取体积平均后得到：

(21)

式中，κ为孔隙介质的渗透率，可分别根据1.1节和1.2节的计算方法确定.

第四步，推导出流固相互作用的流体和固体的本构关系.

固体和流体的本构关系为

(22)

式中，弹性常数aij可根据固体和流体的孔隙度、体积模量来进行计算.具体表示为

(23)

这里：

(24)

等效体积模量Kb的表达式为(Mavko and Nur，1978)：

(25)

式中，v是固体骨架矿物的泊松比；Nc为单位体内的裂缝数；R1i是第i条裂缝的椭圆截面长轴半径长度；di是裂缝延伸长度；V是单元体的体积.

最后，将本构方程和体积平均后的动量守恒方程结合，得到最终形式的波动方程:

(26)

2.3 频散/衰减表征

对波动方程两端取散度，可得到纵波方程.对于各向同性孔隙介质，假设弹性波沿各方向传播一致，利用平面波分析方法可得到纵波频率-波数方程.将纵波频率-波数方程的解记为v，根据如下定义，即可计算纵波频散/衰减曲线：

(27)

3 频散/衰减特征分析

波场在含流体的孔隙介质中传播时会产生频散和衰减现象，即波速随着频率发生改变，同时波的振幅也会发生变化.

下面基于三维裂缝/软孔隙网络模型和三维裂缝-孔隙网络模型分别研究纵波频散与衰减特征.

3.1 基于三维裂缝/软孔隙网络模型的频散衰减特征

利用前面所述频散衰减曲线的表达式，根据岩石物理参数可计算得到基于三维裂缝/软孔隙网络模型的频率依赖的纵波速度和逆品质因子曲线.数值模拟所选取的模型、矿物和流体参数具体见表1.

表1 模型、矿物及流体参数表Table 1 Parameters of model, mineral and pore fluid

3.1.1 裂缝孔隙度的影响

图3、图4给出了不同裂缝孔隙度条件下，快纵波的频散和衰减曲线.本算例中除裂缝孔隙度外，其余参数见表1.

图3、图4中四条曲线的裂缝孔隙度分别为0.001、0.005、0.02和0.06，对应的渗透率分别为7.00×10-17m2、1.57×10-16m2、3.14×10-16m2和5.44×10-16m2.从图3可以看出，随着裂缝孔隙度的减小，快纵波速度增大，起始速度由3116 m·s-1增大到4863 m·s-1.而且随着裂缝孔隙度的增大，快纵波速度表现出更明显的频散现象.从图4可以看出，随着裂缝孔隙度的减小，逆品质因子曲线的峰值下降，特征频率(逆品质因子曲线的峰值对应的频率)向低频方向移动.

图3 快纵波频散曲线随裂缝孔隙度的变化Fig.3 Velocity dispersion changing with fracture porosity

图4 快纵波衰减曲线随裂缝孔隙度的变化Fig.4 Inverse quality factor versus fracture porosity

3.1.2 裂缝纵横比的影响

图5、图6给出了不同裂缝纵横比条件下，快纵波的频散和衰减曲线.本算例中除裂缝纵横比外，其余参数见表1.

图5、图6中四条曲线的裂缝纵横比分别为0.02、0.04、0.08和0.16，对应的渗透率分别为5.44×10-16m2、3.07×10-15m2、1.73×10-14m2和9.60×10-14m2.通过计算可知，随着裂缝纵横比的下降，渗透率下降.这是一个合理的事实，因为随着裂缝纵横比的减小，裂缝趋于闭合.图5显示，随着裂缝纵横比的增大，快纵波速度增大.这主要是因为随着纵横比的增大，干岩石的可压缩性变差，干岩石的体积模量Kb增大，从而引起快纵波速度的增大.图6显示，随着裂缝纵横比的增大，特征频率向低频方向移动.

图5 快纵波频散曲线随裂缝纵横比的变化Fig.5 Velocity dispersion versus aspect ratio of fracture

图6 快纵波衰减曲线随裂缝纵横比的变化Fig.6 Attenuation changing with aspect ratio of fracture

3.1.3 裂缝数密度的影响

图7、图8给出了不同裂缝数密度条件下，快纵波的频散和衰减曲线.本算例中除网络节点数外，其余参数见表1.

图7 快纵波频散曲线随裂缝数密度的变化Fig.7 Velocity dispersion varying with fracture density

图8 快纵波衰减曲线随裂缝数密度的变化Fig.8 Attenuation varying with fracture density

沿X、Y、Z方向上分别排布有3×3×3、4×4×4、5×5×5和6×6×6个孔隙，即裂缝数密度分别为54、144、300和540.计算得到对应的渗透率分别为1.40×10-16m2、3.11×10-16m2、5.44×10-16m2和8.39×10-16m2.通过计算可知，渗透率随裂缝数密度的增加而增大.这表明随着裂缝之间的连通程度升高，岩石渗透率逐步增大.图7显示，频散幅值(速度极大值与极小值的差)对裂缝数密度的变化不敏感，但频散曲线随裂缝数密度的变化发生移动.在相同的频率上，快纵波速度最大相差仅10 m·s-1.图8显示，随着裂缝数密度的增大，特征频率向低频方向移动，而几乎不影响逆品质因子曲线的峰值.

3.1.4 孔隙流体黏度的影响

下面研究孔隙流体的黏度对快纵波频散衰减特征的影响.孔隙中的流体参数见表2，其余参数同表1.

表2 孔隙流体参数表Table 2 Parameters of pore fluid

从图9、图10可以看出，在此模型下的快纵波速度：油饱和砂岩明显大于水饱和砂岩，水饱和砂岩明显大于气饱和砂岩.这主要是因为孔隙流体黏度增大时，孔隙流体与岩石之间的耦合作用加强，相当于提高了岩石的总体刚度，从而提高了快纵波的速度.但油饱和时的速度变化率相对于水饱和时变小.相较于孔隙充填其他流体，气饱和时的快纵波频散现象最不明显，且逆品质因子曲线的峰值最小.

图9 不同孔隙流体黏度下的快纵波频散曲线Fig.9 Velocity dispersion curve as function of fluid viscosity

图10 不同孔隙流体黏度下的快纵波衰减曲线Fig.10 Attenuation curve under different fluid viscosity

模型内不仅存在快纵波，还存在慢纵波.将表1中的岩石物理参数代入公式(27)，得到慢纵波的频散衰减曲线，如图11、图12所示.

图11 慢纵波频散曲线Fig.11 Velocity dispersion curve of slow P-waves

图12 慢纵波衰减曲线Fig.12 Attenuation curve of slow P-waves

3.1.5 慢纵波频散衰减曲线

从图11可以看出，在低频范围内(0～102Hz)慢纵波速度基本为0，且在0～108Hz频率范围内发生频散，速度随频率的增加而增大，在高频范围内稳定在大约740 m·s-1.从图12可以看出，慢纵波逆品质因子曲线在0～106Hz的频率范围内为水平线，后随着频率的增加而减小，呈直线下降趋势.

3.2 基于三维裂缝-孔隙网络模型的频散衰减特征

根据式(3)知，裂缝椭圆截面长轴半径与裂缝两端的孔隙半径有关.因此，裂缝的孔隙度φf和孔隙空间的孔隙度φpore不是互相独立的.

在实际的计算中，给定φpore，首先根据式(5)计算得到孔隙半径，后根据式(3)计算得到裂缝截面的长轴半径，再根据裂缝纵横比得到裂缝的短轴半径.将上述计算得到的裂缝截面的长轴半径代入式(4)，计算得到φf.

为了满足预设的总孔隙度φ，可能需要对孔隙半径和裂缝半径进行重新处理.因此，在三维裂缝-孔隙网络模型中，总孔隙度φ与预设结果相同，但是φpore与φf可能与预设结果不同.

3.2.1 总孔隙度的影响

除裂缝孔隙度外，仍采用表1中的岩石物理参数.图13、图14给出了总孔隙度分别为0.001、0.005、0.02和0.06时的快纵波频散与衰减曲线.通过计算得到对应的渗透率为3.12×10-17m2、7.00×10-17m2，1.40×10-16m2和2.22×10-16m2.从图中可以看出，随着总孔隙度的减小，快纵波速度增大，起始速度由3559 m·s-1增大到4367 m·s-1，逆品质因子曲线的峰值下降，特征频率向低频方向移动.

图13 快纵波频散曲线随总孔隙度的变化Fig.13 Velocity dispersion curve as function of total porosity

图14 快纵波衰减曲线随总孔隙度的变化Fig.14 Attenuation curve as function of total porosity

3.2.2 裂缝参数的影响

从建模过程中可以发现，三维裂缝/软孔隙网络模型可以认为是三维裂缝-孔隙网络模型的特例，即当孔隙空间所占孔隙度为0时，三维裂缝-孔隙网络模型退化为三维裂缝/软孔隙网络模型.若固定孔隙空间所占孔隙度为0，分析裂缝参数对快纵波频散衰减特征的影响，即为文章3.1节的内容.

3.2.3 慢纵波频散衰减曲线

将表1中的岩石物理参数代入公式(27)，得到慢纵波的频散衰减曲线.

从图15、图16可以看出，在该模型下慢纵波频散与衰减的趋势与三维裂缝/软孔隙网络模型相似，但频散的幅值更大，发生频散的频率范围更宽，在高频范围内慢纵波的速度稳定在大约1214 m·s-1.

图15 慢纵波频散曲线Fig.15 Velocity dispersion curve of slow P-waves

图16 慢纵波衰减曲线Fig.16 Attenuation curve of slow P-waves

4 结论

本文基于三维裂缝/软孔隙网络模型和三维裂缝-孔隙网络模型，详细研究了总孔隙度、裂缝孔隙度、裂缝纵横比、裂缝数密度、流体黏度对快纵波频散衰减特征的影响，以及分析了慢纵波的频散衰减特征.

基于三维裂缝/软孔隙网络模型，可以发现：

(1)随着裂缝孔隙度的增大，快纵波速度表现出更明显的频散现象.随着裂缝孔隙度的减小，快纵波速度增大，逆品质因子曲线的峰值下降，特征频率(逆品质因子曲线的峰值对应的频率)向低频方向移动.

(2)随着裂缝纵横比的增大，快纵波速度增大，特征频率向低频方向移动，而频散幅值(速度极大值与极小值的差)和逆品质因子曲线的峰值的大小未表现出规律性的变化.

(3)裂缝数密度的变化几乎不会影响快纵波速度、频散幅值和逆品质因子曲线的峰值，但对特征频率影响显著.随裂缝数密度的增大，特征频率向低频方向移动.

(4)当孔隙充填油、水、气时，油饱和砂岩的快纵波速度明显大于水饱和砂岩，水饱和砂岩明显大于气饱和砂岩.气饱和时的快纵波频散现象最不明显，且逆品质因子曲线的峰值最小.

基于三维裂缝-孔隙网络模型，可以发现：

在此模型下，总孔隙度、裂缝参数等对快纵波频散衰减特征的影响以及慢纵波的频散衰减趋势与三维裂缝/软孔隙网络模型相似.

总的来说，孔隙度的变化主要影响逆品质因子曲线的峰值的大小；裂缝数密度主要控制速度显著变化的范围；裂缝纵横比的变化对纵波速度和特征频率都影响显著.