CN103760600A - 一种含气饱和度反演方法 - Google Patents

Abstract

本发明涉及地质勘探技术领域，具体的讲是一种含气饱和度反演方法。本发明方法通过加权统计的概念以及解析化的目标反演函数表达式，利用测井资料统计分析与统计岩石物理建模相结合的思路，推导了较传统概率统计类方法效率更高、不受弹性参数精度及含气敏感性差异影响的目标反演函数及高效求解策略，方法具有较高的实用性，能有效地运用于三维实际资料，识别出高含气饱和度气藏的位置。

Description

一种含气饱和度反演方法

技术领域

本发明涉及地质勘探技术领域，具体的讲是一种含气饱和度反演方法。

背景技术

近年来我国天然气勘探的重点逐渐由构造气藏向岩性气藏勘探转移。不同于构造气藏，这些新型储层受构造和储层非均质性的影响，成藏条件复杂，识别难度大，投资风险高。因此，开发一套能够有效识别有利储层，分析岩性、流体及储层物性条件，尤其是含气饱和度的地球物理解决方案显得尤为迫切。含气饱和度的空间分布，直接影响油气储层的品质，是地质、地球物理工作者进行储层评价、估算油气储量、确定开发井位的重要依据。

天然气储层在地震资料上常表现为强烈的反射特征，如“亮点”。这种地震响应特征可能是极具商业开采价值的高含气饱和度储层引起，也可能是没有商业开采价值的低含气饱和度储层引起。事实上，不同含气饱和度储层的地震响应特征十分相似，如何利用地震资料进行储层含气饱和度预测，一直以来是地球物理勘探工作者的重要研究内容之一。

对于含气饱和度的反演最为常见的做法是先利用地震振幅信息，通过各种叠前地震反演方法得到弹性参数，然后在弹性参数的基础上进一步反演得到含气饱和度。目前，由地震振幅信息到弹性参数的反演技术在实际应用中已经成熟。难点在于弹性参数到含气饱和度的反演，应用比较广泛的方法主要基于多元统计类技术（Doyen,1988;Fournier,1989）以及确定性岩石物理建模类（Blangy1992;Marion,1997）技术。多元统计类技术主要利用数理统计的方式得到弹性参数与含气饱和度间的统计关系，进行弹性参数到含气饱和度的转化。然而这种统计关系基于纯粹的数学关系式，物理意义不明确且这种统计关系的准确性依赖于训练样本的数量，这就造成物性参数的反演结果具有极大的不确定性。确定性岩石物理建模类技术主要利用流体替代技术结合各种等效介质理论建立起弹性参数与含气饱和度间的确定性转化关系。事实上，对于复杂的储层环境而言，确定性岩石模型往往过于简单，误差较大。Mukerji等2001年首次将统计岩石物理模型介绍到储层预测领域用来评价预测结果的不确定性。Bachrach等2006年基于统计岩石物理模型，提出了含气饱和度与孔隙度的联合概率统计反演方法，认为该方法具有较高的反演精度。

发明人在实现本发明的过程中发现，常规含气饱和度方法在利用弹性参数反演含气饱和度时，一方面没有考虑到输入的不同弹性参数变量间存在有精度差异。这种精度差异主要是由叠前地震反演造成的，一般来说，由叠前地震反演得到的纵、横波速度及密度的精度依次降低。另一方面，不同的弹性参数对含气饱和度变化的敏感度不同，一般来说密度、纵波速度、横波速度对含气饱和度变化的敏感度依次降低。如果在含气饱和度反演过程中同等对待不同精度、敏感度的弹性参数变量，势必对反演结果造成影响。另外，常规概率统计类方法由于无解析化的目标函数表达式且在反演过程中涉及到大量的样本统计过程，导致反演十分耗时，使得这类方法较难应用于实际三维资料。

发明内容

本发明的目的在于，克服现有技术的不足，提供一种含气饱和度反演方法，以克服不同精度、不同含气敏感性的弹性参数对含气饱和度反演结果的影响且具有较高的反演效率。

本发明实施例提供了一种含气饱和度反演方法，包括，

根据含气饱和度参数R_i先验分布特征，得到含气饱和度先验概率密度函数

$P\left(R_i\right)=\underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n);$

上式中，表示高斯分布概率密度函数，和分别为与饱和度参数Ri有关的混合高斯分布部件数、权值、均值与方差；

设各弹性参数变量m^k间相互条件独立，则似然函数P(m_j|R_i)为：

$P(m_j^1,m_j^2,...m_j^{\mathrm{Nm}}|R_i)=\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Π}}}P\left(m_j^k\right|R_i);$

根据训练样本统计出所述似然函数P(m^k|R_i)的分布信息，得到所述似然函数P(m^k|R_i)的解析表达式：

$P\left(m^k\right|R_i)=\underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_i}^{\mathrm{\π}}),$

其中，N(·)表示高斯分布概率密度函数，N_c,α_π,

和

分别为与m^k、R_i有关的混合高斯分布部件数、权系数、均值与方差；

根据所述含气饱和度参数先验分布特征和所述似然函数P(m^k|R_i)的解析表达式以及引入权重系数W_k取对数运算得到反演目标函数的解析化表达式：

$\tilde{R}=\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k{R}_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k{R}_{i1}}^{\mathrm{\π}});$

上式中，ln表示取自然对数；

输入弹性参数数据，取所述反演目标函数最大后验概率密度解作为反演结果

${\tilde{R}}_{\mathrm{final}}=\arg \underset{R_i\∈R}{\mathrm{Max}}\{\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_i}^{\mathrm{\π}})\}$

上式中，

表示返回极大值对应的Ri。

通过上述实施例的方法，通过加入了权重系数的手段使得反演原理更加严谨、反演过程更加合理；推导出了目标函数的解析化表达式及高效的统计策略，使得大规模三维实际资料的含气饱和度定量反演成为现实；通过实际资料含气饱和度反演，精准地识别出高饱和度气藏的有利位置。

附图说明

结合以下附图阅读对实施例的详细描述，本发明的上述特征和优点，以及额外的特征和优点，将会更加清楚。

图1所示为本发明实施例一种含气饱和度反演方法的流程图；

图2所示为本发明实施例一种含气饱和度反演方法的具体流程图；

图3为现有技术的含气饱和度先验分布特征统计示意图；

图4为现有技术与本发明实施例中训练样本总数量随弹性参数取值个数的变化比较图；

图5为发明实施例的似然函数解析化表达式的求取示意图；

图6a为现有技术对纵、横波阻抗及密赋予相同权重时的反演结果；

图6b为本发明实施例对纵、横波阻抗及密度赋予不同权重时的反演结果；

图7为发明实施例针对某研究区过A井剖面含气饱和度反演结果图；

图8为发明实施例针对某研究区A井的井旁道反演结果与A井测井结果对比图；

图9为发明实施例针对某研究区过验证井B井剖面含气饱和度反演结果图；

图10为发明实施例针对某研究区验证井B井的井旁道反演结果与B井测井结果对比图；

图11为发明实施例针对某三维研究区含气饱和度反演结果沿层切片图。

具体实施方式

下面的描述可以使任何本领域技术人员利用本发明。具体实施例和应用中所提供的描述信息仅为示例。这里所描述的实施例的各种延伸和组合对于本领域的技术人员是显而易见的，在不脱离本发明的实质和范围的情况下，本发明定义的一般原则可以应用到其他实施例和应用中。因此，本发明不只限于所示的实施例，本发明涵盖与本文所示原理和特征相一致的最大范围。

如图1所示为本发明实施例一种含气饱和度反演方法的流程图。

包括步骤101，根据含气饱和度参数R_i先验分布特征，例如可以利用现有期望最大估计技术（EM）对含气饱和度先验分布特征进行高斯混合分布（GMM）函数逼近，得到含气饱和度先验概率密度函数 $P\left(R_i\right)=\underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n).$

步骤102，设各弹性参数变量m^k间相互条件独立，则似然函数P(m_j|R_i)为：

$P(m_j^1,m_j^2,...m_j^{\mathrm{Nm}}|R_i)=\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Π}}}P\left(m_j^k\right|R_i).$

步骤103，根据训练样本统计出所述似然函数P(m^k|R_i)的分布信息，例如可以利用现有期望最大估计技术（EM）对似然函数分布信息进行高斯混合分布（GMM）函数逼近，得到所述似然函数P(m^k|R_i)的解析表达式：

$P\left(m^k\right|R_i)=\underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_i}^{\mathrm{\π}}),$

其中，N(·)表示高斯分布概率密度函数，N_c,α_π,

和

分别为与m^k、R_i有关的混合高斯分布部件数、权系数、均值与方差。

步骤104，根据所述含气饱和度参数先验概率密度函数和所述似然函数P(m^k|R_i)的解析表达式以及引入权重系数W_k取对数运算得到反演目标函数的解析化表达式：

$\tilde{R}=\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_{i1}}^{\mathrm{\π}})$

步骤105，输入弹性参数数据，取所述反演目标函数最大后验概率密度解作为反演结果

${\tilde{R}}_{\mathrm{final}}=\arg \underset{R_i\∈R}{\mathrm{Max}}\{\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k{R}_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k{R}_i}^{\mathrm{\π}})\}$

作为本发明的一个实施例，在步骤101中，综合测井资料、地质资料、叠前地震反演资料等信息，采用直方图统计的方式得到所述含气饱和度参数先验分布特征。利用现有期望最大估计技术（EM）对含气饱和度参数先验分布特征进行高斯混合分布函数（GMM）逼近，使得含气饱和度度先验分布特征能用具体函数表达式来表示。在本发明中还可以使用其它现有技术手段，例如可以根据先验信息直方图统计特征，选用均匀分布函数、指数分布、柯西分布函数、高斯分布函数等已知表达式的分布函数对含气饱和度参数先验分布特征进行逼近，使得含气饱和度度先验分布特征能用具体函数表达式来表示。本发明选择高斯混合分布函数（GMM）的原因在于，如果期望最大估计技术(EM)得到的参数足够准确，高斯混合分布函数能够逼近任何形态的直方图统计特征。期望最大估计技术(EM)具体参见：李昌利,沈玉利，2008，期望最大算法及其应用：计算机工程与应用,44(29):61-64.

作为本发明的一个实施例，在所述步骤103中，根据训练样本统计出所述似然函数P(m^k|R_i)的分布信息进一步包括：

基于等效介质理论，根据储层地质特点选择确定性岩石物理模型m=f(R)；

在确定性岩石物理模型的基础上，根据模型与实际测井数据间的差异程度，确定随机误差ε，构成统计岩石物理模型m=f(R)+ε；

根据统计岩石物理模型m=f(R)+ε以及结合步骤101所述含气饱和度先验分布概率密度函数，采用现有基于马尔科夫链的蒙特卡洛随机模拟技术产生Ns个含气饱和度、弹性参数数据，构成含气饱和度与弹性参数数据联合训练样本空间{(m_k,R_k)}_{k＝1,2,...Ns}。

作为本发明的一个实施例，在所述步骤103中，根据训练样本统计出所述似然函数的分布信息，利用现有期望最大估计技术（EM）对似然函数分布信息进行高斯混合分布（GMM）函数逼近，得到所述似然函数的解析表达式，具体包括：

根据所述含气饱和度与弹性参数数据联合训练样本空间，统计出似然函数的直方分布结果，利用现有期望最大（EM）参数估计对似然函数分布信息进行高斯混合分布（GMM）函数逼近，得到似然函数的解析化表达式 $P\left(m^k\right|R_i)=\underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k{R}_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k{R}_i}^{\mathrm{\π}});$

作为本发明的一个实施例，在所述步骤104中，引入权重系数W_k来调节弹性参数数据信息的采用量，所述反演目标函数的解析化表达式为：

$\tilde{R}=\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_{i1}}^{\mathrm{\π}})$

其中，所述权重系数W_k由两部分组成：精度差异权重W_e与含气敏感性权重W_s。精度差异权重W_e由测井数据与弹性参数数据体中抽取的伪井曲线计算相关系数得到，敏感性权重W_s由Dillon敏感指示系数公式计算得到。所述精度差异权重W_e与含气敏感性权重W_s二者相乘，结果间进行归一化即可得到最终的权重系数W_k。

作为本发明的另一个实施例，所述的权重系数W_k还可以为根据经验、常识，而预先定义的数值。

作为本发明的一个实施例，在步骤105中，所述输入的弹性参数数据来自于叠前地震反演成果数据，主要是指纵波速度、横波速度、密度、纵波阻抗、横波阻抗、泊松比、弹性阻抗、拉梅参数等弹性参数数据；

所述反演结果指的是，在含气饱和度值区间R＝[R ₁R₂ ... R_N]内,给定任一含气饱和度值R_i，计算相应的反演目标函数值的大小（表示含气饱和度的后验概率密度），取反演目标函数最大值所对应的含气饱和度值为最终反演结果：

${\tilde{R}}_{\mathrm{final}}=\arg \underset{R_i\∈R}{\mathrm{Max}}\{\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k{R}_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k{R}_i}^{\mathrm{\π}})\}.$

如图2所示为本发明实施例一种含气饱和度反演方法的具体流程图。

记弹性参数数据为m，含气饱和度为R，将含气饱和度按照一定的取值区间分为N类，即R=[R₁ R₂ ... R_N]。依据贝叶斯反演框架，目标参数

为在给定弹性参数数据m的条件下最大后验概率密度所对应的一类含气饱和度取值，表示为：

$\begin{array}{c}{\tilde{R}}_{\mathrm{dinal}}=\arg \underset{R_i\∈R}{\mathrm{Max}}P\left(R_i\right|m)\\ =\arg \underset{R_i\∈R}{\mathrm{Msx}}\left\{\frac{P\left(m\right|R_i)\×P\left(R_I\right)}{P\left(m\right)}\right\}\end{array}$

其中P(R_i)为第i类含气饱和度参数取值的先验概率密度函数；P(m|R_i)为在给定含气饱和度参数值为R_i的前提下，m的条件概率密度函数，称似然函数；P(m)为常数，可以略去。

包括步骤201，综合测井资料、地质资料、叠前地震反演资料等信息，采用直方图统计的方式得到所述含气饱和度参数先验分布特征。含气饱和度先验分布特征的直方图示意图如图3所示。

步骤202，所述先验概率密度函数P(R_i)可根据含气饱和度先验分布特征，利用现有期望最大估计技术（EM）对含气饱和度先验分布信息进行高斯混合分布（GMM）函数逼近，得到含气饱和度先验概率密度函数。

在上述步骤后得到先验概率密度函数P(R_i)，然后的重点在于似然函数P(m|R_i)的确定。现有技术中的做法是利用统计岩石物理模型建立起弹性参数与含气饱和度间的统计关系，结合现有基于马尔科夫链的蒙特卡洛随机模拟技术产生含气饱和度与弹性参数的联合分布样本空间{(m_k,R_k)}_k＝_1,2,...Ns作为训练样本，进而对似然函数进行统计学意义上的求取：

$P\left(m_j\right|R_i)=\frac{n\left(m_j\right)}{n\left(R_j\right)}$

其中n(R_i)表示物性参数值等于R_i的样本个数；n(m_j)表示含气饱和度值等于R_i且弹性参数值为m_j的样本个数。通常情况下，为减少多解性，往往利用多个弹性参数变量进行反演，即对于Nm个弹性参数有：

$P\left(m_j\right|R_i)=P(m_j^1,m_j^2,...m_J^{\mathrm{Nm}}|R_j)=\frac{n(m_j^1,m_j^2,...m_j^{\mathrm{Nm}})}{n\left(R_j\right)}$

要得到合理的统计结果，则对于弹性参数值为

的样本，其在训练样本集中必须出现多次，这就要求训练样本集{(m_k,R_k)}_{k＝1,2,...Ns}非常巨大，因而统计过程相当耗时。

在本发明实施例中的步骤203，根据储层地质特点得到确定性岩石物理模型m=f(R)。

步骤204，在确定性岩石物理模型的基础上，根据模型与实际测井数据间的差异程度，确定随机误差ε，构成统计岩石物理模型m=f(R)+ε。

步骤205，根据统计岩石物理模型m=f(R)+ε以及步骤202中所述含气饱和度先验分布函数P(R)，采用现有基于马尔科夫链的蒙特卡洛随机模拟技术产生Ns个含气饱和度、弹性参数数据，构成含气饱和度与弹性参数联合样本空间{(m_k,R_k)}_{k＝1,2,...Ns}。

步骤206，根据所述含气饱和度与弹性参数联合样本空间，统计出所述似然函数的直方分布结果，利用现有期望最大估计技术（EM）对似然函数分布信息进行高斯混合分布（GMM）函数逼近，得到所述似然函数的解析表达：

$P\left(m^k\right|R_i)=\underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k{R}_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k{R}_i}^{\mathrm{\π}})\text{.}$

具体的，假设各弹性参数变量间相互条件独立，则有：

$P(m_j^1,m_j^2,...m_j^{\mathrm{Nm}}|R_i)=\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Π}}}P\left(m_j^k\right|R_i)$

即多个弹性参数变量联合条件概率密度变为单个弹性参数变量条件概率密度间的乘积，只需弹性参数变量

在样本集中单独出现多次即可。如图4所示，假设储层物性参数取值分为10类，有3个弹性参数变量，每个弹性参数变量取值个数均为j，如果每个样本数据在训练样本集中出现100次，则可以看到新策略所需训练样本总数量(虚线)随弹性参数变量取值个数j的变化较未作独立性假设前（实线）大大减少，因此统计效率也将大幅提高。

根据训练样本可以先统计出P(m^k|R_i)的分布信息，利用现有期望最大估计技术（EM）对似然函数分布信息进行高斯混合分布（GMM）函数逼近，得到P(m^k|R_i)的解析表达式（如图5所示为本发明实施例似然函数解析化表达式的求取示意图）：

$P\left(m^k\right|R_i)=\underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_i}^{\mathrm{\π}})$

其中N(·)表示高斯分布概率密度函数，N_c,α_π,

和

分别为与m^k、R有关的混合高斯分布部件数、权系数、均值与方差。

图5中的直方图代表某区纵波阻抗Ip在含气饱和度R_i为0.55时的统计特征，其中的3条高斯分布状曲线代表利用期望最大估计技术（EM）得到的三个高斯分布函数，这样似然函数便可由这三个高斯分布函数进行混合得到。对似然函数进行解析化能有效解决现有技术对似然函数的求取依赖反复进行样本统计，速度较慢的问题。

事实上，由叠前地震反演得到的各弹性参数变量间存在有精度差异且不同的弹性参数变量对含气饱和度变化的敏感度也不同。如果在反演过程中同等对待不同精度、不同含气敏感性的弹性参数变量，势必对含气饱和度反演结果造成影响。因此不同的弹性参数变量在目标函数中的权重也应不一样。

所以还包括步骤207，引入权重系数W_k来调节弹性信息的采用量，所述权重系数W_k由两部分组成：精度差异权重W_e与含气敏感性权重W_s。利用相关系数公式计算测井数据与弹性参数数据体中抽取的伪井曲线的相关系数即可得到精度差异权重，敏感性权重W_s由现有Dillon（2003）提出的敏感指示系数公式计算得到（具体出处请见Dillon,L.Schwedersky,G.,Vasquez,G.,Velloso,R.,and Nunes,C.,2003,A multiscale DHI elastic attributes evaluation:The Leading Edge,22no.10,1024-1029）。所述精度差异权重W_e与含气敏感性权重W_s二者相乘、进行归一化即可得到最终的权重系数W_k。

步骤208，根据所述似然函数和权重系数W_k得到最终的反演目标函数。

引入权重系数W_k来调节弹性信息的采用量，从而解决弹性参数精度与敏感度差异对含气饱和度反演的影响，取对数运算即可得到反演目标函数：

$\tilde{R}=\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k{R}_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k{R}_{i1}}^{\mathrm{\π}})$

步骤209，向所述反演目标函数输入的弹性参数数据来自于叠前地震反演成果数据，主要是指纵波速度、横波速度、密度、纵波阻抗、横波阻抗、泊松比、弹性阻抗、拉梅参数等弹性参数数据。

步骤210，所述反演结果指的是在含气饱和度值区间R＝[R₁ R₂ ... R_N]内,给定任一含气饱和度值R_i，计算相应的反演目标函数的大小（表示含气饱和度的后验概率密度），取反演目标函数最大值所对应的含气饱和度值为最终反演结果：

${\tilde{R}}_{\mathrm{final}}=\arg \underset{R_i\∈R}{\mathrm{Max}}\{\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_i}^{\mathrm{\π}})\}.$

以输入弹性参数为纵波阻抗、横波阻抗、密度等叠前地震反演成果数据为例。图6a为现有技术未考虑纵波阻抗、横波阻抗、密度等弹性参数间的精度差异以及对含气饱和度的敏感度差异，相当于对纵、横波阻抗及密赋予相同权重时的反演结果，图6b为本发明实施例对纵、横波阻抗及密度分别赋予权重1、0.6、0.8时的反演结果；虚线、实线分别代表含气饱和度模型实际值与反演值。通过反演结果对比可以看到，引入权重系数后，反演结果精度有了较大的提高。

图7、图8分为本发明实施例针对某研究区过A井含气饱和度反演剖面以及井旁道反演结果与测井结果对比图。综合剖面反演结果与测井曲线，可以看到，本发明实施例反演得到的含气饱和度反演结果（图8中虚线）与测井结果（图8中实线）吻合度较高，不同含气饱和度层段均得到正确显示。

图9、图10分别为过该研究区内未参与反演的一口验证井，B井的含气饱和度剖面及井旁道与B井曲线对比结果。可以看到，B井在未参与反演的情况下，B井处仍有较高的反演精度，从而验证了本发明实例含气饱和度预测的准确性。

图11所示为该研究区三维含气饱和度沿层切片，可以看到A井与B井均位于高含气饱和度区域，这与实际试气结果一致。另外，更为重要的是，地球物理工作者可以通过该沿层切片对未探明区域进行含气性预测并根据高含气饱和度值区域确定开发井位。

本发明的目的在于针对传统概率统计类方法易受弹性参数精度及含气敏感性差异影响以及反演效率较低无法应用于三维实际资料的问题，实现高精度、高效的含气饱和度定量反演方法，进而应用于三维实际资料。该发明提出了加权统计的概念以及解析化的目标反演函数表达式，利用测井资料统计分析与统计岩石物理建模相结合的思路，推导了较传统概率统计类方法效率更高、不受弹性参数精度及含气敏感性差异影响的目标反演函数及高效求解策略，方法具有较高的实用性，能有效地运用于三维实际资料，识别出高含气饱和度气藏的位置。

本发明实施例的优点在于：1）本发明实施例考虑到了弹性参数精度及敏感度对反演结果的影响，充分利用了地质、测井、叠前地震反演资料等信息，建立了弹性参数与含气饱和度之间的加权统计关系，理论更加严谨、更加符合实际，反演结果更加准确；2）高效的统计策略以及解析化的目标反演函数表达式，使得反演效率更高，能够应用于大规模三维实际资料，为下一步分析奠定了定量的基础；3）通过三维实际资料的实际应用，精准的描述了高含气饱和度储层的位置。

本发明可以以任何适当的形式实现，包括硬件、软件、固件或它们的任意组合。本发明可以根据情况有选择的部分实现，比如计算机软件执行于一个或多个数据处理器以及数字信号处理器。本文的每个实施例的元素和组件可以在物理上、功能上、逻辑上以任何适当的方式实现。事实上，一个功能可以在独立单元中、在一组单元中、或作为其他功能单元的一部分来实现。因此，该系统和方法既可以在独立单元中实现，也可以在物理上和功能上分布于不同的单元和处理器之间。

在相关领域中的技术人员将会认识到，本发明的实施例有许多可能的修改和组合，虽然形式略有不同，仍采用相同的基本机制和方法。为了解释的目的，前述描述参考了几个特定的实施例。然而，上述的说明性讨论不旨在穷举或限制本文所发明的精确形式。前文所示，许多修改和变化是可能的。所选和所描述的实施例，用以解释本发明的原理及其实际应用，用以使本领域技术人员能够最好地利用本发明和各个实施例的针对特定应用的修改、变形。

Claims (6)

1.一种含气饱和度反演方法，其特征在于包括，

根据含气饱和度参数R_i先验分布特征，得到含气饱和度先验概率密度函数

$P\left(R_i\right)=\underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n);$

上式中，N(·)表示高斯分布概率密度函数，N_r,α_n,

和

分别为与饱和度参数R_i有关的混合高斯分布部件数、权值、均值与方差；

设各弹性参数变量m^k间相互条件独立，则似然函数P(m_j|R_i)为：

$P(m_j^1,m_j^2,...m_j^{\mathrm{Nm}}|R_i)=\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Π}}}P\left(m_j^k\right|R_i);$

根据训练样本统计出所述似然函数P(m^k|R_i)的分布信息，得到所述似然函数P(m^k|R_i)的解析表达式：

$P\left(m^k\right|R_i)=\underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_i}^{\mathrm{\π}}),$

上式中，N(·)表示高斯分布概率密度函数，Nc,απ,

和分别为与m^k、R_i有关的混合高斯分布部件数、权值、均值与方差；

根据所述含气饱和度参数先验分布特征和所述似然函数P(m^k|R_i)的解析表达式以及引入权重系数W_k取对数运算得到反演目标函数的解析化表达式：

$\tilde{R}=\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_{i1}}^{\mathrm{\π}});$

上式中，ln表示取自然对数；

输入弹性参数数据，取所述反演目标函数最大后验概率密度解作为反演结果

${\tilde{R}}_{\mathrm{final}}=\arg \underset{R_i\∈R}{\mathrm{Max}}\{\ln \underset{n=1}{\overset{\mathrm{Nr}}{\mathrm{\Σ}}}{\mathrm{\α}}_{n}N(R_i;{\mathrm{\μ}}_{R_i}^n,{\mathrm{\Σ}}_{R_i}^n)+\underset{k=1}{\overset{\mathrm{Nm}}{\mathrm{\Σ}}}{W}_k\ln \underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k{R}_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k{R}_i}^{\mathrm{\π}})\},$

上式中，

表示返回极大值对应的R_i。

2.根据权利要求1所述的一种含气饱和度反演方法，其特征在于，在根据含气饱和度参数R_i先验分布特征，得到含气饱和度先验概率密度函数之中还包括，综合测井资料、地质资料、叠前地震反演成果资料，采用直方图统计的方式得到所述含气饱和度参数先验分布特征，利用期望最大估计技术对含气饱和度先验分布特征进行高斯混合分布函数逼近，得到含气饱和度先验概率密度函数。

3.根据权利要求1所述的一种含气饱和度反演方法，其特征在于，在根据训练样本统计出所述似然函数P(m^k|R_i)的分布信息进一步包括，

基于等效介质理论，根据储层地质特点选择确定性岩石物理模型m=f(R)；

根据确定性岩石物理模型与实际测井数据间的差异程度，确定随机误差ε，构成统计岩石物理模型m=f(R)+ε；

根据统计岩石物理模型m=f(R)+ε以及结合步骤101所述含气饱和度先验概率密度函数，采用基于马尔科夫链的蒙特卡洛随机模拟技术产生Ns个含气饱和度、弹性参数数据，构成含气饱和度与弹性参数数据联合训练样本空间{(m_k,R_k)}_{k＝1,2,...Ns}。

4.根据权利要求3所述的一种含气饱和度反演方法，其特征在于，在根据训练样本统计出所述似然函数P(m^k|R_i)的分布信息，得到所述似然函数的解析表达式进一步包括，根据所述含气饱和度与弹性参数数据联合训练样本空间，统计出似然函数的直方分布结果，利用期望最大参数估计对似然函数分布信息进行高斯混合分布函数逼近，得到似然函数的解析化表达式 $P\left(m^k\right|R_i)=\underset{\mathrm{\π}=1}{\overset{\mathrm{Nc}}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{\π}}N(R_i;{\mathrm{\μ}}_{m^k|R_i}^{\mathrm{\π}},{\mathrm{\Σ}}_{m^k|R_i}^{\mathrm{\π}}).$

5.根据权利要求1所述的一种含气饱和度反演方法，其特征在于，所述权重系数W_k由两部分组成：精度差异权重W_e与含气敏感性权重W_s；

所述精度差异权重W_e由测井数据与弹性参数数据体中抽取的伪井曲线计算相关系数得到；

所述敏感性权重W_s由Dillon敏感指示系数公式计算得到。

6.根据权利要求1所述的一种含气饱和度反演方法，其特征在于，所述输入的弹性参数数据包括纵波速度、横波速度、密度、纵波阻抗、横波阻抗、泊松比、弹性阻抗、拉梅参数等弹性参数数据。

Images