CN102004263A - 一种基于角度部分叠加地震数据的流体识别的方法 - Google Patents

Abstract

本发明提供一种基于角度部分叠加地震数据的流体识别的方法，分为以下三个步骤：1)快速的估计阻抗反射系数的方法：从Zoeppritz方程的Shuey近似(Shuey，1985)出发，推导得到准确快速估算纵横波阻抗反射系数的公式，计算用于流体识别的地震属性；2)基于角道集部分叠加地震数据的新流体识别属性提取：从Smith和Gidlow(Smith，2003年)提出的流体因子公式出发，得到应用角道集部分叠加地震数据的组合进行流体识别和储层预测；3)基于曲波变换和贝叶斯理论的流体识别方法：以1)和2)为基础，用子波均衡法，对各个角度叠加道集数据进行均衡，以消除在不同角度上的差异，保留由流体引起的不同角度道集上数据的差异；然后，利用曲波域的贝叶斯波场分离技术，得到目的层处的异常，命名为角度流体道集，再求取各种属性得到目的层处的异常；应用角道集部分叠加地震数据的各种组合得到的属性可以比较好的区别含气(油)和含水储层，特别是角度道集的三次幂所得到的属性剖面，可以较好地的区分含气(油)和含水储层。

Description

一种基于角度部分叠加地震数据的流体识别的方法

所属技术领域

本发明属于地震资料处理领域，是一种利用地震数据识别流体的技术。

背景技术

随着勘探难度的加大，对储层预测的精度要求越来越高，需要更有效的预测储层的方法和技术。常规地震属性分析是建立在叠后地震数据基础上的，叠后地震资料经过多次叠加后，信噪比有了较大的提高，但同时损失了大量的振幅信息。叠前地震数据较叠后地震数据包含着更加丰富的振幅和旅行时信息，能更灵敏的反映地下油气藏的变化，但是受噪音的影响较大，信噪比不高。角度部分叠加地震数据对两种数据的优点进行了折中，既有相对较高的信噪比又保留了丰富的信息量。因此，从叠前角道集部分叠加地震数据出发研究识别流体和预测储层的技术具有很高的实际应用价值。

发明内容

本发明的目的在于提供一种基于角度部分叠加地震数据的流体识别的方法。

本发明所采用的技术方案有：一种基于角度部分叠加地震数据的流体识别的方法，分为以下三个步骤：

1)快速的估计阻抗反射系数的方法：从Zoeppritz方程的Shuey近似(Shuey，1985)出发，推导得到准确快速估算纵横波阻抗反射系数的公式，计算用于流体识别的地震属性；

2)基于角道集部分叠加地震数据的新流体识别属性提取：从Smi th和Gidlow(Smith，2003年)提出的流体因子公式出发，得到应用角道集部分叠加地震数据的组合进行流体识别和储层预测；

3)基于曲波变换和贝叶斯理论的流体识别方法：以1)和2)为基础，用子波均衡法，对各个角度叠加道集数据进行均衡，以消除在不同角度上的差异，保留由流体引起的不同角度道集上数据的差异；然后，利用曲波域的贝叶斯波场分离技术，得到目的层处的异常，命名为角度流体道集，再求取各种属性得到目的层处的异常。

步骤1)的具体运算为：

由Zoeppritz方程的shuey简化式的假设条件可知，当入射角不大时(0°＜θ≤30°)Zoeppritz方程可以简化为两项的形式，即

R(θ)≈P+Gsin²θ(1)

其中第一项P为θ＝0°时的反射系数，称为AVO的截距；G为与岩石纵、横波速度和密度有关的项，称为AVO的梯度

在角度叠加道集中，任意选择两个道集数据，组成方程组，可以求出截距P和梯度G

$\left\{\begin{array}{c}{R}_1=P+{G\sin }^2{\mathrm{\θ}}_1\\ R_2=P+G{\sin }^1{\mathrm{\θ}}_2\end{array}\right.---\left(2\right)$

这里θ₁和θ₂为入射角，R₁，R₂分别为相应的反射系数

或者：

$\left\{\begin{array}{c}{R}_1=(1+{\sin }^2{\mathrm{\θ}}_1)\mathrm{Rp}-2{\sin }^2{\mathrm{\θ}}_1\mathrm{Rs}\\ R_2=(1+s{\mathrm{in}}^2{\mathrm{\θ}}_2)\mathrm{Rp}-2{\sin }^2{\mathrm{\θ}}_2\mathrm{Rs}\end{array}\right.---\left(3\right)$

其中，Rp为纵波阻抗反射系数，Rs为横波阻抗反射系数直接推导出纵，横波阻抗反射系数表达式：

$R_p=\frac{{\sin }^2{\mathrm{\θ}}_2}{{\sin }^2{\mathrm{\θ}}_2-{\sin }^2{\mathrm{\θ}}_1}(R_1-\frac{{\sin }^2{\mathrm{\θ}}_1}{{\sin }^2{\mathrm{\θ}}_2}{R}_2)$

$R_s=\frac{1+{\sin }^2{\mathrm{\θ}}_2}{2(s{\mathrm{in}}^2{\mathrm{\θ}}_2-{\sin }^2{\mathrm{\θ}}_1)}(R_1-\frac{1+{\sin }^2{\mathrm{\θ}}_1}{1+{\sin }^2{\mathrm{\θ}}_2}{R}_2)$

角度一般小于30°，由sinθ的泰勒展开式可以得到

sinθ≈θ

R_p和R_s可近似表示为：

$R_p\≈\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_2^2}{R}_2)---\left(4a\right)$

$R_s\≈\frac{1+{\mathrm{\θ}}_2^2}{2({\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_2^2}{R}_2)---\left(4b\right)$

计算R_p，R_s时，采用三次计算取平均的算法，可同时利用三个角度的信息，使计算结果更接近实际值；假设小、中、大三个角度道集的角度θ₁、θ₂和θ₃满足θ₂＝2θ₁，θ₃＝3θ₁，则得到如下估算公式：

$R_p=\frac{1}{3}\{\frac{{\mathrm{\θ}}_3^2}{{\mathrm{\θ}}_3^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_3^2}{R}_3)+\frac{{\mathrm{\θ}}_3^2}{{\mathrm{\θ}}_3^2-{\mathrm{\θ}}_2^2}(R_2-\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_3^2}{R}_3)+\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_2^2}{R}_2)\}$

$=\frac{1}{3}\{\frac{9}{8}\·(R_1-\frac{{\mathrm{\θ}}_1^2}{9{\mathrm{\θ}}_1^2}{R}_3)+\frac{9}{5}\·(R_2-\frac{4{\mathrm{\θ}}_1^2}{9{\mathrm{\θ}}_1^2}{R}_3)+\frac{4}{3}\·(R_1-\frac{{\mathrm{\θ}}_1^2}{4{\mathrm{\θ}}_1^2}{R}_2)\}$

$=\frac{59}{65}(R_1+0.6R_2-0.376R_3)---\left(5a\right)$

$R_s=\frac{1}{3}\{\frac{1+{\mathrm{\θ}}_3^2}{2({\mathrm{\θ}}_3^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_3^2}{R}_3)+\frac{1+{\mathrm{\θ}}_3^2}{2({\mathrm{\θ}}_3^2-{\mathrm{\θ}}_2^2)}(R_2-\frac{1+{\mathrm{\θ}}_2^2}{1+{\mathrm{\θ}}_3^2}{R}_3)+\frac{1+{\mathrm{\θ}}_2^2}{2({\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_2^2}{R}_2)\}$

$=\frac{1}{3}\{\frac{1+9{\mathrm{\θ}}_1^2}{8{\mathrm{\θ}}_1^2}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+9{\mathrm{\θ}}_1^2}{R}_3)+\frac{1+{9\mathrm{\θ}}_1^2}{5{\mathrm{\θ}}_1^2}(R_2-\frac{1+4{\mathrm{\θ}}_1^2}{1+{9\mathrm{\θ}}_1^2}{R}_3)+\frac{1+4{\mathrm{\θ}}_1^2}{3{\mathrm{\θ}}_1^2}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+4{\mathrm{\θ}}_1^2}{R}_2)\}$

$=\frac{11+{59\mathrm{\θ}}_1^2}{72{\mathrm{\θ}}_1^2}(R_1-\frac{8}{5}\frac{2-22{\mathrm{\θ}}_1^2}{11+59{\mathrm{\θ}}_1^2}{R}_2-\frac{3}{5}\frac{13+{37\mathrm{\θ}}_1^2}{11+59{\mathrm{\θ}}_1^2}{R}_3)---\left(5b\right)$

上式可以用R_p和R_s的相对大小表示(令R₁的系数为1)为如下形式：

R_p＝R₁+0.6R₂-0.376R₃(6a)

R_s＝R₁+αR₂-βR₃(6b)

公式(5a)和(5b)是分别用三次平均得到Rp和Rs的估算结果；其中公式(6b)中的系数α和β是随着小角度的入射角θ₁变化的。

步骤2)的具体运算为：

流体因子为含油气情况下当入射角为流体因子角度θ_f时的反射系数值，可以用流体因子角和截距P、梯度G的函数表示为：

流体因子R_{Fluid_Factor}＝P+Gsin²θ_f

(7)

上式显示储层含水时，流体因子为零。

步骤3)的具体运算：

a)角度流体道集

将公式(2)中的两式相减得到下式：

R₂-R₁＝G□(sin²θ₂-sin²θ₁)  (8)

命名为角度流体道集

b)角度部分叠加数据体预处理

首先从每个角度部分叠加道集中的非储层中开时窗提取出一个等效地震子波w_θ，计算每个等效地震子波的反子波

设原始地震数据为s_θ(t)，

由s_θ(t)＝r_θ(t)*w_θ得

r(t，θ)＝s_θ(t)*w(θ)^-1

令不同角度的等效子波的平均子波为

可得预处理后的角道集：

$s_{\mathrm{\θ}}^l\left(t\right)=r_{\mathrm{\θ}}\left(t\right)*\stackrel{\‾}{w}\left(t\right)=s_{\mathrm{\θ}}\left(t\right)*w{\left(\mathrm{\θ}\right)}^{-1}*\stackrel{\‾}{w}\left(t\right),$

c)曲波域的贝叶斯波场分离技术

设地震记录数据S由下式给出：

s＝s₁+s₂+n

s是大角度地震道集，s₁是角度流体道集，s₂是均衡后的小角度道集，n是白噪声

地震信号s_i可以表示成为曲波的叠加，

s₁＝Ax₁+n  i＝1，2，，

其中A是曲波合成矩阵(Candes，1998)，并获得系统方程，即：

s₁＝Ax₁+n₁

s₂＝Ax₂+n₂

其中x₁是角度流体道集的曲波系数；x₂是均衡后小角度数据体的曲波系数，n₁和n₂是白噪声

由贝叶斯理论可以得到：

$P(x_1,x_2|b_1,b_2)=\frac{P(x_1,x_2)P\left(b_1\right|x_1,x_2)P\left(b_2\right|b_1,x_1,x_2)}{P(b_1,b_2)}$

$\∝P(x_1,x_2)P\left(n\right)P\left(n_2\right)$

(9)

$\∝\exp (-A\left(l_n\right)-\frac{{\left|\right|{\mathrm{Ax}}_2-b_2\left|\right|}_2^2}{{\mathrm{\σ}}_2^2}-\mathrm{\η}\frac{{\left|\right|A(x_1+x_2)-(b_1+b_2)\left|\right|}_2^2}{{\mathrm{\σ}}^2})$

得到所求曲波系数的迭代公式

$x_1^{n+1}=[A^T{b}_2-A^{T}A{x}_2^n+A^T{b}_1-A^{T}A{x}_1^n+x_1^n]$

$x_2^{n+1}=[A^T{b}_2-A^{T}A{x}_2^n+x_2^n+\frac{\mathrm{\η}}{\mathrm{\η}+1}(A^T{b}_1-A^{T}A{x}_1^n)]---\left(10\right)$

继而得到角度流体道集的估计值是

${\tilde{s}}_1=A{\tilde{x}}_1.$

基于角道集部分叠加数据，提取流体属性分为以下几步：

1)应用公式(5)和公式(6)估算纵横波阻抗反射系数，将得到的估算值应用到流体因子属性(Fatti，1997)的计算中；

2)将不同的角道集进行相乘运算，得到反映储层流体特征的属性剖面，进行储层预测和流体识别；

3)角度流体道集的提取以及流体识别

A)角度道集预处理，应用前面提出的角度子波互均衡法，以消除非储层在不同角度上的差异，保留由流体引起的不同角度道集上数据的差异。

B)分离算法构建，应用曲波域非线性匹配滤波。

C)流体属性的提取，对分离后的流体道集或者均衡后的角度部分叠加数据进行运算提取叠向地震属性以及反映流体特征的属性，R_p和R_s属性，G属性以及流体因子属性。

本发明的效果为：按照上面的方法原理，首先对Marmousi II模型进行了分析，给出具体分析结果。为了能够说明问题而计算简单只选取了其中含有两个比较有代表性的储层的一部分，如图1所示，储层分别为含气砂岩层与含油砂岩层。

模型中共用到了三个角度的数据，分别为8度、16度、24度。

1)快速估算纵横波阻抗反射系数的方法及其应用

图2是从MarmousiII模型中抽取的第356道数据，此道分别含有气层，油层，水层。模型的纵、横波速度，密度如图所示，油、气、水所在的位置也已经在图中标出。

图3是应用不同方法估算的Rs和Rp值的比较图，图中圆圈是用上面提出的简单公式估算的结果，星号是应用两项shuey近似式计算的结果，黑线是应用纵、横波阻抗差计算的精确结果。从图中可以看出应用本发明的简单公式对Rp，Rs的估算是相当准确的，不论是在储层位置还是非储层位置，以及在气、水层位置和油层位置，都是十分准确的，能够满足生产的要求。

图4是应用Marmousi II模型对本发明进行测试的结果。从图3和图4中可以看出本发明能够比较准确的估算出纵、横波阻抗反射系数。

图5是利用三种方法计算的流体因子属性的对比。图5(a)是利用快速估计波阻抗反射系数法得到纵横波阻抗反射系数，然后将估算值代入流体因子属性计算公式(5)得到的属性剖面。图5(b)是应用角度子波均衡法将角道集进行预处理后，利用曲波域的贝叶斯分离技术将角道集进行分离，得到的角度流体道集剖面。图5(c)和(d)是应用角道集部分叠加地震数据的各种组合计算得到的流体识别属性剖面。从图5的对比中，可以看出三种方法都能比较清楚地刻画储层的位置。从图5(a)中，可以看到含气(油)和含水储层的差值稍大，因此，第一种方法区别含气(油)和含水储层的能力稍强；第二种方法不能够很好的区分含气(油)和含水储层，但是可以很好的压制不含流体地层的信息，特别是压制了含油储层下方的标志层，如图5(b)所示。而应用角道集部分叠加地震数据的各种组合得到的属性可以比较好的区别含气(油)和含水储层，特别是角度道集的三次幂所得到的属性剖面，可以较好地的区分含气(油)和含水储层。综合这三种方法的计算结果就可以预测储层的位置以及储层内所含流体的性质了。

附图说明：

图1为角道集示意图。其中(a)为小角度道集；(b)为中角度道集；(c)为大角度道集。

图2为Marmousi II模型数据。

图3为本发明简单近似结果与精确结果的比较。其中(a)为不同方法估计的Rp比较；(b)为不同方法估计的Rs比较。

图4为本发明的简单计算方法与精确结果的比较。其中(a)为阻抗差计算Rp；(b)为快速估计Rp；(c)为阻抗差计算Rs；(d)为快速估计Rs。

图5为三种方法计算流体属性比较。其中(a)为快速估计波阻抗反射系数法得到的流体因子属性；(b)为基于曲波变换得到的流体因子属性；(c)为三个角度道集相乘；(d)为大角度道集的三次方。

图6为Rp和Rs随着角度增大在角道集数据中的比重变化的示意图。

图7为贝叶斯参数估计流程图。

图8为实施例1的叠前角度道集。其中(a)为某地区8°角道集；(b)为某地区16°角道集；(c)为某地区24°角道集。

图9为实施例1的流体识别属性的比较。其中(a)为方法1计算的流体因子属性，(b)为方法2得到的流体属性(角道集相乘)；(c)为角度流体道集。

图10为实施例2的叠前角度道集。其中(a)为某地区8°角道集；(b)为某地区16°角道集；(c)为某地区24°角道集。

图11为不同方法计算的角度流体因子属性的对比。其中(a)为方法一计算得到的流体因子属性；(b)为方法二计算得到的流体因子属性(角度道集相乘)；(c)为基于曲波变换的角度流体因子属性。

具体实施方式

角道集数据是由叠前CMP道集数据转换而来的，它们反映的是地下同一位置处的信息。只是由于入射波角度的不同而在不同的角道集部分叠加剖面上存在一定的差别，因此不同角度的角道集部分叠加数据之间的地震属性也存在一定的差别。所以本发明从角道集部分叠加数据出发，研究识别流体和预测储层的有效方法和技术。本发明的研究思路是直接将角道集部分叠加数据认为是流体属性(这可以由Smith和Gidlow(Smith，2003年)的流体因子公式得出)，由其直接进行各种运算，以突出目的层处的异常，达到流体识别的目的。(也就是将角度流体道集代入本发明推导出的简单公式中进行计算，得到各种反映流体的流体属性)。据此，可以准确地判定储层的位置及范围，为油藏精细描述提供了一种较为可靠的研究思路。具体分为下面三部分：

(1)快速的估计波阻抗反射系数的方法从Zoeppritz方程的Shuey近似(Shuey，1985)出发，推导得到准确快速估算纵横波阻抗反射系数的公式，计算用于流体识别的地震属性。

(2)基于角道集部分叠加地震数据的新流体识别属性提取从Smith和Gidlow(Smith，Gidlow，2003年)提出的流体因子公式^[2]出发，经过深入研究得到应用角道集部分叠加地震数据的组合进行流体识别和储层预测的新属性技术。

(3)基于曲波变换和贝叶斯理论的流体识别技术以上面的研究方法为基础，提出子波均衡法，对各个角道集数据进行均衡，以消除地层骨架在不同角度上的差异，保留由流体引起的不同角度道集上数据的差异。然后，利用曲波域的贝叶斯波场分离技术，得到目的层处的异常，命名为角度流体道集，再求取各种属性。据此，可以准确地判定储层的位置及范围，为油藏精细描述提供可靠的依据。

下面详细论述方法技术原理：

第一、该技术的基本原理

1)快速的估计波阻抗反射系数的方法

在AVO研究中，纵、横波阻抗反射系数是流体识别和岩性预测的基础。叠前参数都可以由这两个参数直接或者间接计算得到，并且许多流体因子公式都是这两个参数的函数。因此，在岩性预测和流体识别中，这两个参数的求解就变得尤为重要。

由Zoeppritz方程的shuey简化式的假设条件可知，当入射角为中等时(0°＜θ≤30°)Zoeppritz方程可以简化为两项的形式，即

R(θ)≈P+Gsin²θ(1)

其中第一项P为θ＝0°时的反射系数，称为AVO的截距；G为与岩石纵、横波速度和密度有关的项，称为AVO的梯度。

这样在角度叠加道集中，任意选择两个道集数据，组成方程组，就可以求出截距P和梯度G。例如，

$\left\{\begin{array}{c}{R}_1=P+{G\sin }^2{\mathrm{\θ}}_1\\ R_2=P+G{\sin }^1{\mathrm{\θ}}_2\end{array}\right.---\left(2\right)$

这里θ₁和θ₂为入射角，R₁，R₂分别为相应的反射系数，。

公式(2)也可以写成以下形式：

$\left\{\begin{array}{c}{R}_1=(1+{\sin }^2{\mathrm{\θ}}_1)\mathrm{Rp}-2{\sin }^2{\mathrm{\θ}}_1\mathrm{Rs}\\ R_2=(1+s{\mathrm{in}}^2{\mathrm{\θ}}_2)\mathrm{Rp}-2{\sin }^2{\mathrm{\θ}}_2\mathrm{Rs}\end{array}\right.---\left(3\right)$

其中，Rp为纵波阻抗反射系数，Rs为横波阻抗反射系数。

可以看出，不同角度的反射系数可以看成是Rp和Rs的线性组合，随着角度的变化，Rp和Rs所占的比重也在变化，如图6所示：

在此图中，‘1’代表小角度；‘2’代表中角度；‘3’代表大角度，黑色表示Rp在角道集中所占的百分比，灰色表示Rs在角道集中所占的百分比。

从图6中可以看出小角度道集包含的Rp成分较多，随着角度的增大，道集数据中包含的Rs成分越来越多。

由公式(3)可以直接推导出纵，横波阻抗反射系数表达式：

$R_p=\frac{{\sin }^2{\mathrm{\θ}}_2}{{\sin }^2{\mathrm{\θ}}_2-{\sin }^2{\mathrm{\θ}}_1}(R_1-\frac{{\sin }^2{\mathrm{\θ}}_1}{{\sin }^2{\mathrm{\θ}}_2}{R}_2)$

$R_s=\frac{1+{\sin }^2{\mathrm{\θ}}_2}{2(s{\mathrm{in}}^2{\mathrm{\θ}}_2-{\sin }^2{\mathrm{\θ}}_1)}(R_1-\frac{1+{\sin }^2{\mathrm{\θ}}_1}{1+{\sin }^2{\mathrm{\θ}}_2}{R}_2)$

因为角度比较小(一般小于30°)，所以由sinθ的泰勒展开式可以得到

sinθ≈θ

此时相对误差小于5％。

这样R_p和R_s就可以近似表示为：

$R_p\≈\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_2^2}{R}_2)---\left(4a\right)$

$R_s\≈\frac{1+{\mathrm{\θ}}_2^2}{2({\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_2^2}{R}_2)---\left(4b\right)$

当计算R_p，R_s时，可以采用三次计算取平均的算法，这样可以同时利用三个角度的信息，使计算结果更接近实际值。假设小、中、大三个角度道集的角度θ₁、θ₂和θ₃满足θ₂＝2θ₁，θ₃＝3θ₁，则得到如下估算公式：

$R_p=\frac{1}{3}\{\frac{{\mathrm{\θ}}_3^2}{{\mathrm{\θ}}_3^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_3^2}{R}_3)+\frac{{\mathrm{\θ}}_3^2}{{\mathrm{\θ}}_3^2-{\mathrm{\θ}}_2^2}(R_2-\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_3^2}{R}_3)+\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_2^2}{R}_2)\}$

$=\frac{1}{3}\{\frac{9}{8}\·(R_1-\frac{{\mathrm{\θ}}_1^2}{9{\mathrm{\θ}}_1^2}{R}_3)+\frac{9}{5}\·(R_2-\frac{4{\mathrm{\θ}}_1^2}{9{\mathrm{\θ}}_1^2}{R}_3)+\frac{4}{3}\·(R_1-\frac{{\mathrm{\θ}}_1^2}{4{\mathrm{\θ}}_1^2}{R}_2)\}$

$=\frac{59}{65}(R_1+0.6R_2-0.376R_3)---\left(5a\right)$

$R_s=\frac{1}{3}\{\frac{1+{\mathrm{\θ}}_3^2}{2({\mathrm{\θ}}_3^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_3^2}{R}_3)+\frac{1+{\mathrm{\θ}}_3^2}{2({\mathrm{\θ}}_3^2-{\mathrm{\θ}}_2^2)}(R_2-\frac{1+{\mathrm{\θ}}_2^2}{1+{\mathrm{\θ}}_3^2}{R}_3)+\frac{1+{\mathrm{\θ}}_2^2}{2({\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_2^2}{R}_2)\}$

$=\frac{1}{3}\{\frac{1+9{\mathrm{\θ}}_1^2}{8{\mathrm{\θ}}_1^2}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+9{\mathrm{\θ}}_1^2}{R}_3)+\frac{1+{9\mathrm{\θ}}_1^2}{5{\mathrm{\θ}}_1^2}(R_2-\frac{1+4{\mathrm{\θ}}_1^2}{1+{9\mathrm{\θ}}_1^2}{R}_3)+\frac{1+4{\mathrm{\θ}}_1^2}{3{\mathrm{\θ}}_1^2}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+4{\mathrm{\θ}}_1^2}{R}_2)\}$

$=\frac{11+{59\mathrm{\θ}}_1^2}{72{\mathrm{\θ}}_1^2}(R_1-\frac{8}{5}\frac{2-22{\mathrm{\θ}}_1^2}{11+59{\mathrm{\θ}}_1^2}{R}_2-\frac{3}{5}\frac{13+{37\mathrm{\θ}}_1^2}{11+59{\mathrm{\θ}}_1^2}{R}_3)---\left(5b\right)$

为方便计算，上式可以用R_p和R_s的相对大小表示(令R₁的系数为1)为如下形式：

R_p＝R₁+0.6R₂-0.376R₃(6a)

R_s＝R₁+αR₂-βR₃(6b)

公式(5a)和(5b)是分别用三次平均得到Rp和Rs的估算结果。其中公式(6b)中的系数α和β是随着小角度的入射角θ₁变化的，其变化规律如表1所示。

表1：公式(5a)中α和β随角度的变化

系数\角度	5°	6°	7°	8°	9°	10°
							α	0.256	0.2416	0.2251	0.2069	0.1872	0.1663
β	0.696	0.6906	0.6844	0.6776	0.6702	0.6623

2)基于角道集部分叠加地震数据的新流体识别属性提取

Smith和Gidlow(Smith，2003年)在对Shuey(Shuey，1985年)两项近似式(1)分析之后，提出了应用流体因子角识别流体的方法。流体因子角(fluidfactor angle)θ_f就是当储层为含水砂岩时，Shuey(Shuey，1985)两项近似式的反射系数为零时的入射角。此时，流体因子为含油气情况下当入射角为流体因子角度θ_f时的反射系数值，可以用流体因子角和截距P、梯度G的函数表示为：

流体因子R_{Fluid_Factor}＝P+Gsin²θ_f

(7)

上式显示储层含水时，流体因子为零。从公式(1)和(7)的对比可知，Smith和Gidlow(Smith，2003年)提出的以流体因子角表示的流体因子(R_{Fulid_Factor})是shuey近似的特殊情况。公式(7)将流体因子属性与入射角，即角道集剖面联系了起来。鉴于此，提出了直接应用常规角道集的各种组合进行储层预测和流体识别的方法。比较发现，角道集部分叠加数据直接相乘具有与上述流体因子相似的作用，并且可以增大不同流体储层之间的差距从而更好地预测储层。

3)基于曲波变换和贝叶斯理论的流体识别技术

A)角度流体道集

将公式(2)中的两式相减可以得到下式：

R₂-R₁＝G□(sin²θ₂-sin²θ₁)(8)

这就是直接应用角度道集进行流体识别的理论依据。当消除子波随角度变化的影响之后，任意两个角度道集的差，是与AVO梯度有关的量，反映了振幅随入射角的变化量，此道集主要反映流体的变化，所以我们命名为角度流体道集。

B)角度部分叠加数据体预处理

不同角度的道集数据反映的是同一位置处的储层信息，它们只是由于叠加角度的不同才引起角道集部分叠加数据之间存在一定的差异。消除非储层位置处的差异，保留储层处的差异是我们要解决的主要问题。根据角道集自身的特点，本发明采用子波法对其进行互均衡处理。

假设子波随着入射角的变化而变化，同时也把不同角度道集上引起非储层差异的影响归于子波，这种子波是一种随角度变化的等效子波。因此，本发明首先提取不同角度道集的等效子波，然后将其替换为平均子波，这种方法我们称为角度子波均衡法。具体做法是：

首先从每个角度部分叠加道集中的非储层中开时窗提取出一个等效地震子波w_θ，计算每个等效地震子波的反子波

设原始地震数据为s_θ(t)，

由s_θ(t)＝r_θ(t)*w_θ得

r(t，θ)＝s_θ(t)*w(θ)^-1

令不同角度的等效子波的平均子波为

可得预处理后的角道集：

$s_{\mathrm{\θ}}^l\left(t\right)=r_{\mathrm{\θ}}\left(t\right)*\stackrel{\‾}{w}\left(t\right)=s_{\mathrm{\θ}}\left(t\right)*w{\left(\mathrm{\θ}\right)}^{-1}*\stackrel{\‾}{w}\left(t\right),$

C)曲波域的贝叶斯波场分离技术

贝叶斯法的基本思想是首先由地质和测井信息得到关于噪声ε和待求参数x的先验分布p(ε)和p(x)，其次用正演模型和p(ε)求出似然函数p(y|x)，再通过贝叶斯公式得到参数的后验分布p(x|y)，最后通过p(x|y)得到关于参数x的信息，流程图如图7所示。

设地震记录数据S由下式给出：

s＝s₁+s₂+n

s是大角度地震道集，s₁是角度流体道集，s₂是均衡后的小角度道集，n是白噪声。

地震信号s_i可以表示成为曲波的叠加，

s_i＝Ax₁+n  i＝1，2，，

其中A是曲波合成矩阵(Candes，1998)，并获得系统方程，即：

s₁＝Ax₁+n₁

s₂＝Ax₂+n₂

其中x₁是角度流体道集的曲波系数；x₂是均衡后小角度数据体的曲波系数，n₁和n₂是白噪声。

这样，就可以采用贝叶斯理论的后验概率分布函数(PPDF)建立目标函数，在曲波域用迭代方法求解模型参数估计值。在这里假设x1和x2具有的先验分布为Huber分布。由贝叶斯理论可以得到：

$P(x_1,x_2|b_1,b_2)=\frac{P(x_1,x_2)P\left(b_1\right|x_1,x_2)P\left(b_2\right|b_1,x_1,x_2)}{P(b_1,b_2)}$

$\∝P(x_1,x_2)P\left(n\right)P\left(n_2\right)$

(9)

$\∝\exp (-A\left(l_n\right)-\frac{{\left|\right|{\mathrm{Ax}}_2-b_2\left|\right|}_2^2}{{\mathrm{\σ}}_2^2}-\mathrm{\η}\frac{{\left|\right|A(x_1+x_2)-(b_1+b_2)\left|\right|}_2^2}{{\mathrm{\σ}}^2})$

应用反复重加权迭代阈值算法解此最优化问题。得到所求曲波系数的迭代公式

$x_1^{n+1}=[A^T{b}_2-A^{T}A{x}_2^n+A^T{b}_1-A^{T}A{x}_1^n+x_1^n]$

$x_2^{n+1}=[A^T{b}_2-A^{T}A{x}_2^n+x_2^n+\frac{\mathrm{\η}}{\mathrm{\η}+1}(A^T{b}_1-A^{T}A{x}_1^n)]---\left(10\right)$

继而得到角度流体道集的估计值是

${\tilde{s}}_1=A{\tilde{x}}_1.$

第二、该技术的实现过程

基于角道集部分叠加数据，提取流体属性分为以下几步：

1)应用公式(5)和公式(6)快速准确的估算纵横波阻抗反射系数，将得到的估算值应用到流体因子属性的计算中。例如将估算值代入Fatti(Fatti，1997年)提出的流体因子公式或者Smith和Gidlow(Smith，2003年)提出的以流体因子角表示的流体因子属性公式中计算流体因子，得到流体识别的属性值，进行流体识别和储层预测。

2)将不同的角道集进行相乘运算，得到反映储层流体特征的属性剖面，进行储层预测和流体识别。

3)角度流体道集的提取以及流体识别

a角度道集预处理，应用前面提出的角度子波互均衡法，以消除非储层在不同角度上的差异，保留由流体引起的不同角度道集上数据的差异。

b分离算法构建，这一步关键是构建迭代的分离算法，此方法是稳健的并且能够减小角度子波均衡法中的误差。传统的在匹配滤波之后再相减的分离方法，会产生具有残留振幅，调谐，倾角，相位和子波误差的信号预测值。应用曲波域非线性匹配滤波能够改进这种方法的分离值，如果预测信号是位置和倾角的函数，此种匹配法能处理重要的振幅误差。

c流体属性的提取。对分离后的流体道集(或者均衡后的角度部分叠加数据)进行各种运算提取叠向地震属性以及反映流体特征的属性，如R_p和R_s属性，G属性以及流体因子属性(应用上面提出的简单公式推导计算得到)，研究它们之间的差异。据此，可以准确地判定储层的位置及范围，为油藏精细描述提供较为可靠的依据。

4)将三种方法计算结果进行对比分析，预测有利的储层发育带和储层所含流体性质。

根据上述理论，对实际地震数据进行了叠前地震属性提取及分析。首先对角道集数据进行角度子波法预处理，然后从均衡后的角道集数据中提取流体因子。

图8和图10为两个实例的三个角道集剖面。图9和图11为应用三种方法提取的叠前地震属性。从图9a和图11a中可以看出应用快速估算波阻抗反射系数的方法得到纵横波阻抗反射系数后，再应用于流体因子公式(7)得到的流体因子属性与原始的大角度道集有点相似，这与前面的分析是十分一致的，但是还是可以突出储层的位置。应用不同角道集相乘的方法得到的新流体因子属性能够比较的好的突出储层的位置，并且对骨架信息有比较好的压制(如图9b和图11b所示)。经过角度子波法预处理和曲波波场分离之后非目的层处的差异减小，因此利用曲波波场分离计算出的属性能更明显的突出储层(如图9c和图11c所示)。从结果中也可以看出曲波波场分离能够突出储层位置。

从对以上模型与实际资料的处理结果中可以看出，从角道集部分叠加数据中提取的地震属性在模型的分析中应用效果较好。从对实际地震数据的分析结果中可以看出，应用方法2和方法3从角道集部分叠加数据中提取出的流体因子属性能够较清晰的刻画出油气储层的边界。

Claims (5)

1.一种基于角度部分叠加地震数据的流体识别的方法，其特征在于分为以下三个步骤：

1)快速的估计波阻抗反射系数的方法：从Zoeppritz方程的Shuey近似(Shuey，1985)出发，推导得到准确快速估算纵横波阻抗反射系数的公式，计算用于流体识别的地震属性；

2)基于角道集部分叠加地震数据的新流体识别属性提取：从Smith和Gidlow(Smith，2003年)提出的流体因子公式出发，得到应用角道集部分叠加地震数据的组合进行流体识别和储层预测；

3)基于曲波变换和贝叶斯理论的流体识别方法：以1)和2)为基础，用子波均衡法，对各个角度叠加道集数据进行均衡，以消除在不同角度上的差异，保留由流体引起的不同角度叠加道集上数据的差异；然后，利用曲波域的贝叶斯波场分离技术，得到目的层处的异常，命名为角度流体道集，再求取各种属性，以进行储层预测和流体识别。

2.根据权利要求1所述的一种基于角度部分叠加地震数据的流体识别的方法，其特征在于步骤1)的具体运算为：

由Zoeppritz方程的shuey简化式的假设条件可知，当入射角不大时(0°＜θ≤30°)Zoeppritz方程可以简化为两项的形式，即

R(θ)≈P+Gsin²θ(1)

其中第一项P为θ＝0°时的反射系数，称为AVO的截距；G为与岩石纵、横波速度和密度有关的项，称为AVO的梯度

在角度叠加道集中，任意选择两个道集数据，组成方程组，可以求出截距P和梯度G

$\left\{\begin{array}{c}{R}_1=P+{G\sin }^2{\mathrm{\θ}}_1\\ R_2=P+G{\sin }^1{\mathrm{\θ}}_2\end{array}\right.---\left(2\right)$

这里θ₁和θ₂为入射角，R₁，R₂分别为相应的反射系数

或者：

$\left\{\begin{array}{c}{R}_1=(1+{\sin }^2{\mathrm{\θ}}_1)\mathrm{Rp}-2{\sin }^2{\mathrm{\θ}}_1\mathrm{Rs}\\ R_2=(1+s{\mathrm{in}}^2{\mathrm{\θ}}_2)\mathrm{Rp}-2{\sin }^2{\mathrm{\θ}}_2\mathrm{Rs}\end{array}\right.---\left(3\right)$

其中，Rp为纵波阻抗反射系数，Rs为横波阻抗反射系数

直接推导出纵，横波阻抗反射系数表达式：

$R_p=\frac{{\sin }^2{\mathrm{\θ}}_2}{{\sin }^2{\mathrm{\θ}}_2-{\sin }^2{\mathrm{\θ}}_1}(R_1-\frac{{\sin }^2{\mathrm{\θ}}_1}{{\sin }^2{\mathrm{\θ}}_2}{R}_2)$

$R_s=\frac{1+{\sin }^2{\mathrm{\θ}}_2}{2(s{\mathrm{in}}^2{\mathrm{\θ}}_2-{\sin }^2{\mathrm{\θ}}_1)}(R_1-\frac{1+{\sin }^2{\mathrm{\θ}}_1}{1+{\sin }^2{\mathrm{\θ}}_2}{R}_2)$

因为上面假设角度小于30°，所以由sinθ的泰勒展开式可以得到

sinθ≈θ

R_p和R_s可近似表示为：

$R_p\≈\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_2^2}{R}_2)---\left(4a\right)$

$R_s\≈\frac{1+{\mathrm{\θ}}_2^2}{2({\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_2^2}{R}_2)---\left(4b\right)$

计算R_p，R_s时，采用三次计算取平均的算法，可同时利用三个角度的信息，使计算结果更接近实际值；假设小、中、大三个角度道集的角度θ₁、θ₂和θ₃满足θ₂＝2θ₁，θ₃＝3θ₁，则得到如下估算公式：

$R_p=\frac{1}{3}\{\frac{{\mathrm{\θ}}_3^2}{{\mathrm{\θ}}_3^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_3^2}{R}_3)+\frac{{\mathrm{\θ}}_3^2}{{\mathrm{\θ}}_3^2-{\mathrm{\θ}}_2^2}(R_2-\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_3^2}{R}_3)+\frac{{\mathrm{\θ}}_2^2}{{\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2}(R_1-\frac{{\mathrm{\θ}}_1^2}{{\mathrm{\θ}}_2^2}{R}_2)\}$

$=\frac{1}{3}\{\frac{9}{8}\·(R_1-\frac{{\mathrm{\θ}}_1^2}{9{\mathrm{\θ}}_1^2}{R}_3)+\frac{9}{5}\·(R_2-\frac{4{\mathrm{\θ}}_1^2}{9{\mathrm{\θ}}_1^2}{R}_3)+\frac{4}{3}\·(R_1-\frac{{\mathrm{\θ}}_1^2}{4{\mathrm{\θ}}_1^2}{R}_2)\}$

$=\frac{59}{65}(R_1+0.6R_2-0.376R_3)---\left(5a\right)$

$R_s=\frac{1}{3}\{\frac{1+{\mathrm{\θ}}_3^2}{2({\mathrm{\θ}}_3^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_3^2}{R}_3)+\frac{1+{\mathrm{\θ}}_3^2}{2({\mathrm{\θ}}_3^2-{\mathrm{\θ}}_2^2)}(R_2-\frac{1+{\mathrm{\θ}}_2^2}{1+{\mathrm{\θ}}_3^2}{R}_3)+\frac{1+{\mathrm{\θ}}_2^2}{2({\mathrm{\θ}}_2^2-{\mathrm{\θ}}_1^2)}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+{\mathrm{\θ}}_2^2}{R}_2)\}$

$=\frac{1}{3}\{\frac{1+9{\mathrm{\θ}}_1^2}{8{\mathrm{\θ}}_1^2}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+9{\mathrm{\θ}}_1^2}{R}_3)+\frac{1+{9\mathrm{\θ}}_1^2}{5{\mathrm{\θ}}_1^2}(R_2-\frac{1+4{\mathrm{\θ}}_1^2}{1+{9\mathrm{\θ}}_1^2}{R}_3)+\frac{1+4{\mathrm{\θ}}_1^2}{3{\mathrm{\θ}}_1^2}(R_1-\frac{1+{\mathrm{\θ}}_1^2}{1+4{\mathrm{\θ}}_1^2}{R}_2)\}$

$=\frac{11+{59\mathrm{\θ}}_1^2}{72{\mathrm{\θ}}_1^2}(R_1-\frac{8}{5}\frac{2-22{\mathrm{\θ}}_1^2}{11+59{\mathrm{\θ}}_1^2}{R}_2-\frac{3}{5}\frac{13+{37\mathrm{\θ}}_1^2}{11+59{\mathrm{\θ}}_1^2}{R}_3)---\left(5b\right)$

上式可以用R_p和R_s的相对大小表示(令R₁的系数为1)为如下形式：

R_p＝R₁+0.6R₂-0.376R₃(6a)

R_s＝R₁+αR₂-βR₃(6b)

公式(5a)和(5b)是分别用三次平均得到Rp和Rs的估算结果；其中公式(6b)中的系数α和β是随着小角度的入射角θ₁变化的。

3.根据权利要求1所述的一种基于角度部分叠加地震数据的流体识别的方法，其特征在于步骤2)的具体运算为：

流体因子为含油气情况下当入射角为流体因子角度θ_f(Smith，2003)时的反射系数值，可以用流体因子角和截距P、梯度G的函数表示为：

流体因子R_{flui_Factor}＝P+Gsin²θ_f

(7)

上式显示储层含水时，流体因子为零。

4.根据权利要求1所述的一种基于角度部分叠加地震数据的流体识别的方法，其特征在于步骤3)的具体运算：

a)角度流体道集

将公式(2)中的两式相减得到下式：

R₂-R₁＝G□(sin²θ₂-sin²θ₁)(8)

命名为角度流体道集

b)角度部分叠加数据体预处理

首先从每个角度部分叠加道集中的非储层中开时窗提取出一个等效地震子波w_θ，计算每个等效地震子波的反子波

设原始地震数据为s_θ(t)，

由s_θ(t)＝r_θ(t)*w_θ得

r(t，θ)＝s_θ(t)*w(θ)^-1

令不同角度的等效子波的平均子波为可得预处理后的角度叠加道集：

$s_{\mathrm{\θ}}^l\left(t\right)=r_{\mathrm{\θ}}\left(t\right)*\stackrel{\‾}{w}\left(t\right)=s_{\mathrm{\θ}}\left(t\right)*w{\left(\mathrm{\θ}\right)}^{-1}*\stackrel{\‾}{w}\left(t\right),$

c)曲波域的贝叶斯波场分离技术

设地震记录数据S由下式给出：

s＝s₁+s₂+n

s是大角度地震道集，s₁是角度流体道集，s₂是均衡后的小角度道集，n是白噪声

地震信号s₁可以表示成为曲波的叠加，

s_i＝Ax₁+n    i＝1，2，，

其中A是曲波合成矩阵(Candes，1998)，并获得系统方程，即：

s₁＝Ax₁+n₁

s₂＝Ax₂+n₂

其中x₁是角度流体道集的曲波系数；x₂是均衡后小角度数据体的曲波系数，n₁和n₂是白噪声

由贝叶斯理论可以得到：

$P(x_1,x_2|b_1,b_2)=\frac{P(x_1,x_2)P\left(b_1\right|x_1,x_2)P\left(b_2\right|b_1,x_1,x_2)}{P(b_1,b_2)}$

$\∝P(x_1,x_2)P\left(n\right)P\left(n_2\right)$

$\∝\exp (-A\left(l_n\right)-\frac{{\left|\right|{\mathrm{Ax}}_2-b_2\left|\right|}_2^2}{{\mathrm{\σ}}_2^2}-\mathrm{\η}\frac{{\left|\right|A(x_1+x_2)-(b_1+b_2)\left|\right|}_2^2}{{\mathrm{\σ}}^2})$

得到所求曲波系数的迭代公式

$x_1^{n+1}=[A^T{b}_2-A^{T}A{x}_2^n+A^T{b}_1-A^{T}A{x}_1^n+x_1^n]$

$x_2^{n+1}=[A^T{b}_2-A^{T}A{x}_2^n+x_2^n+\frac{\mathrm{\η}}{\mathrm{\η}+1}(A^T{b}_1-A^{T}A{x}_1^n)]---\left(10\right)$

继而得到角度流体道集的估计值是

${\tilde{s}}_1=A{\tilde{x}}_1.$

5.根据权利要求1所述的一种基于角度部分叠加地震数据的流体识别的方法，其特征在于基于角道集部分叠加数据，提取流体属性分为以下几步：

1)应用公式(5)和公式(6)估算纵横波阻抗反射系数，将得到的估算值应用到流体因子属性(Fatti，1997)的计算中；

2)将不同的角道集进行相乘运算，得到反映储层流体特征的属性剖面，进行储层预测和流体识别；

3)角度流体道集的提取以及流体识别

A)角度道集预处理，应用前面提出的角度子波互均衡法，以消除非储层在不同角度上的差异，保留由流体引起的不同角度道集上数据的差异。

B)分离算法构建，应用曲波域非线性匹配滤波。

C)流体属性的提取，对分离后的流体道集或者均衡后的角度部分叠加数据进行运算提取叠向地震属性以及反映流体特征的属性，R_p和R_s属性，G属性以及流体因子属性。

Images