在含流体储层的岩性识别和油气预测中, 孔隙度和孔隙结构对储层有效弹性性质的作用和影响十分重要^[1-7]。SUN等^[3]的研究表明:岩石孔隙扁率(结构)的变化, 可使岩石的波速变化在2 000m/s以上。但目前, 关于孔隙结构参数本身的研究较少, 孔隙结构与孔隙度对储层弹性性质影响的定量化研究更少。究其原因, 一方面是孔隙结构十分复杂, 难以像孔隙度那样用简单明确的参数表示; 二是在目前应用十分广泛的含流体双相介质的Gassmann流体方程中未包含孔隙结构参数, 使得利用该方程来定量化研究孔隙结构的影响变得复杂和困难; 三是孔隙结构的量化表示不统一, 例如, ESHELBY等用孔、洞、缝的纵横比(扁率)表示孔、缝的形状, 并结合与岩石压实性质有关的基质模量和泊松比等参数, 形成岩石结构参数^[8-11]; PRIDE等^[12]在构建结构参数时, 以固结度表示孔、缝密度和连通性; NUR^[13]以临界孔隙度的倒数表示孔隙度大小和岩石性态; 还有很多学者用电性和放射性等以及某些定性方法来表示孔隙结构^{[1, 14-15]}。这些关于孔隙结构的表示方法均有各自的适用范围, 在应用时如何正确选择和评价十分困难。

为克服上述困难, 我们将重点研究以下几方面的问题:探讨孔隙结构参数的含义和数学表述; 引入带有孔隙度和孔隙结构的岩石骨架统一模型和公式; 将岩石骨架统一公式代入Gassmann流体方程中, 形成可定量计算孔隙度和孔隙结构的扩展Gassmann方程; 在此基础上推导出孔隙度和孔隙结构对储层等效弹性属性影响的关系式, 作为定量比较二者对储层属性影响的理论基础; 通过数值模拟计算, 验证方法的正确性和适用性。

1 孔隙结构参数及其表述 1.1 孔隙结构的定义

地质上, 孔隙结构是指岩石所具有的孔、洞、缝的大小、形状和相互连通关系, 是岩石物理的重要参数。但相对岩石的孔隙度、渗透率和饱和度而言, 孔隙结构参数的定义很不统一, 表述也不一致:有的学者主要考虑孔、缝的大小、形状和岩石基质属性; 有的学者除了考虑孔、洞、缝的大小、形状和基质属性外, 还要研究与储层界面现象有关的表面张力、吸附作用、湿润度以及孔隙吼道和毛细管作用等的影响。前一类学者, 主要是地质、地球物理勘探工作者。他们研究的问题相对宏观一些, 且易于孔隙结构的量化表示。如前所述, 他们用孔、洞、缝的纵横比(扁率)表示孔、缝的形状, 以固结度表示孔、缝密度与连通性, 以临界孔隙度表示孔隙度大小和岩石性态。后一类学者主要是地球物理测井和油藏开发工作者, 他们对岩石孔喉等微观特征更为关注^{[1, 14-15]}, 涉及的问题更复杂。我们根据地震、地质问题的综合需求和易于定量表示的原则, 认为孔隙结构是指孔、喉的形状、大小及其连通性。

1.2 孔隙结构参数的统一表示

HE等^[16]根据现有文献中应用广泛的岩石骨架参数模型和相应的数学公式, 主要包括Eshelby-Walsh公式、Pride公式、Geertsma公式^[17]、Nur公式、Keys-Xu公式^[18]和Krief公式^[19]等, 综合提炼出带有孔隙结构参数的岩石骨架参数模型:

$ \left\{ \begin{array}{l} {\beta _D} = \frac{{{\beta _s}\left( {1 + q\eta } \right)}}{{\left( {1 - p\eta } \right)}}\\ {\mu _D} = \frac{{{\mu _s}\left( {1 - p\eta } \right)}}{{\left( {1 + {q^ * }\eta } \right)}} \end{array} \right. $

(1)

式中:p和q(或q^*)表示孔隙结构参数, 在岩石骨架模型中作为孔隙度η的系数; β_s为岩石基质压缩系数; β_D为骨架压缩系数; μ_D表示剪切模量。公式(1)有如下优点:①对p, q(或q^*)赋与不同的变量和数值可得到完全等价于Eshelby-Walsh, Pride, Geertsma, Nur, Keys-Xu和Krief的岩石骨架模型的表示式, 具有通用性(见表 1); ②可获得不同于现有模型的新关系式, 具有更强的适用性; ③孔隙结构参数p, q综合了孔隙度大小与范围、孔隙形状(扁率)、岩石基质属性、孔隙连通性、岩石固结度等属性, 使丰富的孔隙结构信息可用两个参数来表述; ④将(1)式代入Gassmann流体方程, 可获得既有孔隙度, 又含有孔隙结构的新方程。

2 基于Gassmann流体方程的孔隙结构参数的计算 2.1 Gassmann流体方程的扩展

Gassmann流体方程是流体替换和预测的基本公式, 适用于低频波传播问题的研究, 广泛用于基于地震资料的孔隙流体识别研究^{[2, 13, 19-26]}。问题在于该方程以岩石体积模量的倒数—压缩系数β为基本参数, 未包含孔隙结构参数:

$ \left\{ \begin{array}{l} \frac{1}{{{\beta _E} - {\beta _s}}} = \frac{1}{{{\beta _D} - {\beta _s}}} + \frac{1}{{\left( {{\beta _p} - {\beta _s}} \right)\eta }}\\ {\mu _D} = {\mu _E} \end{array} \right. $

(2)

式中:β_E为岩石饱和流体的有效压缩系数; β_p为岩石孔隙流体的压缩系数; μ_E为有效剪切模量。在不考虑剪切模量μ的情况下, (2)式有5个变量, 已知其中的4个, 可求出另一个变量(参数)。例如, 利用预先给定的β_E, β_D, β_p和β_s, 可方便地求出孔隙度η。为使(2)式能够包含孔隙结构参数, 我们将(1)式代入(2)式, 得到新方程:

$ \left\{ \begin{array}{l} \frac{1}{{{\beta _E} - {\beta _s}}} = \frac{1}{{\frac{{{\beta _s}\left( {1 + q\eta } \right)}}{{\left( {1 + p\eta } \right)}} - {\beta _s}}} + \frac{1}{{\left( {{\beta _p} - {\beta _s}} \right)\eta }}\\ {\mu _E} = \frac{{{\mu _s}\left( {1 - p\eta } \right)}}{{\left( {1 - {q^ * }\eta } \right)}} = {\mu _D} \end{array} \right. $

(3)

(3) 式包含了孔隙结构参数, 我们称其为Gassmann扩展流体方程。

2.2 Gassmann扩展流体方程的简化和实用化

考虑到流体压缩系数β_p往往远大于岩石基质压缩系数^[22], 即β_p≫β_s, β_p-β_s≈β_p, 将其代入(3)式, 且暂不考虑剪切模量, 可获得含有孔隙结构参数的计算孔隙度的新公式:

$ \eta = \frac{{\left( {p + q} \right)\left( {{\beta _E}{\beta _s} - \beta _s^2} \right) - {\beta _p}{\beta _s} + {\beta _p}{\beta _E}}}{{q{\beta _p}{\beta _s} + p{\beta _p}{\beta _E}}} $

(4)

为了方便计算孔隙结构参数, 并易于同孔隙度对岩石(或储层)有效弹性性质影响进行比较, 实际应用中, 需要合并p和q成一个参数。我们对(4)式右边的分母添加辅助项, 得:

$ \begin{array}{l} \eta = \frac{{\left( {p + q} \right)\left( {{\beta _E}{\beta _s} - \beta _s^2} \right) - {\beta _p}{\beta _s} + {\beta _p}{\beta _E}}}{{p{\beta _p}{\beta _E} + q{\beta _p}{\beta _s} + q{\beta _p}{\beta _E} - q{\beta _p}{\beta _E}}}\\ \;\;\; = \frac{{\left[ {\left( {p + q} \right){\beta _s} + {\beta _p}} \right]\left( {{\beta _E} - {\beta _s}} \right)}}{{\left( {p + q} \right){\beta _p}{\beta _E}\left[ {1 - \frac{{q\left( {{\beta _E} - {\beta _s}} \right)}}{{\left( {p + q} \right){\beta _E}}}} \right]}} \end{array} $

(5)

设 $\frac{{q({\beta _E} - {\beta _S})}}{{\left( {p + q} \right){\beta _E}}} = x$ 和S=p+q, 再设:

$ \begin{array}{l} C = \frac{{\left[ {\left( {p + q} \right){\beta _s} + {\beta _p}} \right]\left( {{\beta _E} - {\beta _s}} \right)}}{{\left( {q + p} \right){\beta _p}{\beta _E}}}\\ \;\;\; = \frac{{\left( {S{\beta _s} + {\beta _p}} \right)\left( {{\beta _E} - {\beta _s}} \right)}}{{S{\beta _p}{\beta _E}}} \end{array} $

(5) 式变为:

$ \eta = C\frac{1}{{1 - x}} $

(6)

在|x|＜1的情况下, 对1/(1-x)进行二项式展开, 得:

$ \frac{1}{{1 - x}} = 1 + x - {x^2} + {x^3} \cdots $

(7)

舍去(7)式中x的一次方以上项, 就得到了含单一孔隙结构参数S的孔隙度计算公式:

$ \begin{array}{l} \eta = \frac{{S\left( {{\beta _E}{\beta _s} - \beta _s^2} \right) - {\beta _p}{\beta _s} + {\beta _p}{\beta _E}}}{{S{\beta _p}{\beta _E}}}\\ \;\;\; = \frac{{\left( {S{\beta _s} + {\beta _p}} \right)\left( {{\beta _E} - {\beta _s}} \right)}}{{S{\beta _p}{\beta _E}}} \end{array} $

(8)

由(8)式可得孔隙结构参数的计算公式:

$ S = \frac{{{\beta _p}{\beta _E} - {\beta _p}{\beta _s}}}{{{\beta _E}{\beta _p}\eta + \beta _s^2 - {\beta _s}{\beta _E}}} $

(9)

饱和岩石有效弹性的计算公式:

$ {\beta _E} = \frac{{S\beta _s^2 + {\beta _p}{\beta _s}}}{{{\beta _p} + S{\beta _s} - S\eta {\beta _p}}} $

(10)

如果舍去(7)中x二次方以上项, 则得:

$ \begin{array}{l} \eta = C\left( {1 + x} \right) = \frac{{\left[ {\left( {p + q} \right){\beta _s} + {\beta _p}} \right]\left( {{\beta _E} - {\beta _S}} \right)}}{{\left( {q + p} \right){\beta _p}{\beta _E}}} \cdot \\ \;\;\;\;\;\;\left[ {1 + \frac{{q\left( {{\beta _E} - {\beta _S}} \right)}}{{\left( {p + q} \right){\beta _E}}}} \right] = \frac{{\left( {S{\beta _s} + {\beta _p}} \right)\left( {{\beta _E} - {\beta _S}} \right)}}{{S{\beta _p}{\beta _E}}} \cdot \\ \;\;\;\;\;\;\left[ {1 + \frac{{q\left( {{\beta _E} - {\beta _S}} \right)}}{{\left( {p + q} \right){\beta _E}}}} \right] \end{array} $

(11)

考虑到在大多数情况下, q/(p+q)≈1或者0, 见表 1。当其值为0时, (11)式退化成(8)式; 当其值为1时, 得:

$ \eta = \frac{{\left( {S{\beta _s} + {\beta _p}} \right)\left( {{\beta _E} - {\beta _S}} \right)}}{{S{\beta _p}{\beta _E}}}\left( {1 + \frac{{{\beta _E} - {\beta _S}}}{{{\beta _E}}}} \right) $

(12)

由于β_E的值恒大于β_s^[22], 所以(12)式的精度满足大多数实际问题的要求。

2.3 误差分析

计算孔隙度的公式(8)和公式(12)以及在此基础上推导的计算孔隙结构和储层饱和弹性参数的公式均存在一定的误差。该误差取决于公式(7)的截断误差。由于(7)的变量, 并且其值总是小于1, 则截断误差取决于x小于1的程度, 因为q/(p+q)≈1。可按实际问题中β_E和β_s的值确定x。若x≤0.1, 则舍去x²以上项, 其误差不超过1%; 若x≥0.1, 则需要保留x的更高次方项。

3 孔隙度和孔隙结构参数对岩石有效弹性性质影响的比较

孔隙度和孔隙结构对饱和岩石等效弹性性质的影响是一个倍受关注的问题。如果孔隙结构相对孔隙度对岩石等效弹性的影响很小, 则无须考虑孔隙结构的影响, 否则, 应当仔细估计其作用。下面我们通过定量分析得出孔隙度和孔隙结构参数在什么情况下对岩石弹性性质有什么样的影响。

3.1 定量比较公式的推导

如前所述, 在实际问题中存在:β_p≫β_s, β_s＜β_E, (8)式可进一步简化为:

$ \begin{array}{*{20}{c}} {\eta = \frac{{S\left( {{\beta _E}{\beta _s} - \beta _s^2} \right) - {\beta _p}{\beta _s} + {\beta _p}{\beta _E}}}{{S{\beta _p}{\beta _E}}} = \frac{{{\beta _s}}}{{{\beta _p}}} - }\\ {\frac{{\beta _s^2}}{{{\beta _p}{\beta _E}}} + \frac{1}{S} - \frac{1}{S}\frac{{{\beta _s}}}{{{\beta _E}}} \approx \frac{1}{S}\left( {1 - \frac{{{\beta _s}}}{{{\beta _E}}}} \right)} \end{array} $

(13)

由(13)式得:

$ {\beta _E} = \frac{{{\beta _S}}}{{1 - S\eta }} $

(14)

从(14)式出发, 分别求取孔隙结构S和孔隙度η对饱和流体有效模量β_E的偏导数。得:

$ \frac{{\partial {\beta _E}}}{{\partial S}} = \frac{{\eta {\beta _s}}}{{{{\left( {1 - S\eta } \right)}^2}}} $

(15)

$ \frac{{\partial {\beta _E}}}{{\partial \eta }} = \frac{{S{\beta _s}}}{{{{\left( {1 - S\eta } \right)}^2}}} $

(16)

设比值 $R = \frac{{\partial {\beta _E}}}{{\partial S}}/\frac{{\partial {\beta _E}}}{{\partial \eta }}$ 得:

$ R = \frac{\eta }{S} $

(17)

我们称R为孔隙结构相对孔隙度对储层有效弹性性质的影响因子。此因子是从简化公式(13)出发推导出来的, 但是仍然具有一般性, 如果从精度更高的公式(12)出发进行比值R的推导, 我们可得到与(17)式完全一样的结果, 只是过程复杂一些, 见附录A。

3.2 影响因子R值的大小比较和数值模拟分析

(17)式简单、明了, 但是内容十分丰富, 能够用它全面评价孔隙度和孔隙结构对储层弹性的不同影响。按R值的大小, 得到了下述3种基本情况:

1) R值较大时, 孔隙结构对岩石饱和流体有效压缩系数β_E的影响大于孔隙度的影响。这对应孔隙度比较大, 孔隙结构比较小的情况(参看公式(17))。图 1是孔隙结构S与岩石有效压缩系数β_E的交会图, 孔隙度η作为参量。该图是按附录中的公式(A4)计算的, 式中的基质压缩系数β_s一般按目标地层的平均值给出, 也可根据岩石的主要矿物成分, 通过查表获得^[10], 我们在计算中采用的β_s值为14×10^－3GPa^－1, 这是一种通用值。由图 1可见, η=0.2的曲线是3条曲线中η值最大的, 而该曲线横坐标所对应的孔隙结构S是3条曲线中最小的(横坐标范围为1~8), 曲线陡度最大, 岩石有效压缩系数β_E随孔隙结构的增加, 变化值达400×10^－3GPa^－1以上。说明孔隙结构的微小变化, 会引起β_E的很大变化。在此种情况下, 孔隙结构对岩石弹性性质的影响起主要作用。作为对比, 还可看到, η=0.1曲线的横坐标范围为1~18, 该曲线的陡度中等; η=0.05曲线的陡度较缓, 孔隙结构变化量ΔS达到35, 引起的压缩系数变化量Δβ_E相对于前者, 减小了一倍(变化值达200×10^－3GPa^－1)。说明在不同条件下, 孔隙结构对岩石弹性性质的影响可能成倍增加。

2) R值较小时, 孔隙度对岩石饱和流体有效压缩系数β_E的影响大于孔隙结构的影响。这对应孔隙结构参数比较大, 孔隙度比较小的情况。图 2是孔隙度η与岩石有效压缩系数β_E的交会图, 孔隙结构S作为参量, β_E的取值与第1种情况一样。同样, 在图 2中, S=15的那条曲线是4条曲线中S最大的, 该曲线横坐标的孔隙度最小, 而且曲线最陡, 孔隙度的微小变化会引起β_E的很大变化。此时孔隙度对岩石弹性性质的影响起关键作用。然而这种情况会发生变化, 在S值逐渐减小时, 所对应曲线的斜率随之减小。S=1时, 随孔隙度增大, β_E几乎不变, 孔隙度对岩石弹性性质没有影响, 转化成第1种情况。

3) R值适中时, 孔隙结构与孔隙度对饱和流体储层等效模量β_E的影响相当。图 1和图 2是利用本文推导的公式计算的。为验证公式的正确性和上述结论的合理性, 我们利用广泛应用的Keys-Xu模型^[18]计算了孔隙度和地震纵、横波速度交会图, 见图 3。图中的参数为孔隙的扁率, 它与本文的结构参数呈反比:较小的扁率对应较大的结构参数。由于岩石纵横波速度与弹性参数存在确定的转换关系, 因此孔隙度与地震速度交会图等价于孔隙度与岩石弹性参数的交会图。由图 3可见, 在结构参数较大(扁率为0.01)和孔隙度较小(小于0.1)时, 曲线斜率大, 纵波速度由5 000m/s以上, 减小到2 000m/s以上, 减小量为2 000~3 000m/s, 横波速度由3 000m/s以上降低到不到1 000m/s, 减小量也为2 000~3 000m/s。说明孔隙度对弹性的影响起主要作用, 类似于前面的第2种情况。同时还可看到, 随着孔隙度增大, 在扁率固定时, 曲线迅速变缓, 孔隙度引起的速度变化很小, 起主要作用的是扁率(孔隙结构), 与第1种情况类似。但无论那种情况, 孔隙结构和孔隙度均可引起地震纵、横波速度产生2 000~3 000m/s变化。这与我们的研究结果一致, 也与其它学者的研究相符^[3]。

需要说明的是, 对于前面的第3种情况(R值适中时, 孔隙结构与孔隙度对饱和流体储层有效模量β_E的影响相当), 需要判断什么是“适中”。为此我们对S进行归一化处理, 使其数值像孔隙度一样, 也分布在0和1之间。这样在R=1时, 孔隙结构与孔隙度对饱和流体岩石有效模量β_E的影响相当。同样, R＜1, 孔隙度的贡献大于孔隙结构; R＞1, 则相反, 孔隙结构的作用大于孔隙度的作用。S归一化处理后, R的作用很好记忆和理解, 便于问题的分析。但在制作模版时, 还需要利用未归一化的S值。

4 结束语

利用孔隙结构参数的数值化表示以及将岩石骨架统一模型引入Gassmann方程, 我们获得了能够同时计算孔隙度和孔隙结构的扩展Gassmann方程。基于此方程, 推导了饱和岩石弹性模量随孔隙度和孔隙结构等参数变化的简化公式, 并对简化公式进行了适用性和精度分析, 为定量解释和反演孔隙度和孔隙结构提供了可能。导出了全面评价和定量比较孔隙度和孔隙结构对储层弹性性质的影响因子R。根据R值的大小, 得出3种基本情况:影响因子较大时, 孔隙结构对饱和岩石有效压缩系数的影响起主要作用, 孔隙结构引起的有效压缩系数的变化, 可达一倍以上; 影响因子较小时, 孔隙度对饱和岩石有效压缩系数的影响起主要作用; 影响因子适中时, 二者的影响相当。作为对比, 采用常见的Keys-Xu模型计算了储层纵、横波速度随孔隙结构和孔隙度的变化, 经分析比较, 无论是孔隙度还是孔隙结构均可使岩石的纵、横波速度的最大变化量达到2 000~3 000m/s, 变化规律与所推导公式的结果一致。由于孔隙结构的复杂性, 其丰富的内涵仍需在实践中继续挖掘。

附录A 影响因子R的进一步推导

重写正文中的公式(12):

$ \eta = \frac{{\left( {S{\beta _s} + {\beta _p}} \right)\left( {{\beta _E} - {\beta _S}} \right)}}{{S{\beta _p}{\beta _E}}}\left( {1 + \frac{{{\beta _E} - {\beta _S}}}{{{\beta _E}}}} \right) $

按此式求解饱和等效弹性参数β_E。

将文中(13)式引入(12)式, 得:

$ \eta = \frac{1}{S}\left( {1 - \frac{{{\beta _S}}}{{{\beta _E}}}} \right)\left( {1 + \frac{{{\beta _E} - {\beta _S}}}{{{\beta _E}}}} \right) $

(A1)

对(A1)式整理后得:

$ \left( {2 - S\eta } \right) - 3\frac{{{\beta _S}}}{{{\beta _E}}} + {\left( {\frac{{{\beta _S}}}{{{\beta _E}}}} \right)^2} = 0 $

(A2)

此式有两解, 但β_S/β_E不能大于1, 其值也不可为负。取其有物理意义的解:

$ \frac{{{\beta _S}}}{{{\beta _E}}} = \frac{{3 - \sqrt {9 - 4\left( {2 - S\eta } \right)} }}{2} = \frac{{3 - \sqrt {1 + 4S\eta } }}{2} $

(A3)

于是得:

$ {\beta _E} = \frac{{2{\beta _S}}}{{3 - \sqrt {1 + 4S\eta } }} $

(A4)

式中, $\sqrt {1 + 4S\eta } < 3$ 。将(A4)分别对S和η求偏导数, 得:

$ \frac{{\partial {\beta _E}}}{{\partial S}} = \frac{{{\beta _S}{{\left( {1 + 4S\eta } \right)}^{ - \frac{1}{2}}}4\eta }}{{{{\left( {3 - \sqrt {1 + 4S\eta } } \right)}^2}}} $

(A5)

$ \frac{{\partial {\beta _E}}}{{\partial \eta }} = \frac{{{\beta _S}{{\left( {1 + 4S\eta } \right)}^{ - \frac{1}{2}}}4S}}{{{{\left( {3 - \sqrt {1 + 4S\eta } } \right)}^2}}} $

(A6)

于是得:

$ R = \frac{{\partial {\beta _E}}}{{\partial S}}/\frac{{\partial {\beta _E}}}{{\partial \eta }} = \frac{\eta }{S} $

(A7)

此结果与正文公式(17)完全一致。