CN104375182A - 一种裂缝型储层流体的识别方法及装置 - Google Patents

Abstract

本申请公开了一种裂缝型储层流体的识别方法及装置，该识别方法包括以下步骤：将三维地震纵波资料进行处理，生成多方位的道集数据；利用所生成的道集数据进行反演得到储层的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度和各向同性梯度的比值利用裂缝横向宽度与纵向高度的比值d、储层流体的第二拉梅常数μ'以及背景介质拉梅常数λ,μ对进行修正，获得流体因子I_fluid；利用该流体因子I_fluid对储层流体进行识别。识别装置相应的包括：道集数据生成单元、第一级流体因子计算单元、第二级流体因子获取单元、流体因子获取单元以及流体识别单元。通过本申请中的识别方法及装置能减小表征储层流体异常特性的流体因子的相对误差，准确识别裂缝型储层流体的类型。

Description

一种裂缝型储层流体的识别方法及装置

技术领域

本申请涉及石油勘探领域，尤其涉及一种裂缝型储层流体的识别方法及装置。

背景技术

裂缝型油气藏的产量占全世界石油天然气总产量的一半以上，在国内的油气勘探中，裂缝型油气藏是一个重要的勘探领域。碳酸盐岩、泥岩、砂岩、火山岩和变质岩中均有裂缝发育。裂缝既是油、气、水等地下流体的运移通道，也是流体的储集空间，裂缝能极大地提升储层的渗透率。因此，识别裂缝中所贮存的流体，对于寻找油气储层具有重大意义。

利用地震资料研究岩石裂缝所贮存流体的类型和特性，是根据岩石物理等效介质理论，将与储层流体有关的异常特性表征为流体因子，再依托流体因子实现储层流体的类型判识。流体识别的质量主要取决于两方面，一是反演计算的弹性参数是否可靠；二是构建的流体因子对流体类型是否敏感且稳定。储层内发育裂缝时对外表现为各向异性特征，所以用各向异性参数来构建流体因子。所谓各向异性是指地震波的物理性质(振幅、速度、频率等)随观测方位的变化而变化。对于垂直排列的平行裂缝，可以将其等效为HTI介质(HorizontalTransverse Isotropy，即具有水平对称轴的横向各向同性介质)模型来研究。若入射角为θ，测线方位角为φ，Ruger认为HTI介质分界面处的纵波反射系数是各向异性参数、入射角和方位角的函数：

$\begin{array}{c}R(\mathrm{\θ},\mathrm{\φ})=\frac{1}{2}\frac{\mathrm{\ΔZ}}{\stackrel{\‾}{Z}}+\frac{1}{2}\{\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\ΔG}}{\stackrel{\‾}{G}}+[\mathrm{\Δ}{\mathrm{\δ}}^V+2{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\mathrm{\Δ\γ}\left]{\cos }^2(\mathrm{\φ}-{\mathrm{\φ}}_s)\right\}{\sin }^2\mathrm{\θ}\\ +\frac{1}{2}(\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}+{\mathrm{\Δ\ϵ}}^V{\cos }^4(\mathrm{\φ}-{\mathrm{\φ}}_s)+{\mathrm{\Δ\δ}}^V{\sin }^2(\mathrm{\φ}-{\mathrm{\φ}}_s){\cos }^2(\mathrm{\φ}-{\mathrm{\φ}}_s)){\sin }^2\mathrm{\θ}{\tan }^2\mathrm{\θ}\end{array}---\left(1\right)$

上式中：φ_s为裂缝平面的对称轴方向，即裂缝面法向，(φ-φ_s)为测线与裂缝面法向的夹角；α,β,Z,G分别是纵波速度、横波速度、纵波阻抗和剪切模量；ε^V,δ^V,γ为各向异性参数。分别表示上下两层介质间纵波速度的平均值、横波速度的平均值、纵波阻抗的平均值和剪切模量的平均值，Δα,ΔZ,ΔG,Δε^V,Δδ^V,Δγ分别表示上下两层介质间纵波速度的差值、纵波阻抗的差值、剪切模量的差值以及三个各向异性参数的差值。

在石油工业领域，目前常用下面公式(2)进行AVAZ反演，该算法忽略了上述Ruger公式(1)的四次项会带来的系统误差，并且计算各向异性梯度参数B^ani时需要开平方，这将导致反演结果存在二义性，由此计算的各向同性梯度参数B^iso和各向异性梯度参数B^ani存在较大相对误差，最大相对误差高达10％，从而使得后续利用此两者构建的流体因子也会存在较大误差，并且现有算法中并没有对流体因子进行修正，不能很好的达到准确识别储层流体类型的目的。

R(θ,φ)≈R₀+[B^iso+B^ani cos²(φ-φ_s)]sin²θ         (2)

其中: $R_0=\frac{1}{2}\frac{\mathrm{\ΔZ}}{\stackrel{\‾}{Z}},B^{\mathrm{iso}}=\frac{1}{2}[\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\ΔG}}{\stackrel{\‾}{G}}],B^{\mathrm{ani}}=\frac{1}{2}[\mathrm{\Δ}{\mathrm{\δ}}^V+{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\mathrm{\Δ\γ}].$

发明内容

本申请的目的是提供一种通过进行修正来减小所构建的流体因子的相对误差的裂缝型储层流体的识别方法及装置，以准确判断地下裂缝介质是否充填油气。

为解决上述技术问题，本申请通过以下技术方案实现：

本申请提供了一种裂缝型储层流体的识别方法，该识别方法包括以下步骤：

S1，将获取的整个研究区域的三维地震纵波资料进行各向异性偏移处理，生成道集数据，所述道集数据包括三维地震纵波射线的方位角、入射角和反射振幅的信息；

S2，利用所生成的道集数据中反射振幅R与入射角θ和方位角φ之间的关系，并考虑大于30度的入射角θ对反射振幅R的影响，即考虑大偏移距对反射振幅R的影响，进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度B^ani和各向同性梯度B^iso的比值，将所述比值记为第一级流体因子；

S3，利用裂缝横向宽度与纵向高度的比值d对所述第一级流体因子进行修正，获得第二级流体因子，所述裂缝横向宽度与纵向高度的比值d是通过运用背景介质裂缝结构的各向异性反演计算所得到的；

S4，利用所获得的第二级流体因子对储层流体类型进行区分。

优选的，所述利用所生成的道集数据中反射振幅R与入射角θ和方位角φ之间的关系，并考虑入射角θ大于30度时对反射振幅R的影响，进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，包括：

利用反射振幅R与入射角θ和方位角φ之间的下述关系式进行AVAZ反演，得到储层的各向同性梯度B^iso、各向异性梯度B^ani以及其它四个参数R₀,C^iso,C^ani1和C^ani2。

R(θ,φ)＝R₀+(B^iso+B^anicos²φ')sin²θ+(C^iso+C^ani1cos⁴φ'+C^ani2sin²φ'cos²φ')sin²θtan²θ

进一步的，所述AVAZ反演包括：

将所述反射振幅R与入射角θ和方位角φ之间的关系式简化为：

R(θ,φ)＝AR₀+BB^iso+CB^ani+DC^iso+EC^ani1+FC^ani2，其中A、B、C、D、E和F是入射角和方位角的组合：

A＝1

B＝sin²θ

C＝sin²θcos²φ'

D＝sin²θtan²θ

E＝sin²θtan²θcos⁴φ'

F＝sin²θtan²θsin²φ'cos²φ'

将上述公式构建超定方程组：M_n×6X_6×1＝R_n×1

其中n为覆盖次数， $M={\left[\begin{array}{cccccc}{A}_1& B_1& C_1& D_1& E_1& F_1\\ A_2& B_1& C_2& D_2& E_1& F_2\\ ...& ...& ...& ...& ...& ...\\ A_n& B_1& C_n& D_n& E_n& F_n\end{array}\right]}_{n\×6},X={\left[\begin{array}{c}{R}_0\\ B^{\mathrm{iso}}\\ B^{\mathrm{ani}}\\ C^{\mathrm{iso}}\\ C^{\mathrm{ani}1}\\ C^{\mathrm{ani}2}\end{array}\right]}_{6\×1},R={\left(\begin{array}{c}{R}_1\\ R_2\\ ...\\ R_n\end{array}\right)}_{n\×1}$

通过求解所述超定方程组，得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani。

优选的，所述步骤S2包括在获得所述第二级流体因子后，利用储层流体的第二拉梅常数μ'以及背景介质的拉梅常数λ,μ对第二级流体因子进行修正，获得流体因子I_fluid，所述背景介质的拉梅常数λ,μ是通过利用所述道集数据中的反射振幅反演计算得到的；

所述步骤S3为利用流体因子I_fluid对储层流体进行区分。

优选的，所述利用流体因子I_fluid对储层流体进行区分包括：

将所获得的流体因子I_fluid与预设的阈值进行对比来判断储层流体的类型，

或，

将所获得的流体因子I_fluid与现有的流体因子谱库中的数据进行对比来判断储层流体的类型。

本申请还提供了一种裂缝型储层流体的识别装置，该装置包括：

道集数据生成单元，所述道集数据生成单元用于将获取的整个研究区域的三维地震纵波资料进行各向异性偏移处理，生成道集数据，所述道集数据包括三维地震纵波射线的方位角、入射角和反射振幅的信息；

第一级流体因子计算单元，所述第一级流体因子计算单元用于利用所生成的道集数据，通过对方位角和入射角进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度B^ani和各向同性梯度B^iso的比值，将所述比值记为第一级流体因子

第二级流体因子获取单元，所述第二级流体因子获取单元用于利用裂缝横向宽度与纵向高度的比值d对所述第一级流体因子进行修正，获得第二级流体因子裂缝横向宽度与纵向高度的比值d是通过运用背景介质裂缝结构的各向异性反演计算得到的；

流体识别单元，所述流体识别单元用于利用所获得的第二级流体因子对储层流体进行识别。

优选的，所述第一级流体因子计算单元包括第一子单元，所述第一子单元用于利用反射振幅R与入射角θ和方位角φ之间的下述关系式进行AVAZ反演，得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani以及其它四个参数R₀,C^iso,C^ani1和C^ani2。

R(θ,φ)＝R₀+(B^iso+B^ani cos²φ')sin²θ+(C^iso+C^ani1cos⁴φ'+C^ani2sin²φ'cos²φ')sin²θtan²θ

优选的，所述第一子单元包括：

关系式简化单元，所述关系式简化单元用于将所述反射振幅R与入射角θ和方位角φ之间的关系式简化为：

R(θ,φ)＝AR₀+BB^iso+CB^ani+DC^iso+EC^ani1+FC^ani2，其中A、B、C、D、E和F是入射角和方位角的组合：

A＝1

B＝sin²θ

C＝sin²θcos²φ'

D＝sin²θtan²θ

E＝sin²θtan²θcos⁴φ'

F＝sin²θtan²θsin²φ'cos²φ'

超定方程组构建单元，所述超定方程组构建单元用于将关系式简化单元简化后的关系式构建超定方程组：M_n×6X_6×1＝R_n×1

其中n为覆盖次数， $M={\left[\begin{array}{cccccc}{A}_1& B_1& C_1& D_1& E_1& F_1\\ A_2& B_1& C_2& D_2& E_1& F_2\\ ...& ...& ...& ...& ...& ...\\ A_n& B_1& C_n& D_n& E_n& F_n\end{array}\right]}_{n\×6},X={\left[\begin{array}{c}{R}_0\\ B^{\mathrm{iso}}\\ B^{\mathrm{ani}}\\ C^{\mathrm{iso}}\\ C^{\mathrm{ani}1}\\ C^{\mathrm{ani}2}\end{array}\right]}_{6\×1},R={\left(\begin{array}{c}{R}_1\\ R_2\\ ...\\ R_n\end{array}\right)}_{n\×1}$

超定方程组求解单元，所述超定方程组求解单元用于求解所述超定方程组，得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani。

优选的，所述第二级流体因子获取单元包括流体因子获取单元，所述流体因子获取单元用于利用储层流体的第二拉梅常数μ'以及背景介质的拉梅常数λ,μ对第二级流体因子进行修正，获得流体因子I_fluid，所述背景介质的拉梅常数λ,μ是通过利用所述道集数据中的反射振幅反演计算得到的；

所述流体识别单元用于利用流体因子I_fluid对储层流体进行区分。

进一步的，所述流体识别单元包括第二子单元或第三子单元，所述第二子单元用于将所获得的流体因子I_fluid与预设的阈值进行对比来判断储层流体的类型，所述第三子单元用于将所获得的流体因子I_fluid与现有的流体因子谱库中的数据进行对比来判断储层流体的类型。

与现有技术相比，本申请具有以下有益效果：

通过本申请中的技术方案能减小所构建的流体因子的相对误差，以达到准确判断地下裂缝介质是否充填油气的目的。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本申请中裂缝型储层流体的识别方法的流程示意图。

图2为地震波射线在裂缝介质中的传播示意图。

图3为含裂隙灰岩的方位角道集示意图。

图4为本申请实施例一中裂缝型储层流体的识别方法的流程示意图。

图5为入射角θ、方位角φ以及反射振幅R之间的关系示意图。

图6为sin²θ和sin²θtan²θ与入射角θ的关系曲线。

图7为利用现有技术中的公式(2)与利用本申请中的公式(4)进行反演得到的R₀、各向同性梯度B^iso和各向异性梯度B^ani的相对误差的柱状图。

图8为新疆塔里木盆地哈拉哈塘地区三维地震资料宽方位采集时的偏移距与方位角的关系示意图。

图9为油气勘探目的层流体因子切片示意图。

图10为本申请实施例二中裂缝型储层流体的识别装置的结构示意图。

具体实施方式

为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。

本申请提供一种裂缝型储层流体的识别方法，如图1所示，包括如下步骤：

S1，将获取的整个研究区域的三维地震纵波资料进行各向异性偏移处理，生成道集数据，所述道集数据包括三维地震纵波的方位角、入射角和反射振幅的信息；

S2，利用所生成的道集数据中反射振幅R与入射角θ和方位角φ之间的关系，并考虑大于30度的入射角θ对反射振幅R的影响，进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度B^ani和各向同性梯度B^iso的比值，将该比值记为第一级流体因子

S3，利用裂缝横纵比(即裂缝的横向宽度与纵向高度的比值)d对第一级流体因子进行修正，获得第二级流体因子，该横纵比d是通过运用背景介质裂缝结构的各向异性反演计算所得到的，在此声明，后面所提到的横纵比d均为裂缝的横向宽度与纵向高度的比值d；

S4，利用所获得的第二级流体因子对储层流体进行识别。

与现有的获取储层流体的流体因子的方法相比，通过本申请中的技术方案能减小所构建的流体因子的相对误差，以达到准确判断地下裂缝介质是否充填油气的目的。

下面结合附图，对本发明的具体实施方式作进一步的详细说明。

实施例一

图2为地震波射线在裂缝介质中的传播示意图，图中入射角θ为地震波射线与反射面法向的夹角，方位角φ为地震波的激发点和接收点的连线(即测线)与x轴的夹角。储层内发育裂缝时多对外表现为各向异性特征，其反射系数(即反射振幅)不仅随入射角变化，也随方位角变化，利用入射角θ、方位角φ以及反射振幅R三者的关系可判断地下介质裂缝中是否含有油气，以图3中所示的含裂隙灰岩的方位角道集为例，图3中入射角θ＝20°，横坐标为方位角φ，纵坐标为时间t，(a)为含水，(b)为含油，(c)为含气。图中曲线上黑色部分填充代表了反射振幅的大小，黑色部分填充的越多，代表反射振幅越大，从图中可以看出该方位角道集存在明显的方位各向异性，含气时不但反射振幅最大，其振幅变化程度也比含油或水时大。本申请利用入射角θ、方位角φ以及反射振幅R三者的关系，通过高精度的AVAZ(Amplitude variation with Angle and Azimuth振幅随入射角、方位角变化)反演来获取表征储层流体特性的流体因子，以对裂缝型储层流体进行准确识别。

参考图4所示，本申请实施例所提供裂缝型储层流体的识别方法包括以下步骤：

步骤S101,将获取的三维地震纵波资料进行各向异性偏移处理，生成道集数据。具体的，

对所获取的整个研究区域的三维地震纵波资料进行各向异性偏移处理，生成道集数据，该道集数据中包含有三维地震纵波射线的入射角θ、方位角φ以及反射振幅R的相关信息，入射角θ、方位角φ以及反射振幅R三者之间具有一定的关系，如图5所示。该图为反射振幅的等值线图，横轴表示方位角的大小，纵轴表示入射角的大小，图中曲线上的数值表示反射振幅的大小，当地震波沿着不同的射线传播时，振幅不仅随着入射角变化，还随着方位角而变化。例如，当方位角φ为0°时，随着入射角θ的增大，反射振幅R增大；当方位角φ为90°时，随着入射角θ的增大，反射振幅R减小；对于同一个入射角θ，反射振幅R关于90°的方位角对称。

步骤S102，利用所生成的道集数据中反射振幅R与入射角θ和方位角φ之间的关系，并考虑大于30度的入射角对反射振幅的影响，进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度B^ani和各向同性梯度B^iso的比值，将该比值记为第一级流体因子。具体的，

各向异性裂缝介质模型中含有一系列垂直排列且相互平行的裂缝。对于垂直排列的平行裂缝，可以将其等效为HTI介质(Horizontal Transverse Isotropy，即具有水平对称轴的横向各向同性介质)模型来研究。Ruger认为HTI介质分界面处的纵波反射系数(即反射振幅)是各向异性参数、入射角和方位角的函数：

$\begin{array}{c}R(\mathrm{\θ},\mathrm{\φ})=\frac{1}{2}\frac{\mathrm{\ΔZ}}{\stackrel{\‾}{Z}}+\frac{1}{2}\{\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\ΔG}}{\stackrel{\‾}{G}}+[\mathrm{\Δ}{\mathrm{\δ}}^V+2{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\mathrm{\Δ\γ}\left]{\cos }^2(\mathrm{\φ}-{\mathrm{\φ}}_s)\right\}{\sin }^2\mathrm{\θ}\\ +\frac{1}{2}(\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}+{\mathrm{\Δ\ϵ}}^V{\cos }^4(\mathrm{\φ}-{\mathrm{\φ}}_s)+{\mathrm{\Δ\δ}}^V{\sin }^2(\mathrm{\φ}-{\mathrm{\φ}}_s){\cos }^2(\mathrm{\φ}-{\mathrm{\φ}}_s)){\sin }^2\mathrm{\θ}{\tan }^2\mathrm{\θ}\end{array}---\left(1\right)$

上式中：R为反射振幅，θ为入射角、φ为方位角，φ_s为裂缝平面的对称轴方向，即裂缝面法向，(φ-φ_s)为测线与裂缝面法向的夹角；α,β,Z,G分别是纵波速度、横波速度、纵波阻抗和剪切模量；ε^V,δ^V,γ为各向异性参数。分别表示上下两层介质间纵波速度的平均值、横波速度的平均值、纵波阻抗的平均值和剪切模量的平均值，Δα,ΔZ,ΔG,Δε^V,Δδ^V,Δγ分别表示上下两层介质间纵波速度的差值、纵波阻抗的差值、剪切模量的差值以及三个各向异性参数的差值。

Ruger反射系数公式(1)包含入射角三角函数的二次项sin²θ和四次项sin²θtan²θ，两者取值大小随入射角θ变化，如图6所示。从图中可以看出当θ≤30°时，sin²θtan²θ≤1/12＜＜1，此时可忽略该四次项对反射振幅的影响；而当30°≤θ≤45°时，sin²θtan²θ随着入射角θ的增大而迅速增大，此时不能忽略该四次项对反射振幅的影响；当θ＞45°时，sin²θtan²θ＞sin²θ＞1/2，该四次项对反射振幅的影响大过二次项，更不能忽略。

所以为了提高反演精度，考虑入射角θ大于30度时对反射振幅R的影响，本申请采用保留Ruger公式中四次项sin²θtan²θ的高精度AVAZ算法，且在四次项中引入三个参数C^iso，C^ani1和C^ani2。通过NMOZ(方位动校正速度椭圆)等方法计算裂缝面法向角度φ_s，记测线与裂缝法向的夹角为φ'＝φ-φ_s，公式(1)可简化为：

R(θ,φ)＝R₀+(B^iso+B^anicos²φ')sin²θ

+(C^iso+C^ani1cos⁴φ'+C^ani2sin²φ'cos²φ')sin²θtan²θ     (3)

其中： $R_0=\frac{1}{2}\frac{\mathrm{\ΔZ}}{\stackrel{\‾}{Z}},B^{\mathrm{iso}}=\frac{1}{2}[\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\ΔG}}{\stackrel{\‾}{G}}],B^{\mathrm{ani}}=\frac{1}{2}[\mathrm{\Δ}{\mathrm{\δ}}^V+{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\mathrm{\Δ\γ}],C^{\mathrm{iso}}=\frac{1}{2}\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}},$ $C^{\mathrm{ani}1}=\frac{1}{2}{\mathrm{\Δ\ϵ}}^V,C^{\mathrm{ani}2}=\frac{1}{2}{\mathrm{\Δ\δ}}^V.$

将上述公式(3)简化为：

R(θ,φ)＝AR₀+BB^iso+CB^ani+DC^iso+EC^ani1+FC^ani2     (4)

其中A、B、C、D、E和F是入射角和方位角的组合：

A＝1

B＝sin²θ

C＝sin²θcos²φ'

D＝sin²θtan²θ

E＝sin²θtan²θcos⁴φ'

F＝sin²θtan²θsin²φ'cos²φ'

假设覆盖次数为n，将公式(4)构建成如下超定方程组：

M_n×6X_6×1＝R_n×1       (5)

其中， $M={\left[\begin{array}{cccccc}{A}_1& B_1& C_1& D_1& E_1& F_1\\ A_2& B_1& C_2& D_2& E_1& F_2\\ ...& ...& ...& ...& ...& ...\\ A_n& B_1& C_n& D_n& E_n& F_n\end{array}\right]}_{n\×6},X={\left[\begin{array}{c}{R}_0\\ B^{\mathrm{iso}}\\ B^{\mathrm{ani}}\\ C^{\mathrm{iso}}\\ C^{\mathrm{ani}1}\\ C^{\mathrm{ani}2}\end{array}\right]}_{6\×1},R={\left(\begin{array}{c}{R}_1\\ R_2\\ ...\\ R_n\end{array}\right)}_{n\×1}$

利用奇异值分解法求解方程组(5)，求解结果为：

X＝M⁺R        (6)

上式中M⁺为矩阵M的广义逆，由式(6)即可获得储层流体的各向同性梯度B^iso和各向异性梯度B^ani，同时也可以得到R₀，C^iso，C^ani1和C^ani2这四个参数，其中各向异性梯度B^ani可用来刻画裂缝的发育程度。

图7为利用现有技术中的公式(2)与利用本申请中的公式(4)进行反演得到的R₀、各向同性梯度B^iso和各向异性梯度B^ani的相对误差的柱状图。从图中可以看出利用本申请所提供的公式(4)进行反演所得到的R₀、各向同性梯度B^iso和各向异性梯度B^ani的相对误差均小于0.5％，这远远小于利用现有技术中的公式(2)进行反演所得结果的相对误差。由此可见，利用本申请中的算法进行反演计算所得到的结果的精度相比常规AVAZ反演有大幅度提高。

由于从求解结果式(6)中无法得到各向同性梯度B^iso和各向异性梯度B^ani的表达式，所以利用公式(3)中的各向同性梯度B^iso和各向异性梯度B^ani的表达式并对其进行简化后来计算各向异性梯度B^ani与各向同性梯度B^iso的比值，在此声明一下，虽然利用的是公式(3)中的表达式，但各向同性梯度B^iso和各向异性梯度B^ani是从式(6)中得到的，各向同性梯度B^iso和各向异性梯度B^ani的表达式如下：

$\begin{array}{c}{B}^{\mathrm{iso}}=\frac{1}{2}[\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\ΔG}}{\stackrel{\‾}{G}}]\\ B^{\mathrm{ani}}=\frac{1}{2}[{\mathrm{\Δ\δ}}^V+{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\mathrm{\Δ\γ}]\end{array}---\left(7\right)$

将相关数值代入到公式(7)中并对其进行简化，简化结果如下：

$B^{\mathrm{iso}}=\frac{{(\sqrt{A_5}-\sqrt{A_5-A_{3}e})}^2}{A_{3}e}$

$\begin{array}{c}{B}^{\mathrm{ani}}=\frac{{(A_1+A_{2}e-A_{3}e)}^2-{(A_1-A_{2}e-A_{4}e)}^2}{4(A_5-A_{3}e)(A_1+A_{2}e-A_{3}e)}\\ -\frac{8A_{2}e}{{(\sqrt{A_5}+\sqrt{A_5-A_{3}e})}^2}\end{array}---\left(8\right)$

其中：A₁＝λ+μ，

$A_2=\frac{16\mathrm{\πd}(\mathrm{\λ}+2\mathrm{\μ}){\mathrm{\μ}}^2}{(12(\mathrm{\λ}+2\mathrm{\μ}){\mathrm{\μ}}^{\′}+3\mathrm{\πd\μ}(3\mathrm{\λ}+4\mathrm{\μ}))},$

$A_3=\frac{4\mathrm{\πd}{\mathrm{\λ}}^2(\mathrm{\λ}+2\mathrm{\μ})\mathrm{\μ}}{(12(\mathrm{\λ}+2\mathrm{\μ})\mathrm{\μ}({\mathrm{\λ}}^{\′}+{2\mathrm{\μ}}^{\′})+3\mathrm{\πd}{\mathrm{\μ}}^2(\mathrm{\λ}+\mathrm{\μ}))},$

$A_4=\frac{4\mathrm{\πd\λ}{(\mathrm{\λ}+2\mathrm{\μ})}^2\mathrm{\μ}}{(3(\mathrm{\λ}+2\mathrm{\μ})\mathrm{\μ}({\mathrm{\λ}}^{\′}+2{\mathrm{\μ}}^{\′})+3\mathrm{\πd}{\mathrm{\μ}}^2(\mathrm{\λ}+\mathrm{\μ}))},$

A₅＝λ+2μ。

上式中λ,μ为背景介质的拉梅常数，λ',μ'为储层流体的拉梅常数，d为裂缝横纵比。

在裂缝密度e＝0处对公式(8)进行Taylor级数展开，可得：

$\begin{array}{c}{B}^{\mathrm{iso}}=\frac{A_3}{4A_5}e+O\left(e^2\right)\≈k^{\mathrm{iso}}e\\ B^{\mathrm{ani}}=\frac{6A_2-A_3+A_4}{2A_5}e+O\left(e^2\right)\≈k^{\mathrm{ani}}e\end{array}---\left(9\right)$

从式(9)可以看出各向同性梯度B^iso和各向异性梯度B^ani近似与裂缝密度e成正比。

选取碳酸盐岩为背景介质，其纵横波速度及密度分别为6000m/s、3100m/s、2670kg/m³，模拟不同裂缝密度和充填物(即储层流体，例如油、气、水)。由公式(9)计算AVO梯度，即各向同性梯B^iso和各向异性梯度B^ani，计算所得的各向同性梯度B^iso的最大相对误差小于5％，各向异性梯度B^ani的最大相对误差小于3％。

将背景介质设为碎屑岩，其纵横波速度及密度分别为4000m/s、2100m/s、2530kg/m³，模拟不同裂缝密度和充填物(即储层流体，例如油、气、水)。根据公式(8)计算AVO梯度，计算所得的各向同性梯度B^iso的最大相对误差小于4％，各向异性梯度B^ani的最大相对误差小于3％。

从上述两个实例可以看出本实施例中的各向同/异性AVO梯度线性近似公式(9)适用于多种岩性，且最大相对误差小于5％，由此可见本申请中的上述算法适合于构建不同岩性的裂缝型储层的流体因子。

然后，计算各向异性梯度B^ani与各向同性梯度B^iso的比值，即第一级流体因子，即为：

$I_{\mathrm{fluid}}^{\left(1\right)}=\frac{B^{\mathrm{ani}}}{B^{\mathrm{iso}}}=\frac{2(6A_2-A_3+A_4)}{A_3}---\left(10\right)$

将上述A2，A3和A4的计算公式代入式(10)中化简后可得：

$\begin{array}{c}{I}_{\mathrm{fluid}}^{\left(1\right)}=\frac{48(\mathrm{\λ}+2\mathrm{\μ})\mathrm{\μ}}{\mathrm{\π}{\mathrm{\λ}}^2(3\mathrm{\λ}+4\mathrm{\μ})}({\mathrm{\λ}}^{\′}+\frac{2(\mathrm{\λ}+2\mathrm{\μ})}{(3\mathrm{\λ}+4\mathrm{\μ})}{\mathrm{\μ}}^{\′})\frac{1}{d}\\ +\frac{4\mathrm{\μ}}{\mathrm{\λ}}(\frac{12\mathrm{\μ}(\mathrm{\λ}+\mathrm{\μ})}{\mathrm{\λ}(3\mathrm{\λ}+4\mathrm{\μ})}+1)\end{array}---\left(11\right)$

在此声明，虽然第一级流体因子的计算利用到式(9)中的表达式，但其中各向异性梯度B^ani与各向同性梯度B^iso的具体数值是从式(6)中得到的。

经数值模型实验证实，利用本申请中的算法计算所得的各向异性梯度B^ani与各向同性梯度B^iso的比值约等于常数，不随裂缝密度变化。对于水、油、气这三种流体，相对应的值之比为20:10:1，区分度高，可识别裂缝型储层流体的类型。

步骤S103，利用裂缝横纵比d对第一级流体因子进行修正，获得第二级流体因子然后利用储层流体的第二拉梅常数μ'以及背景介质拉梅常数λ,μ对第二级流体因子进行修正，获得流体因子I_fluid。具体的，

为了使得第一级流体因子更加接近真实值，现利用裂缝横纵比d对公式(11)中的第一级流体因子进行修正，得到不受裂缝横纵比d影响的第二级流体因子，即：

$\begin{array}{c}{I}_{\mathrm{fluid}}^{\left(2\right)}=(I_{\mathrm{fluid}}-g)d\\ =\frac{48(\mathrm{\λ}+2\mathrm{\μ})\mathrm{\μ}}{\mathrm{\π}{\mathrm{\λ}}^2(3\mathrm{\λ}+4\mathrm{\μ})}({\mathrm{\λ}}^{\′}+\frac{2(\mathrm{\λ}+2\mathrm{\μ})}{(3\mathrm{\λ}+4\mathrm{\μ})}{\mathrm{\μ}}^{\′})\end{array}---\left(12\right)$

其中为由背景介质唯一决定的截距，裂缝横纵比d可通过现有的运用背景介质裂缝结构的各向异性反演计算得到，其具体计算方法可参考杨凤英、印兴耀和刘博发表于CPS/SEG 2014北京国际地球物理会议上的文章“基于岩石物理的碳酸盐岩储层参数自适应预测”中有关孔隙纵横比的计算，本申请中所提到的裂缝横纵比d等同于那篇文章中的孔隙纵横比。

可以将通过式(12)计算得到的第二级流体因子作为最终的流体因子，但为了计算更加准确，还需要利用储层流体的第二拉梅常数μ'以及背景介质拉梅常数λ,μ来对此流体因子进行修正，但该修正并非必须。

由于储层流体受压力和温度的影响，储层流体的第二拉梅常数μ'不是理想流体所对应的零值，因此公式(10)中保留了μ'。但实际上μ'远小于λ'，因此在计算中可以忽略第二拉梅常数μ'的影响，在公式(12)中忽略第二拉梅常数μ'的影响后得到第三级流体因子，即：

$I_{\mathrm{fluid}}^{\left(3\right)}=\frac{\mathrm{\π}{\mathrm{\λ}}^2(3\mathrm{\λ}+4\mathrm{\μ})}{48(\mathrm{\λ}+2\mathrm{\μ})\mathrm{\μ}}{\mathrm{\λ}}^{\′}---\left(13\right)$

式(13)表明第三级流体因子近似与储层流体的拉梅常数λ'成正比，并且其受背景介质拉梅常数λ,μ的影响。为了消除背景介质拉梅常数λ,μ对第三级流体因子的影响，使得计算结果更加的准确，先利用背景介质拉梅常数λ,μ对第三级流体因子进行修正，背景介质的拉梅常数λ,μ可采用现有技术中的计算方法，即通过反射振幅进行反演得到背景介质的纵波速度、横波速度以及密度，再利用所获得的背景介质的纵波速度、横波速度以及密度进行计算得到。

利用背景介质拉梅常数λ,μ对第三级流体因子进行修正后得到的流体因子为：

$I_{\mathrm{fluid}}=I_{\mathrm{fluid}}^{\left(3\right)}\frac{\mathrm{\π}{\mathrm{\λ}}^2(3\mathrm{\λ}+4\mathrm{\μ})}{48(\mathrm{\λ}+2\mathrm{\μ})\mathrm{\μ}}\≈{\mathrm{\λ}}^{\′}---\left(14\right)$

步骤S104，利用所获得的流体因子I_fluid对储层流体进行识别。具体的，

从式(14)中可以看出，流体因子I_fluid可近似等于储层流体的第一拉梅常数λ'，所以只需将计算得到的流体因子I_fluid与预设的有关油、气、水的流体因子阈值进行对比，即可判断出整个研究区域中背景介质的裂缝中是否含有油、气、水，也可以将计算所得到的流体因子I_fluid与现有的流体因子谱库中储存的有关油、气、水的第一拉梅常数λ'进行对比来判断储层流体的类别。

以下以该方法在新疆某区块实际资料上的应用为例子具体说明该方法的实际应用效果。

图8是新疆塔里木盆地哈拉哈塘地区三维地震资料宽方位采集时的偏移距与方位角的关系示意图，其中纵轴代表偏移距，横轴代表方位角。从图中可以看出，当偏移距小于3600米时，测线方位角范围为0°～180°，属于全方位覆盖；当偏移距大于6400米时，测线方位角范围为0°～35°、145°～180°，属于局部方位覆盖。对于这两种方位角范围，利用本申请中的裂缝型储层流体的识别方法都能正确识别储层流体的类型。

图9是油气勘探目的层流体因子切片，图中白色圆圈表示井点位置，颜色灰度越深代表含油概率越大。从图中可以看出北部的XK9-3井试油出水，正好落在白色附近的浅灰区，即含水程度较高，与试油结果吻合。综合研究后该井向西侧钻即XK9-3C井，试油获0.2531万吨/天的高产油流。此后钻探的XK9-1、XK9-2、XK902、XK9004等4口井均获得日产超过1万吨的高产油流，这证明了利用本申请中所提供的算法预测流体类型结果的可靠性。

实施例二

图10为本实施例中裂缝型储层流体的识别装置的结构示意图。该识别装置包括道集数据生成单元901、第一级流体因子计算单元902、第二级流体因子获取单元903以及流体识别单元905，其中第二级流体因子获取单元903中包括流体因子获取单元904，道集数据生成单元901用于将获取的整个研究区域的三维地震纵波资料进行各向异性偏移处理，生成道集数据，所生成的道集数据包括三维地震纵波射线的方位角、入射角和反射振幅的信息。

第一级流体因子计算单元902用于利用所生成的道集数据，通过对方位角和入射角进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度B^ani和各向同性梯度B^iso的比值，即第一级流体因子，第一级流体因子计算单元902包括第一子单元，该第一子单元用于利用反射振幅R与入射角θ和方位角φ之间的关系式(3)(具体见实施例一中步骤S102)进行高精度的AVAZ反演来得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，该第一子单元包括：关系式简化单元9021、超定方程组构建单元9022以及超定方程组求解单元9023，其中关系式简化单元9021用于将反射振幅R与入射角θ和方位角φ之间的关系式(3)简化成实施例一中步骤S102中的公式(4)，超定方程组构建单元9022用于将关系式简化单元9021简化后的关系式构建超定方程组，该超定方程组的具体形式见实施例一中步骤S102中的公式(5)，超定方程组求解单元9023用于求解超定方程组(5)，得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani。

第二级流体因子获取单元903用于利用通过运用背景介质裂缝结构的各向异性反演计算所得的裂缝横纵比d对第一级流体因子进行修正，获得第二级流体因子。

流体因子获取单元904用于利用储层流体的第二拉梅常数μ'以及背景介质的拉梅常数λ,μ对第二级流体因子进行修正，获得流体因子I_fluid。

体识别单元905用于利用所获得的流体因子I_fluid对储层流体进行识别，其包括第二子单元或第三子单元，其中第二子单元用于将所获得的流体因子I_fluid与预设的阈值进行对比来判断储层流体的类型，第三子单元用于将所获得的流体因子I_fluid与现有的流体因子谱库中的数据进行对比来判断储层流体的类型。如果该识别装置不包括流体因子获取单元904，那么该流体识别单元905就用于利用第二级流体因子对储层流体进行识别，其第二子单元或第三子单元分别用于将第二级流体因子与预设的阈值进行对比或与现有的流体因子谱库中的数据进行对比来判断储层流体的类型。

Claims (10)

1.一种裂缝型储层流体的识别方法，其特征在于，包括以下步骤：

S1，将获取的整个研究区域的三维地震纵波资料进行各向异性偏移处理，生成道集数据，所述道集数据包括三维地震纵波射线的方位角、入射角和反射振幅的信息；

S2，利用所生成的道集数据中反射振幅R与入射角θ和方位角φ之间的关系，并考虑大于30度的入射角θ对反射振幅R的影响，进行反演得到储层的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度B^ani和各向同性梯度B^iso的比值，将所述比值记为第一级流体因子

S3，利用裂缝横向宽度与纵向高度的比值d对所述第一级流体因子进行修正，获得第二级流体因子所述裂缝横向宽度与纵向高度的比值d是通过运用背景介质裂缝结构的各向异性反演计算所得到的；

S4，利用所获得的第二级流体因子对储层流体类型进行区分。

2.根据权利要求1所述的裂缝型储层流体的识别方法，其特征在于，所述利用所生成的道集数据中反射振幅R与入射角θ和方位角φ之间的关系，并考虑入射角θ大于30度时对反射振幅R的影响，进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，包括：

利用反射振幅R与入射角θ和方位角φ之间的下述关系式进行AVAZ反演，得到储层的各向同性梯度B^iso、各向异性梯度B^ani以及其它四个参数R₀,C^iso,C^ani1和C^ani2：

R(θ,φ)＝R₀+(B^iso+B^anicos²φ')sin²θ+(C^iso+C^ani1cos⁴φ'+C^ani2sin²φ'cos²φ')sin²θtan²θ。

3.根据权利要求2所述的裂缝型储层流体的识别方法，其特征在于，所述AVAZ反演包括：

将所述反射振幅R与入射角θ和方位角φ之间的关系式简化为：

R(θ,φ)＝AR₀+BB^iso+CB^ani+DC^iso+EC^ani1+FC^ani2，其中A、B、C、D、E和F是入射角和方位角的组合：

A＝1

B＝sin²θ

C＝sin²θcos²φ'

D＝sin²θtan²θ

E＝sin²θtan²θcos⁴φ'

F＝sin²θtan²θsin²φ'cos²φ'

将上述公式构建超定方程组：M_n×6X_6×1＝R_n×1

其中n为覆盖次数， $M={\left[\begin{array}{cccccc}{A}_1& B_1& C_1& D_1& E_1& F_1\\ A_2& B_1& C_2& D_2& E_1& F_2\\ ...& ...& ...& ...& ...& ...\\ A_n& B_1& C_n& D_n& E_n& F_n\end{array}\right]}_{n\×6},$ $X={\left[\begin{array}{c}{R}_0\\ B^{\mathrm{iso}}\\ B^{\mathrm{ani}}\\ C^{\mathrm{iso}}\\ C^{\mathrm{ani}1}\\ C^{\mathrm{ani}2}\end{array}\right]}_{6\×1},$ $R={\left(\begin{array}{c}{R}_1\\ R_2\\ ...\\ R_n\end{array}\right)}_{n\×1}$

通过求解所述超定方程组，得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani。

4.根据权利要求1-3任一项所述的裂缝型储层流体的识别方法，其特征在于，所述步骤S2还包括：在获得所述第二级流体因子后，利用储层流体的第二拉梅常数μ'以及背景介质的拉梅常数λ,μ对第二级流体因子进行修正，获得流体因子I_fluid，所述背景介质的拉梅常数λ,μ是通过利用所述道集数据中的反射振幅信息反演计算得到的；

所述步骤S3为利用流体因子I_fluid对储层流体进行区分。

5.根据权利要求4所述的裂缝型储层流体的识别方法，其特征在于，所述利用流体因子I_fluid对储层流体进行区分包括：

将所获得的流体因子I_fluid与预设的阈值进行对比来判断储层流体的类型，

或，

将所获得的流体因子I_fluid与现有的流体因子谱库中的数据进行对比来判断储层流体的类型。

6.一种裂缝型储层流体的识别装置，其特征在于，包括：

道集数据生成单元，所述道集数据生成单元用于将获取的整个研究区域的三维地震纵波资料进行各向异性偏移处理，生成道集数据，所述道集数据包括三维地震纵波射线的方位角、入射角和反射振幅的信息；

第一级流体因子计算单元，所述第一级流体因子计算单元用于利用所生成的道集数据，通过对方位角和入射角进行反演得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani，并计算各向异性梯度B^ani和各向同性梯度B^iso的比值，将所述比值记为第一级流体因子

第二级流体因子获取单元，所述第二级流体因子获取单元用于利用裂缝横向宽度与纵向高度的比值d对所述第一级流体因子进行修正，获得第二级流体因子裂缝横向宽度与纵向高度的比值d是通过运用背景介质裂缝结构的各向异性反演计算得到的；

流体识别单元，所述流体识别单元用于利用所获得的第二级流体因子对储层流体进行识别。

7.根据权利要求6所述的裂缝型储层流体的识别装置，其特征在于，所述第一级流体因子计算单元包括第一子单元，所述第一子单元用于利用反射振幅R与入射角θ和方位角φ之间的下述关系式进行AVAZ反演，得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani以及其它四个参数R₀,C^iso,C^ani1和C^ani2：

R(θ,φ)＝R₀+(B^iso+B^anicos²φ')sin²θ+(C^iso+C^ani1cos⁴φ'+C^ani2sin²φ'cos²φ')sin²θtan²θ。

8.根据权利要求7所述的裂缝型储层流体的识别装置，其特征在于，所述第一子单元包括：

关系式简化单元，所述关系式简化单元用于将所述反射振幅R与入射角θ和方位角φ之间的关系式简化为：

R(θ,φ)＝AR₀+BB^iso+CB^ani+DC^iso+EC^ani1+FC^ani2，其中A、B、C、D、E和F是入射角和方位角的组合：

A＝1

B＝sin²θ

C＝sin²θcos²φ'

D＝sin²θtan²θ

E＝sin²θtan²θcos⁴φ'

F＝sin²θtan²θsin²φ'cos²φ'

超定方程组构建单元，所述超定方程组构建单元用于将关系式简化单元简化后的关系式构建超定方程组：M_n×6X_6×1＝R_n×1

其中n为覆盖次数， $M={\left[\begin{array}{cccccc}{A}_1& B_1& C_1& D_1& E_1& F_1\\ A_2& B_1& C_2& D_2& E_1& F_2\\ ...& ...& ...& ...& ...& ...\\ A_n& B_1& C_n& D_n& E_n& F_n\end{array}\right]}_{n\×6},$ $X={\left[\begin{array}{c}{R}_0\\ B^{\mathrm{iso}}\\ B^{\mathrm{ani}}\\ C^{\mathrm{iso}}\\ C^{\mathrm{ani}1}\\ C^{\mathrm{ani}2}\end{array}\right]}_{6\×1},$ $R={\left(\begin{array}{c}{R}_1\\ R_2\\ ...\\ R_n\end{array}\right)}_{n\×1}$

超定方程组求解单元，所述超定方程组求解单元用于求解所述超定方程组，得到储层流体的各向同性梯度B^iso和各向异性梯度B^ani。

9.根据权利要求6-8任一项所述的裂缝型储层流体的识别装置，其特征在于，所述第二级流体因子获取单元包括流体因子获取单元，所述流体因子获取单元用于利用储层流体的第二拉梅常数μ'以及背景介质的拉梅常数λ,μ对第二级流体因子进行修正，获得流体因子I_fluid，所述背景介质的拉梅常数λ,μ是通过利用所述道集数据中的反射振幅反演计算得到的；

所述流体识别单元用于利用流体因子I_fluid对储层流体进行区分。

10.根据权利要求9所述的裂缝型储层流体的识别装置，其特征在于，所述流体识别单元包括第二子单元或第三子单元，所述第二子单元用于将所获得的流体因子I_fluid与预设的阈值进行对比来判断储层流体的类型，所述第三子单元用于将所获得的流体因子I_fluid与现有的流体因子谱库中的数据进行对比来判断储层流体的类型。