2017, Vol. 60
地下固体介质弹性模量是固体地球定量表征的重要参数，固体介质弹性模量越大，地震波速度越大，地层的抗压、抗形变能力越强，越不容易遭到破坏.因此，地下固体介质弹性模量反演对于地球圈层的认知、工程建设或煤田安全生产方案的制定都是非常重要的.在油气工程中，地下固体介质的弹性模量越大，地层就越坚固，这也是钻井方案设计的重要依据之一.
近年来随着油气资源发现难度的增大，储集空间品质变差及钻探成本的提高，地下流体识别及横波测井预测越来越受到人们的重视，发展了多种方法(Murphy et al., 1993；Yin and Zhang, 2014；石玉梅等, 2014, 2016；Han and Batzle, 2004；熊晓军等，2012；林凯等，2013；Grana，2016)，这些方法中准确的固体介质弹性模量非常重要.
由于很难通过实验室测量，目前的地下介质固体基质弹性模量基本上都是采用简化方法获得.这些方法可分为两大类：矿物组分法和统计分析法.矿物组分法通常根据组成固体介质各组分的模量及所占的体积比例，利用不同的平均模式进行计算(Mavko et al., 2009).实际上，地层固体基质弹性模量除与其矿物组分有关外，成岩作用(如压实、胶结等)也是不可忽视的重要因素，同样的矿物组分，不同的固结程度，其固体基质的弹性模量具有很大的差别，这也是组分类似的同类岩石，它们的纵、横波速度具有很宽变化范围的根本原因.统计法利用速度与孔隙度的经验关系外推至孔隙度为零时估计基质固体弹性模量.统计法一方面需要大量的数据进行统计拟合，获得经验公式，另一方面，外推的结果有时存在非常大的误差(Han and Batzle, 2004).近年来，基于岩石物理模型的固体基质弹性模量反演取得了很大的进展(孙成禹等，2015；Lin et al., 2014；熊晓军等，2012；Dupuy et al., 2016a, 2016b；Grana，2016).孙成禹等(2015)按岩性划分固体基质(如砂岩、泥岩等)，然后利用Gassmann方程和Voigt-Reuss-Hill平均，采用自适应模拟退火算法从测井数据中反演单一岩性固体基质的弹性模量，如砂岩的固体基质弹性模量.Lin等(2014)和熊晓军等(2012)基于Gassmann-Biot-Geertsma方程和Russell流体因子公式分别计算流体项，并通过使这两种方法计算结果差值的极小化，迭代反演骨架固体体积模量和干骨架泊松比.但Gassmann-Biot-Geertsma方程和Russell流体因子公式两种方法本身就存在一定的差异，这种差异对反演结果的影响需要进一步讨论.Dupuy等(2016a, 2016b)基于广义动态孔隙弹性模型，采用邻域优化算法从粘弹性地震属性参数(如速度，密度，品质因子)中反演孔隙介质弹性参数(包括固体基质弹性模量、孔隙度、渗透率、流体参数等)，并通过对不同组合粘弹性参数反演结果的误差分析，详细讨论了孔隙弹性参数反演的敏感性或可行性.由于采用的模型没有涉及表征孔隙结构特征的参数以及文中实际地震数据处理(合成或野外采集)时假设固体基质弹性模量已知，因此，反演中忽略了孔隙形态对结果的影响.
本文将地层固体部分统称为固体基质，例如，含有泥质的砂岩地层，其固体基质就是由其中的泥质和砂质两部分组成.由于地下地层多含有孔隙，尤其是含流体的储集层，因此，本文首先简要分析常用的流体饱和孔隙介质岩石物理型，在此基础上，通过数值计算分析流体饱和多孔介质等效弹性模量与固体基质、孔隙度、孔隙结构和孔隙流体的关系，探讨利用等效介质弹性模量反演固体基质弹性模量的可行性，研究基于岩石物理模型的地下固体介质弹性模量反演方法，并通过对三类岩石样本实验室测量数据的反演和结果分析，验证方法的有效性，同时通过不同条件下对实验数据的反演，进一步讨论影响方法效果的因素，为方法的合理应用提供依据.
2 流体饱和孔隙介质弹性模量地球物理测量获得的是地下介质(即地层)的等效弹性模量，它主要取决于固体基质弹性模量、孔隙度、孔隙形态和孔隙流体这四个因素.这里，我们首先基于岩石物理模型讨论这四种因素对等效介质弹性模量的贡献.
2.1 流体饱和孔隙介质岩石物理模型目前，描述流体饱和孔隙介质弹性特征的岩石物理模型有多种，常用的有Gassmann模型、Xu-White模型、微分等效介质模型、Kuster-Toksöz模型等.Gassmann方程(Gassmann，1951)形式简单，容易使用，但由于只涉及到地层的总孔隙度，没有表征孔隙形态的参数，且假设饱和流体的等效介质剪切模量等于干骨架的剪切模量，这就隐含了球形孔隙的假设，因此，很难描述含有不同形态孔隙的实际地层.Xu-White模型(Xu and White, 1995)通过引入孔隙纵横比使研究饱和流体的泥质砂岩纵、横波速度特征变得容易，但其中砂和泥的孔隙纵横比固定，即各自只采用一个纵横比孔隙，而实际地层中的孔隙一般都存在多种形态(Timur et al., 1971；Brace et al., 1972；Nur and Simmons, 1969；Brown and Korringa, 1975；Yan et al., 2002；Seymour，2005)，采用单一纵横比的孔隙可能引起较大的预测偏差(Kuster and Toksöz，1974；Toksöz et al., 1976；白俊雨等，2012).近些年人们关注比较多的微分等效模型(Norris，1985；Zimmerman，1984；Berryman，2002；Cho et a1.，2013)考虑了多重孔隙情况，能更好地表征实际地层，但由于缺乏解析形式，不易于使用，目前，主要用于研究干燥或饱和流体情况下的等效介质弹性特征.Kuster和Toksöz(1974)基于散射理论，导出了长波长一阶近似条件下流体饱和孔隙介质弹性模量方程，即Kuster-Toksöz方程.该方程以“孔隙纵横比谱”的形式表征多重孔隙，能更好地描述实际地层，且相对于微分等效介质模型，直接表征了地层的等效弹性模量与固体基质模量、孔隙流体模量以及孔隙(孔隙度和孔隙形态)的关系，易于使用.本文，我们采用Kuster-Toksöz方程进行数值计算和反演方法研究.
Kuster-Toksöz方程(Kuster and Toksöz，1974；Toksöz et al., 1976)如下：
|
(1a) |
|
(1b) |
其中，Ke、Ks和K′分别为等效介质、固体基质和孔隙填充物的体积模量；μe、μs和μ′分别为等效介质、固体基质和孔隙填充物的剪切模量；α为孔隙纵横比；c(α)为纵横比等于α的孔隙所占的比例，
四阶张量Tiijj和Tijij为Ks、μs和K′、μ′以及α的函数，按下列方法计算(Toksöz et al., 1976)，
|
(2a) |
|
(2b) |
其中
|
下面，我们基于方程(1)，通过计算分别讨论固体基质、孔隙度及孔隙形态、孔隙流体对地下等效介质弹性模量的贡献.由于固体基质弹性模量或孔隙度或孔隙纵横比或饱和度具有不同的量纲，为便于研究等效介质弹性模量对这些参数的敏感性，我们定义了等效介质弹性模量相对变化率.
记Me=(Ke，μe)为等效介质弹性模量，x为归一化后的固体基质弹性模量或孔隙度或孔隙纵横比或饱和度，Δx为归一化后参数的变化量，ΔMe为参数x变化Δx时引起的等效介质弹性模量的变化量，按下式计算等效介质弹性模量的相对变化率，
|
(3) |
由上式可知，相对变化率越大，敏感性越大，从等效弹性模量中提取相应参数的可靠性越大.
(1) 等效介质弹性模量与固体基质
方程(1)中保持孔隙度、孔隙纵横比和孔隙流体(即孔隙填充物)不变, 并取ϕ=15%，K′=2.25 GPa，μ′=0，先后变化基质体积模量Ks和剪切模量μs，计算相应的等效介质体积模量Ke和剪切模量μe，结果如图 1所示.图中显示了孔隙纵横比分别为1、0.5和0.1三种情况下等效介质弹性模量随基质弹性模量变化曲线(图 1a为保持μs不变(μs=44 GPa)，Ks从28 GPa变化到48 GPa时Ke和μe的变化曲线；图 1b为保持Ks不变(Ks=38 GPa)，μs从34 GPa变化到54 GPa时Ke和μe的变化曲线).可以看出，Ke和μe随Ks或μs的变化基本都呈线性变化，且基质体积模量对等效介质体积模量的影响比较大，对等效剪切模量影响小；同样，基质剪切模量对等效介质剪切模量的影响比较大，而对其体积模量影响小.Ks变化引起的Ke和μe的相对变化率分别为：①α=1时，dKe≈38.3，dμe≈1.2；②α=0.5时，dKe≈36.8，dμe≈1.5；③α=0.1时，dKe≈15.4，dμe≈3.2.μs变化引起的Ke和μe的相对变化率分别为：①α=1时，dKe≈3.1，dμe≈42.7；②α=0.5时，dKe≈3.6，dμe≈41.5；③α=0.1时，dKe≈6.5，dμe≈23.9.可见，固体基质的弹性模量对等效介质模量的影响比较大，且其程度与介质中孔隙的形态有一定的关系，孔隙越接近球形，影响越大.
|
图 1 等效介质弹性模量(Ke和μe)随基质弹性模量(Ks和μs)变化曲线 (a) s=44 GPa时，Ke和μe随Ks变化曲线；(b) Ks=38 GPa时，Ke和μe随μs变化曲线. Fig. 1 Effective elastic modulus (Ke and μe) versus matrix elastic modulus (Ks and μs) (a) Ke and μe versus Ks at μs=44 GPa; (b) Ke and μe versus μs at Ks=38 GPa. |
(2) 等效介质弹性模量与孔隙度
方程(1)中保持固体基质、孔隙纵横比和孔隙流体不变(分别取Ks=38 GPa，μs=44 GPa，K′=2.25 GPa，μ′=0)，计算孔隙度从0%变化到35%时等效介质弹性模量，结果如图 2所示(类似于固体基质情况，图 2中也显示了3个纵横比孔隙情况下等效介质体积模量(a)与剪切模量(b)随孔隙度的变化曲线.图中可见，①等效介质体积模量和剪切模量随孔隙度的增大都呈比较明显的减小趋势；②含有小纵横比孔隙介质的等效体积模量和剪切模量随孔隙度的增大减小的速度快.孔隙度变化引起的等效介质体积模量和剪切模量的相对变化率分别为：①α=1时，dKe≈-47.1，dμe≈-66.5；②α=0.5时，dKe≈50.3，dμe≈70.1；③α=0.1时，dKe≈115.4，dμe≈140.2.可以看出，在其它因素不变的情况下，孔隙度对等效介质弹性模量的影响比较大，且其影响程度受孔隙形态的影响也非常明显，尤其是介质中含有小纵横比孔隙.
|
图 2 不同孔隙纵横比情况下等效介质弹性模量随孔隙度变化曲线 (a)等效体积模量；(b)等效剪切模量. Fig. 2 Effective elastic modulus versus porosity for different aspect ratios (a) Effective bulk modulus; (b) Effectve shear modulus. |
(3) 等效介质弹性模量与孔隙形态
方程(1)中保持固体基质、孔隙度和孔隙流体不变(分别取Ks=38 GPa，μs=44GPa，=15%，K′=2.25 GPa，μ′=0)，我们计算了8种纵横比孔隙情况下等效介质的弹性模量，如图 3所示.图中可见，①等效介质体积模量和剪切模量受孔隙形态的影响程度基本相同；②对于含有小纵横比(α=0.05~0.125)孔隙介质，其等效弹性模量受孔隙形态的影响比较大，相对变化率分别为dKe≈130.4和dμe≈ 155.4；随着孔隙纵横比的增大，等效弹性模量的变化变得比较平缓，孔隙纵横比在0.125~0.5之间时相对变化率分别为dKe≈20.8和dμe≈22.1；孔隙纵横比大于0.5时，等效弹性模量随孔隙纵横比的变化变得非常平缓，相对变化率只有dKe≈1.2，dμe≈1.5，基本可以忽略.可见，根据等效介质弹性模量可以很好地预测裂隙介质中裂隙的纵横比，但很难预测α>0.5介质中孔隙的纵横比.
|
图 3 等效介质弹性模量随孔隙纵横比变化 Fig. 3 Effective elastic modulus versus aspect ratios of pore |
(4) 等效介质弹性模量与孔隙流体
保持固体基质和孔隙度及孔隙纵横比不变(取Ks=38 GPa，μs=44 GPa，ϕ=15%)，孔隙流体为气或水或气水混合物(气的体积模量Kg=0.00183 GPa，水的体积模量Kw=2.25 GPa)，由于水和气的粘滞系数非常小，计算时取μ′=0.混合流体体积模量用Brie经验公式计算(取e=3).图 4显示了3种纵横比孔隙介质不同饱和度对应的等效弹性模量.图中可以看出，①饱和度对等效介质弹性模量的影响比较小，尤其是对于含有大纵横比的孔隙介质；②饱和度对等效介质体积模量的影响比剪切模量略大；③饱和度对等效介质弹性模量的影响程度与孔隙形态有关，大纵横比孔隙介质(如α=0.5，1.0)，等效弹性模量随饱和度的变化非常平缓，说明饱和度对等效弹性模量影响小，含有小纵横比的孔隙介质(如α=0.1)，等效体积模量随饱和度变化比较明显，且其剪切模量随饱和度变化也比前两者明显.三种纵横比孔隙介质等效弹性模量的相对变化率分别为：①α=1时，dKe≈0.74，dμe=0；②α=0.5时，dKe≈0.87，dKe≈0.02；③α=0.1时，dKe≈4.71，dμe≈1.08.可见，存在小纵横比孔隙的介质，其等效体积模量对孔隙流体的变化比较敏感.
|
图 4 不同纵横比孔隙情况下等效介质弹性模量随饱和度变化曲线 Fig. 4 Effective elastic modulus versus saturation for different aspect ratios of pore |
综合比较上述四种因素变化对应的等效介质弹性模量相对变化率，总体来说，孔隙度对等效介质弹性模量的影响最大，固体基质次之，孔隙形态的影响受孔隙纵横比控制，小纵横比孔隙影响大，大纵横比孔隙影响比较小，孔隙流体的影响最小，因此，从等效介质弹性模量中能很好地预测介质中的孔隙度、固体基质弹性模量和裂隙介质中孔隙的纵横比，但对于流体饱和度和近球形孔介质中孔隙的纵横比，利用等效介质弹性模量预测难度比较大.
3 固体基质弹性模量反演Kuster-Toksöz方程中没有干骨架参数，从而避免了多个未知数求解的困难，但Kuster-Toksöz方程中固体基质体积模量Ks和剪切模量μs不仅显式地表现在方程(1a)和(1b)中，张量Tiijj和Tijij中也含有Ks、μs信息，因此，直接利用方程(1a)和(1b)从测量的地层纵、横波速度及密度或体积模量和剪切模量中获取固体基质弹性模量是非常困难的.
3.1 固体基质弹性模量的非线性方程系统及数值求解
记
|
在孔隙填充物为流体，且其剪切模量可忽略不计(即μ′=0)的情况下，方程(1a)和(1b)可改写为
|
(4a) |
|
(4b) |
为简洁起见，分别用K和μ表示Ks和μs，重新整理方程(4a)和(4b)，得
|
(5a) |
|
(5b) |
其中，
方程(5a)和(5b)组成了一个二维非线性方程系统.该方程系统中系数ai、bi、ci和di(i=1, 2)都与固体基质弹性模量Ks和μs有关，因此不易直接求解.这里，我们采用Newton-Raphson迭代求解方法.
分别用函数f和g表示方程系统(5)中(5a)和(5b)的左端项，
|
(6a) |
|
(6b) |
并记m =(K, μ)(m为以K和μ为分量的向量)，F =(f, g)(F为以函数f和g为分量的函数向量).定义目标函数，
|
(7) |
式中上标“t”表示转置.
非线性方程系统(5)的解就是要求使目标函数S极小，且同时满足F=0.
用泰勒级数展开F，
|
(8) |
其中，δm =(δK, δμ)(δm为m的牛顿修改步长)；J为雅可比矩阵，
|
(9) |
忽略二次项，并置F (m +δ m)=0，由式(8)得
|
(10) |
在每步迭代中适当调整(10)式获得的δm，使其满足
|
(11) |
并用下式进行迭代修改，
|
(12) |
即沿目标函数下降的方向对方程系统(5)的根进行修改，直至S < ε(ε为在一定精确度下给定的适当小的正数).
该方法与函数极小化法类似，不同之处即采用使F =0的牛顿步长求目标函数S极小化，而不是采用使▽S=0的牛顿步长求S的极小化.
取K′=μ′=0，这时，方程(1a)和(1b)中的等效弹性模量Ke和μe即为干骨架弹性模量Kd和μd.重新整理方程(1a)和(1b)，可得干骨架弹性模量公式:
|
(13a) |
|
(13b) |
反演出固体基质弹性模量后，利用方程(13a)和(13b)就可以计算出干骨架的体积模量和剪切模量.
基于Kuster-Toksöz方程反演固体基质弹性模量需要已知等效介质弹性模量、孔隙填充物弹性模量(如为流体即为流体体积模量)和孔隙参数(即孔隙的体积分数和孔隙纵横比谱).实际工作中，等效介质弹性模量和孔隙度可由实验室测量或测井数据或地震反演获得，孔隙纵横比谱一般由薄片的电镜扫描获得.
3.2 初始值的确定方程系统(5)迭代求解时初始值的选择非常重要，这里，我们采用下列方法估计固体基质弹性模量的初始值m0.
(1) 根据实验室或地球物理测量获得的纵波速度(VP)、横波速度(VS)和密度(ρ)计算等效介质弹性模量：
|
(14a) |
|
(14b) |
(2) 利用下式计算初始值m0：
|
(15a) |
|
(15b) |
|
(15c) |
其中，me=(Ke, μe)为地层的等效弹性模量；m ′=(K′, μ′)为孔隙填充物的弹性模量；ϕ为孔隙度. m ′和ϕ可由先验信息获得.
由(15)式可以获得比较好的固体基质弹性模量初始估计值(未引Mavko et al., 1998)，从而保证Newton-Raphson方法有效快速收敛.
4 实验数据反演这里，我们利用Toksöz等(1976)实验室测量数据进行反演，验证方法的可行性和有效性.
4.1 实验数据准备表 1—4数据来源于文献(Toksöz et al，1976)中的实验室测量数据.表 1—3列出了三种岩石(Biose砂岩、Bedford灰岩、Troy花岗岩)的基础数据，即固体基质弹性模量、密度、孔隙度、孔隙谱和孔隙流体密度、体积模量、声波速度，表 4列出了三种岩石分别饱和水或盐水时实验室测量的纵、横波速度(由该文献中图 3、图 4和图 5拾取获得).三种岩石的孔隙度和孔隙结构特征对三类岩石(砂岩、灰岩、花岗岩)具有一定的代表性.
|
表 1 固体基质弹性模量、密度和孔隙度(Toksöz等, 1976) Table 1 Matrix properties of rocks and porosities (Toksöz et al., 1976) |
|
|
表 2 孔隙谱(压力P=1 bar(105Pa))(Toksöz等, 1976) Table 2 Pore Spectra at P=1 bar(105Pa)(Toksöz et al., 1976) |
|
|
表 3 饱和流体性质(压力：0.1 MPa，20 ℃) Table 3 Properties of saturating fluids at 20 ℃ and pressure of 0.1 MPa |
|
|
表 4 不同压力下测量的三种岩石纵、横波速度(单位：m·s-1)(Toksöz等, 1976) Table 4 Compressional and shear wave velocities at differential pressure (Toksöz et al., 1976) |
|
图 5 不同压力下Troy花岗岩固体基质反演结果 (a)体积模量; (b)剪切模量. Fig. 5 Inversion results of matrix elastic modulus for Trov granite at different pressure (a) Bulk modulus; (b) Shear modulus. |
根据前面所述的方法，我们对表 4所列的三种岩石不同压力下测量数据进行了反演，结果如图 5—7所示(图中散点为反演结果，虚线为真实值).图 5为Troy花岗岩反演结果.可以看出，反演的剪切模量相对比较集中，分布在真实值附近，反演的体积模量整体偏低，但接近真实值，表明，反演获得了比较好的Troy花岗岩固体基质弹性模量估计值.图 6显示了Boise砂岩的反演结果.与Troy花岗岩结果类似，Boise砂岩反演的剪切模量分布在真实值附近，反演的体积模量在小压力时偏低，高压力下偏高，最大反演误差约为8%.图 7为Bedford灰岩反演结果.反演的Bedford灰岩体积模量值比较分散，整体偏高，但剪切模量与前两岩石结果类似，与真实值很接近.
|
图 6 不同压力下Boise砂岩固体基质反演结果 (a)体积模量; (b)剪切模量. Fig. 6 Inversion results of matrix elastic modulus for Boise sandstone at different pressure (a) Bulk modulus; (b) Shear modulus. |
|
图 7 不同压力下Bedford灰岩固体基质反演结果 (a)体积模量; (b)剪切模量. Fig. 7 Inversion results of matrix elastic modulus for Bedford limestone at different pressure (a) Bulk modulus; (b) Shear modulus. |
我们对不同压力下测量值的反演结果进行平均，结果列于表 5.可以看出，三种岩石不同压力下测量数据反演结果的平均能够更好地估计岩石基质的弹性模量.
|
|
表 5 不同压力下三种岩石固体基质弹性模量反演结果平均值 Table 5 Mean value of inverted matrix modulus at differential pressure |
利用公式(13a)和(13b)和表 5列出的基质固体弹性模量平均值，我们计算了三种岩石干骨架弹性模量，Troy花岗岩(ϕ=0.33%)：Kd=57.0886 GPa，μd=30.6909 GPa；Boise砂岩(ϕ=25%)：Kd= 10.94982 GPa，μd=9.59221 GPa；Bedford灰岩(ϕ=12.6%)：Kd=13.50419 GPa，μd=8.73914 GPa.根据文献中图 2、图 3和图 4实际测量的干骨架纵波速度和横波速度(Toksöz et al., 1976)，我们计算了三种岩石相应的干骨架体积模量和剪切模量，分别为Troy花岗岩：Kd=58.39487 GPa，μd=30.25094 GPa；Boise砂岩：Kd=11.68965 GPa，μd=9.92651 GPa；Bedford灰岩：Kd=14.28898 GPa，μd=8.41845 GPa.与实测值比较，预测的Troy花岗岩干骨架的体积模量和剪切模量误差分别为-1.30627 GPa和0.43996 GPa，Boise砂岩干骨架的体积模量和剪切模量误差分别为-0.73983 GPa和-0.3343 GPa，Bedford灰岩干骨架的体积模量和剪切模量误差分别为-0.78479 GPa和0.32069 GPa，基本上都反映了实际干骨架的弹性模量值.
三类岩石样本不同压力下的测量数据进行反演，都能获得比较可靠的固体介质弹性模量，表明上述基于Kuster-Toksöz方程的固体基质弹性模量反演方法具有比较好的稳定性.
4.3 误差分析记Ms=(Ks，μs)，为反演获得的固体基质弹性模量，Msr=(Ksr，μsr)为列于表 1的实际固体基质弹性模量，按下式计算反演的相对误差:
|
(16) |
图 8显示了这三种岩石样本测量数据反演结果的相对误差直方图.总体来看，剪切模量的误差比体积模量小.图中还可以看出，对于体积模量，Troy花岗岩反演结果相对误差较小，Boise砂岩和Bedford灰岩反演结果的相对误差相当；对于剪切模量，Troy花岗岩和Boise砂岩反演结果要好于Bedford灰岩.综合体积模量和剪切模量反演结果，相对于Bedford灰岩，Troy花岗岩和Boise砂岩反演的结果更好地反映了实际值.
|
图 8 不同压力下三种岩石固体基质弹性模量反演误差直方图 (a)体积模量; (b)剪切模量. Fig. 8 Inversion errors of matrix elastic modulus for the three rock at different pressure (a) Error of bulk modulus (ΔKs); (b) Error of shear modulus (Δμs). |
Kuster-Toksöz方程中含有四个参量，即固体基质弹性模量Ms、等效介质弹性模量Me、孔隙填充物弹性模量M′和孔隙纵横比谱.在利用地球物理测量获得的等效介质弹性参数反演固体基质弹性模量时，孔隙结构和孔隙填充物(在油气地震勘探中为流体)是两个主要的影响因素.由于描述孔隙结构和孔隙填充物两因数的参量以比较复杂的形式体现在方程系统的系数中，因此，这里我们通过不同条件下对实验测量的数据进行反演来讨论这两因数对反演结果的影响.
(1) 孔隙纵横比
孔隙纵横比在一定程度上反映孔隙形态，它是压力的函数.表 2列出了压力等于1 bar(105Pa)时Troy花岗岩、Boise砂岩和Bedford灰岩的孔隙纵横比谱数据，表 4数据是在不同压力下测得的纵、横波速度.我们以表 2所列孔隙谱为基础，对表 4的测量数据进行反演.
首先对不同压力下的测量数据采用同一孔隙谱数据进行反演，然后考虑不同压力下孔隙的可能闭合情况进行反演，例如，Boise砂岩，测量数据的压力范围为0.03125~0.5375×108Pa，反演时采用了3个孔隙谱：①直接采用表 2中所列的孔隙谱(简称“闭0”)，②表 2所列的孔隙谱中最小纵横比的孔隙(即纵横比为0.1×10-3的孔隙)闭合(简称“闭1”)，③表 2所列的孔隙谱中最小两个纵横比的孔隙(即纵横比为0.1×10-3和0.5×10-3的孔隙)闭合(简称“闭2”)，…….表 6列出了其中某一压力下不同孔隙谱的反演结果.表中可以看出，不同孔隙结构对固体基质体积模量反演结果的影响小，对剪切模量的反演结果影响大，如Boise砂岩，三个孔隙谱情况下孔隙度的最大差异(闭2与闭0)只有0.048%，反演的体积模量最大差异为0.99508 GPa，相对差异为3.3%，而剪切模量的最大差异为1.84551 GPa，相对差异达10.2%.Bedford和Troy结果类似.可见，孔隙结构参数对固体基质剪切模量反演非常重要，且裂隙越发育，这种影响越大.
|
|
表 6 不同孔隙结构情况下反演结果 Table 6 Inverted matrix modulus for differential pore structure |
(2) 孔隙流体
利用Kuster-Toksöz方程从等效介质弹性模量中反演固体介质弹性模量，除了需要已知孔隙纵横比谱外，还需要已知孔隙填充物(及孔隙流体)的弹性参数.为了说明孔隙流体对反演结果的影响，我们以Boise砂岩为例，假定孔隙纵横比谱不变，对盐水饱和情况下的测量数据(表 4所列)分别用盐水、水和煤油三种孔隙流体进行反演，图 9显示了不同压力下测量数据的反演结果.可以看出，孔隙流体为水或盐水时，两者的反演结果接近，尤其是剪切模量，两者几乎重合.相对于前两者，孔隙流体为煤油时反演结果差异较大，体积模量最为明显.这三种孔隙流体的粘度都很小，剪切模量近似为零(赵邦六等，2013)，因此，反演结果的差异主要表现在体积模量上，水和盐水的体积模量差值为0.12 GPa，而煤油和盐水体积模量的差值为1.04 GPa.煤油与盐水或水情况下反演的剪切模量间的微小差异主要来自于密度的影响.
|
图 9 不同压力下孔隙中饱和盐水、水和煤油三种情况Boise砂岩反演结果 (a)体积模量; (b)剪切模量. Fig. 9 Inversion results of matrix elastic modulus of Boise sandstone saturated with brine, water and kerosene at different pressure (a) Bulk modulus; (b) Shear modulus. |
流体对反演结果的影响主要是孔隙流体的弹性模量.当孔隙中流体为单一成分时，孔隙流体弹性模量即为该成分流体的弹性模量，如孔隙中只填充水，则孔隙流体的弹性模量即为孔隙中水的弹性模量；当孔隙中流体为混合流体时，则孔隙流体弹性模量取决于各组分的弹性模量及各自的体积百分数(即饱和度)，不同的饱和度，混合流体的模量不同.这里，我们虽然只采用了盐水、水和煤油三种流体完全饱和情况，但三者的流体体积模量是不同，因此也反映了混合流体情况下饱和度对反演结果的影响.
(3) 等效介质弹性模量
本文是从等效介质弹性模量中反演固体介质弹性模量.等效介质弹性模量可以由测量的纵、横波速度和密度计算获得，也可以由地震反演获得.一般情况下，测量或反演都会存在一定的误差，那么等效介质弹性模量误差如何影响反演结果？这里，我们仍然以Boise砂岩为例.
Boise砂岩基质弹性模量真实值(表 1)：Ksr=30 GPa，μsr=18 GPa，孔隙度ϕ=25%，孔隙谱如表 2所列，完全水饱和情况下的等效弹性模量：Ker=15.0940 GPa，μer=8.3039 GPa.假设测量或反演获得的等效体积模量和剪切模量误差均为5%，即Ke=15.848721 GPa，μe=8.719116 GPa.我们分三种情况进行反演：①只有等效体积模量存在误差，②只有等效剪切模量存误差，和③等效体积模量和剪切模量同时存在误差.反演结果用Ksi和μsi表示.①当只有体积模量存在误差时，获得的Ksi=32.67 GPa，μsi=17.96 GPa，误差分别为：dKsi=8.9%，dμsi= -0.23%；②当只有剪切模量存在误差时，获得的Ksi=29.50 GPa，μsi=18.96 GPa，误差为：dKsi= -1.6%，dμsi=5.3%；③当体积模量和剪切模量同时存在5%的误差时，反演结果Ksi=32.08 GPa，μsi =18.91 GPa，反演误差为dKs=6.9%，dμs=5.0%.可见，对于Boise砂岩，①等效介质弹性模量误差引起的反演结果误差基本在同一数量级，且等效介质弹性模量误差对基质体积模量的影响比剪切模量大；②等效介质体积模量的测不准主要影响基质体积模量的反演结果，等效介质剪切模量的测不准主要影响基质剪切模量的反演结果.
6 结论地下固体介质弹性模量反演对于地球内部探测和定量表征具有重要的意义.基于Kuster-Toksöz方程的固体介质弹性模量反演充分考虑了地层中孔隙的体积百分比、孔隙纵横比和孔隙流体，能获得比较可靠的地下固体介质弹性模量反演结果.
利用Kuster-Toksöz方程从地层等效介质弹性模量中反演固体介质弹性模量时，需要的孔隙结构参数(即孔隙纵横比谱数据)可通过岩芯薄片电镜扫描获得，但在实际工作中需要考虑地层的地下压力状况.孔隙流体对反演结果有一定的影响，在流体粘滞性可以忽略的情况下，如气、水或轻质石油，孔隙流体对固体介质弹性模量反演的影响主要体现在体积模量上，因此，利用测井数据反演地下固体介质体积模量时，测井饱和度非常重要.另一方面，等效介质弹性模量是反演的基础数据，其数值越准确，反演结果的可靠性越大.
地下固体介质弹性模量反演，一方面为后续地震的定量表征(如流体识别、饱和度预测等)提供了基础数据，同时，通过地下固体介质弹性模量，我们也可以很好地认识地层(如硬度、固结程度等)，有助于设计合理的工程方案.
Kuster-Toksöz方程中除了等效介质弹性模量、孔隙体积百分比(即孔隙度)、孔隙填充物和孔隙纵横比谱外，没有涉及介质的其它参数，如各向异性、粘弹性等，因此，对于更为复杂的地下地层固体介质弹性模量反演，需要开展进一步的研究.
Bai J Y, Song Z X, Su L, et al.
2012. Error analysis of shear-velocity prediction by the Xu-White model. Chinese J. Geophys., 55(2): 589-595.
DOI:10.6038/j.issn.0001-5733.2012.02.021 |
|
Berryman J G. 2002. A differential scheme for elastic properties of rocks with dry or saturated cracks.//Proc. 15th ASCE Engineering Mechanics Conference. New York:Columbia University.
|
|
Brace W F, Silver E, Hadley K, et al.
1972. Cracks and pores:A closer look. Science, 178(4057): 162-164.
DOI:10.1126/science.178.4057.162 |
|
Brown R J S, Korringa J.
1975. On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics, 40(4): 608-616.
DOI:10.1190/1.1440551 |
|
Cho Y J, Lee W J, Park S K, et al.
2013. Effect of pore morphology on deformation behaviors in porous Al by FEM simulations. Advanced Engineering Materials, 15(3): 166-169.
DOI:10.1002/adem.201200145 |
|
Dupuy B, Garambois S, Virieux J.
2016a. Estimation of rock physics properties from seismic attributes-Part 1:Strategy and sensitivity analysis. Geophysics, 81(3): M35-M53.
DOI:10.1190/geo2015-0239.1 |
|
Dupuy B, Garambois S, Asnaashari A, et al.
2016b. Estimation of rock physics properties from seismic attributes-Part 2:Applications. Geophysics, 81(4): M55-M69.
DOI:10.1190/geo2015-0492.1 |
|
Gassmann F.
1951. Vber die elastizität poröser medien. Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich, 96(1): 1-23.
|
|
Grana D.
2016. Bayesian linearized rock-physics inversion. Geophysics, 81(6): D625-D641.
DOI:10.1190/geo2016-0161.1 |
|
Han D H, Batzle M L.
2004. Gassmann's equation and fluid-saturation effects on seismic velocities. Geophysics, 69(2): 398-405.
DOI:10.1190/1.1707059 |
|
Kuster G T, Toksöz M N.
1974. Velocity and attenuation of seismic waves in two-phase media:Part 1. Theoretical formulations. Geophysics, 39(5): 587-606.
DOI:10.1190/1.1440450 |
|
Lin K, He Z H, Xiong X J, et al.
2013. S-wave velocity inversion based on adaptive extraction of matrix mineral modulus. Oil Geophysical Prospecting, 48(2): 262-267.
|
|
Lin K, He Z H, Xiong X J, et al.
2014. AVO forwarding modeling in two-phase media:multiconstrained matrix mineral modulus inversion. Applied Geophysics, 11(4): 395-404.
DOI:10.1007/s11770.014-0457-x |
|
Mavko G, Mukerji T, Dvorkin J. 2009.
The Rock Physics Handbook:Tools for Seismic Analysis of Porous Media. Cambridge: Cambridge University Press.
|
|
Murphy W, Reischer A, Hsu K.
1993. Modulus decomposition of compressional and shear velocities in sand bodies. Geophysics, 58(2): 227-239.
DOI:10.1190/1.1443408 |
|
Norris A N.
1985. A differential scheme for the effective moduli of composites. Mech. Mater., 4(1): 1-16.
DOI:10.1016/0167-6636(85)90002-X |
|
Nur A, Simmons G.
1969. The effect of viscosity of a fluid phase on velocity in low-porosity rocks. Earth Planet. Sci. Lett., 7(2): 99-108.
DOI:10.1016/0012-821X(69)90021-1 |
|
Seymour R H. 2005. The Effects of Stress on the Seismic Velocities of Reservoir Rocks[Ph. D. thesis]. London:The University of London.
|
|
Shi Y M, Zhang Y, Yao F C, et al.
2014. Methodology of seismic imaging for hydrocarbon reservoirs based on acoustic full waveform inversion. Chinese J. Geophys., 57(2): 607-617.
DOI:10.6038/cjg20140224 |
|
Shi Y M, Cao H, Sun X P, et al.
2016. Formation modulus decomposition and its application to discrimination of pore fluid. Chinese J. Geophys., 59(11): 4278-4286.
DOI:10.6038/cjg20161128 |
|
Sun C Y, Li J J, Tang J, et al.
2015. Reservoir moduli inversion and seismic reflection chracteristics analysis based on logging data. Geophysical Prospecting for Petroleum, 54(3): 350-358.
DOI:10.3969/j.issn.1000-1441.2015.03.015 |
|
Timur A, Hempkins W B, Weinbrandt R M.
1971. Scanning electron microscope study of pore systems in rocks. J. Geophys. Res., 76(20): 4932-4948.
DOI:10.1029/JB076i020p04932 |
|
Toksöz M N, Cheng C H, Timur A.
1976. Velocities of seismic waves in porous rocks. Geophysics, 41(4): 621-645.
DOI:10.1190/1.1440639 |
|
Xiong X J, Lin K, He Z H.
2012. A method for S-wave velocity estimation based on equivalent elastic modulus inversion. Oil Geophysical Prospecting, 47(5): 723-727.
|
|
Xu S Y, White R E.
1995. A new velocity model for clay-sand mixtures. Geophysical Prospecting, 43(1): 91-118.
DOI:10.1111/gpr.1995.43.issue-1 |
|
Yan J, Li X Y, Liu E R.
2002. Effects of pore aspect ratios on velocity prediction from well-log data. Geophysical Prospecting, 50(3): 289-300.
DOI:10.1046/j.1365-2478.2002.00313.x |
|
Yin X Y, Zhang S X.
2014. Bayesian inversion for effective pore-fluid bulk modulus based on fluid-matrix decoupled amplitude variation with offset approximation. Geophysics, 79(5): R221-R232.
DOI:10.1190/geo2013-0372.1 |
|
Zhao B L, Shi Y M, Yao F C, et al.
2013. Prediction of the remaining oil distribution using multi-component seismic full waveform elastic inversion. Acta Petrolei Sinica, 34(2): 328-333.
DOI:10.7623/syxb201302015 |
|
Zimmerman R W.
1984. Elastic moduli of a solid with spherical pores:new self-consistent method. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 21(6): 339-343.
|
|
白俊雨, 宋志翔, 苏凌, 等.
2012. 基于Xu-White模型横波速度预测的误差分析. 地球物理学报, 55(2): 589–595.
DOI:10.6038/j.issn.0001-5733.2012.02.021 |
|
林凯, 贺振华, 熊晓军, 等.
2013. 基于基质矿物模量自适应提取横波速度反演方法. 石油地球物理勘探, 48(2): 262–267.
|
|
石玉梅, 张研, 姚逢昌, 等.
2014. 基于声学全波形反演的油气藏地震成像方法. 地球物理学报, 57(2): 607–617.
DOI:10.6038/cjg20140224 |
|
石玉梅, 曹宏, 孙夕平, 等.
2016. 地层模量分解及在流体识别中的应用. 地球物理学报, 59(11): 4278–4286.
DOI:10.6038/cjg20161128 |
|
孙成禹, 李晶晶, 唐杰, 等.
2015. 基于测井资料的储层模量反演与地震反射特征分析. 石油物探, 54(3): 350–358.
DOI:10.3969/j.issn.1000-1441.2015.03.015 |
|
熊晓军, 林凯, 贺振华.
2012. 基于等效弹性模量反演的横波速度预测方法. 石油地球物理勘探, 47(5): 723–727.
|
|
赵邦六, 石玉梅, 姚逢昌, 等.
2013. 多分量地震全波形弹性反演预测砂岩油藏剩余油分布. 石油学报, 34(2): 328–333.
DOI:10.7623/syxb201302015 |
|