CN106556861A - 一种基于全方位地震资料的方位avo反演方法 - Google Patents

Abstract

本发明提供了一种基于全方位地震资料的方位AVO反演方法，包括：基于方位角度道集形式的全方位地震资料，计算不同方位角及入射角下的反射系数；将所述反射系数与实测的方位角度道集建立目标函数；推导目标函数的敏感度矩阵；基于目标函数和敏感度矩阵，利用最优化算法求解方位AVO反演。实现了从全方位角度道集反演纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向，为全方位地震资料的应用及预测裂缝分布提供重要资料。

Description

一种基于全方位地震资料的方位AVO反演方法

技术领域

本发明涉及石油地球物理勘探技术领域，具体地涉及基于全方位地震资料的方位AVO反演。

背景技术

裂缝分布和发育信息对石油和天然气的勘探及开发有重要指导作用，特别是碳酸盐岩储层的裂缝、孔隙和孔洞是油气的主要富集部位和运移通道，是碳酸盐岩油气成藏的必要条件之一。为了有效的识别裂缝，地球物理工作者提出了各种裂缝检测方法，比如基于叠后地震资料的应力场分析与不连续性检测、基于叠前地震资料的AVO与各向异性分析。在叠后资料应用方面，基于应力场分析的裂缝预测适用于构造引起的裂缝，基于不连续检测的各种裂缝预测方法，如相干分析、广义希尔伯特变换和基于图像边缘检测等，能够有效的圈定不规律的缝洞发育带与均匀介质的地震波场之间的界限，但受噪音影响较大。随着宽方位(或者全方位)地震勘探的应用，如何从全方位地震资料中挖掘更加丰富的地下介质信息成为一个重要的研究热点。

由于裂缝性储层常存在明显的HTI各向异性的影响，主要表现为振幅、速度、反射波形及相位随测线方位变化而变化，因此在宽方位叠前资料应用中可以利用纵波信息随方位的变化来预测裂缝。比如利用方位旅行时、方位衰减、方位速度(VVAZ，速度随方位角变化)和方位振幅变化(AVAZ，振幅随方位角变化)等来预测裂缝。在利用纵波方位各向异性特征预测裂缝方面，我国学者对此开展了较多的应用研究，如程冰洁等基于宽方位转换波地震资料利用转换波AVAZ裂缝检测方法对川西深层裂缝性气藏进行了裂缝检测；孙炜对用于纵波方位各向异性裂缝预测的地震属性优选进行了研究，并利用成像测井及岩心资料上的裂缝信息验证了预测结果的可靠性和准确性。目前的方法是分别对不同方位角下的地震数据进行属性提取或者地震反演，然后对不同方位角度下的属性进行各向异性椭圆拟合，由椭圆的长、短轴的比值和长短轴的方向来表征裂缝发育强度和方向，从而预测储层裂缝的分布特征。其中用于椭圆拟合的地震属性可以是AVO属性，也可是波阻抗等。这类基于椭圆拟合的方法受到椭圆拟合算法精度的影响很大，并且抗噪性低。

发明内容

本发明直接从全方位地震资料出发，以全方位地震资料为基础，综合利用各个方位角下的地震资料，利用P波方位各向异性分析原理进行方位AVO(AVAZ)反演，从全方位地震资料中反演纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向，为裂缝型储层预测提供重要分析技术手段。反演得到的各向异性梯度及裂缝方向可直接表征裂缝发育强度和方向。

本发明提供了一种基于全方位地震资料的方位AVO反演方法，包括：基于方位角度道集形式的全方位地震资料，计算不同方位角及入射角下的反射系数；将所述反射系数与实测的方位角度道集建立目标函数；推导目标函数的敏感度矩阵；基于目标函数和敏感度矩阵，利用最优化算法求解方位AVO反演。

进一步地，由P波方位各向异性分析原理中的Corrigan-Ruger近似方程推导出方位各向异性线性近似方程(1)，计算不同方位角及入射角下的反射系数：

所述公式(1)中，θ是入射角、φ是方位角、表示方位角为φ入射角为θ时的反射系数、A、B、C、D是Corrigan-Ruger近似方程中的临时变量，四个临时变量与纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向的关系为：

所述公式(2)中，A_iso是纵波阻抗反射系数、B_iso是各项同性梯度、B_ani是各向异性梯度、φ_sym是HTI介质对称轴方向。

进一步地，上述公式(1)为线性的，对于不同方位角和入射角的情况简写为矩阵形式：

d＝Gm (3)

其中，

所述公式(4)中，m是由临时变量构成的向量，G是由方位角和入射角构成的正演矩阵，d是不同方位角和入射角下的振幅值，R(θ_i,φ_j)为方位角为φ_j入射角为θ_i时的振幅值，M为方位角个数，N为入射角个数。

进一步地，所述将所述反射系数与实测的方位角度道集建立目标函数的步骤，利用公式(3)正演计算出不同方位角和入射角下的反射系数，并与相应的方位角和入射角下的观测数据求取误差来构建反演问题的目标函数：

式中，d_s为不同方位角下的观测数据，即由实测全方位地震资料经过处理得到的方位角度道集，R(φ_i,θ_j)是方位角为φ_j，入射角为θ_i时的振幅值。

进一步地，所述推导目标函数的敏感度矩阵的步骤，包括对目标函数J(m)求取偏导数：

所述公式(6)中，为目标函数J(m)对临时变量A、B、C、D的偏导数，由此组成目标函数的敏感度矩阵即

进一步地，利用最优化算法中的共轭梯度法求解方位AVO反演。利用最优化算法中的共轭梯度法迭代求解，直到满足终止条件为止。

优选地，所述终止条件是目标函数误差达到预定值或目标函数敏感度矩阵范数达到预定值。

优选地，目标函数误差的预定值为0.1，目标函数敏感度矩阵范数的预定值为0.1。

本发明利用P波方位各向异性分析原理进行AVAZ反演，从全方位地震资料中反演纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向，反演得到的各向异性梯度及裂缝方向可直接表征裂缝发育强度和方向，为裂缝型储层预测提供重要分析技术手段，提高了裂缝型储层预测的效率和准确性。

本发明的方法是基于P波方位各向异性分析原理中的Corrigan-Ruger近似方程，综合的利用全方位地震资料，从全方位地震资料中反演得到反映裂缝发育强度和方向的数据，当某一方位较其他方位数据存在较大误差时算法也能够得到较为稳定的结果，抗噪能力强。

附图说明

通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。

图1是本发明方法的流程图。

图2是实施例1中的方位角度道集，作为输入数据给定。

图3a-d是实施例1中方位AVO(AVAZ)反演结果，依次为P波反射系数、各向同性梯度、各向异性梯度、方向。

图4是实施例1中方位AVO(AVAZ)反演结果的沿层切片，A_iso为P波反射系数、B_iso为各向同性梯度、B_ani为各向异性梯度。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

本发明公开一种基于全方位地震资料的方位AVO反演方法，属于石油地球物理勘探技术领域。首先利用Corrigan-Ruger近似方程计算不同方位角及入射角下的反射系数，继而与方位角度道集建立目标函数，然后推导了目标函数的敏感度矩阵，最后利用最优化算法求解方位AVO反演(AVAZ)反演问题。实现了从全方位角度道集反演纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向，为全方位地震资料的应用及预测裂缝分布提供重要资料。

具体地，如图1所示，本发明提供了一种基于全方位地震资料的方位AVO反演方法，包括：基于方位角度道集形式的全方位地震资料，计算不同方位角及入射角下的反射系数；将所述反射系数与实测的方位角度道集建立目标函数；推导目标函数的敏感度矩阵；基于目标函数和敏感度矩阵，利用最优化算法求解方位AVO反演。

首先，基于方位角度道集形式的全方位地震资料，计算不同方位角及入射角下的反射系数。具体地，可由P波方位各向异性分析原理中的Corrigan-Ruger近似方程推导出方位各向异性线性近似方程，方位各向异性线性近似方程具体为公式(1)所示：

所述公式(1)中，方位各向异性线性近似方程能够计算出不同方位角及入射角下的反射系数，θ是入射角、φ是方位角、表示方位角为φ入射角为θ时的反射系数、A、B、C、D是Corrigan-Ruger近似方程中的临时变量，这四个临时变量与纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向的关系为：

所述公式(2)中，A_iso是纵波阻抗反射系数、B_iso是各项同性梯度、B_ani是各向异性梯度、φ_sym是HTI介质对称轴方向。

优选地，上述公式(1)为线性的，对于不同方位角和入射角的情况可简写为矩阵形式：

d＝Gm (3)

其中，

所述公式(4)中，m是由临时变量构成的向量，G是由方位角和入射角构成的正演矩阵，d是不同方位角和入射角下的振幅值，R(θ_i,φ_j)为方位角为φ_j入射角为θ_i时的振幅值，方位角个数为M，入射角个数为N。

接下来，将所述反射系数与实测的方位角度道集建立目标函数。具体地，利用公式(3)正演计算出不同方位角和入射角下的反射系数，并与相应的方位角和入射角下的观测数据求取误差来构建反演问题的目标函数，具体由公式(5)：

式中，d_s为不同方位角下的观测数据，即由实测全方位地震资料经过处理得到的方位角度道集，m是由临时变量A、B、C、D构成的向量，即m＝(A,B,C,D)，G是由方位角和入射角构成的正演矩阵，具体形式见上部步骤中公式(4)。R(φ_i,θ_j)是方位角为φ_j，入射角为θ_i时的振幅值。

接下来，推导目标函数的敏感度矩阵。对反演问题的目标函数J(m)求取偏导数，详见公式(6)：

所述公式(6)中，为目标函数J(m)对临时变量A、B、C、D的偏导数，由此组成目标函数的敏感度矩阵即

最后，基于目标函数和敏感度矩阵，利用最优化算法求解方位AVO反演。可利用最优化算法中的共轭梯度法求解。具体地，基于目标函数和敏感度矩阵，利用最优化算法中的共轭梯度法迭代求解，直到满足终止条件为止，终止条件是目标函数误差达到较小的值(例如0.1)或目标函数敏感度矩阵范数达到较小的值(例如0.1)。

为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。

下面以某地区的实际地震资料为例，利用本发明方法从全方位地震资料中反演纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向，进而说明本发明的效果。

图2是实施例1中的方位角度道集，作为输入数据给定。图中黑色方框内为一个方位角度道集，由若干角度道集组成，每个角度道集所对应的方位角不同。图中方框的方位角道集由5个方位角(分别为0、36、72、108、144度)，每个方位角对应一组角度道集(入射角度道集，图中入射角为0度到45度)。

图3a-d是实施例1中方位AVO(AVAZ)反演结果，以图2中方位角度道集为输入数据，进行方位AVO(AVAZ)反演，得到图3a中P波反射系数、图3b中各向同性梯度、图3c中各向异性梯度、图3d中方向。

图4是实施例1中方位AVO(AVAZ)反演结果的沿层切片，A_iso为P波反射系数、B_iso为各向同性梯度、B_ani为各向异性梯度。

实施例1

本实施例是本发明用于从全方位地震资料中反演纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向的实施实例，一种基于全方位地震资料的方位AVO反演方法的步骤如下：

本发明中以方位角度道集的形式给定全方位地震资料，给定临时变量(A、B、C、D)的初始模型，然后利用方位各向异性线性近似方程计算不同方位角及入射角下的反射系数；

将所述反射系数与实测的方位角度道集建立目标函数；

推导目标函数的敏感度矩阵；

基于目标函数和敏感度矩阵，利用最优化算法中的共轭梯度法迭代求解上述步骤，直到满足终止条件为止，终止条件是目标函数误差达到较小的值(如0.1)或目标函数敏感度矩阵范数达到较小的值(如0.1)。

本发明直接从全方位地震资料出发，以全方位地震资料为基础，综合利用各个方位角下的地震资料，利用P波方位各向异性分析原理进行AVAZ反演，从全方位地震资料中反演纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向，为裂缝型储层预测提供重要分析技术手段。反演得到的各向异性梯度及裂缝方向可直接表征裂缝发育强度和方向。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

Claims (10)

1.一种基于全方位地震资料的方位AVO反演方法，包括：

基于方位角度道集形式的全方位地震资料，计算不同方位角及入射角下的反射系数；

将所述反射系数与实测的方位角度道集建立目标函数；

推导目标函数的敏感度矩阵；

基于目标函数和敏感度矩阵，利用最优化算法求解方位AVO反演。

2.根据权利要求1所述的方位AVO反演方法，其中由P波方位各向异性分析原理中的Corrigan-Ruger近似方程推导出方位各向异性线性近似方程(1)，计算不同方位角及入射角下的反射系数：

$R_P^{HTI}(\mathrm{\φ},\mathrm{\θ})=A+\[B+Ccos\left(2\mathrm{\φ}\right)+Dsin\left(2\mathrm{\φ}\right)\]{\sin }^2\mathrm{\θ}---\left(1\right)$

所述公式(1)中，θ是入射角、φ是方位角、表示方位角为φ入射角为θ时的反射系数、A、B、C、D是Corrigan-Ruger近似方程中的临时变量，四个临时变量与纵波反射系数、各向同性梯度、各向异性梯度及裂缝发育方向的关系为：

$\left\{\begin{array}{c}{A}_{iso}=A\\ B_{iso}=B-\sqrt{C^2+D^2}\\ B_{ani}=2\sqrt{C^2+D^2}\\ {\mathrm{\φ}}_{sym}=arctan(D/C)/2\end{array}\right.;---\left(2\right)$

所述公式(2)中，A_iso是纵波阻抗反射系数、B_iso是各项同性梯度、B_ani是各向异性梯度、φ_sym是HTI介质对称轴方向。

3.根据权利要求2所述的方位AVO反演方法，其中上述公式(1)为线性的，对于不同方位角和入射角的情况简写为矩阵形式：

d＝Gm (3)

其中，

$d=\left[\begin{array}{c}R\left({\mathrm{\θ}}_1,{\mathrm{\φ}}_1\right)\\ R\left({\mathrm{\θ}}_2,{\mathrm{\φ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ R\left({\mathrm{\θ}}_i,{\mathrm{\φ}}_j\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ R\left({\mathrm{\θ}}_N,{\mathrm{\φ}}_M\right)\end{array}\right],G=\left[\begin{array}{cccc}1& {\sin }^2{\mathrm{\θ}}_1& \cos 2{\mathrm{\φ}}_1{\sin }^2{\mathrm{\θ}}_1& \sin 2{\mathrm{\φ}}_1{\sin }^2{\mathrm{\θ}}_1\\ 1& {\sin }^2{\mathrm{\θ}}_2& \cos 2{\mathrm{\φ}}_1{\sin }^2{\mathrm{\θ}}_2& \sin 2{\mathrm{\φ}}_1{\sin }^2{\mathrm{\θ}}_2\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 1& {\sin }^2{\mathrm{\θ}}_i& \cos 2{\mathrm{\φ}}_j{\sin }^2{\mathrm{\θ}}_i& \sin 2{\mathrm{\φ}}_j{\sin }^2{\mathrm{\θ}}_i\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 1& {\sin }^2{\mathrm{\θ}}_N& \cos 2{\mathrm{\φ}}_M{\sin }^2{\mathrm{\θ}}_N& \sin 2{\mathrm{\φ}}_M{\sin }^2{\mathrm{\θ}}_N\end{array}\right],m=\left[\begin{array}{c}A\\ B\\ C\\ D\end{array}\right]---\left(4\right)$

所述公式(4)中，m是由临时变量构成的向量，G是由方位角和入射角构成的正演矩阵，d是不同方位角和入射角下的振幅值，R(θ_i,φ_j)为方位角为φ_j入射角为θ_i时的振幅值，M为方位角个数，N为入射角个数。

4.根据权利要求3所述的方位AVO反演方法，其中所述将所述反射系数与实测的方位角度道集建立目标函数的步骤，利用公式(3)正演计算出不同方位角和入射角下的反射系数，并与相应的方位角和入射角下的观测数据求取误差来构建反演问题的目标函数：

$J\left(m\right)=\frac{1}{2}{(Gm-d_s)}^T(Gm-d_s)=\frac{1}{2}\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(R\right({\mathrm{\φ}}_i,{\mathrm{\θ}}_j)-{\mathrm{ds}}_{ij})---\left(5\right)$

式中，d_s为不同方位角下的观测数据，即由实测全方位地震资料经过处理得到的方位角度道集，R(φ_i,θ_j)是方位角为φ_j，入射角为θ_i时的振幅值。

5.根据权利要求4所述的方位AVO反演方法，其中所述推导目标函数的敏感度矩阵的步骤，包括对目标函数J(m)求取偏导数：

$\left\{\begin{array}{c}\∂J\left(m\right)/\∂A=\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(R\right({\mathrm{\φ}}_i,{\mathrm{\θ}}_j)-{\mathrm{ds}}_{ij})\\ \∂J\left(m\right)/\∂B=\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(R\right({\mathrm{\φ}}_i,{\mathrm{\θ}}_j)-{\mathrm{ds}}_{ij}){\sin }^2{\mathrm{\θ}}_j\\ \∂J\left(m\right)/\∂C=\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(R\left({\mathrm{\φ}}_i,{\mathrm{\θ}}_j\right)-{\mathrm{ds}}_{ij}\right)\cos \left(2{\mathrm{\φ}}_i\right){\sin }^2{\mathrm{\θ}}_j\\ \∂J\left(m\right)/\∂D=\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(R\right({\mathrm{\φ}}_i,{\mathrm{\θ}}_j)-{\mathrm{ds}}_{ij})sin\left(2{\mathrm{\φ}}_i\right){\sin }^2{\mathrm{\θ}}_j\end{array}\right.---\left(6\right)$

所述公式(6)中，为目标函数J(m)对临时变量A、B、C、D的偏导数，由此组成目标函数的敏感度矩阵即 $\▿J\left(m\right)=(\∂J(m)/\∂A,\∂J(m)/\∂B,\∂J(m)/\∂C,\∂J(m)/\∂D).$

6.根据权利要求5所述的方位AVO反演方法，其中利用最优化算法中的共轭梯度法求解方位AVO反演。

7.根据权利要求6所述的方位AVO反演方法，其中利用最优化算法中的共轭梯度法迭代求解，直到满足终止条件为止。

8.根据权利要求7所述的方位AVO反演方法，其中所述终止条件是目标函数误差达到预定值或目标函数敏感度矩阵范数达到预定值。

9.根据权利要求8所述的方位AVO反演方法，其中目标函数误差的预定值为0.1，目标函数敏感度矩阵范数的预定值为0.1。

10.根据权利要求1所述的方位AVO反演方法，其中反演得到的各向异性梯度及裂缝方向表征裂缝发育强度和方向。