0 引言

孔隙度是指示储层物性的关键参数之一，可用于精细油藏研究。岩石物理反演是定量预测孔隙度的主要手段，其思路是以岩石物理模型为基础，利用纵、横波速度等弹性参数反演孔隙度、泥质含量等物性参数^[1-5]。现今的岩石物理反演方法主要以非线性岩石物理模型为基础，构建反演目标函数，并采用非线性算法求解，如蒙特卡洛法、模拟退火法等。但此类算法多解性强，且计算量巨大^[6-8]。为此，近年来人们研究了线性反演方法。Grana^[9]提出了对岩石物理模型线性近似的方法，并以此为基础，实现了三参数(储层孔隙度、泥质含量与饱和度)线性化反演，大幅降低了计算量与反演多解性，为在实际生产中推广应用奠定了基础。凌东明等^[10]针对远离均值点模型近似结果偏差较大的问题，采用岩性约束相分段近似Xu-White模型，提高了储层物性预测精度。张佳佳等^[11]以前人研究成果为基础，利用阻尼最小二乘算法直接求解线性化岩石物理反演问题，提高了三参数反演的计算效率。

现有的岩石物理线性化反演方法均为三参数反演，是将叠前弹性反演得到的纵、横波速度及密度作为已知量，求取孔隙度、泥质含量与饱和度三个未知量^[9-11]。与纵、横波速度相比，密度对反射系数贡献较小，反演密度需要更大的叠前地震道集角度范围^[12-13]。在储层深埋条件下，由于叠前地震道集反射角较小，导致密度反演结果不可靠；在忽略密度项的情况下，岩石物理三参数反演方程为欠定形式，无法得到唯一解，从而限制了现有岩石物理线性化反演方法对深埋储层的预测。

在前人研究基础之上，本文提出一种基于迭代算法的储层孔隙度反演方法。该方法是在岩石物理模型线性近似的基础上，推导纵、横波速度与孔隙度、泥质含量的线性关系式，然后基于贝叶斯理论构建孔隙度迭代反演目标函数，再利用二分法迭代求解。该方法不依赖于密度项，适用于储层深埋条件的孔隙度反演。

1 孔隙度迭代反演方法 1.1 线性反演的局限性分析

岩石物理模型是岩石物理反演的基础，建立了储层物性参数与弹性参数之间的关系，可表示为

$ \boldsymbol{d}=\boldsymbol{f}(\boldsymbol{m})+\boldsymbol{e} $

(1)

其中

$ \left\{\begin{array}{l} \boldsymbol{d}=\left(\begin{array}{c} V_{\mathrm{P}} \\ V_{\mathrm{S}} \\ \rho \end{array}\right) \\ \boldsymbol{m}=\left(\begin{array}{c} \varphi \\ \phi \\ S_{\mathrm{W}} \end{array}\right) \end{array}\right. $

(2)

式中：d为由岩石纵波速度V_P、横波速度V_S和密度ρ组成的向量；f(m)为岩石物理模型函数，m为由孔隙度φ、泥质含量ϕ和含水饱和度S_W组成的向量；e为误差向量。岩石物理反演以式(1)为基础求解模型参数m。在一般情况下，f(m)为非线性函数，需采用复杂的迭代算法求解，而此类算法多解性强，且计算量巨大，难以用于生产。为此，采用泰勒级数展开f(m)，保留一阶项得到线性近似式

$ \boldsymbol{d} \cong \boldsymbol{f}\left(\boldsymbol{m}_{0}\right)+\boldsymbol{J}_{\boldsymbol{m}_{0}} \cdot\left(\boldsymbol{m}-\boldsymbol{m}_{0}\right)+\boldsymbol{\delta} $

(3)

其中

$ \boldsymbol{J}=\left(\begin{array}{lll} \frac{\partial V_{\mathrm{P}}}{\partial \varphi} & \frac{\partial V_{\mathrm{P}}}{\partial \phi} & \frac{\partial V_{\mathrm{P}}}{\partial S_{\mathrm{W}}} \\ \frac{\partial V_{\mathrm{S}}}{\partial \varphi} & \frac{\partial V_{\mathrm{S}}}{\partial \phi} & \frac{\partial V_{\mathrm{S}}}{\partial S_{\mathrm{W}}} \\ \frac{\partial \rho}{\partial \varphi} & \frac{\partial \rho}{\partial \phi} & \frac{\partial \rho}{\partial S_{\mathrm{W}}} \end{array}\right) $

(4)

式中：J_m0为f(m)在已知点m₀处的雅可比矩阵，在反演过程中可将m₀设置为常值模型；δ为e的一阶近似。在实际应用中需根据目的层特征选择合适的岩石物理模型建模，得到适合的模型参数，进而计算f(m₀)与J_m0。

令b=f(m₀)－J_m0·m₀，对式(3)重新整理，得到储层物性参数与弹性参数之间的线性关系式

$ \left(\begin{array}{c} V_{\mathrm{P}}-b^{(\mathrm{P})} \\ V_{\mathrm{S}}-b^{(\mathrm{S})} \\ \rho-b^{(\rho)} \end{array}\right)=\left(\begin{array}{lll} \frac{\partial V_{\mathrm{P}}}{\partial \varphi} & \frac{\partial V_{\mathrm{P}}}{\partial \phi} & \frac{\partial V_{\mathrm{P}}}{\partial S_{\mathrm{W}}} \\ \frac{\partial V_{\mathrm{S}}}{\partial \varphi} & \frac{\partial V_{\mathrm{S}}}{\partial \phi} & \frac{\partial V_{\mathrm{S}}}{\partial S_{\mathrm{W}}} \\ \frac{\partial \rho}{\partial \varphi} & \frac{\partial \rho}{\partial \phi} & \frac{\partial \rho}{\partial S_{\mathrm{W}}} \end{array}\right) \cdot\left(\begin{array}{c} \varphi \\ \phi \\ S_{\mathrm{W}} \end{array}\right)+\left(\begin{array}{l} \delta_{\mathrm{P}} \\ \delta_{\mathrm{S}} \\ \delta_{\rho} \end{array}\right) $

(5)

式中：b^(P)、b^(S)和b^(ρ)分别为泰勒级数展开后纵、横波速度和密度对应的系数项；δ_P、δ_S、δ_ρ分别为泰勒级数展开后纵、横波速度和密度对应的误差项。现有的岩石物理线性反演方法均以式(5)为基础，构建目标函数求解^[9-11]。式(5)为正定方程组，将叠前弹性反演得到的纵、横波速度及密度作为已知量，可以得到唯一的孔隙度、泥质含量与含水饱和度参数。但在储层深埋条件下，叠前道集反射角较小，无法达到反演密度需要的最小角度范围。若忽略密度项，则式(5)改为

$ \left(\begin{array}{l} V_{\mathrm{P}}-b^{(\mathrm{P})} \\ V_{\mathrm{S}}-b^{(\mathrm{S})} \end{array}\right)=\left(\begin{array}{lll} \frac{\partial V_{\mathrm{P}}}{\partial \varphi} & \frac{\partial V_{\mathrm{P}}}{\partial \phi} & \frac{\partial V_{\mathrm{P}}}{\partial S_{\mathrm{W}}} \\ \frac{\partial V_{\mathrm{S}}}{\partial \varphi} & \frac{\partial V_{\mathrm{S}}}{\partial \phi} & \frac{\partial V_{\mathrm{S}}}{\partial S_{\mathrm{W}}} \end{array}\right) \cdot\left(\begin{array}{c} \varphi \\ \phi \\ S_{\mathrm{W}} \end{array}\right)+\left(\begin{array}{c} \delta_{\mathrm{P}} \\ \delta_{\mathrm{S}} \end{array}\right) $

(6)

式(6)包含两个方程、三个未知量，为欠定方程组，以此为基础进行岩石物理反演无法得到唯一解。因此，在储层深埋条件下无法得到可靠密度反演结果，导致常规岩石物理线性反演方法不适用。

1.2 孔隙度反演方程推导

式(6)为储层深埋条件下的纵、横波速度与孔隙度、泥质含量及饱和度的线性关系，雅可比矩阵系数项反映了不同物性参数对速度的贡献程度。在油、水两相饱和岩石中，孔隙度变化对速度影响最大，泥质含量次之，饱和度最小^[14]，即纵、横波速度对饱和度变化不敏感。为此，将饱和度作为扰动项，则式(6)变为

$ \begin{aligned} &\left(\begin{array}{l} V_{\mathrm{P}}-b^{(\mathrm{P})}-\frac{\partial V_{\mathrm{P}}}{\partial S_{\mathrm{W}}} \cdot S_{\mathrm{W}} \\ V_{\mathrm{S}}-b^{(\mathrm{S})}-\frac{\partial V_{\mathrm{S}}}{\partial S_{\mathrm{W}}} \cdot S_{\mathrm{W}} \end{array}\right)= \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left(\begin{array}{cc} \frac{\partial V_{\mathrm{P}}}{\partial \varphi} & \frac{\partial V_{\mathrm{P}}}{\partial \phi} \\ \frac{\partial V_{\mathrm{S}}}{\partial \varphi} & \frac{\partial V_{\mathrm{S}}}{\partial \phi} \end{array}\right) \cdot\left(\begin{array}{l} \varphi \\ \phi \end{array}\right)+\left(\begin{array}{l} \delta_{\mathrm{P}} \\ \delta_{\mathrm{S}} \end{array}\right) \end{aligned} $

(7)

由此可见，当给定初始饱和度时，式(7)为正定方程组，利用纵、横波速度即可反演唯一的孔隙度与泥质含量。因此，可以式(7)为基础，基于贝叶斯理论构建反演目标函数。

为公式推导方便起见，令

$ \boldsymbol{d}^{\prime}=\left(\begin{array}{l} V_{\mathrm{P}}-b^{(\mathrm{P})}-\frac{\partial V_{\mathrm{P}}}{\partial S_{\mathrm{W}}} \cdot S_{\mathrm{W}} \\ V_{\mathrm{S}}-b^{(\mathrm{S})}-\frac{\partial V_{\mathrm{S}}}{\partial S_{\mathrm{W}}} \cdot S_{\mathrm{W}} \end{array}\right) $

$ \boldsymbol{G}=\left(\begin{array}{ll} \frac{\partial V_{\mathrm{P}}}{\partial \varphi} & \frac{\partial V_{\mathrm{P}}}{\partial \phi} \\ \frac{\partial V_{\mathrm{S}}}{\partial \varphi} & \frac{\partial V_{\mathrm{S}}}{\partial \phi} \end{array}\right) $

$ \boldsymbol{m}^{\prime}=\left(\begin{array}{l} \varphi \\ \phi \end{array}\right) $

则式(7)简化为

$ \boldsymbol{d}^{\prime}=\boldsymbol{G} \cdot \boldsymbol{m}^{\prime}+\boldsymbol{\delta} $

(8)

依据贝叶斯理论，后验概率密度分布函数正比于似然函数与先验分布的乘积^[15]。假设误差项δ与先验模型参数m′均服从高斯分布，基于式(8)可构建后验概率密度分布函数近似式

$ \begin{aligned} &P\left(\boldsymbol{m}^{\prime} \mid \boldsymbol{d}^{\prime}\right) \cong \lambda \times \\ &\exp \left[-\frac{1}{2}\left(\boldsymbol{d}^{\prime}-\boldsymbol{G} \cdot \boldsymbol{m}^{\prime}\right)^{\mathrm{T}} \cdot \boldsymbol{C}_{\boldsymbol{\delta}}^{-1} \cdot\left(\boldsymbol{d}^{\prime}-\boldsymbol{G} \cdot \boldsymbol{m}^{\prime}\right)\right] \times \\ &\exp \left[-\frac{1}{2}\left(\boldsymbol{m}^{\prime}-\boldsymbol{m}^{\prime}{}_{0}\right)^{\mathrm{T}} \cdot \boldsymbol{C}_{\boldsymbol{m}^{\prime}}^{-1} \cdot\left(\boldsymbol{m}^{\prime}-\boldsymbol{m}^{\prime}{}_{0}\right)\right] \end{aligned} $

(9)

式中：λ为常数；C_δ、C_m′分别为δ、m′的协方差矩阵，与常值模型m′₀共同作为先验信息(均由统计测井数据得到)约束求解过程。求解式(9)概率最大值等效于求解目标函数

$ \begin{aligned} L=&\left(\boldsymbol{d}^{\prime}-\boldsymbol{G} \cdot \boldsymbol{m}^{\prime}\right)^{\mathrm{T}} \cdot \boldsymbol{C}_{\boldsymbol{\delta}}^{-1} \cdot\left(\boldsymbol{d}^{\prime}-\boldsymbol{G} \cdot \boldsymbol{m}^{\prime}\right)+\\ &\left(\boldsymbol{m}^{\prime}-\boldsymbol{m}^{\prime}{}_{0}\right)^{\mathrm{T}} \cdot \boldsymbol{C}_{\boldsymbol{m}^{\prime}}^{-1} \cdot\left(\boldsymbol{m}^{\prime}-\boldsymbol{m}^{\prime}{}_{0}\right) \end{aligned} $

(10)

的最小值。为此，令∂L/∂m′=0，得到反演方程

$ \begin{aligned} \boldsymbol{m}^{\prime}=&\left(\boldsymbol{G}^{\mathrm{T}} \cdot \boldsymbol{C}_{\boldsymbol{\delta}}^{-1} \cdot \boldsymbol{G}+\boldsymbol{C}_{\boldsymbol{m}^{\prime}}^{-1}\right)^{-1} \cdot \\ &\left(\boldsymbol{G}^{\mathrm{T}} \cdot \boldsymbol{C}_{\boldsymbol{\delta}}^{-1} \cdot \boldsymbol{d}^{\prime}+\boldsymbol{C}_{\boldsymbol{m}^{\prime}}^{-1} \cdot \boldsymbol{m}^{\prime}{}_{0}\right) \\ =& \boldsymbol{m}^{\prime}{}_{0}+\boldsymbol{C}_{\boldsymbol{m}^{\prime}} \cdot \boldsymbol{G}^{\mathrm{T}} \cdot\left(\boldsymbol{G} \cdot \boldsymbol{C}_{\boldsymbol{m}^{\prime}} \cdot \boldsymbol{G}^{\mathrm{T}}+\boldsymbol{C}_{\boldsymbol{\delta}}\right)^{-1} \cdot \\ &\left(\boldsymbol{d}^{\prime}-\boldsymbol{G} \cdot \boldsymbol{m}^{\prime}{}_{0}\right) \end{aligned} $

(11)

1.3 饱和度敏感性分析

利用式(11)求解孔隙度与泥质含量时，需给定饱和度项，即反演精度与饱和度准确程度相关。为验证反演结果对饱和度的敏感性，应用理论模型进行分析。

已知油、水两相流体饱和岩石孔隙度、泥质含量与含水饱和度曲线(图 1a)，通过岩石物理正演得到纵波速度与横波速度(图 1b)。在正演过程中，选用适用砂、泥岩储层的Xu-White模型^[16]作为岩石物理模型。为进行极限测试，给定不同常值饱和度模型，利用纵、横波速度反演孔隙度与泥质含量，进而分析饱和度对反演精度的影响。

图 2为不同饱和度模型的反演与正演结果。由图可见：当含水饱和度为1.00时，反演、正演结果均与模型曲线存在较大偏差(图 2a)；当含水饱和度为0时(与饱和度均值较接近)，反演、正演效果得到明显改善(图 2b)；当含水饱和度为0.50时(进一步逼近饱和度均值)，反演、正演效果进一步提升(图 2c)；当含水饱和度为模型均值(0.29)时，反演、正演结果与模型具有较好的一致性(图 2d)。

图 3为图 2对应的反演结果绝对误差。由图可见，设定饱和度越接近模型饱和度均值，绝对误差越小，具体表现为：当饱和度为1.00时，孔隙度平均绝对误差为0.034，泥质含量平均绝对误差为0.106(图 3a)；当饱和度为0时，孔隙度平均绝对误差为0.013，泥质含量平均绝对误差为0.042(图 3b)；当饱和度为0.50时，孔隙度平均绝对误差为0.01，泥质含量平均绝对误差为0.033(图 3c)；当饱和度为模型饱和度均值(0.29)时，孔隙度平均绝对误差为0.003，泥质含量平均绝对误差为0.011(图 3d)。因此，对于油、水两相的流体饱和岩石模型，反演结果对饱和度模型的精度要求低，给定合理的常值饱和度模型，即可保证相对反演结果的可靠性，说明反演结果对饱和度不敏感。

1.4 孔隙度迭代反演流程

由上述分析可知，反演精度取决于对真实饱和度均值的搜索，常值模型的饱和度越接近真实饱和度均值，反演精度越高。为此，本文采用二分法迭代求解式(11)，思路如下：

(1) 分别设置初始模型的含水饱和度为0、1.00进行第一次反演，得到孔隙度与泥质含量，并基于岩石物理模型，结合反演结果正演纵、横波速度，并计算正演结果与实测结果的平均相关系数γ₀、γ₁；

(2) 按二分法，对初始饱和度模型求平均，得到第二次反演所需的饱和度常值模型，即为0.50，并计算正演相关系数γ₂，对比γ₀，γ₁，γ₂，保留两项最大相关系数及对应饱和度进入下一轮反演；

(3) 如果反演结果的精度满足要求，则终止迭代，输出最终反演结果，否则，求取平均孔隙度继续迭代，直到反演结果满足要求或达到最大迭代次数为止。

上述迭代反演方法可忽略密度项影响，仅利用纵、横波速度即可得到相对可靠的孔隙度(图 4)。

2 模型试算

为验证孔隙度迭代反演方法的可行性与反演结果的有效性，利用一维正演模型(图 1)进行测试。以正演过程采用的岩石物理模型为基础，利用纵、横波速度迭代反演孔隙度与泥质含量。

图 5为一维模型孔隙度迭代反演结果，图 6为图 5对应的反演结果绝对误差。由图可见：随着迭代次数增加，饱和度均值趋于模型饱和度均值0.29，反演结果逐渐收敛(图 5)；孔隙度、泥质含量反演结果绝对误差接近0(图 6)，并且收敛速度较快，2~4次迭代即可得到最优解。上述结果说明所提迭代反演方法可行、有效。

3 实际应用

M区位于塔里木盆地北部凹陷，研究目标为底水块状油藏，埋深大于5000m，目的层为砂岩储层，厚度约为100m，以滨岸相沉积为主。储层物性相对均匀，且地震资料品质较好，适合孔隙度迭代反演测试。该区W1井测井序列相对完整，包含纵、横波速度及密度数据，通过计算可得到孔隙度、泥质含量及饱和度。以W1井数据为基础进行岩石物理建模，选用适合砂岩储层的Xu-White模型^[16]。

图 7为M区岩石物理模型正演结果，可见实测曲线与正演结果较吻合，证明所建岩石物理模型有效。通过岩石物理建模，确定M区砂、泥岩组分孔隙横纵比分别为0.190、0.008，岩石物理基本建模参数如表 1所示。

基于岩石物理建模结果，结合叠前同步反演得到的纵、横波速度等弹性参数，并将饱和度分别为1和0的初始含水饱和度模型作为输入，进行迭代反演，求取孔隙度与泥质含量(图 8)。

图 8为M区孔隙度与泥质含量迭代反演剖面。由反演结果可见一条高泥质含量条带(图 8b)，且对应位置的反演孔隙度出现明显低值区(图 8a)，稳定泥岩段上覆于巨厚块状砂体之上，与钻井成果及地质认识相吻合。因此，反演结果较好地表征了目的层段储盖组合。为了进一步验证反演效果，提取井点位置的反演曲线进行分析。

图 9为W1、W2井井位处目的层段反演结果。由图可见，孔隙度、泥质含量反演结果与实测曲线的相对变化规律基本一致，但二者的绝对值在局部存在差异，其原因是迭代得到的常值模型饱和度与实际饱和度的差异所致。因此，反演结果对目的层段储层物性具有一定指示作用，对于纯油层或纯水层，反演精度高，但当流体性质存在纵向分异时，反演精度有所下降。

4 结论

本文以岩石物理线性近似模型为基础，推导了纵、横波速度与孔隙度、泥质含量的线性关系式，并基于贝叶斯理论构建反演目标函数，利用二分法求解，迭代反演了储层孔隙度、泥质含量。该方法不依赖于密度项，适用于储层深埋条件下的孔隙度预测，拓宽了岩石物理线性反演方法的适用性。理论模型与实际数据测试结果表明，所提方法具有一定可行性，且应用效果较明显。

针对油、水两相流体饱和岩石，本文方法反演结果对饱和度模型不敏感，应用迭代反演方法可以得到相对可靠的孔隙度预测结果。对于纯油层或纯水层，反演精度高，但当流体性质存在纵向分异时，反演精度有所下降。对于气层，反演结果对饱和度模型的精度要求较高，可能会降低孔隙度迭代反演方法的适用性。关于含气储层的孔隙度反演方法有待进一步研究和完善。