Gassmann方程是利用地震、测井和岩石物理实验室资料进行横波速度预测、流体替代及速度评价的主要工具^[1].但是，在Gassmann方程中干岩石骨架的弹性模量求取一直是岩石物理界还没有完全解决的问题.因此，很多地球物理学家给出了各种干岩石骨架弹性模量的计算方法.目前比较流行的干岩石骨架模型，有经验模型和等效介质模型两种.

经验模型通常简单假设弹性模量与孔隙度有关，不考虑岩石孔隙结构.最常见的经验模型是临界孔隙度模型和Pride模型等^[2-3]，这些模型对干岩石骨架的模量比都有各自不同的假设.临界孔隙度模型认为干岩石骨架的模量比为常数，与孔隙度无关，而Pride模型则认为干岩石骨架的模量比随孔隙度增大.但这些假设并不完全合理，实验室数据证明，干岩石骨架的模量比既有可能随孔隙度增大，也有可能随孔隙度减少^[4-7].因此，这些模型在用于研究岩石的弹性性质时往往具有局限性.

等效介质模型，如Kuster-Toksöz理论^[8]和微分等效介质理论^[9]，建立了干岩石骨架的弹性模量与孔隙度和孔隙形状之间的关系.岩石物理理论认为，岩石速度强烈受孔隙形状影响，它是造成速度-孔隙度关系发散的主要因素^[10-11]. Li和Zhang系统地由微分等效介质理论推导出干岩石骨架弹性模量公式，这里称之为DEM解析模型，这些公式可以客观反映骨架模量比及其弹性模量随孔隙度和孔隙形状参数的变化特征^{[7, 12-13]}.等效介质模型目前已广泛用于横波速度预测中^{[10, 14]}.等效介质模型中通常用理想的椭球来描述岩石孔隙的形状^{[8, 10, 15]}，然而表征孔隙形状的孔隙纵横比参数一般难以从实验室或野外测量中获得，流行的Xu-White模型假定岩石的孔隙由硬孔和软孔两部分组成，它们的体积由砂岩和粘土矿物的含量确定，且纵横比为常数^[10].这些假设条件影响了Xu-White模型的横波速度预测精度，因此，有学者提出在使用Xu-White模型时通过经验关系得到变化的孔隙纵横比来提高精度^[16-17]，但是经验关系基于统计规律，其应用效果依赖于统计样品，不具普遍性.

本文基于DEM解析模型^[13]和Gassmann方程提出一种根据纵波速度反演岩石的等效孔隙纵横比并进行横波速度预测和储层孔隙结构评价的方法.首先介绍Gassmann方程的原理，然后引入基于微分等效介质理论推得的包含孔隙纵横比的干岩石骨架模量解析公式，接着建立基于Gassmann方程应用非线性全局寻优算法来反演等效孔隙纵横比算法，最后用实验室测量的花岗岩和砂岩数据研究孔隙纵横比随有效压力的变化特征，同时对实验室测量的砂岩数据^[18]和实际测井数据进行等效孔隙纵横比反演和横波速度预测.

在低频假设下，Gassmann ^{[1, 19]}给出了流体饱和岩石的体积模量K_sat和剪切模量G_sat的计算公式

(1)

和

(2)

式中，K_dry、K_m和K_f分别为干岩石骨架、岩石基质和孔隙流体的体积模量，G_dry为干岩石骨架剪切模量，φ为总孔隙度.

由Gassmann方程获得饱和流体岩石的体积模量、剪切模量后，饱和流体岩石的纵、横波速度可表示为

(3)

(4)

式中，V_p为纵波速度，V_s为横波速度，ρ_sat为饱和流体岩石的密度.

利用Gassmann方程(式(1)和(2))研究饱和岩石的弹性性质时，需要预先提供有关干岩石骨架的弹性模量，但干岩石骨架的弹性模量并没有统一公式.通常有如下几种方法确定：(1)岩石样品的实验室直接测量^{[18, 20]}，该方法获得的干岩石骨架模量精度高，但一般情况很难得到.(2)经验关系^[2-3]或等效介质理论^{[8-10, 12-13]}预测，无论经验关系或等效介质理论，它们的预测精度依赖于模型的假设^{[12, 21]}.

基于经验关系的干岩石骨架模型简单假设模量与孔隙度有关，而不考虑岩石的孔隙结构.最常见的经验模型是临界孔隙度模型，它假设在岩石的孔隙度小于某个临界孔隙度下干岩石骨架的体积和剪切模量都随孔隙度线性变化，Nur等(1998)建议用下面公式建立干岩石骨架与基质的体积和剪切模量之间的关系^[2]：

(5)

(6)

式中，φ_c为临界孔隙度，G_m为岩石基质的剪切模量.对于砂岩，临界孔隙度φ_c≈0.40.

Pride等(2004)建议固结的干岩石骨架的体积和剪切模量用下面简单的形式^[3]：

(7)

(8)

式中，c是固结系数，表示岩石的固结程度，对于砂岩，通常取2＜c＜20.

无论是临界孔隙度模型还是Pride模型都具有局限性，它们都隐含了对干岩石骨架模量比的不适当假设.由式(5)-(8)知，临界孔隙度模型隐含定义了干岩石骨架模量比是常数，与孔隙度无关，也就是干岩石骨架的纵横波速度比与基质的纵横波速度比是相等的，但是这种假设有时与实验室测量结果不相符^{[4-5, 21]}.而Pride模型由于要求固结系数为正值，因此，Pride模型隐含定义了干岩石骨架模量比随孔隙度是单调增加的，而实际砂岩的骨架模量比随孔隙度的变化既可能是增加的(如石英砂岩)也可能是下降的(如长石砂岩)^[6-7].

包裹体理论，如微分等效介质理论，认为干岩石骨架的模量比不仅与孔隙度有关，而且与孔隙的几何形态有关^{[9, 15]}.微分等效介质理论的精确数值解可以很好地模拟理想的椭球孔和币形缝的干岩石骨架模量比随孔隙度的变化特征. Li和Zhang通过对微分等效介质理论中关于体积模量和剪切模量的两个相互耦合的常微分方程进行解耦，系统地得到关于干岩石骨架的模量比以及干岩石骨架的体积模量和剪切模量的DEM解析模型^{[7, 12-13]}.对于椭球孔，

(9)

(10)

式中，，，，，α为包裹体的纵横比.常数a和b分别为满足P^*i-Q^*i=a+bK^*(y)/G^*(y)的梯度和截距^{[7, 13]}，常数b可以简单设置为激化因子之差P^*i-Q^*i的一阶导数在K^*/G^*=K_m/G_m时所取的值，b已知后常数a就能够由公式P^*i-Q^*i=a+bK^*(y)/G^*(y)来计算得到^{[7, 12-13]}，椭球包裹体的激化因子之差P^*i-Q^*i及其一阶导数见文献[13]中的附录A.P^*i和Q^*i为所谓的等效体积模量K^*和等效剪切模量G^*的激化因子，它们通常根据Eshelby的成果和Wu的成果来计算^{[1, 15]}.

式(9)和(10)隐含定义了干岩石骨架的模量比为

(11)

图 1a是利用式(11)计算的孔隙纵横比α=1，0.2，0.1，0.05，0.01，0.001时石英砂岩的干岩石骨架模量比随孔隙度变化特征.这里将石英作为背景基质，其体积和剪切模量分别取K_m=37 GPa和G_m=44 GPa.可以看到，当α＞0.1时，干岩石骨架的模量比随着孔隙度的增加而增大.当α=0.1时，干岩石骨架的模量比几乎为常数，与孔隙度无关.当α＜0.1时，干岩石骨架的模量比随孔隙度的增加而减小.可见，假设干岩石骨架模量比为常数的临界孔隙度模型和模量比随孔隙度增加而增大的Pride模型仅为式(9)-(11)的特例.图 1b是利用式(11)计算的孔隙纵横比α=1时石英砂岩和长石砂岩的干岩石骨架模量比随孔隙度变化情况，长石砂岩的弹性模量取为K_m=37.5 GPa和G_m=15 GPa.可以看出，对于相同孔隙形状的砂岩，长石砂岩的干岩石骨架模量比随孔隙度的增加而减小，与石英砂岩明显不同. Knackstedt等(2003)运用有限元法来估计模型系统的弹性性质得到类似的结果^[6].图 1表明，式(11)可以很好地描述不同岩性干岩石骨架模量比随孔隙度的变化特征，克服了经验模型存在的局限性，比常规的骨架模型如临界孔隙度模型等适用范围更广.

图 1

Fig. 1

图 1 干燥岩石的体积和剪切模量比随孔隙度变化特征(a)不同孔隙纵横比的石英砂岩; (b)含球孔的石英砂岩与长石砂岩. Fig. 1 Volume and shear modulus ratios versus porosity for dry rock (a) Quartz sandstone with various pore aspect ratios; (b) Quartz and feldspar sandstone with spherical pore.

式(9)和(10)表明：干岩石骨架的弹性模量不仅是矿物颗粒的体积和剪切模量以及孔隙度的函数，而且与岩石的孔隙纵横比密切相关.把式(9)-(10)分别代入到式(1)-(4)中，我们就建立了纵、横波速度与孔隙纵横比之间的非线性关系式，利用这样的关系式，我们既可以根据岩石基质矿物组分的体积和剪切模量以及孔隙度、流体饱和度和孔隙纵横比正演得到岩石的纵、横波速度，也可以根据岩石的纵波(和横波)速度来反演等效孔隙纵横比.为了利用式(1)-(4)和(9)-(10)反演岩石孔隙的纵横比，这里我们把岩石的孔隙简单等效为具有单一纵横比的理想椭球孔，反演中采用的目标函数OF为

(12)

式中，K_sat^cal和G_sat^cal分别为预测的饱和岩石的体积和剪切模量，W_p，W_s为加权因子，W_p+W_s=1.只有纵波资料时，W_p=1，W_s=0；同时有纵、横波资料时，可取W_p=0.5，W_s=0.5.我们采用具有非线性全局寻优的模拟退火算法来求解式(12).当没有横波测井资料时，可由纵波反演等效孔隙纵横比来构建横波实现横波速度预测，而利用纵、横波同时反演的等效孔隙纵横比所构建得到的纵、横波速度则可用于速度评价，因为实际纵、横波速度测井数据常存在噪音，由其计算的泊松比有时会出现负值异常，特别是在井孔扩径段，声波时差受周波跳跃的影响会出现异常值，横波受到的影响更大，这时的时差数据往往都不可靠，需要进行必要的校正.

我们首先采用纵、横波两种资料同时反演来研究岩石孔隙纵横比随有效压力的变化情况.图 2a是使用Coyner测量的Westerly花岗岩数据^[22]进行等效孔隙纵横比反演得到的结果. Coyner(1977)给出了干燥花岗岩样品在有效压力为5至100MPa时的弹性模量和孔隙度数据.岩石的孔隙度值在0MPa时为0.8%，随着压力的增大孔隙度缓慢降低，在100MPa时降至0.71%左右.反演结果表明，花岗岩样品的孔隙纵横比随有效压力的增加而增大，当有效压力从5增加到100MPa时，花岗岩样品的孔隙纵横比从0.009增加到0.023.图 2b是使用Coyner测量的Navajo石英砂岩数据^[22]进行等效孔隙纵横比反演得到的结果. Coyner(1977)测量了干燥砂岩样品在有效压力为0至100MPa时的纵横波速度，并给出了不同有效压力下的孔隙度值，变化范围在11.5%~11.8%之间^[22].反演结果表明，砂岩样品的孔隙纵横比随有效压力的增加而增大，当有效压力从0增加到100MPa时，砂岩样品的孔隙纵横比从0.08增加到0.15. 图 2表明，本文的方法既可以反演得到像裂缝型储层如花岗岩的低孔隙纵横比，又可以反演得到像孔隙型储层如砂岩的高孔隙纵横比，利用这些孔隙纵横比数据就可以对储层的孔隙结构进行评价.根据测量的速度估计孔隙纵横比分布最早是由Cheng和Toksöz(1979)提出^[23]，他们是针对一系列不同有效压力情形下测量的速度数据来估计岩石的孔隙结构，由于是基于Kuster-Toksöz理论通过迭代的方法逐步求出，计算过程复杂，而且纵横比分布的反演结果具有非唯一性，依赖于预先给定的纵横比，特别是大孔隙的纵横比^[24]；另一方面，该方法仅限于实验室测量的岩石物理数据，而难以有效应用到测井数据中，因为无法同时获得目的层段不同有效压力下的测量数据.此外，Kuster-Toksöz理论本身适合于低孔隙度岩石的等效弹性模量估计，其应用效果受到理论假设限制^{[8, 25]}.

图 2

Fig. 2

图 2 由纵、横波同时反演得到的岩石孔隙纵横比随有效压力的变化特征(a)花岗岩; (b)砂岩. Fig. 2 Pore aspect ratios inverted from P and S-wave data versus pressure (a) Granite; (b) Sandstone.

其次，我们采用Han^[18]测量的75个饱和砂岩的岩石物理数据来分析等效孔隙纵横比反演在横波速度预测方面的应用效果.这些数据是在有效压力40MPa下测量的，砂岩样品的孔隙度分布范围为5%~30%，粘土含量分布范围为0~50%.用于分析的矿物石英的体积和剪切模量分别为41.42GPa和37.87MPa，粘土的体积和剪切模量分别为23.98MPa和6.96 MPa，孔隙流体的体积模量根据Xu和White提供的声波时差和密度计算得到^[10].对于泥质砂岩，其基质体积和剪切模量使用Voigt-Reuss-Hill平均方法计算^[1].图 3a为反演的等效孔隙纵横比与孔隙度之间的关系，实心点为只用纵波速度资料反演得到的等效孔隙纵横比，而空心点为用纵、横波速度同时反演得到的等效孔隙纵横比. 75个数据中有3个样品的孔隙纵横比大于0.4，其孔隙结构接近于球型，图中没有显示.反演结果表明，空心点与实心点比较接近，说明单纯利用纵波资料反演等效孔隙纵横比是可信的，进一步证明我们提出的由孔隙纵横比预测横波的方法是正确的.对于砂岩来说，总体上，在孔隙度大于15%时，孔隙纵横比分布比较稳定，在0.15上下变化；而在孔隙度小于15%时，孔隙纵横比的变化范围较大.泥质含量低，岩石的孔隙纵横比低，但是孔隙纵横比与泥质含量相关性差.根据反演的等效孔隙纵横比和式(9)-(10)以及式(1)-(4)可以正演得到纵波和横波速度.图 3b是预测与测量的横波速度之间的对比，可以看到，绝大多数样品的预测值与实测值吻合良好，仅有两个样品的横波速度预测结果与实测存在较大偏差，这充分说明碎屑砂岩的孔隙结构可以用椭球形状来表征.我们这里对岩石孔隙具有单一纵横比的假设与Xu-White模型不同，Xu-White模型使用两个孔隙纵横比，一个用于砂岩，另一个用于泥岩，但这两个孔隙纵横比值在横波速度预测中必须固定，图 3c是Xu和White采用同样的砂岩样品所预测的横波速度与测量的横波速度之间对比^[10].可以看到，本文的结果与Xu-White模型的预测结果相比精度更高，Xu-White模型预测的横波速度整体偏高，这说明我们提出的基于孔隙纵横比反演的横波速度预测方法更加合理可信.

图 3

Fig. 3

图 3 砂岩实验室测量数据的等效孔隙纵横比反演与横波速度预测效果分析(a)反演的等效孔隙纵横比与孔隙度关系; (b)本文DEM解析模型预测与测量的横波速度对比; (c) Xu-White模型预测与测量的横波速度对比[10]. Fig. 3 Inverted effective pore aspect ratios and predicted shear velocity for sandstone measured in laboratory (a) Relation between inverted effective pore aspect ratio and porosity; (b) Comparison of predicted and measured shear velocities for DEM analytical model; (c) Comparison of predicted and measured shear-velocities for Xu-Whitemodel[10].

我们采用渤海湾一口气井的测井数据来评价等效孔隙纵横比反演和横波速度预测结果，该井储层为砂岩，储层段在2712~2880m之间，其中2712~2778m为气层，2778~2880m则为水层.

图 4a从左到右的第一栏至第三栏分别显示了目的层段的纵波速度、密度、泥质含量、孔隙度和含水饱和度曲线.第四栏为反演得到的等效孔隙纵横比数据，图中红色显示的为只用纵波速度反演的等效孔隙纵横比，而黑色则是用纵、横波速度同时参与反演的等效孔隙纵横比.两种反演结果几乎完全吻合.图 4b为等效孔隙纵横比分布图，孔隙纵横比主要分布在0.18~0.32之间，大部分位于0.2附近，说明该段储层为孔隙型储层.此外，还有一部分储层的等效孔隙纵横比达到1，接近于球孔，这部分对应的纵波速度相对都很高.图 4a中的第五栏为实测横波速度与预测横波速度对比，其中红色曲线为由第四栏中只由纵波速度反演的等效孔隙纵横比正演得到的横波速度曲线，而黑色曲线则为实际测量的横波速度.可以看到，由纵波中提取的横波与实测横波吻合良好，几乎一致，说明本方法从实测声波资料中提取的横波信息是可靠的，可以得到准确的结果.

图 4

Fig. 4

图 4 某气井的等效孔隙纵横比反演与横波速度预测效果分析(a)测量及反演的测井曲线; (b)等效孔隙纵横比频数直方图. Fig. 4 Inverted pore aspect ratios and predicted shear velocity for a gas well (a) Measured and inverted log curves; (b) Histogram of the effective pore aspect ratios.

本文提出了根据纵波(和横波)速度反演岩石等效孔隙纵横比进行储层孔隙结构评价和横波速度预测的方法.首先，根据已知岩石的纵波(和横波)速度、密度、孔隙度、饱和度和泥质含量等参数，运用Gassmann方程和DEM解析模型建立各参数之间的关系，应用非线性全局寻优算法来反演等效孔隙纵横比；其次，将反演得到的等效孔隙纵横比代入到Gassmann方程和干岩石骨架模型中计算横波速度，以及利用孔隙纵横比进行储层孔隙结构评价.

与现有的经验模型相比，由微分等效介质理论建立的DEM解析模型，用孔隙纵横比描述了孔隙的几何特征，弥补了经验模型无法描述孔隙形状的不足，突破了经验模型如临界孔隙度模型所假设的模量比为常数，与孔隙度无关的局限，可以准确反映不同岩性的弹性模量及其模量比随孔隙度的变化特征.

实验室和井孔测量数据应用表明，反演得到的等效孔隙纵横比可以准确对储层的孔隙结构进行评价.对于裂缝型储层如花岗岩，其孔隙纵横比通常小于0.025，而对于孔隙型储层如砂岩，其孔隙纵横比通常大于0.08.单纯利用纵波与同时利用纵、横波反演等效孔隙纵横比的结果几乎完全一致.并且，由纵波中提取的横波与实测横波吻合良好，这说明本文提出的等效孔隙纵横比反演及其横波速度预测方法是有效的.

本文通过将岩石孔隙等效为具有单一纵横比的理想椭球孔，从纵波数据中来反演等效孔隙纵横比用于储层的孔隙结构评价，这种评价方法比较适合于实际储层孔隙类型为单一的孔隙型或裂缝型储层.对于多重孔隙类型的复杂储层，依然可以利用本文的方法反演得到一个等效的孔隙纵横比用于横波速度预测和评价，至于如何从声波数据中反演出各自的孔隙纵横比进行储层孔隙评价将是下一步亟需解决的问题.

感谢匿名评阅人对本文提出的宝贵修改意见.