CN104018826B - 基于非线性算法的裂缝及裂缝内流体识别方法 - Google Patents

Abstract

一种基于非线性算法的裂缝及裂缝内流体识别方法。本发明基于三项Rüger反射系数方程，通过不同方位相减的办法推导出包含Rüger方程第三项(sin²θtan²θ项)和Thomsen各向异性参数的非线性方程。利用非线性算法－粒子群算法，从所推导的方程中求解出裂缝密度和Thomsen各向异性参数等。利用求解出的裂缝密度、各向异性参数推导出了裂缝流体指示因子表达式。本发明在实际应用中，基于所提出的非线性方程得到的裂缝密度远远比常规线性AVOZ反演得到的裂缝密度接近理论值，反演的各向异性参数的值也非常接近理论值；利用非线性方程的解得到的裂缝流体指示因子也很好地识别出了裂缝流体类型，充分说明了该方法的可行性。

Description

基于非线性算法的裂缝及裂缝内流体识别方法

技术领域

本发明属于石油勘探领域，涉及一种裂缝储层综合评价中裂缝以及流体类型预测的方法，尤其涉及一种利用三项Rüger方程的裂缝参数非线性反演方法。

背景技术

裂缝既是油气的储集空间，也是流体的运移通道，能够将储层中分散的孔洞储集空间连通起来。因此，识别裂缝和裂缝内流体对油气勘探具有重大意义。储层内发育裂缝时对外表现为各向异性特征，对于垂直排列的平行裂缝，可以将其等效为HTI介质(HorizontalTransverseIsotropy，即具有水平对称轴的横向各向同性介质)模型来研究，常用Thomsen参数表征HTI介质的各向异性特征(Thomsen，1995)。描述这类垂直排列的平行裂缝的参数主要有裂缝密度和裂缝走向等。目前工业界大多用常规AVOZ方法进行裂缝探测，即基于两项Rüger反射系数公式(Rüger，1998)。实际上，三项Rüger反射系数公式更能表征裂缝特征参数，并且，该公式中的第三项包含Thomsen各向异性参数。下面给出Rüger（1998）推导的反射系数公式：

R_PP(θ,φ)=A+(Biso+Banicos²φ)sin²θ+Csin²θtan²θ(1)

其中，

$A=\frac{1}{2}\frac{\mathrm{\ΔZ}}{\stackrel{\‾}{Z}};\mathrm{Biso}=\frac{1}{2}[\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\ΔG}}{\stackrel{\‾}{G}}]$

$\mathrm{Bani}=\frac{1}{2}[{\mathrm{\δ}}^v+2{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\mathrm{\γ}]{\cos }^2\mathrm{\φ};C=\frac{1}{2}[\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}+{\mathrm{\ϵ}}^v{\cos }^4\mathrm{\φ}+{\mathrm{\δ}}^v{\cos }^2\mathrm{\φ}{\sin }^2\mathrm{\φ}]$

式中：α、β分别为纵波速度和横波速度；ε^v、δ^v、γ是Thomsen各向异性参数；为采集方位角；为裂缝走向；θ为入射角；Z为纵波垂直入射时的波阻抗；G为横波的切向模量。在（1）式中的第二项（即与sin²i相乘的项）是AVO梯度，梯度项分为两部分：各向同性梯度Biso和各向异性梯度Bani（该参数能够指示裂缝引起的各向异性程度即裂缝密度）。

基于Rüger近似公式的AVOZ反演方法一般只用其前两项，不用第三项(含sin²θtan²θ项)。实际上，在用Rüger反射系数近似公式描述地下某个成像点的各向异性特征时非常受第三项的影响，尤其是对入射角相对较大的情况。目前随着地震采集技术的提高，可以获得大偏移距数据，相应地，反射波入射角范围也较大。当入射角范围较大(大于30°)时，Rüger公式中的第三项(sin²θtan²θ)对纵波反射系数的贡献也很大，不能忽略。因此，为了获得较精确的地下裂缝特征参数，如何充分利用Rüger方程中的第三项是需要解决的一个问题。

同时，裂缝流体也是目前地球物理工作者研究的热点和难点。根据Schoenberg和Sayers(1995)的线性滑动理论，裂缝流体指示因子(fluidindicator)可以写成如下形式：

$\mathrm{Fluidindicator}\≈g\frac{\mathrm{\ΔN}}{\mathrm{\ΔT}}---\left(2\right)$

Bakulin等(2000)推导出各向异性参数和线性滑动理论中法向柔量△N、切向柔量△T的关系，如下公式：

$\left\{\begin{array}{c}{\mathrm{\ϵ}}^{\left(v\right)}=-2g(1-g)\mathrm{\ΔN}\\ {\mathrm{\δ}}^{\left(v\right)}=-2g[(1-2g)\mathrm{\ΔN}+\mathrm{\ΔT}]\\ {\mathrm{\γ}}^{\left(v\right)}=-\mathrm{\ΔT}/2\end{array}\right.---\left(3\right)$

同时，根据Rüger（1998）明确γ^(v)和γ的关系：γ^(v)=-γ。

据Schoenberg和Douma(1988)，在裂缝饱含不同流体类型时，法向柔量△N和切向柔量△T可以表示为：

(1)对于干裂缝(干的或者充填气)，k'=μ'=0，则：

$\left\{\begin{array}{c}\mathrm{\ΔN}=\frac{4e}{3g(1-g)}\\ \mathrm{\ΔT}=\frac{16e}{3(3-2g)}\end{array}\right.---\left(4\right)$

(2)对于湿裂缝(充填油或者水)，则有：

$\left\{\begin{array}{c}\mathrm{\ΔN}=0\\ \mathrm{\ΔT}=\frac{16e}{3(3-2g)}\end{array}\right.---\left(5\right)$

其中即横纵波速度比的平方；e是Hudson理论中的裂缝体密度。

基于以上理论，虽然有很多人也对如何识别裂缝内流体类型做了研究，但是大多只是定性描述。

发明内容

本发明提出一种基于非线性算法的裂缝及裂缝内流体识别方法，它是基于三项Rüger方程，该方法能够用于实际资料的裂缝以及裂缝内流体类型的识别。通过本发明提供的方法，能够识别出裂缝发育区，并能识别出裂缝内充填的流体类型，为更精细的裂缝储层描述提供依据。

本发明实现上述目的的具体方案如下：

一种基于非线性算法的裂缝及裂缝内流体识别方法，其特征在于，所述方法包括如下步骤：

步骤1：把分方位共反射点道集分别记为：方位1、方位2、方位3…方位n，接着把n-1个方位的共反射点道集分别减去某一个方位的共反射点道集，得到n-1个方位的共反射点振幅差异道集，为非线性反演准备数据；

步骤2：利用n个方位的共反射点道集进行裂缝反演，算出每个反射点上每个时间采样点上的裂缝方位；

步骤3：把步骤2中得到的裂缝方位输入求初值的线性方程，得到求解非线性方程所需要的初值，这些初值包括裂缝密度、各向异性参数；

步骤4：利用步骤1所述的振幅差异道集，步骤2得到的裂缝方位、步骤3得到的裂缝密度和各向异性参数作为初始值，基于所述非线性方程进行反演，得到裂缝密度、各向异性参数；

步骤5：利用步骤4得到的裂缝密度、各向异性参数求得裂缝流体指示因子的大小，进行流体类型识别。

所述的步骤2中所用的裂缝反演方法是基于两项Rüger方程的常规线性AVOZ反演。

所述步骤2中的每个反射点上每个时间采样点上的裂缝方位取自世界著名的Rüger反射系数方程。

所述的步骤3的初值的求取是基于一个线性方程，而不是随机地给定非线性反演所需要的初始值。

所述的步骤4的实施是基于正演信息的指导，是通过数值模型得到：从世界著名的三项Rüger方程出发，通过不同方位相减的办法推导出一个新的非线性方程，通过求解该方程得到裂缝密度和各向异性参数。数值模型的应用表明：基于非线性方程得到的裂缝密度远远比常规线性AVOZ反演得到的裂缝密度接近理论值，另外，反演的各向异性参数的值也非常接近理论值。

所述的步骤4中非线性方程的求解利用的是粒子群算法。

所述的步骤5中裂缝流体指示因子，是从1995年Schoenberg和Sayers推导的裂缝流体指示因子g△N/△T出发，结合各向异性参数和法向柔量△N、切向柔量△T的关系，推导出裂缝流体指示因子一个表达形式，而且数值模型的应用表明利用非线性方程的解得到的裂缝流体指示因子也很好地识别出了裂缝流体类型。

本发明与已有技术相比产生的有益效果是：

在反射波入射角较大的情况下(>30°)，基于两项Rüger反射系数方程的线性AVOZ反演不能够精确描述地下裂缝诱导的各向异性特征。然而，目前随着野外采集技术的提高，地震资料的偏移距越来越大，也意味着入射角越来越大。

1、本发明基于三项Rüger反射系数方程，通过不同方位相减的办法推导出包含Rüger方程第三项(sin²θtan²θ项)和Thomsen各向异性参数的非线性方程。利用非线性算法－粒子群算法，从所推导的方程中求解出裂缝密度和Thomsen各向异性参数等。利用求解出的裂缝密度、各向异性参数推导出了裂缝流体指示因子的新的表达形式。

2、数值模型的应用表明：基于所提出的非线性方程得到的裂缝密度远远比常规线性AVOZ反演得到的裂缝密度接近理论值，反演的各向异性参数的值也非常接近理论值；

3、利用非线性方程的解得到的裂缝流体指示因子也很好地识别出了裂缝流体类型。

4、将本发明应用到塔里木哈7井区奥陶系碳酸盐岩储层中，反演求出的裂缝参数与FMI结果吻合，且比常规线性AVOZ反演结果更精确，得到的裂缝流体因子也和井生产情况吻合，充分说明了该方法的可行性。

附图说明

图1是基于非线性算法的裂缝及裂缝内流体识别方法流程图；

图2（a）是方位角道集图；

图2（b）不同方位做差后的道集图；

图3（a）是无噪音情况下非线性反演和常规AVOZ反演的裂缝密度Bani对比图；

图3（b）是无噪音情况下非线性反演和常规AVOZ反演的各向异性参数ε^v对比图；

图3（c）是无噪音情况下非线性反演和常规AVOZ反演的各向异性参数δ^v对比图；

图4表示无噪音情况下第19个CDP点的目标函数值随迭代次数的变化情况图；

图5（a）是含15%噪音情况下非线性反演和常规AVOZ反演的裂缝密度Bani对比图；

图5（b）是含15%噪音情况下非线性反演和常规AVOZ反演的各向异性参数ε^v对比图；

图5（c）是含15%噪音情况下非线性反演和常规AVOZ反演的各向异性参数δ^v对比图；

图6表示含15%噪音情况下第19个CDP点的目标函数值随迭代次数的变化情况图。

图7（a）是无噪音情况下裂缝流体因子的理论值、初始值、非线性反演结果的对比图；

图7（b）是含15%噪音情况下裂缝流体因子的理论值、初始值、非线性反演结果的对比图；

图8是不同方位的偏移距道集图；

图9是不同方位做差后的偏移距道集图；

图10是测井储层参数综合评价图；

图11是过井叠前时间偏移剖面图；

图12是基于非线性反演的裂缝密度剖面图；

图13是裂缝流体指示因子图。

具体实施方式

以下结合实施例与附图说明对本发明做进一步说明。

图1是基于非线性算法的裂缝及裂缝内流体识别方法流程图。

一种基于非线性算法的裂缝及裂缝内流体识别方法，包括如下步骤：

步骤1：把分方位共反射点道集分别记为：方位1、方位2、方位3…方位6，接着把5个方位的共反射点道集分别减去某一个方位的共反射点道集，得到5个方位的共反射点振幅差异道集，为非线性反演准备数据；

步骤2：利用6个方位的共反射点道集进行裂缝反演，算出每个反射点上每个时间采样点上的裂缝方位（裂缝方位取自世界著名的Rüger反射系数方程（1））；

步骤3：把步骤2中得到的裂缝方位输入本发明所推导的求初值的线性方程（8），得到求解非线性方程所需要的初值，这些初值包括裂缝密度Bani、各向异性参数ε^v和δ^v；

步骤4：利用步骤1的振幅差异道集，步骤2得到的裂缝方位、步骤3得到的裂缝密度和各向异性参数作为初始值，基于本发明推导的非线性方程（6）进行反演，得到裂缝密度Bani、各向异性参数ε^v和δ^v；

步骤5：利用步骤4得到的裂缝密度、各向异性参数求得裂缝流体指示因子的大小，进行流体类型识别。

所述的步骤2中所用的裂缝反演方法是基于两项Rüger方程的常规线性AVOZ反演。

所述步骤2中的每个反射点上每个时间采样点上的裂缝方位取自世界著名的Rüger反射系数方程。

所述的步骤3的初值的求取是基于一个线性方程，而不是随机地给定非线性反演所需要的初始值。

所述的步骤4的实施是基于正演信息的指导，是通过数值模型得到：从世界著名的三项Rüger方程出发，通过不同方位相减的办法推导出一个新的非线性方程，通过求解该方程得到裂缝密度和各向异性参数。数值模型的应用表明：基于非线性方程得到的裂缝密度远远比常规线性AVOZ反演得到的裂缝密度接近理论值，另外，反演的各向异性参数的值也非常接近理论值。

所述的步骤4中非线性方程的求解利用的是粒子群算法。

所述的步骤5中裂缝流体指示因子，是从1995年Schoenberg和Sayers推导的裂缝流体指示因子g△N/△T出发，结合各向异性参数和法向柔量△N、切向柔量△T的关系，推导出裂缝流体指示因子一个表达形式，而且数值模型的应用表明利用非线性方程的解得到的裂缝流体指示因子也很好地识别出了裂缝流体类型。

本发明所提供的基于非线性算法的裂缝及裂缝内流体识别方法的基本原理如下：

从公式(1)出发，对于两个不同方位的反射地震数据，分别记为R_pp(θ,φ₁)和R_pp(θ,φ₂)，对两者做差，把这个差值记为△R_pp(θ,△φ)，并且△φ=φ₁-φ₂，通过一系列的化简，可以得到公式(6)，

$\begin{array}{c}{\mathrm{\ΔR}}_{\mathrm{pp}(\mathrm{\θ},\mathrm{\Δ\φ})}=\frac{B}{4}({a\cos 2\mathrm{\φ}}_s+{b\sin 2\mathrm{\φ}}_s){\sin }^2\mathrm{\θ}+\frac{E}{4}({a\cos 2\mathrm{\φ}}_s+{b\sin 2\mathrm{\φ}}_s){\sin }^2\mathrm{\θ}{\tan }^2\mathrm{\θ}\\ +\frac{E-D}{16}[{2\mathrm{ac}\cos 4\mathrm{\φ}}_s+(\mathrm{ad}+\mathrm{bc}){\sin 4\mathrm{\φ}}_s]{\sin }^2\mathrm{\θ}{\tan }^2\mathrm{\θ}\end{array}---\left(6\right)$

式中，

a=cos2φ_i-cos2φ₀；b=sin2φ_i-sin2φ₀

c=cos2φ_i+cos2φ₀；d=sin2φ_i+sin2φ₀

B=δ^v+8g△r；E=ε^v；D=δ^v

φ_i，φ₀均表示野外采集测线方位，φ₀同时还表示被减去的方位；φ_s为裂缝方位；表示垂直横纵波速度比的平方；ε^v，δ^v，r是Thomsen各向异性参数。非线性方程(6)中的B表示常规线性AVOZ反演Bani，E，D分别是Thomsen各向异性参数ε^v和δ^v。方程(6)中B，E，D和φ_s是要求解的未知量。

要求解方程(6)就要建立目标函数，本发明中将不同方位的叠前共反射点道集数据作为非线性反演的输入，利用公式(6)求解计算出的△R(θ,φ)与实际观测的误差E作为目标函数：

$E\left(t\right)=\underset{i}{\overset{N}{\mathrm{\Σ}}}\underset{j}{\overset{M}{\mathrm{\Σ}}}{[{\mathrm{\ΔR}}_{\mathrm{obs}}({\mathrm{\θ}}_i,{\mathrm{\φ}}_j)-{\mathrm{\ΔR}}_{\mathrm{cal}}({\mathrm{\θ}}_i,{\mathrm{\φ}}_j)]}^2---\left(7\right)$

式中，△R_obs(θ,φ)为实际观测的两个方位反射振幅之差；△R_cal(θ,φ)是由公式(6)计算得到的理论值；θ_i为第i个入射角和透射角的平均角度；φ_j指第j个方位数据的平均方位角；t为迭代次数。通过求解非线性方程(7)可以得到描述裂缝特征的参数裂缝密度Bani(各向异性强度)、裂缝方位φ_s，以及各向异性参数ε^v和δ^v。

本发明使用一种新型的非线性反演方法－粒子群优化算法(PSO)方法来求解方程(7)，这就需要B，E，D和φ_s的初值。如果已知裂缝方位φ_s(本文中使用的φ_s来自于AVOZ反演)，把其代入方程(6)，那么可以得到线性方程(8)，如下：

△R_pp(i,△φ)=BM₀sin²i+5EM₁sin²itan²i+DM₂sin²itan²i(8)

式中，

$\left\{\begin{array}{c}{M}_0=\frac{a\cos {2\mathrm{\φ}}_s+{b\sin 2\mathrm{\φ}}_s}{4}\\ M_1=\frac{{4a\cos 2\mathrm{\φ}}_s+{4b\sin 2\mathrm{\φ}}_s+{2\mathrm{ac}\cos 4\mathrm{\φ}}_s+(\mathrm{ad}+\mathrm{bc})\sin {4\mathrm{\φ}}_s}{16}\\ M_2=-\frac{2\mathrm{ac}\cos {4\mathrm{\φ}}_s+(\mathrm{ad}+\mathrm{bc}){\sin 4\mathrm{\φ}}_s}{16}\end{array}\right.$

式(8)可以写为，

AX=△R(9)

这样，通过方程(8)求得的B，E，D作为求解目标函数(7)的初值。

接着，结合方程(3)和方程(6)中各向异性参数的表示形式，可以把式(2)中裂缝流体指示因子简化成如下形式：

$\mathrm{Fluidindicator}\≈g\frac{\mathrm{\ΔN}}{\mathrm{\ΔT}}=\frac{{2\mathrm{\ϵ}}^v-{\mathrm{\δ}}^v-B}{{\mathrm{\δ}}^v-B}=\frac{2E-D-B}{D-B}=2\frac{E-D}{D-B}+1---\left(10\right)$

公式(10)中的E，B，D均可以从方程(7)中求出。

本发明中裂缝探测的适用范围如下：对于无裂缝和裂缝走向呈现不同方向的井（多组裂缝），各向异性程度弱，建立在Rüger方程基础上裂缝参数非线性反演方法无法预测出井附近裂缝发育区和井中的流体类型；只有对于裂缝走向呈现某一优势方向（单组裂缝）的井，其各向异性程度较强时，才能预测其中的裂缝和裂缝流体类型。

图2用数值模型正演模拟说明不同方位道集做差后能量差异更明显。首先利用Rüger公式（包含第三项）进行模型正演。设计一双层介质模型，上层为均匀各向同性介质层，下层为HTI介质层，模型参数如表1，并且HTI介质中充填气。同时，设定裂缝纵横比为0.001，裂缝密度范围0.005～0.15，间隔是0.005，即共有30个CDP点。利用Hudson(1981，1986)等效介质理论和Thomsen理论(1986)，求得各个CDP点的Thomsen各向异性参数。设定入射角为15°、25°和35°，方位角为0°、30°、60°、90°、120°、150°。利用公式(1)得到反射系数，再将30个CDP点的反射系数序列分别与40Hz的零相位雷克子波进行褶积，得到的合成地震记录。图1(a)是第30个CDP点的方位角道集，图1(b)是第30个CDP点的0°、30°、60°、120°、150°方位分别减去90°方位后的道集，从图中可以看出做差后的不同道集之间的能量差异明显增大，这为下面进行非线性反演得到比较好的结果奠定了数学基础。

表1各向同性/HTI双层单界面模型参数

图3是不含噪音的情况下非线性反演和常规AVOZ反演的结果对比情况，非线性反演结果中包括裂缝密度Bani、各向异性参数ε^v和δ^v，充分说明了非线性反演比常规AVOZ反演更精确。求解计算误差函数也是目标函数E(t)最小的具体步骤如下：①每个CDP点都以90°方位作为被减方位，得到5个方位做差后的道集，作为目标函数的观测值；②以合成的地震记录为输入，利用AVOZ反演求得裂缝方位φ_s，求解过程详见Sun等(2012a)；③把步骤②得到的φ_s输入方程(8)求出用于非线性反演的初始值E0，B0，D0；④把通过步骤②得到的φ_s和步骤③得到的E0，B0，D0作为目标函数(7)的初值进行PSO反演得到E，B，D和φ_s，PSO反演的基本原理详见Sun等(2012b)。本发明所展示的均是反射界面处的结果。图3(a)是通过不同方法求得的裂缝密度Bani对比结果，黑色实线表示理论值，带圈实线表示线性反演的结果，即AVOZ反演的结果，绿色虚线表示初始值，蓝色点表示非线性反演的结果。可以看出，在裂缝内充填流体的情况下，用线性AVOZ反演的Bani值偏离理论值很多，尤其是裂缝密度较大的情况，在比较接近实际值的初始值的约束下，通过非线性反演得到的Bani值比较接近于理论值。图3(b)和图3(c)分别是通过非线性反演得到的各向异性参数ε^v和δ^v与理论值、以及初始值的结果，实线表示理论值，虚线表示初始值，蓝色点表示非线性反演的结果。从图3(b)和图3(c)中可以看出，非线性反演得到的ε^v和δ^v都比较接近于理论值，非线性反演结果相对于理论值的偏离程度小于初始值相对于理论值的偏离程度。

图4表示无噪音情况下第19个CDP点的目标函数值随迭代次数的变化情况，这说明随着迭代次数的增加，观测值和计算值之间的误差越来越小，所求的解越来越接近理论值。

图5是含15%随机噪音的情况下非线性反演和常规AVOZ反演的结果对比情况，非线性反演结果中包括裂缝密度Bani、各向异性参数ε^v和δ^v，这不仅说明了非线性反演比常规AVOZ反演更精确，更说明了该方法有一定的抗噪音能力。其具体求解步骤同不含噪音的情况相同。图5(a)是通过不同方法求得的裂缝密度Bani对比结果，黑色实线表示理论值，带圈实线表示线性反演的结果，即AVOZ反演的结果，绿色虚线表示初始值，蓝色点表示非线性反演的结果。可以看出，所求得的裂缝密度Bani的整体变化特征和无噪音时候相同，由于有噪音的存在AVOZ反演结果、非线性反演结果和初始值均表现不同程度的抖动，相对于初始Bani，通过非线性反演得到的Bani更接近于理论值。图11和图12分别是通过非线性反演得到的各向异性参数ε^v和δ^v与理论值、以及初始值的结果，实线表示理论值，虚线表示初始值，蓝色点表示非线性反演的结果。从这两张图中可以看出，非线性反演得到的ε^v和δ^v都比较接近于理论值，并且非线性反演结果相对于理论值的偏离程度小于初始值相对于理论值的偏离程度。

图6是第19个CDP点的目标函数值随迭代次数的变化情况，可见用本发明提出的非线性反演算法，在噪音存在的情况下目标函数仍然能够快速收敛，误差值也相对较小且合理，通过图4和图6对比可以看出由于噪音的存在，目标函数收敛所需要的迭代次数会相应地增大。本发明只显示了第19个CDP点的目标函数随迭代次数的变化情况，实际上不同的CDP点，其收敛情况是不同的，从数值模型的反演结果上看，用PSO方法求解所推导的非线性方程是可行的，结果也比较接近于理论值。

图7同样的，在其他模型参数不变的条件下按着上面的步骤求出了裂缝内充填水(包括无噪音和有噪音两种情况)的裂缝密度Bani(方程(6)中的B)，各向异性参数ε^v(方程(6)中的E)和δ^v(方程(6)中的D)。通过公式(10)可以知道如果已知B、E、D就可以求出裂缝流体指示因子。通过公式(10)和(11)不难发现当裂缝内充填流体时裂缝流体指示因子近似为0，当裂缝内充填气时裂缝流体指示因子不为0，且其值远大于0，根据此可以用流体指示因子判断裂缝内的流体类型。图7(a)和图7(b)均是用非线性方程的解B、E、D求得的裂缝流体指示因子，其中图7(a)显示的是无噪音的结果，图7(b)是添加15%的随机噪音的结果。可以看出当裂缝内饱含水时，流体指示因子较小接近于0，当裂缝充填气时，流体指示因子的值较大，通过裂缝流体指示因子二者可以明显被区分开来。通过图7(b)可看出，虽然受噪音的影响，初始值和非线性反演的结果均在一定程度上偏离理论值，但是这不影响利用流体因子区分裂缝流体类型。

以下以该方法在新疆某区实际三维数据的一条测线为例具体说明该方法的实际应用效果。

图8是不同方位的偏移距道集。其中每个玫瑰图中红色表示对应的采集方位。

图9不同方位做差后的道集。它是由图8中方位2～方位6分别减去方位1后得到的5个振幅差异道集。可以看出，相减后的振幅差异道集的信噪比明显降低，随机噪音比重大幅增加，但是在目标层(3950ms)附近，依然表现出一定的特征，经过相减后的道集之间的差异更明显。把不同方位的振幅差异道集按不同的入射角范围叠加作为非线性反演的输入。

图10是研究井的测井储层参数综合评价图。本发明选择一条过井测线进行研究，目标层在一间房组和鹰山组一段。该井钻遇串珠群储层，它是由裂缝连通的溶蚀孔而成的。玫瑰图是由FMI资料统计，由玫瑰图显示该井张开缝发育，并且FMI显示张开缝倾角较陡，各向异性程度强，声成像解释结论表明该井段裂缝发育，测井评价结果也显示目标层内裂缝孔隙度比较高。

图11表示过井的叠前时间偏移剖面。地震表现为串珠响应特征。

图12表示利用本发明求出的过井裂缝密度。可以看出，利用本发明求出的裂缝密度在井的周围表现高值，即强各向异性特征，与井资料吻合。

图13表示利用本发明求出的裂缝流体指示因子。利用非线性反演得到的各向异性参数ε^v、δ^v和裂缝密度Bani值，根据公式(13)得到流体指示因子。根据公式(4)和(5)可知，当裂缝内充填油或水时，流体指示因子的值较低，而充填气时，流体指示因子值较高。生产数据表明该井是口油井，井累计产液52197t，产油51392t，流体指示因子在井位置表现低值，与生产情况吻合，这表明利用本发明求的裂缝流体因子可以较好地识别裂缝流体类型。

Claims (3)

1.一种基于非线性算法的裂缝及裂缝内流体识别方法，其特征在于，所述方法包括如下步骤：

步骤1：把分方位共反射点道集分别记为：方位1、方位2、方位3…方位n，接着把n-1个方位的共反射点道集分别减去某一个方位的共反射点道集，得到n-1个方位的共反射点振幅差异道集，为非线性反演准备数据；

步骤2：利用n个方位的共反射点道集进行裂缝反演，算出每个反射点上每个时间采样点上的裂缝方位；

步骤3：把步骤2中得到的裂缝方位输入求初值的线性方程，得到求解非线性方程所需要的初值，所述初值包括裂缝密度、各向异性参数；

步骤4：利用步骤1所述的振幅差异道集，步骤2得到的裂缝方位、步骤3得到的裂缝密度和各向异性参数作为初始值，基于所述非线性方程进行反演，得到裂缝密度、各向异性参数；

步骤5：利用步骤4得到的裂缝密度、各向异性参数求得裂缝流体指示因子的大小，进行流体类型识别；

所述的步骤2中所用的裂缝反演方法是基于两项Rüger方程的常规线性AVOZ反演；

所述步骤2中的每个反射点上每个时间采样点上的裂缝方位取自Rüger反射系数方程；

所述的步骤3的初值的求取是基于一个线性方程，而不是随机地给定非线性反演所需要的初始值；

步骤3中所述求初值的线性方程是：

△R_pp(i,△φ)＝BM₀sin²i+5EM₁sin²itan²i+DM₂sin²itan²i(8)

式中，

步骤4中所述非线性方程是：

式中，

a＝cos2φ_i-cos2φ₀；b＝sin2φ_i-sin2φ₀

c＝cos2φ_i+cos2φ₀；d＝sin2φ_i+sin2φ₀

B＝δ^v+8g△r；E＝ε^v；D＝δ^v；

其中：φ_i，φ₀均表示野外采集测线方位，φ₀同时还表示被减去的方位；φ_s为裂缝方位；表示垂直横纵波速度比的平方；ε^v，δ^v，r是Thomsen各向异性参数，非线性方程(6)中的B表示常规线性AVOZ反演Bani，E，D分别是Thomsen各向异性参数ε^v和δ^v，方程(6)中B，E，D和φ_s是要求解的未知量；

所述的步骤5中裂缝流体指示因子的大小，进行流体类型识别，是从1995年Schoenberg和Sayers推导的裂缝流体指示因子出发，结合各向异性参数和法向柔量△N、切向柔量△T的关系，推导出裂缝流体指示因子一个表达形式，而且数值模型的应用表明利用非线性方程的解得到的裂缝流体指示因子也很好地识别出了裂缝流体类型。

2.根据权利要求1所述的基于非线性算法的裂缝及裂缝内流体识别方法，其特征在于，所述的步骤4的实施是基于正演信息的指导，是通过数值模型得到：从三项Rüger方程出发，通过不同方位相减的办法推导出一个非线性方程，通过求解该方程得到裂缝密度和各向异性参数。

3.根据权利要求1所述的基于非线性算法的裂缝及裂缝内流体识别方法，其特征在于，所述的步骤4中非线性方程的求解利用的是粒子群算法。