0 引言

岩石物理模型主要是指一些有关岩石弹性、黏弹性或者孔隙弹性的理论模型，其中比较著名的是Biot孔隙弹性模型^[1]和低频Gassmann模型^[2]。岩石物理模型一般分为包含物模型和接触胶结模型，前者将岩石的孔隙度、矿物基质、孔隙形状、孔隙流体与岩石的弹性模量联系起来建立模型^[3-10]，后者将岩石近似为分散的矿物颗粒集合体，弹性参数则由矿物颗粒与颗粒接触的变形性和刚度来确定，这种模型大多是基于Hertz-Mindlin模型建立的理论公式^[11-14]。Dvorkin等^[15]描述了在球形颗粒接触时加入少量矿物胶结后产生的影响。在包含物模型和接触胶结模型中均涉及到很多岩石物理参数，如岩石的体积模量、剪切模量、孔隙度、密度、含水饱和度，矿物基质的体积模量、剪切模量以及孔隙形状等。这些岩石物理参数一般都是未知的，有的可以从测井资料和实验室测量数据中直接或者间接获得，有的则需要利用已知信息通过反演得到。

测井资料中包含了岩石的矿物含量、泥质含量、密度、含水饱和度以及纵横波速度等信息，利用这些信息可以求取岩石的有效体积模量和剪切模量，但是有些时候，由于受到井口条件、岩石流体饱和度以及岩性变化的影响^[16]，从测井资料中得到的岩石物理参数不一定准确。相比之下，在实验室条件下测量的数据很精确，可以通过测量得到岩石在不同状态（干燥、饱水和饱油）下的纵波速度、横波速度以及密度，还可以通过岩石薄片鉴定确定岩石中矿物的组分及其含量，但是从实验室数据中无法直接获取矿物基质的弹性模量。如果构成岩石的矿物是单一的，可以直接参考Mavko等^[17]给出的纯矿物的弹性模量，如果构成岩石的矿物是多种的，可以采用基于等效介质理论的平均模型方法获取各种矿物的弹性模量。

Voigt ^[18]和Reuss ^[19]先后提出了对矿物性质进行体积平均的空间平均模型，后来Hill ^[20]将二者的理论模型进行了一个算术平均，称之为Voigt-Reuss-Hill平均模型，该模型比较适合于计算矿物成分的有效体积模量和剪切模量及其可能的上下限值。Hashin等^[21]提出了一种新的平均模型的概念，它用“平均矿物”模量表示一种混合矿物模型，该模型说明对于任意岩石给定各成分的弹性模量后，其等效模量将位于界限之间，但精确值将取决于孔隙的几何信息。基于等效介质理论的平均模型方法虽然可以求取矿物基质的体积模量和剪切模量，但是需要准确地知道组成岩石的各种矿物所占的百分比以及单矿物的弹性模量，而且每个测井深度点对应的矿物基质弹性模量是有差别的，如果取值不正确，计算的结果往往误差较大。因此，一些学者研究并提出了更加有效的岩石物理参数反演方法，可以在已知测井数据的情况下直接反演得到矿物基质弹性模量以及其他岩石物理参数。云美厚等^[22]利用测井数据，拟合了关于矿物基质弹性模量、孔隙度、泥质含量、有效压力和温度的多元函数关系式，利用这一关系式可以反演得到未知的岩石物理参数。张金强等^[23]根据测井资料，通过选取岩石物性较均匀的地层进行岩石物理参数的统计、分析，获取了均匀岩性对应的矿物基质弹性模量。Lin等^[24]提出了一种基于Gassmann理论的自适应矿物基质体积模量反演方法（简称为林凯方法），可较好地解决流体替换过程中矿物基质体积模量的取值问题，能够自适应逐点反演每个测井深度点对应的高精度矿物基质体积模量，但是该方法针对的岩石孔隙介质模型只能是满足Gassmann理论的球形孔隙模型。对于岩石孔隙结构，也有学者做了一些研究。胡晓丽等^[25]利用KT模型研究了孔隙纵横比对AVO响应的影响；刘航宇等^[26]利用分形理论研究了碳酸盐岩储层孔隙结构的分形特征；葛小波等^[27]利用分形理论研究了致密砂岩储层微观孔隙结构的分形特征；闫建平等^[28]利用压汞、核磁共振、物性、薄片及CT分析等方法研究了低渗透滩坝砂岩的孔隙结构。

Kuster等^[7]根据散射理论，得到了一个双相介质包含物模型，即KT模型。该模型中涉及到的岩石物理参数比Gassmann模型更多，不仅包括矿物基质体积模量，而且还包括矿物基质剪切模量以及岩石中孔隙纵横比的大小和不同孔隙纵横比孔隙的体积分数，同时不需要求取干燥岩石的体积模量。由于利用前人提出的反演方法不能够同时获取矿物基质剪切模量和孔隙纵横比谱参数，因此，在KT模型理论公式的基础上，提出一种基于KT模型流体替换的岩石物理参数反演方法，针对实验室条件下测量的细砂岩样品，同时反演所有样品的矿物基质体积模量、剪切模量和不同孔隙纵横比孔隙的体积分数，并检验利用KT模型进行流体替换对实验室测量数据的适应性，再利用经反演获取的多种岩石物理参数分析砂岩储层的孔隙结构特征，以及不同孔隙纵横比孔隙的体积分数与饱水岩石弹性模量之间的敏感性关系。

1 理论模型和反演方法 1.1 包含3种孔隙纵横比的KT模型

KT模型将岩石孔隙度和孔隙形状与纵横波速度联系在了一起。假设孔隙的形状为椭圆形，那么孔隙的纵横比可定义为椭圆体的短轴与长轴之比。饱和多孔介质等效模量的KT模型表达式推广后可写为

$ \left( {k_{{\rm{sat}}}^* - {k_{{\rm{ma}}}}} \right)\frac{{\left( {{k_{{\rm{ma}}}} + \frac{4}{3}{\mu _{{\rm{ma}}}}} \right)}}{{\left( {k_{{\rm{sat}}}^* + \frac{4}{3}{\mu _{{\rm{ma}}}}} \right)}} = \sum\limits_{i = 1}^N {{c_i}\left( {{k_i} - {k_{{\rm{ma}}}}} \right){P^{{\rm{mi}}}}} $

(1)

$ \left( {\mu _{{\rm{sat}}}^* - {\mu _{{\rm{ma}}}}} \right)\frac{{\left( {{\mu _{{\rm{ma}}}} + \frac{4}{3}{\zeta _{{\rm{ma}}}}} \right)}}{{\left( {\mu _{{\rm{sat}}}^* + \frac{4}{3}{\zeta _{{\rm{ma}}}}} \right)}} = \sum\limits_{i = 1}^N {{c_i}\left( {{\mu _i} - {\mu _{{\rm{ma}}}}} \right){Q^{{\rm{mi}}}}} $

(2)

$ {\zeta _{{\rm{ma}}}} = \frac{{{\mu _{{\rm{ma}}}}}}{6}\frac{{9{k_{{\rm{ma}}}} + 8{\mu _{{\rm{ma}}}}}}{{{k_{{\rm{ma}}}} + 2{\mu _{{\rm{ma}}}}}} $

(3)

$ \varphi = \sum\limits_{i = 1}^N {{c_i}} $

(4)

式（1）~（4）中：k_sat^*为岩石的有效体积模量（饱水或饱油时岩石的体积模量），GPa；μ_sat^*为岩石的有效剪切模量（饱水或饱油时岩石的剪切模量），GPa；k_ma和μ_ma分别为矿物基质的体积模量和剪切模量，GPa；k_i和μ_i分别为第i种孔隙类型中孔隙流体的体积模量和剪切模量，GPa；P^mi和Q^mi为在矿物基质中加入第i种孔隙类型时的系数；ζ_ma为岩石中矿物的压缩系数；φ为孔隙度，%；c_i为第i种孔隙类型的体积占岩石总孔隙体积的百分比，%。

P^mi系数和Q^mi系数是矿物基质弹性模量k_ma和μ_ma、孔隙流体弹性模量k_fl以及孔隙纵横比的α函数，Wu ^[6]提出的具有任意高宽比的椭球状包含物公式以及Berryman ^[29]提出的具有球状、针状和硬币状的3D形态的包含物公式都可以用来计算这2个系数。

Cheng等^[30]根据扫描电镜成像结果，并利用数值计算方法近似得到了一些国际上典型砂岩样品的孔隙纵横比谱数据（表 1）。孔隙纵横比谱数据采用科学计数法表示。

从表 1可知，常见的砂岩孔隙结构由1.00和0.10以及一系列小于0.01的孔隙纵横比反映出来。因此，可以利用Wu ^[6]提出的包含物公式求取P^mi和Q^mi，并利用3种孔隙纵横比（1.00，0.10和0.01）来描述所有砂岩的孔隙结构特征。

图 1为3种孔隙纵横比的示意图。其中，α₁ = 1.00代表球形孔隙，相当于地层中在围压下难以闭合的溶蚀孔、洞；α₂ = 0.10代表粒间孔隙，相当于岩石的背景孔隙；α₃ = 0.01代表裂缝孔隙，这种孔隙在围压下很容易闭合。

1.2 矿物基质弹性模量的求取及反演流程

在确定矿物基质弹性模量的取值时，引用了Krief等^[31]总结出的干燥岩石体积模量与矿物基质体积模量之间的Biot系数关系式，即

$ {k_{{\rm{dry}}}} = {k_{{\rm{ma}}}}{\left( {1 - \varphi } \right)^{\frac{3}{{1 - \varphi }}}} $

(5)

式中：k_dry为干燥岩石的体积模量，GPa。

式（5）变换后可得

$ {k_{{\rm{ma}}}} = {k_{{\rm{dry}}}}{\left( {1 - \varphi } \right)^{\frac{3}{{\varphi - 1}}}} $

(6)

在岩石物理学中，体积模量的相对关系式为

$ {k_{{\rm{ma}}}} > {k_{{\rm{sat}}}} > {k_{{\rm{dry}}}} $

(7)

林凯等^[32]根据Biot系数建立了矿物基质体积模量与岩石有效体积模量、孔隙度之间的关系式为

$ {k_{{\rm{sat}}}} < {k_{{\rm{ma}}}} < {k_{{\rm{sat}}}}{\left( {1 - \varphi } \right)^{\frac{3}{{\varphi - 1}}}} $

(8)

另外，在岩石物理学中，岩石的矿物基质体积模量和剪切模量与岩石的有效剪切模量之间也应该存在一定的相互关系，即

${k_{{\rm{ma}}}} > {\mu _{{\rm{ma}}}} > {\mu _{{\rm{sat}}}} $

(9)

由此可见，式（8）和式（9）分别确定了矿物基质体积模量和矿物基质剪切模量的变化区间。为了同时精确地获取多种岩石物理参数，结合模拟退火组合优化算法，提出了基于KT模型流体替换的岩石物理参数反演方法，其实现步骤为：①输入实测数据（测井数据或实验室测量数据）中的纵波速度、横波速度、含水饱和度、密度及孔隙度，然后利用速度公式计算已知饱和度下岩石的体积模量和剪切模量。②利用式（8）确定矿物基质体积模量的变化区间，并在该区间内随机产生一个矿物基质体积模量，再利用式（9）确定矿物基质剪切模量的变化区间，并在该区间内随机产生一个矿物基质剪切模量；随机产生1.00，0.10和0.01等3种孔隙纵横比孔隙的体积分数，并满足三者之和与孔隙度大小相等；将产生的岩石物理参数代入到KT模型中，利用KT模型公式进行流体替换，替换为已知含水饱和度下岩石的体积模量和剪切模量，并与利用测井数据或实验室测量数据计算得到的岩石弹性模量进行对比，如果两者存在较大误差，则不断地利用模拟退火优化算法对岩石物理参数进行优选，直到满足利用KT模型计算得到的岩石弹性模量与通过实测数据计算得到的岩石弹性模量之间的误差达到最小为止。③输出误差达到最小时获得的矿物基质体积模量和剪切模量以及3种孔隙纵横比孔隙的体积分数。

2 反演的多种岩石物理参数分析

为了检验本文方法的应用效果，在实验室条件下，针对来自高邮凹陷戴南组埋深2 000~4 000 m的42块细砂岩岩心样品，进行了岩石物理测试。首先，利用超声波脉冲透射法（横波频率和纵波频率分别取200 kHz和700 kHz）测量了岩样饱水和饱油状态下的纵波速度和横波速度；然后，采用电子天平测量了岩样饱水和饱油状态下的质量，并用电子游标卡尺测量了岩样的长度、直径，每块岩样的每项参数均测量3次，取平均值，计算岩样的密度；最后，对每块岩样均进行了覆压孔隙度测量。在岩石物理测试之后，利用本文方法反演出了42块细砂岩样品的孔隙纵横比谱和矿物基质弹性模量，并对样品的多种岩石物理参数进行了分析。

2.1 孔隙纵横比谱和矿物基质弹性模量

图 2为利用本文方法反演出的42块细砂岩样品的孔隙纵横比谱。从图 2可以看出，该砂岩储层以发育球形孔隙为主（约占总孔隙体积的75%），其次是粒间孔隙（约占总孔隙体积的20%），裂缝孔隙发育较少（约占总孔隙体积的5%）。

随着地层埋深增大，孔隙结构同时受到压实作用和胶结作用的影响，前者会使裂缝孔隙增加，后者会使裂缝孔隙减少，因此，裂缝孔隙在岩石总孔隙中的含量取决于地层中压实作用占主导地位还是胶结作用占主导地位。从本次实验反演得到的孔隙纵横比谱（图 2）可知，大多数样品裂缝孔隙的体积分数都很低，可见在该砂岩储层中，岩石组构主要是受到了胶结作用的影响。

图 3为分别利用林凯方法和本文方法反演出的42块细砂岩样品的矿物基质弹性模量。从图 3可以看出，2种方法反演出的矿物基质体积模量整体上趋势一致，大多数样品为20~45 GPa；利用本文方法反演出的所有样品的矿物基质剪切模量（15~ 40 GPa）均较其矿物基质体积模量稍低。

2.2 理论数据与测量数据的对比

为了检验利用本文方法反演得到的多种岩石物理参数在KT模型中进行流体替换对实验室数据的适应性，首先分别利用林凯方法和本文方法反演获取了42块细砂岩样品的多种岩石物理参数，然后将它们分别代入到Gassmann模型和KT模型中，求取了所有样品在饱油状态下的纵波速度，最后与实验室测量得到的饱油岩石纵波速度进行了对比。

图 4为利用理论方法计算的饱油岩石纵波速度与实验室测量数据的对比，图中反映了利用林凯方法反演出所有样品的岩石物理参数后，再利用Gassmann模型计算的饱油岩石纵波速度与实验室测量数据的对比，以及利用本文方法反演出所有样品的岩石物理参数后，再利用KT模型计算的饱油岩石纵波速度与实验室测量数据的对比。从图 4可以看出，利用林凯方法和本文方法计算的饱油岩石纵波速度与实验室测量数据的吻合程度均较高，并且本文方法的计算结果与实验室测量数据的吻合程度更高。在考虑相对误差小于3%的情况下，利用本文方法计算的数据有38个（占样品总数的90%）与实验室测量数据相吻合，而利用林凯方法计算的数据有31个（占样品总数的74%）与实验室测量数据相吻合。由此可见，对于存在多种孔隙结构的砂岩样品来说，利用KT模型进行流体替换模拟的适应性要比Gassmann模型更好。

2.3 孔隙结构对岩石弹性模量的敏感性

岩石中矿物基质的体积模量、孔隙度和孔隙纵横比的变化都会对岩石的弹性性质及纵横波速度产生影响，而孔隙度和孔隙纵横比的影响相对更大，远大于岩性变化产生的影响^[33]。因此，针对42块细砂岩样品，有必要分析并找到它们的孔隙结构与饱水岩石弹性模量之间的敏感性关系。

图 5和图 6分别为42块细砂岩样品的饱水岩石体积模量和饱水岩石剪切模量与孔隙度的交会图。图 5和图 6中每个圆形片都代表一个岩石样品，样品中球形孔隙、粒间孔隙和裂缝孔隙的体积分数均为利用本文方法反演获取。从图 5和图 6可以明显看出，饱水岩石弹性模量与孔隙度之间存在一定的负相关关系；在孔隙度小于10%的样品中，裂缝孔隙的体积分数较高，随着孔隙度的增大，裂缝孔隙的体积分数逐渐降低；在孔隙度为10%~ 20%的样品中，粒间孔隙的体积分数较高；在所有样品中，球形孔隙的体积分数都非常高。另外，从图 5和图 6还可以看出，随着孔隙度的增大，饱水岩石的体积模量和剪切模量均呈减小的趋势，并且当孔隙度小于10%时，减小的趋势更明显，而当孔隙度大于10%时，减小的趋势逐渐变缓，说明孔隙度对岩石体积模量和剪切模量的敏感性随着孔隙度的增大而减小。

对于某些细砂岩样品来说，它们的孔隙度差别很小，但是饱水岩石体积模量却差别很大。这是因为孔隙纵横比对岩石弹性模量的影响要大于孔隙度产生的影响。下面针对一些细砂岩样品，分析不同孔隙纵横比孔隙的体积分数与饱水岩石弹性模量之间的敏感性关系。

表 2中列出了从42块细砂岩样品中挑选出来的8块典型样品的岩石物理参数数据。从表 2可以看出，26号样品和40号样品的孔隙度都约为2%，二者球形孔隙的体积分数几乎相同，但粒间孔隙和裂缝孔隙的体积分数却有较大差别。由于任何低孔隙纵横比孔隙的体积分数增大都会使岩石的弹性模量降低^[34]，对比这2块样品的弹性模量可知，裂缝孔隙的体积分数对岩石弹性模量的敏感性要大于粒间孔隙的体积分数。

另外，从表 2还可以看出，16号样品和38号样品的孔隙度都约为7%，二者的孔隙结构特征完全不同，即16号样品以粒间孔隙为主（体积分数为67.8%），而38号样品则以球形孔隙为主（体积分数为91.8%），它们的饱水岩石体积模量差别很大；5号样品和9号样品的孔隙度都约为13%，二者的孔隙结构特征也完全不同，即5号样品以粒间孔隙为主（体积分数为59.6%），而9号样品则以球形孔隙为主（体积分数为94.5%），它们的饱水岩石体积模量却差别很小。这说明，不同孔隙纵横比孔隙的体积分数对岩石弹性模量的敏感性均随着孔隙度的增大逐渐减小。同样，对比1号样品和6号样品的岩石物理参数数据也可以看出此特点。

3 结论和认识

（1）本文提出的基于KT模型流体替换的岩石物理参数反演方法，可同时反演出砂岩样品的矿物基质体积模量、剪切模量以及孔隙纵横比谱，而且将反演出的多种岩石物理参数代入到KT模型中进行流体替换对实验室数据的适应性比利用Gassmann模型的效果更好。

（2）利用本文方法反演出的孔隙纵横比谱信息，可以从多个角度描述砂岩储层的孔隙结构特征。反演出的裂缝孔隙信息有助于了解砂岩储层中裂缝的数量、走向、分布及发育情况，为裂缝预测研究提供理论数据。

（3）岩石中存在的各种矿物成分对岩石的弹性性质均有一定的影响，当某块岩石中存在很多体积模量较大的矿物成分时，该岩石的体积模量和剪切模量也较大。

（4）在所有的岩石物理参数中，裂缝孔隙的体积分数对岩石弹性模量的敏感性最高，但是这种敏感性会随着孔隙度的增大而逐渐减小。

（5）对于砂岩油气储层岩石物理参数的获取，应该更多地选择和利用像KT模型这样考虑了孔隙纵横比、孔隙大小、孔隙几何形状、孔隙分布和孔隙连通性等岩石孔隙结构特征参数的岩石物理模型进行反演。