张佳佳，印兴耀，张广智，谷一鹏，樊祥刚

（1.中国石油大学（华东）地球科学与技术学院，山东青岛 266580；2.深层油气地质与勘探教育部重点实验室，山东青岛 266580；3.海洋国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266071）

0 引言

地震数据中蕴含着丰富的地震弹性参数和储集层物性参数信息[1]，储集层物性参数预测主要分两步：第1步是利用常规地震反演方法由地震数据获取地震弹性参数，现有的地震反演技术已经非常成熟[2-5]，能够反演得到较为准确的弹性参数，例如纵横波速度、密度等；第 2步是利用岩石物理反演由地震弹性参数预测储集层参数，如孔隙度、饱和度、泥质含量等[6-8]。岩石物理反演首先需要构建一个理论的或经验的岩石物理模型，即建立地震弹性参数与储集层物性参数之间的关系[9-10]。其次就是岩石物理反演采用的反演算法，现阶段绝大多数采用了概率统计学反演方法等[11-12]。

岩石物理模型和岩石物理反演算法都直接影响储集层物性参数反演的精度。例如在许多岩石物理模型中，接触胶结模型[13]数学形式相对简单，被广泛应用于岩石物理反演中[14-15]。但是实际生产中弹性参数与物性参数之间的关系非常复杂，例如利用等效介质理论[16]来构建岩石物理模型。如果选择复杂岩石物理模型则需要采用概率统计学反演算法，如贝叶斯方法或蒙特卡洛算法等。国内外很多学者在这两方面开展了大量研究工作，如Bachrach等建立了物性参数与弹性参数之间的统计关系，采用贝叶斯反演方法实现了孔隙度和含水饱和度的联合反演[17]。Spikes等利用接触胶结模型建立了物性参数与弹性参数之间的关系，采用贝叶斯反演方法实现了孔隙度、泥质含量和含水饱和度的预测[18]。Bosch利用Wyllie时间平均方程[14]对岩石进行建模，采用马尔可夫链蒙特卡洛抽样方法和最小二乘迭代反演物性参数[19]。邓继新等提出了基于随机岩石物理模型的贝叶斯反演方法，实现对孔隙度与含水饱和度进行联合反演[20]。Grana和Rossa利用接触胶结模型[13]进行岩石物理建模，在贝叶斯理论框架下进行储集层参数后验概率分布估计，实现储集层参数反演[21]。印兴耀等建立了弹性阻抗与储集层物性参数之间的统计关系，利用贝叶斯理论估计物性参数的后验概率分布，实现储集层物性参数反演[11]。Grana提出了对胶结接触模型线性近似的方法，利用贝叶斯理论求解岩石物理反演问题，估计储集层物性参数[22]。林恬等利用包含物模型构建岩石物理模型，提出了约束最小二乘法与信赖域的储集层参数地震反演方法[12]。Li等利用等效介质理论[23]构建三维岩石物理模版，采用最小二乘原理寻找三维模版网格点的方法预测物性参数。杨培杰建立了砂泥岩储集层物性参数和地震弹性参数之间的定量关系，基于贝叶斯理论同步反演储集层孔隙度和含水饱和度[24]。上述岩石物理反演方法选用的岩石物理模型比较简单，例如统计经验关系和胶结接触模型，另一个问题就是很多反演算法利用的是概率统计学反演算法，需要进行全局寻优或者随机抽样，求解岩石物理反演问题的计算效率较低。因此，本文提出了一种岩石物理模型线性化方法，对岩石物理模型进行泰勒展开，得到岩石物理反演问题的一阶近似表达式。然后再利用阻尼最小二乘算法直接求解岩石物理反演问题，获得岩石物理反演问题的解析解，不需要全局寻优或者随机抽样，而是直接求逆运算，计算效率较高。本文利用 Gassmann方程[25]和临界孔隙度模型[26]对含泥质砂岩储集层进行岩石物理建模，建模过程比较复杂，构建的地震弹性参数和储集层物性参数之间的关系是非线性的。实际测井数据和地震数据验证表明，线性化岩石物理反演方法可以获得较为准确的物性参数，可以用来近似实际岩石物理模型。

1 方法原理

1.1 线性化岩石物理反演

岩石物理反演问题一般可写为：

（1）式中，d为弹性参数，m为待预测的物性参数，f为弹性参数与物性参数之间的关系式，即岩石物理模型。为简单起见，首先假设（1）式中所有变量均为单一参数的标量，然后再把它们推广到多参数，例如d为纵波速度，m为孔隙度，f为Wyllie时间平均方程[19]或者Raymer改进公式[27]。

实际生产中，岩石物理模型f并不像Wyllie时间平均方程或者Raymer改进公式这样简单，往往是非常复杂而且是非线性的。地球物理中经常利用泰勒展开来近似复杂的非线性函数，特别是在地震成像中[28]。为了求解岩石物理反演问题，本文也采用了基于泰勒展开的线性化岩石物理反演方法[22]。

利用泰勒展开对岩石物理反演问题进行线性化，这里使用一阶泰勒展开，那么（1）式可以改写为：

对（2）式中的多项式进行展开和合并，岩石物理反演问题可以改写为线性形式：

再把这种线性化岩石物理反演方法推广到多个参数，那么弹性参数d和物性参数m为包含多参数变量的矢量形式。例如，d为纵波速度vp、横波速度vs和密度ρ等弹性参数组成的向量，m为孔隙度φ、含水饱和度Sw和泥质含量C等物性参数组成的向量。同样岩石物理模型f也写成矢量形式，那么岩石物理反演问题可以写为：

对公式（4）进行一阶泰勒展开，变为

其中，Jm0为函数f在自变量m0处的Jacobian矩阵。再将其改写成线性化形式则变为：

1.2 岩石物理模型建模

本文选用 Gassmann方程[25]和临界孔隙度模型[26]对含泥质砂岩进行岩石物理建模，具体建模过程如下：

①岩石基质的体积模量Km、剪切模量μm和密度ρm可以利用Voigt-Reuss-Hill平均方法[29]计算：

②饱和流体的体积模量Kfl和密度ρfl计算公式为：

这里利用流体均匀混合公式[9]来计算饱和流体的体积模量。

③岩石骨架的体积模量Kdry和剪切模量μdry可以利用临界孔隙度模型[26]计算：

这里临界孔隙度值φc虽然是经验值，譬如砂岩取40%、灰岩取 60%等，但是临界孔隙度值包含了岩石的孔隙类型或者孔隙形状对纵横波速度的影响[30]。

④饱和岩石的体积模量Ksat和剪切模量μsat可以用Gassmann方程[25]计算：

⑤饱和岩石的纵波速度vP、横波速度vs和密度ρ计算公式为：

由于上述岩石物理建模过程即（7）—（18）式是非线性的，所以可以对其进行线性化处理。将弹性参数d（即纵波速度vP、横波速度vS和密度ρ）分别对物性参数m（孔隙度φ、泥质含量C和含水饱和度Sw）求导，计算其对应的Jacobian矩阵：

上式中，纵波速度vP对孔隙度φ的偏导数为：

横波速度vs对孔隙度φ的偏导数为：

密度ρ对孔隙度φ的偏导数为：

纵波速度vP对含水饱和度Sw的偏导数为：

横波速度vs对含水饱和度Sw的偏导数为：

密度ρ对含水饱和度Sw的偏导数为：

纵波速度vP、横波速度vs和密度ρ对泥质含量C的偏导数太复杂，故本文没有给出具体表达式。

再将 Jacobian矩阵代入到线性化岩石物理反演问题（6）式中。求解线性化之后的岩石物理反演问题可以有很多方法，例如最小二乘算法，阻尼最小二乘算法或者贝叶斯方法。这里选用阻尼最小二乘方法来求解，可得：

1.3 线性化岩石物理模型精度分析

为了检验线性化岩石物理反演的精度，分别利用实际岩石物理模型以及线性化岩石物理模型进行含泥质砂岩的岩石物理正演模拟分析。假设含泥质砂岩的矿物成分为石英和黏土，孔隙流体为水和气混合，相应的弹性模量和密度分别为：Kq= 3 6GPa ，μq= 3 6GPa ，ρq= 2 .65g/cm3，Kc= 2 1GPa ，μc=15GPa，ρc=2.45 g/cm3，Khc=0.0208GPa，ρhc= 0 .001g/cm3，Kw=2.25 GPa，ρw= 1 .03g/cm3。线性化岩石物理模型均在孔隙度、泥质含量和饱和度的均值处进行一阶泰勒展开。

图1分别展示了弹性参数（纵波速度、横波速度及密度）与孔隙度之间的变化关系，图中含泥质砂岩的泥质含量C以10%为间隔从0到50%变化，含水饱和度Sw为常数30%，孔隙度φ从0到30%均匀变化。由图1a和1b可以看出，线性化岩石物理模型与实际岩石物理模型在孔隙度为 15%左右处吻合最好，在孔隙度较低和较高时有一定的误差，这是因为在孔隙度均值15%处进行的一阶泰勒展开。由图1c可以看出，线性化岩石物理模型与实际岩石物理模型的密度完全一样，这是因为密度计算公式本身就是线性的，线性化岩石物理模型与实际岩石物理模型是完全相同的。

图1 实际岩石物理模型及线性化岩石物理模型计算的弹性参数与孔隙度之间的变化关系

图2分别展示了弹性参数（纵波速度、横波速度及密度）与泥质含量之间的变化关系。图中含泥质砂岩的孔隙度φ从以10%间隔从0到30%变化，含水饱和度Sw为常数50%，泥质含量C从0到50%均匀变化。由图2a和 2b可以看出，线性化岩石物理模型与实际岩石物理模型在 25%处吻合最好，在泥质含量较低和较高时有一定的误差，这是因为在泥质含量均值处（25%）进行的一阶泰勒展开。由图2c同样可以看出，线性化岩石物理模型与实际岩石物理模型的密度是完全一样的，这是因为密度计算公式本身就是线性的，线性化岩石物理模型与实际岩石物理模型是完全相同的。

图3分别展示了弹性参数（纵波速度、横波速度及密度）与含水饱和度之间的变化关系。图中含泥质砂岩的孔隙度φ以10%间隔从0到30%变化，泥质含量C为常数25%，含水饱和度Sw从0到100%均匀变化。由图3a和3b可以看到，利用均匀饱和方式计算饱和流体的体积模量，实际岩石物理模型和线性化岩石物理模型计算的纵波速度吻合很好，仅在含水饱和度Sw接近 100%左右误差较大，是符合近似要求的。由图3c同样可以看出，线性化岩石物理模型与实际岩石物理模型的密度是完全一样的，这是因为密度计算公式本身就是线性的，线性化岩石物理模型与实际岩石物理模型是完全相同的。

图2 实际岩石物理模型及线性化岩石物理模型计算的弹性参数与泥质含量之间的变化关系

图3 实际岩石物理模型及线性化岩石物理模型计算的弹性参数与含水饱和度之间的变化关系

2 实际应用对比

2.1 测井数据结果对比

首先将线性化岩石物理反演方法应用于实际测井数据。A井位于中国东部海上油田，储集层为含气砂岩，测井曲线包括孔隙度、泥质含量和含水饱和度等物性参数，以及纵波速度、横波速度和密度等弹性参数（见图4）。测井数据经过井震标定后由深度域转换成时间域，并按照2 ms进行重采样。主力产气层为时间深度为2 760～2 800 ms的一套厚砂岩储集层，其上覆2 718 ms左右有一层较薄的泥岩层，2 810～2 822 ms有一套较薄的砂岩储集层，其下覆为一套较厚泥岩层。

图4 A井实际测井数据

在进行线性化岩石物理反演之前，首先进行正演模拟分析。含泥质砂岩的矿物成分为石英和黏土，孔隙流体为水和气混合，相应的弹性模量和密度分别：Kq= 3 8GPa ，μq= 4 0GPa ，ρq= 2 .65g/cm3，Kc=15 GPa，μc= 7 GPa，ρc= 2 .55g/cm3，Khc= 0 .0208GPa ，ρhc= 0 .001g/cm3，Kw=2.25GPa ，ρw= 1 .03g/cm3。

利用精确岩石物理模型以及线性化岩石物理模型，由图4中的孔隙度、泥质含量和含水饱和度等物性参数来计算纵波速度、横波速度和密度，如图5所示。图中蓝色实线代表实际测井数据，黑色实线代表实际岩石物理模型正演计算结果，红色虚线代表线性化岩石物理模型正演计算结果。由图5可见，线性化岩石物理模型与实际岩石物理模型正演计算结果基本一致，与实际测井曲线有一定误差，但整体变化趋势相同。线性化岩石物理模型与实际弹性参数之间的平均误差，纵波速度约为4.0%，横波速度约为3.0%，密度约为1.2%。在2 820～2 860 ms泥岩段误差较大，主要是由于泥岩的弹性模量变化较大，很难获取精确值。

图5 实际岩石物理模型与线性化岩石物理模型正演计算的弹性参数以及实际测井弹性参数

其次，再对测井曲线进行线性化岩石物理反演，即利用图4d—4f中纵波速度、横波速度和密度来反演图4a—4c中的孔隙度、泥质含量和含水饱和度。使用的弹性模量和密度等参数与图5中一致，反演结果如图6所示。由图6可见，线性化岩石物理模型反演结果与实际测井曲线基本一致，虽然幅值没有完全吻合，但整体变化趋势是相同的。这里给出了线性化岩石物理模型反演结果与实际测井曲线之间的相关系数和平均误差，其中孔隙度的相关系数为64.8%，平均误差约为28%，含水饱和度的相关系数为85.4%，平均误差约为17%，泥质含量的相关系数为90.2%，平均误差约为10%。

图6 线性化岩石物理模型反演与实际测井物性参数

2.2 地震数据结果对比

在进行测井曲线的岩石物理反演之后，再对地震数据进行线性化岩石物理反演。选择一条过A井的二维地震剖面，由叠前同时反演[31]分别获得了纵波速度、横波速度和密度等弹性参数的二维剖面（见图7）。由二维弹性参数剖面上可以看到，无论是纵波速度、横波速度还是密度在2 800 ms左右都有一套较厚的低值区域，这是主要的目的层段。再利用图7中的弹性参数进行物性参数反演（见图8），图8a—8c分别展示了线性化岩石物理反演得到的孔隙度、泥质含量和含气饱和度，图8d—8f分别展示了利用弹性参数与物性参数之间的经验关系拟合得到的孔隙度、泥质含量和含气饱和度。由反演结果可见，两种反演结果得到的主要目的层段物性参数与A井上的物性参数吻合较好，并且能够获取主要目的层段的物性参数空间展布情况，特别是线性化岩石物理反演含气饱和度效果要优于经验关系拟合结果。

图7 地震弹性参数反演剖面

3 讨论

图8 线性化岩石物理反演的物性参数

本文提出了一种线性化岩石物理反演方法，由于岩石物理模型一般并不是线性的，所以本文方法较适用于非线性化不是特别强的岩石物理模型，尤其适用于孔隙结构较简单的储集层，如碎屑岩储集层；对于孔隙结构复杂的储集层，如碳酸盐岩储集层并不太适用。另外值得注意的是，不同的岩石类型应该选用不同的岩石物理模型。本文是针对含泥质砂岩储集层进行建模，选择了临界孔隙度模型和Gassmann方程是合理的，如果针对其他类型储集层则需要选用其他的岩石物理模型。文中给出的岩石物理正演模型表明，实际岩石物理模型与线性化岩石物理模型是可以近似的。利用线性化岩石物理模型的优势就是可以获取岩石物理反演问题的解析解。本文中选择阻尼最小二乘方法来求解线性反演问题，求解速度较快，但对岩石物理建模中使用的弹性模量和密度等参数较为敏感。另外，不同于其他岩石物理反演方法，线性化岩石物理反演方法获得的泥质含量和饱和度反演精度要优于孔隙度的反演精度。主要是由于在 Jacobian矩阵的逆矩阵中与泥质含量以及饱和度有关的系数项较大，而与孔隙度有关的系数项较小，所以导致泥质含量与饱和度的反演效果要优于孔隙度反演效果。当然，本文中的弹性参数是由AVO叠前地震反演获得，下一步还可以直接利用地震数据进行岩石物理反演。

4 结论

提出了一种线性化岩石物理反演方法，即通过对实际岩石物理模型进行泰勒级数展开，得到一阶近似的线性化岩石物理模型，再利用阻尼最小二乘方法求解线性化后的岩石物理反演问题，获得精确的解析解。利用了 Gassmann方程和临界孔隙度模型对含泥质砂岩进行实际岩石物理建模，并且给出了其线性化岩石物理模型表达式。本方法适用于线性或轻微非线性的岩石物理模型，对于高度非线性岩石物理模型可能不适应。实际测井数据和地震数据应用表明，线性化岩石物理反演方法能够反演得到比较好的物性参数结果。

符号注释：

C——泥质含量，%；C0——泥质含量均值，%；d——弹性参数；d——弹性参数矢量形式；f——岩石物理模型；f′——岩石物理模型f的一阶导数；f——岩石物理模型矢量形式；f（m0）——岩石物理模型在m0（即φ0、C0和Sw0）处的值；I——单位矩阵；Jm0——Jacobian矩阵在m0（即φ0、C0和Sw0）处的值；Kc——泥质矿物的体积模量，GPa；Kdry——岩石骨架的体积模量，GPa；Kfl——饱和流体的体积模量，GPa；Km——岩石基质的体积模量，GPa；Khc——油或气等碳氢化合物的体积模量，GPa；Kq——砂质矿物的体积模量，GPa；Ksat——饱和岩石的体积模量，GPa；Kw——水的体积模量，GPa；m——待预测的物性参数；m0——物性参数m某一具体值，一般选择m0为待预测模型参数的均值；m——待预测的物性参数矢量形式；m0——物性参数矢量形式m某一具体值，即φ0、C0和Sw0；Sw——含水饱和度，%；Sw0——饱和度均值，%；vP——饱和岩石的纵波速度，km/s；vs——饱和岩石的横波速度，km/s；ρ——饱和岩石的密度，g/cm3；ρc——矿物的密度，g/cm3；ρfl——饱和流体的密度，g/cm3；ρhc——油或气等碳氢化合物的密度，g/cm3；ρm——岩石基质的密度，g/cm3；ρq——砂质矿物的密度，g/cm3；ρw——水的密度，g/cm3；φ——孔隙度，%；φ0——孔隙度均值，%；φc——岩石的临界孔隙度，%；μc——泥质矿物的剪切模量，GPa；μdry——岩石骨架的剪切模量，GPa；μm——岩石基质的剪切模量，GPa；μq——砂质矿物的剪切模量，GPa；μsat——饱和岩石的剪切模量，GPa；ε——阻尼因子。