1 引言

作为一种岩性油气藏，砾岩含油气资源的有效勘探及开发问题近年来受到各大油田所关注.由于砾岩储层具有近物源、快速堆积形成的特征，具体表现为岩性成分多样、矿物颗粒大小不一、孔隙结构复杂，储层非均质性极强，这些因素造成砾岩储层的地震预测难度极大 (昝灵等，2011).国内学者基于岩石结构和孔隙结构特征 (刘敬奎，1983；张守伟，2011)，将砾岩储层划分为单模态 (近似描述为单一孔隙结构)、双模态 (双重孔隙) 以及复模态 (多重孔隙).在复杂的砾岩结构中，建立储层的物性、含油气性 (岩性、孔隙度、渗透率、含烃饱和度) 与地球物理参数 (纵波、横波、密度) 之间的定量关系，是开展砾岩储层地球物理预测与评价的基础 (马丽娟等，2002；印兴耀等，2015).

储层多孔介质中弹性波传播理论的研究，经历了由单孔介质向双重孔隙介质的发展演化过程.所谓“单孔介质”是将地下岩石理想化为只含一种均一孔隙的多孔介质.针对单孔介质波传播的研究最早始于20世纪50年代，Gassmann (1951)讨论了单孔介质中孔隙流体的存在对介质整体弹性的影响，没有考虑孔隙流体与骨架之间的相对运动；Biot (1956)则分析了单孔介质中孔隙流体的Poiseuille流动，讨论了纵、横波的传播特征，这种宏观流体流动可导致超声波频带内波的频散和衰减.Biot-Gassmann理论为多孔介质地震波传播的研究奠定了基础.

随后基于Biot理论，研究者们进一步开展了双重孔隙介质的研究，所谓“双重孔隙介质”是将地下岩石内部具有不同压缩性的区域，抽象为含双重孔隙结构的多孔介质，这种结构可由岩石内部多相流体不均匀分布或孔隙本身的压缩性差异所导致.考虑多相孔隙流体的非均匀分布所形成的含水孔与含气孔，Pride等 (2004)基于双重孔隙介质理论，分析了非饱和岩石中的地震波传播特征.巴晶等 (2012, 2013) 则从Biot-Rayleigh理论出发，研究了非饱和岩石中纵波的频散和衰减特征，并基于岩石物理模型进行了非均质天然气藏预测.一些相关研究通过岩石物理实验，证实部分饱和岩石中气体的非均质分布可造成弹性波的强烈衰减 (Cadoret et al., 1998；Toms-Stewart et al., 2009).

针对孔隙结构非均质性所形成的“双重”孔隙，Berryman和Wang (2000)基于孔-裂隙储层提出了一种双重孔隙介质理论.Dvorkin和Nur (1993)、唐晓明 (2011)、Tang等 (2012)则基于BISQ理论，研究了孔-裂隙之间的“喷射流效应”，分析了裂隙流体对波衰减和频散特征的影响.考虑到实际岩石中软孔隙 (裂隙) 的纵横比及其对应含量在较宽的范围内呈谱分布，邓继新等 (2015)将喷射流模型进行了扩展，以研究岩石复杂孔隙分布特征引起的速度频散.吴建鲁等 (2015)具体分析了硬币型、尖灭型裂隙与球形孔隙之间“喷射流效应”.Brajanovski等 (2010)曾将孔-裂隙结构及部分饱和流体之间的局部流体流动效应简单叠加，以解释纵波的强衰减.

Ba等 (2011)将描述球状流体周期振荡的Rayleigh理论与Biot理论结合，建立了双重孔隙地震波传播理论方程--Biot-Rayleigh理论，并于2014年研究了球状嵌入体内部的流体速度场，对双重孔隙地震波传播方程进行了修正，分析了不同储层环境条件下的纵波传播规律，随后，Ba等 (2015)基于Biot-Rayleigh理论，结合实际储层中非均质分布的孔隙流体与双重孔隙结构叠合形成的“双重-双重孔隙”特征，给出了可确凿描述部分饱和复杂储层的波动方程.

针对单一孔隙、双重孔隙介质中孔隙流体流动对地震波传播特征的影响问题，目前前人已开展系列研究 (Müller et al., 2010；巴晶，2013；孙卫涛等，2015)，而相对更为复杂的孔隙介质模型 (如三重孔隙介质) 中的地震波传播问题，仍有待开展进一步研究.近年来，复杂储层已成为油田重点勘探、开发的研究对象，本文拟将由砾岩、砂岩和泥岩相互充填所组成的三重孔隙结构作为研究对象，推导复杂砾岩储层地震波传播方程.

此处所谓“三重孔隙介质”，是将岩石内部不同压缩性的区域，抽象为含三种均一孔隙骨架成分的复合物.本文首先结合砾岩储层的孔隙结构特征，根据孔隙类型划分建立静态的砾岩骨架模型；其次，分析波激励下砾岩孔隙区域、砂岩孔隙区域及泥岩孔隙区域相互间的局部流体流动机制，建立三重孔隙砾岩储层的势能函数、动能函数和耗散函数，基于哈密顿原理，推导弹性波传播动力学方程组；最后，采用实测砾岩储层参数，通过与前人成果对比检验本文方程的合理性，并通过分析不同储层参数对波响应特征的影响，建立储层波频散、衰减特征与孔隙、流体性质间的定量关系.

2 研究基础

罗明高 (1991)分析了沉积环境对储层孔隙结构的控制作用，认为砂岩中多常见单模态 (单一孔隙组分) 和双模态 (两种孔隙组分) 结构，而砾岩中常见双模态与复模态孔隙结构，砾岩储层的孔隙结构要比砂岩储层复杂得多.从储层宏观地质背景和细观孔隙特征角度看，砾岩储层多杂基支撑，以砾岩为骨架的孔隙空间大多被砂级颗粒充填，而由砂粒组成的孔隙又被黏土颗粒充填 (曲卫光，2010).据罗明高针对碎屑岩结构模态定义的“三模态充填式”及“三模态悬浮式”：前者较为罕见，由于粒间的两次嵌套充填，使孔隙减少，如砂质砾岩，在局部砾石间孔隙充填砂，砂间充泥，即三重孔隙充填式结构，见示意图 1a；后者更为常见，由于砂、泥含量较高，砾悬浮于砂泥中，孔隙发育差，通常为泥质含砾砂岩或泥砂质砾岩，即三重孔隙悬浮式结构，见示意图 1b.

三重孔隙充填式结构中，砾岩作为背景相介质，砂岩作为嵌入体，泥岩作为砂岩的嵌入体，假设泥岩所形成的嵌入体半径明显小于砂岩嵌入体，可简化理论方程，具体充填式三重孔隙 (Filling-mode Triple Porosity, FTP) 介质模型如图 2a所示；对于三重孔隙悬浮式结构，将砂岩作为背景相介质，砾岩和泥岩作为砂岩的嵌入体，假设泥岩与砾岩形成的嵌入体稀疏分布，可忽略嵌入体间的相互影响，则成为悬浮式三重孔隙 (Suspension-mode Triple Porosity，STP) 介质模型 (图 2b).图 2中下标1, 2, 3分别表示砾岩孔隙区、砂岩孔隙区和泥岩孔隙区，砾岩、砂岩和泥岩的局部孔隙度分别为φ₁₀、φ₂₀和φ₃₀，体积比率分别为v₁、v₂和v₃，则绝对孔隙度分别为φ₁=φ₁₀v₁、φ₂=φ₂₀v₂和φ₃=φ₃₀v₃.

本文所讨论的三重孔隙介质模型，需满足以下假设：(1) 三种孔隙区域应分别满足传统孔隙介质理论对单孔介质的基础假定；(2) 砾岩、砂岩以及泥岩的相对组分比率中，嵌入体的体积比率不超过30%；(3) 宏观平均近似假设；(4) 最小力学微元假设；(5) 均一球形嵌入体假设 (具体假设3~5详见文献 (巴晶等，2012) 论述，此处不再赘述).在满足以上基础假设的前提下，基于此前所提出的Biot-Rayleigh理论 (Ba et al., 2014)，分析地震波激励下，砾岩孔隙区、砂岩孔隙区以及泥岩孔隙区之间的孔隙流体局部流动，建立三重孔隙介质模型的动能、势能以及耗散函数，推导三重孔隙介质的波传播动力学方程.

3 理论推导

基于上述讨论，本节分别从充填式与悬浮式三重孔隙介质模型着手，以孔隙介质力学为基础，推导三重孔隙介质的动力学方程.

3.1 充填式三重孔隙介质势能函数

基于Biot (1962)所提出的单孔介质理论，与Ba等 (2011)提出的双重介质理论方程，充填式三重孔隙介质的势能函数可写为:W=W(I₁, I₂, I₃, ξ₁, ξ₂, ξ₃)，其中I₁，I₂，I₃表示固体应变矩阵的三个旋转不变量 (巴晶等，2012)，ξ₁，ξ₂，ξ₃表示三类孔隙中的流体体应变.仅考虑应力应变的线性联系，泰勒展开后舍去二阶以上的高阶项，并考虑纵波激励后，孔隙区域1与2和2与3边界间流体流动所引起的体应变增量ζ₁₂和ζ₂₃，则势能函数形式为:

(1)

其中，N是固体骨架剪切模量，A是固相弹性参数，Q_m(m=1，2，3) 是固体、流体耦合弹性参数，R_m(m=1，2，3) 是流相弹性参数.

分析地震波激励下局域流流动对系统势能的影响.首先，孔隙区域1流入孔隙区域2的流体变化量为-φ₁ζ₁₂，孔隙区域2的孔隙度为φ₂；孔隙区域2流入孔隙区域1的流体变化量为φ₂ζ₁₂，孔隙区域1的孔隙度为φ₁.在周期性胀缩振动的每一个瞬时状态，流体质量守恒, 满足φ₁φ₂ζ₁₂+φ₂(-φ₁ζ₁₂)=0，并考虑到φ₁ζ₁₂=1-V₀⁽¹²⁾/V⁽¹²⁾，这里V₀⁽¹²⁾和V⁽¹²⁾分别表示砂岩嵌入体半径R₁₂、孔隙区域1与2边界的动态流体球半径R_f12所形成的球体体积，因此，可以得到

(2)

其次，孔隙区域2流入孔隙区域3的流体变化量为-φ₂ζ₂₃，孔隙区域3的孔隙度为φ₃；孔隙区域3流入孔隙区域2的流体变化量为φ₃ζ₂₃，孔隙区域2的孔隙度为φ₂.在周期性胀缩振动的每一个瞬时状态，流体质量守恒, 满足φ₂φ₃ζ₂₃+φ₃(-φ₂ζ₂₃)=0，并考虑到φ₂ζ₂₃=1-V₀⁽²³⁾/V⁽²³⁾，这里V₀⁽²³⁾和V⁽²³⁾分别表示泥岩嵌入体半径R₂₃、孔隙区域2与3边界的动态流体球半径R_f23所形成的球体体积，因此，可以得到

(3)

3.2 充填式三重孔隙介质动能函数

基于前人理论 (Ba et al., 2011；Biot，1962)，考虑地震波激励下，孔隙区域1与2以及孔隙区域2与3之间所诱导的局域流，可将充填式三重孔隙介质的动能函数写为

(4)

其中ρ₀₀、ρ₀₁、ρ₀₂、ρ₀₃、ρ₁₁、ρ₂₂、ρ₃₃表示三重孔隙介质中的7个密度参数，可根据岩石中固体颗粒的密度、孔隙流体的密度、孔隙度以及不同孔隙结构的体积含量进行计算，具体表达式在本文3.6节给出.u₁，u₂，u₃与U₁，U₂，U₃分别表示固体与流体三个方向上的位移分量，各位移分量的上标 (1)，(2)，(3) 分别表示砾岩、砂岩和泥岩三类孔隙，动能函数与固体的绝对运动速度矢量及流体的绝对运动速度矢量相关，其中第一项表示三重孔隙介质固体动能函数，第二至第四项分别表示三种孔隙流固耦合动能函数，第五至第七项分别表示三种孔隙流体动能函数，T_LFF12和T_LFF23分别表示孔隙区域1与2和2与3之间局域流振荡的动能函数.

分析地震波激励下局域流流动对系统动能的影响.首先，基于Rayleigh (1917)理论方程，可得到单个嵌入体内孔隙区域1与2之间的局域流动能形式

(5)

这里ρ_f为流体密度.

假设系统内存在N₀个砂岩嵌入体，孔隙区域1与2边界的动态流体球总体积为：4/3πR_f12³N₀=(φ₂φ₃₀+φ₂₀φ₃)/(φ₂₀φ₃₀)，进而得到

(6)

考虑孔隙区域1与2边界的局域流流动所造成的体应变增量ζ₁₂，可得到的表达式

(7)

故整个系统中，孔隙区域1与2之间的局域流振荡动能为

(8)

其次，考虑泥岩嵌入体内外的流体速度场 (Ba et al., 2014)，进而得到整个系统中，孔隙区域2与3之间的局域流振荡动能为

(9)

3.3 充填式三重孔隙介质耗散函数

基于孔隙流体与固体骨架的摩擦耗散机制，可得到充填式三重孔隙介质耗散函数的具体形式:

(10)

式中b₁、b₂和b₃分别表示孔隙区域1、2和3的耗散系数，D_LFF12和D_LFF23分别表示孔隙区域1与2和2与3之间局域流振荡的耗散函数.

与推导动能函数同理，分析地震波激励下局域流流动对系统耗散函数的影响.可得到整个系统中，孔隙区域1与2之间的局域流耗散函数为：

(11)

这里η为流体黏滞系数，κ₁为砾岩骨架的渗透率.

考虑泥岩嵌入体内外的流体速度场 (Ba et al., 2014)，整个系统中，孔隙区域2与3之间的局域流耗散函数为

(12)

这里κ₂与κ₃分别为砂岩骨架与泥岩骨架的渗透率.

3.4 充填式三重孔隙介质波动方程

基于哈密顿原理可以直接推导地震波动力学方程，其中拉格朗日能量密度可表示为：

(13)

带耗散的拉格朗日方程可写为:

(14)

式中x分别代表位移矢量u、U⁽¹⁾、U⁽²⁾和U⁽³⁾.

局域流矢量ζ的拉格朗日方程可表示为：

(15)

这里ζ分别代表局域流矢量ζ₁₂和ζ₂₃.

推导得到充填式三重孔隙 (FTP) 介质的地震波传播动力学方程组 (具体推导过程见附录A)：

(16a)

(16b)

(16c)

(16d)

(16e)

(16f)

3.5 悬浮式三重孔隙介质波动方程

悬浮式三重孔隙 (STP) 介质的动能函数、势能函数以及耗散函数的形式与充填式基本相同，因此其动力学传播方程组形式也基本相同，唯一的差别在于孔隙区域1与2之间的局域流控制方程，故式 (16e) 修改为

(17)

3.6 弹性参数与密度参数

本文采用虚拟实验的方法 (Johnson，1986)，推导三重孔隙介质波动方程中各弹性参数、密度参数与岩石、流体基本性质之间的关系.首先，推导三重孔隙介质的应力应变关系

(18)

这里i, j, m=1, 2, 3，τ为固体应力分量，S为流体应力分量，e为固体的体应变，ξ_m为第m类孔隙区域流体体应变，若i=j，δ_ij=1；若i≠j，δ_ij=0.

(1) 假设岩石样品只受剪切力作用，则e=ξ_m=0，并且孔隙流体不承担剪力，则可得μ_b=N，这里μ_b为岩石骨架的等效剪切模量.

(2) 将岩石用一个橡胶皮套包裹，然后将包裹岩石的皮套置于静水压环境下，且允许流体从中流出，此时三重孔隙介质的应力应变关系可写为：

(19)

式中K_b为岩石骨架的体积模量.

(3) 将岩石置于统一的静水压流体环境中，此时三重孔隙介质的应力应变关系可写为：

(20)

(21)

式中K_se为岩石基质的等效体积模量.

文中K_b、μ_b及K_se均采用Hashin-Shtrikman边界理论 (Hashin and Strikman, 1963) 计算.

假定每一孔隙区域是相对孤立的，针对三重孔隙介质，需引入两个新变量β₁和β₂，以描述孔隙区域1与2以及孔隙区域2与3之间骨架压缩比例，则

(22)

(23)

这里K_g，K_s和K_m分别表示砾岩、砂岩和泥岩的体积模量.

基于上述假设所得到的八条先验性关系，推导八个弹性常数的显式表达式.

(24)

沿用孔隙介质理论对密度参数的定义 (Gassmann，1951)，首先对不同压缩性孔隙区域的弯曲度进行定义：

(25)

进而推导三重孔隙介质的密度参数：

(26)

这里ρ_g，ρ_s和ρ_m分别表示砾岩、砂岩和泥岩的密度.

4 数值算例 4.1 与前人理论的对比分析

理论上，三重孔隙介质模型可退化为双重孔隙介质模型，故基于相同的岩石参数，将两者的预测结果进行对比分析，可验证本文理论的合理性.为确保算例中所选岩石参数的一致性，本文采用张守伟 (2011)所发表的实验数据，具体包括：砾岩基质的体积模量为72 GPa，剪切模量为38 GPa，骨架的体积模量为51.12 GPa，剪切模量为24.42 GPa，密度为2610 kg·m^-3，孔隙度为0.092，渗透率为0.05×10^-12m²；砂岩基质的体积模量58.7 GPa，剪切模量为31 GPa，骨架的体积模量为8.42 GPa，剪切模量为6.26 GPa，密度为2600 kg·m^-3，孔隙度为0.228，渗透率为0.1×10^-12m²；泥岩基质的体积模量为43.4 GPa，剪切模量为20.2 GPa，骨架的体积模量为34.77 GPa，剪切模量为17.751 GPa，密度为2450 kg·m^-3，孔隙度为0.04，渗透率为0.02×10^-12m²；水体积模量为2.5 GPa，密度为1040 kg·m^-3，黏滞系数为0.001 Pa·s.

算例A

在充填式三重孔隙 (FTP) 介质模型算例中，砾岩、砂岩与泥岩的体积含量分别为0.8、0.16与0.04，所采用的嵌入体尺寸R₁₂=1 m，R₂₃=0.05 m.图 3为充填式三重孔隙介质模型与双重孔隙介质模型所预测的纵波速度及衰减曲线.图中FTP表示存在砾岩、砂岩和泥岩三类孔隙，砾岩含砂岩包体，砂岩含泥岩包体；FTP (v₃=0) 表示泥岩含量为0，仅存在砾岩与砂岩两类孔隙，体积含量分别为0.8与0.2；FTP (v₂=0) 表示砂岩含量为0，仅存在砾岩与泥岩两类孔隙，体积含量分别为0.8与0.2.图中同时给出了Berryman理论 (Berryman and Milton, 1991) 双孔模型在低频极限下分别对应于两类双孔情况 (砾岩与砂岩/砾岩与泥岩) 的预测结果.

如图 3a所示，FTP (v₃=0)、FTP (v₂=0) 所预测的纵波速度与Berryman理论低频极限吻合度高，即本文方程与Berryman和Milton (1991)所提出的Gassmann公式在混合孔隙介质中的推广得到了相互印证，Berryman理论在双重孔隙介质低频极限是适用的，然而不能适用于更复杂的三重孔隙介质；FTP (v₃=0) 双孔介质模型所预测速度小于FTP (v₂=0) 所预测速度，原因在于砂岩与泥岩含量相当时，泥岩骨架的体积模量明显高于砂岩；同理，将前者与充填式三孔介质模型 (FTP) 进行对比，发现两者的骨架孔隙度虽相差不大，但预测的纵波速度存在差别，原因仍在于泥岩骨架的影响.

如图 3b所示，充填式三孔介质模型 (FTP) 所预测的衰减曲线，呈现三处“衰减峰”，是由砾岩孔隙区域与砂岩孔隙区域、砂岩孔隙区域与泥岩孔隙区域之间的局域流以及Biot全局流所形成；同理，FTP (v₃=0)、FTP (v₂=0) 双孔模型所预测的两处“频散台阶”和“衰减峰”，分别是由砾岩孔隙区域与砂岩孔隙区域之间的局域流以及Biot全局流、砾岩孔隙区域与泥岩孔隙区域之间的局域流以及Biot全局流所形成.

算例B

悬浮式三重孔隙 (STP) 介质模型算例中，砾岩、砂岩与泥岩的体积含量分别为0.16、0.8与0.04，所采用的嵌入体尺寸与算例A相同.图 4为悬浮式三重孔隙介质模型与双重孔隙介质模型所预测的纵波速度与衰减曲线.图中STP表示存在砾岩、砂岩和泥岩三类孔隙，砾岩孔隙与泥岩孔隙分别充填在砂岩孔隙中；STP (v₃=0) 表示泥岩含量为0，仅存在砾岩与砂岩两类孔隙，体积含量分别为0.2与0.8；STP (v₁=0) 表示砾岩含量为0，仅存在砂岩与泥岩两类孔隙，体积含量分别为0.8与0.2.

在图 4a中，STP (v₃=0) 双孔介质模型所预测的速度整体高于STP (v₁=0) 双孔介质模型以及悬浮式三孔 (STP) 介质模型，原因在于砾岩的体积模量高于泥岩、砂岩.如图 4b所示，悬浮式三孔介质模型 (STP) 所预测的衰减曲线，呈现三处“衰减峰”，是由砂岩孔隙区域与砾岩孔隙区域、砂岩孔隙区域与泥岩孔隙区域之间的局域流以及Biot全局流所形成；同样，STP (v₃=0)、STP (v₁=0) 双孔模型所预测的两处“频散台阶”和“衰减峰”，分别是由砂岩孔隙区域与砾岩孔隙区域之间的局域流以及Biot全局流、砂岩孔隙区域与泥岩孔隙区域之间的局域流以及Biot全局流所形成.

4.2 不同储层环境下纵波的传播特征

本节基于充填式三重孔隙介质模型，从不同储层环境出发，分析不同流体类型以及岩石参数对纵波传播特征的影响.

算例C

采用与算例A相同的岩石参数，仅仅改变的是流体类型，所采用的流体参数包括：水、油、气的体积模量分别为2.5、1.32、0.0001 GPa，密度分别为1040、880、100 kg·m^-3，黏滞系数分别为0.001、0.006、0.00001 Pa·s.

图 5a给出了充填式三重孔隙介质模型不同流体类型情况下预测的纵波速度曲线，由于气体积模量远小于水 (油)，含气孔隙介质预测的纵波速度明显低于含水 (油) 孔隙介质.图 5b为充填式三重孔隙介质模型所预测的纵波衰减曲线，结果显示：砾、砂以及砂、泥孔隙区域之间的局域流所形成衰减峰随黏滞系数的增大向低频方向移动；而Biot全局流所形成衰减峰则随黏滞系数的增大向高频方向移动.

算例D

采用与算例A相同的岩石和流体参数，仅改变嵌入体的半径, 则有：D1: R₁₂=0.5 m, R₂₃=0.05 m；D2:R₁₂=1 m, R₂₃=0.05 m；D3:R₁₂=2 m, R₂₃=0.05 m；D4: R₁₂=1 m, R₂₃=0.001 m；D5:R₁₂=1 m, R₂₃=0.02 m.图 6为充填式三重孔隙介质模型不同嵌入体半径所预测的纵波速度与衰减曲线，结果表明：泥岩嵌入体半径保持不变时，随着砂岩嵌入体半径的增大，对应的第一个纵波速度频散台阶与衰减峰向低频方向移动，第二个频散台阶与衰减峰基本不受影响；同理，砂岩嵌入体半径不变时，随着泥岩嵌入体半径的增大，对应的第二个速度频散台阶与衰减峰也向低频方向移动.嵌入体半径的改变，仅导致纵波速度与衰减曲线沿频率轴左右移动，而不影响其幅值的大小，这与文献 (Ba et al., 2011) 所得结论是一致的.

算例E

采用与算例A相同的岩石参数和流体参数，仅仅改变嵌入体的含量，设置E1：v₂=0.08，v₃=0.02；E2：v₂=0.18，v₃=0.02；E3：v₂=0.28，v₃=0.02；E4：v₂=0.18，v₃=0.04；E5：v₂=0.18，v₃=0.06.针对不同的嵌入体含量，图 7为充填式三重孔隙介质模型所预测的纵波速度与衰减曲线.由图 7a可知，在泥岩嵌入体含量保持不变的前提下，随着砂岩嵌入体的含量增加，纵波速度降低；在砂岩嵌入体含量不变的前提下，随着泥岩嵌入体含量的增加，纵波速度降低，原因在于砂岩、泥岩骨架较砾岩比较“软”，故砂/泥岩嵌入体含量增加，纵波速度降低，而相对于泥岩骨架，砂岩骨架更“软”，故砂岩含量的增加，导致纵波速度减小的幅值较大.如图 7b所示，在泥岩嵌入体含量不变的前提下，随着砂岩嵌入体的含量增加，地震频带内，衰减幅值随之增加；测井频带，衰减幅值则随之减小；实验室频带，衰减幅值几乎没有变化.反之，在砂岩嵌入体含量不变的前提下，随着泥岩嵌入体含量的增加，地震频带，衰减幅值几乎没有变化；测井频带，衰减幅值则随之增加；实验室频带，衰减幅值则没有变化.衰减曲线中的“第一个衰减峰”，主要由砾岩孔隙区域与砂岩孔隙区域之间的局域流导致；曲线中的“第二个衰减峰”，主要由砂岩孔隙区域与泥岩孔隙区域之间的局域流导致；曲线中的“第三个衰减峰”，由Biot全局流导致.

算例F

采用与算例A相同的岩石参数和流体参数，以及张守伟 (2011)所发表的实测砾岩数据.首先，砾岩骨架的体积模量为27.14 GPa，剪切模量为19.51 GPa，孔隙度为0.0836，仅改变砂岩骨架的孔隙度，设置F1：砂岩骨架的体积模量为3.6 GPa，剪切模量为2.94 GPa，孔隙度为0.235；F2：砂岩骨架的体积模量为8.42 GPa，剪切模量为6.26 GPa，孔隙度为0.228；F3：砂岩骨架的体积模量为15 GPa，剪切模量为9.71 GPa，孔隙度为0.184.其次，采用算例F2中的砂岩骨架参数，仅改变砾岩骨架的孔隙度，设置F4：砾岩骨架的体积模量为38.3 GPa，剪切模量为21.8 GPa，孔隙度为0.056；F5：砾岩骨架的体积模量为51 GPa，剪切模量为31 GPa，孔隙度为0.036.

图 8为充填式三重孔隙介质模型不同砂、砾岩孔隙度所预测的纵波速度与衰减曲线，如图 8a所示，在砾岩、泥岩孔隙度保持不变的前提下，随着砂岩孔隙度的降低，纵波速度增加，原因在于砂岩骨架的体积模量增大，而速度频散幅值减小，则是砂岩与砾岩骨架之间弹性参数差异减小的缘故；同理，在砂岩、泥岩孔隙度保持不变的情况下，随着砾岩孔隙度的降低，纵波速度增加，此时砂、砾之间弹性差异反而增大，导致频散幅度增大.如图 8b可知，在砾岩、泥岩孔隙度保持不变的情况下，在10^-2~10⁴ Hz频带内，衰减峰随砂岩孔隙度的降低而减小，而在10⁴~10⁸ Hz频带内，衰减随砂岩孔隙度变化不明显；在砂岩、泥岩孔隙度保持不变的情况下，在10^-2~10⁴ Hz频带内，衰减峰随砾岩孔隙度的降低而增大，在10⁴~10⁸ Hz频带内，衰减峰随砾岩孔隙度变化不明显.

算例G

采用与算例A相同的岩石参数和流体参数，仅仅改变背景相的渗透率.对于充填式三重孔隙介质模型，设置G1：κ₁=0.05×10^-13m²；G2：κ₁=0.25×10^-13m²；G3：κ₁=0.5×10^-13m².图 9为充填式三重孔隙介质模型不同砾岩渗透率所预测的纵波速度与衰减曲线，结果表明：在10^-2~10²频带内，随着渗透率的降低，纵波速度的频散台阶与衰减峰向低频方向移动，这与刘炯等 (2009)得到结论是一致的；而在10²~10⁶ Hz频带内，频散台阶与衰减峰则向高频方向移动，这与传统的Biot理论所得结论是一致的.由于砂、泥岩骨架的渗透率参数不变，砂、泥岩孔隙区域之间的局域流所形成衰减峰则不变；由于砾岩骨架的渗透率发生变化，故由砾、砂孔隙区域之间的局域流以及Biot全局流所形成衰减峰发生平移.

4.3 慢纵波速度的分析

算例H

采用与算例A相同的岩石和流体参数，分析充填式三重孔隙介质模型所预测的三类慢纵波 (P2、P3和P4) 的特征.图 10为充填式三重孔隙介质模型所预测慢纵波速度与衰减曲线.如图 10a所示，慢纵波P2、P3和P4为耗散波，其速度在低频带趋于0，随着频率的增大而增大，而在高频带趋于定值.由图 10b可知，P2慢纵波的逆品质因子低频带趋于1，随着频率的增加而减小.考虑到慢纵波主要由孔隙流体与固体骨架之间的摩擦产生，因此，三重孔隙介质中孔隙流体与砾岩骨架、砂岩骨架和泥岩骨架之间的摩擦，分别形成P2、P3和P4慢纵波.相比较砾岩骨架与泥岩骨架，砂岩骨架的局部孔隙度、渗透率大，所诱导的P3慢纵波速度要高于前两者.

5 结论

本文根据宏观砾岩储层的地质背景与细观孔隙结构的非均质性，基于双重孔隙介质Biot-Rayleigh理论，推导了砾岩储层三重孔隙结构的地震波传播动力学方程组，并通过虚拟实验的方法，推导了弹性参数与密度参数的具体表达式.

通过与双重孔隙介质模型的对比，结果显示砾岩孔隙区域与砂岩孔隙区域、砂岩孔隙区域与泥岩孔隙区域之间的局域流以及Biot全局流，是造成不同频段弹性波频散和衰减的主要原因.

基于相同的岩石骨架参数，分析了不同流体类型对纵波速度和衰减的影响.含气孔隙介质中所预测的纵波速度明显低于含水 (油) 孔隙介质中所预测速度.在频率轴坐标系中，随黏滞系数的增大，砾与砂及砂与泥孔隙区域之间的局域流所形成衰减峰向低频方向移动，Biot全局流形成衰减峰则向高频方向移动.

调查了不同砾岩储层参数对纵波速度和衰减的影响.结果显示：嵌入体尺寸及背景相介质渗透率仅影响纵波速度频散曲线沿频率轴平移，不影响其上下限幅值；嵌入体含量以及孔隙度的变化导致了岩石骨架参数的改变，不仅影响速度曲线沿频率轴平移，而且会影响速度上下限幅值.砾包砂包泥三重孔隙介质模型所预测的衰减曲线中，低频段“第一个衰减峰”主要由砾岩孔隙区域与砂岩孔隙区域之间的局域流导致，中间频段“第二个衰减峰”主要由砂岩孔隙区域与泥岩孔隙区域之间的局域流导致，超声频段“第三个衰减峰”由Biot全局流导致.

给出了三重孔隙介质理论方程所预测的慢纵波特征，结果显示砾包砂包泥介质中孔隙流体与砾岩骨架、砂岩骨架和泥岩骨架之间的摩擦，分别形成P2、P3和P4慢纵波，且砂岩骨架 (局部孔隙度大) 内部的Biot流体流动所造成的耗散要强于另外两个组分 (局部孔隙度较小).

本文通过研究孔隙非均质性所导致的“三重”孔隙介质，提供了砾岩储层中地震波传播的基础方程.文中所建三重孔隙结构可充分描述砾岩储层的孔隙结构特征，但模型参数较多，进一步应用可结合储层地质特征的先验信息，进行储层岩石物理模型的建立，与岩石声学实验、测井以及地震数据进行对比、标定，构建砾岩储层的岩石物理学图版，结合地震反演与解释方法，实现砾岩储层的定量预测.另一方面，也可对理论方程开展正演研究，分析、归纳其波场特征，在不丢失其重要特点、基本规律的前提下，逐步对模型简化，实现“由简入繁，由繁入简”的过程，最后简化的方程，将可能直接适用于储层参数反演.

附录A

基于哈密顿原理可推导地震波动力学方程，其中拉格朗日能量密度可表示为：

(A1)

带耗散的拉格朗日方程可写为:

(A2)

代入固体位移u各分量，则 (A2) 式简化为：

(A3)

等号右边第一项为：

(A4)

等号右边第二项为：

(A5)

等号右边第三项为：

(A6)

整理后，导出式 (16a).

代入流体位移U⁽¹⁾、U⁽²⁾和U⁽³⁾各分量，可得

(A7)

(A8)

(A9)

(A10)

(A11)

(A12)

(A13)

(A14)

(A15)

进一步可得式 (16b)-(16d).

局域流标量ζ的拉格朗日方程可表示为：

(A16)

令ζ=ζ₁₂, (A16) 式整理为

(A17)

令ζ=ζ₂₃，整理为

(A18)

可导出式 (16e)-(16f).