CN104199091A - 基于方位地震数据的纹理属性方法 - Google Patents

Abstract

本发明的基于方位地震数据的纹理属性方法，首先对宽方位角三维地震数据体进行梯度计算，得到每个方位角数据体的三个方向的梯度数据体；其次，根据不同方位角数据体所反应的地层下的地质结构的差异性，利用上一步得到的各个方位角数据体的三个方向的梯度数据体来计算目标点的相关性矩阵；最后对相关性矩阵进行特征值分解，并由得到的特征值计算基于宽方位角数据体的纹理属性体。

Description

基于方位地震数据的纹理属性方法

技术领域

本发明属于地质领域，具体涉及一种方位纹理属性。

背景技术

纹理是按一定的规则或者一定的形式变化而产生的图案。在自然界中，许多植物、动物、矿物都有独特的纹理特征。利用这些纹理特征，人们可以方便的识别和区分事物。地层由于受到构造运动的影响，也会产生断层、裂缝等地质现象，从而形成地底下的纹理。在进行地质勘探中，我们可以得到一系列的由地震数据形成的图像，例如地震数据剖面图。对于这些地震图像，我们可以借助图像处理的方法、手段来计算三维地震数据的纹理属性，以便凸显地底下的内部特征，从而有效地识别断层或裂缝，为找油、找气提供更好的依据。

在如今的地质领域中，纹理属性的提取都是针对常规的叠后数据来处理的，它们以叠后地震数据体为基础输入，能够通过纹理属性的计算来预测地层中尺度较大、数量较少的断裂，但对较小尺度、数量更庞大的小型或微断裂系统难于预测，其根本原因在于叠后地震资料的信息量比较小，缺乏偏移距信息和方位角信息，这样就失去了以各向异性来检测微断裂的理论依据和基础。

发明内容

针对上述问题，本发明提出一种基于方位地震数据的纹理属性方法。

本发明的技术方案为：基于方位地震数据的纹理属性方法，具体包括以下步骤：

S1：计算宽方位角地震数据的梯度；

S2：根据步骤S1得到的梯度，获得各个方位角的三个方向的三维梯度数据体，在得到三维梯度数据体的基础上，获得各个方位角的三维地震数据体在x、y和z方向上的相关性矩阵T；

$T=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right];$

S3：对相关性矩阵T进行特征值分解：Tv＝λv；

其中，λ表示对角矩阵，对角线上的值为特征值，分别为λ₁，λ₂，λ₃，它们是由大到小排列的，即λ_i≥λ_i+1，v表示特征向量矩阵，列向量表示所述特征值λ₁，λ₂，λ₃对应的特征向量，分别为v₁，v₂，v₃；

S4：计算基于宽方位角GST的方位纹理属性体；

所述宽方位角的纹理属性包括：

(1)宽方位角连续属性C_con：

$C_{\mathrm{con}}=\frac{{\mathrm{\λ}}_1}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_2+{\mathrm{\λ}}_3};$

宽方位角连续性属性C_con的取值范围为0到1；

(2)宽方位角平面属性C_plane：

$C_{\mathrm{plane}}=\frac{{\mathrm{\λ}}_1-{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_2};$

(3)宽方位角线条属性C_line：

$C_{\mathrm{line}}=\frac{{\mathrm{\λ}}_2-{\mathrm{\λ}}_3}{{\mathrm{\λ}}_2+{\mathrm{\λ}}_3};$

所述宽方位角线条属性C_line的取值范围为0到1；

(4)宽方位角断裂系统属性C_fault：

$C_{\mathrm{fault}}=\frac{2{\mathrm{\λ}}_2({\mathrm{\λ}}_2-{\mathrm{\λ}}_3)}{({\mathrm{\λ}}_1+{\mathrm{\λ}}_2)({\mathrm{\λ}}_2+{\mathrm{\λ}}_3)};$

所述宽方位角断裂系统属性C_fault的取值范围为0到1；

(5)宽方位角混沌纹理属性C_chaos：

$C_{\mathrm{chaos}}=\frac{2{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_3}-1;$

(6)宽方位角倾角属性dip：

$\mathrm{dip}=a\tan (\sqrt{v_1{\left(1\right)}^2+v_1{\left(2\right)}^2}/v_1\left(3\right));$

其中，atan(·)表示求反正切值，v₁(1)表示特征向量v₁的第一个值，v₁(2)表示特征向量v₁的第二个值，v₁(3)表示特征向量v₁的第二个值；

(7)宽方位角方位属性azimuth：

azimuth＝atan2(v₁(2),v₁(1))；

(8)宽方位角相干属性coherence：

$\mathrm{coherence}=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_s(i,j)}{3({\mathrm{\λ}}_1+{\mathrm{\λ}}_2+{\mathrm{\λ}}_3)}.$

本发明从宽方位角地震资料出发，以宽方位角三维地震资料的波场特征为基础，利用振幅随方位角变化的特征，充分利用宽方位角三维地震资料的充足信息量(如偏移距信息和方位角信息等)来预测尺度较小、数量更庞大的小型或者微断裂系统，从而得到更加稳定、鲁棒性更好的方位纹理属性体。

附图说明

图1是基于宽方位地震数据的纹理属性体的计算流程图。

图2是宽方位地震数据的梯度计算示意图。

图3是宽方位地震数据的梯度计算流程图。

图4是宽方位地震数据的相关性矩阵的计算流程图。

图5是宽方位角数据的Coherence属性。

图6是常规叠后数据的Coherence属性。

图7是常规叠后数据的Dip属性。

图8是宽方位角数据的Dip属性。

图9是宽方位角数据的Dip属性与地质专家提供的溶孔、溶洞分布预测叠合图。

图10是产气井与断层关系图。

图11是常规叠后数据的Continuety属性。

图12是宽方位角数据的Continuety属性。

图13是宽方位角数据的Continuety属性与构造解释断层叠合图。

具体实施方式

总流程图如图1所示，具体方案如下：

S1：首先需要计算输入宽方位角地震数据的梯度，其目的是需要用梯度信息来反映随着方向的不同而检测到的振幅的变化，进而来获得宽方位角的地震纹理属性。其具体思路如下：对每个宽方位角的三维地震数据体，分别在x、y和z方向进行求导运算，从而获得三个方向新的数据体。具体运算采用地震数据与三维差分高斯滤波器作卷积的方式来实现，其目的是降低噪声对求导过程的影响。经典的三维高斯滤波器的表达式为：

$G(x_j,y_j,z_j;\mathrm{\σ})=\exp [(x_j^2+y_j^2+z_j^2)/\left(2{\mathrm{\σ}}^2\right)]---(2-2-1)$

其中，x_j，y_j和z_j分别表示沿着x，y和z轴的以待求导的分析点中心，σ为方差或者这里称为尺度因子。这里用S_θ(x,y,z)表示各个方位角的三维原始振幅数据体，其中θ∈{0°，30°，60°，90°，120°,150°}。每个角度的三个方向的梯度数据体的计算表达式如下：

$G_{(x,\mathrm{\θ})}=S_{\mathrm{\θ}}(x,y,z)*(\∂G/\∂x)---(2-2-2)$

$G_{(y,\mathrm{\θ})}=S_{\mathrm{\θ}}(x,y,z)*(\∂G/\∂y)---(2-2-3)$

$G_{(z,\mathrm{\θ})}=S_{\mathrm{\θ}}(x,y,z)*(\∂G/\∂z)---(2-2-4)$

其中，*为卷积运算符号，和分别为x、y和z方向的三维差分高斯滤波器，G_(x,θ)、G_(y,θ)和G_(z,θ)即是所需要求得的x，y和z方向的三维梯度数据体，它们可以反映在x，y和z方向上原始的三维方位角振幅数据体中振幅值的变化情况。其具体示意图如图2所示：

计算宽方位角地震数据的梯度的具体流程图如图3所示：

S2：本发明在获得各个方位角的三个方向的三维梯度数据体后，紧接着在梯度结构的基础上，获得各个方位角的三维地震数据体在x、y和z方向上的相关性矩阵，该相关性矩阵即可反映出各个方位角数据体在x、y和z方向的差异，以便通过方位角数据体的纹理属性提的计算来预测地层下的地质结构。相关性矩阵是是一个3*3的对称矩阵，其定义式表示如下：

$T=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right]---(2-2-5)$

T代表在计算单点时所得到的相关性矩阵。它的矩阵因子具体表达式如公式2-2-6，其中θ∈{0°，30°，60°，90°，120°,150°}。T_ij是通过在时间上在各个角度的三个方向梯度数据体内取窗口长度为N_t方式来计算获得，分别代表各个方向上的相关系数，反映的是各个方位角的三维地震数据在各个方向上的相似程度。在计算目标点的相关性矩阵的过程中，既利用了单个数据体的三个方向的数据信息，又利用了不同方位角数据体的信息，成功的利用了整个四维数据体给出的所有信息，以便后面计算出的方位角的纹理属性体更加的准确、稳定。

$\left\{\begin{array}{c}{T}_{11}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(x,\mathrm{\θ})}(x,y,t)\×G_{(x,\mathrm{\θ})}(x,y,t)\\ T_{12}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(x,\mathrm{\θ})}(x,y,z)\×G_{(y,\mathrm{\θ})}(x,y,t)\\ T_{13}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(x,\mathrm{\θ})}(x,y,z)\×G_{(z,\mathrm{\θ})}(x,y,z)\\ T_{22}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(y,\mathrm{\θ})}(x,y,z)\×G_{(y,\mathrm{\θ})}(x,y,z)\\ T_{23}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(y,\mathrm{\θ})}(x,y,z)\×G_{(z,\mathrm{\θ})}(x,y,z)\\ T_{33}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(z,\mathrm{\θ})}(x,y,z)\×G_{(z,\mathrm{\θ})}(x,y,z)\end{array}\right.---(2-2-6)$

求相关矩阵的具体流程图如图4所示：

S3：对2-2-5中的相关性矩阵T如下的特征值分解：

Tv＝λv       (2-2-7)

其中，λ表示对角矩阵，对角线上的值为特征值，分别为λ₁，λ₂，λ₃，它们是由大到小排列的，即λ_i≥λ_i+1。v表示特征向量矩阵，列向量表示以上特征值λ₁，λ₂，λ₃对应的特征向量，分别为v₁，v₂，v₃。

对相关性矩阵做特征值分解的目的是得到宽方位角资料中从各个角度的三维地震数据体中所反应的目标点周围的地质结构。

当λ₁＝λ₂＝λ₃＝0时，说明从各个角度的三维地震数据中测得的目标点周围地质结构为没有差异，振幅强度一致。

当λ₁>λ₂＝λ₃＝0时，说明从各个角度的三维地震数据中测得的目标点周围地质结构为平面线性结构，在某一个方向上有变化。

当λ₁>λ₂>λ₃＝0时，说明从各个角度的三维地震数据中测得的目标点周围地质结构为线性结构。但是当λ₁＝λ₂时，说明该线性结构是各向同性的。

当λ₁>λ₂>λ₃>0时，说明从各个角度的三维地震数据中测得的目标点周围地质结构为非线性的。但是当λ₁＝λ₂＝λ₃时，说明在目标点周围的结构为各向同性的。

S4：计算基于宽方位角GST的方位纹理属性体。宽方位角的纹理属性有很多，下面给出主要的几种：

(1)宽方位角连续属性

$C_{\mathrm{con}}=\frac{{\mathrm{\λ}}_1}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_2+{\mathrm{\λ}}_3}---(2-2-8)$

宽方位角连续性属性的取值范围为0到1，值越大，说明在目标点周围各个角度的数据体所获得的地震数据体的连续性越好，而值越小的地方则可能是断层、裂缝所在的地方。

(2)宽方位角平面属性

$C_{\mathrm{plane}}=\frac{{\mathrm{\λ}}_1-{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_2}---(2-2-9)$

(3)宽方位角线条属性

$C_{\mathrm{line}}=\frac{{\mathrm{\λ}}_2-{\mathrm{\λ}}_3}{{\mathrm{\λ}}_2+{\mathrm{\λ}}_3}---(2-2-10)$

宽方位角线条属性的取值范围为0到1，值越大，说明在目标点周围从各个角度的数据体中所获得的地震数据体的线性越强，值越小，则说明目标点周围并不呈现出线条状，则可能与周围的点相比变化不大。

(4)宽方位角断裂系统属性

$C_{\mathrm{fault}}=\frac{2{\mathrm{\λ}}_2({\mathrm{\λ}}_2-{\mathrm{\λ}}_3)}{({\mathrm{\λ}}_1+{\mathrm{\λ}}_2)({\mathrm{\λ}}_2+{\mathrm{\λ}}_3)}---(2-2-11)$

宽方位角断裂系统属性的取值范围为0到1，值越大，说明在目标点周围从各个角度的数据体中所获得的地震数据体的断裂性越强，断层和裂缝越明显，值越小，则目标点周围的地质结构为连续结构。

(5)宽方位角混沌纹理属性

$C_{\mathrm{chaos}}=\frac{2{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_3}-1---(2-2-12)$

(6)宽方位角倾角属性

$\mathrm{dip}=a\tan (\sqrt{v_1{\left(1\right)}^2+v_1{\left(2\right)}^2}/v_1\left(3\right))---(2-2-13)$

(7)宽方位角方位属性

azimuth＝atan2(v₁(2),v₁(1))       (2-2-14)

(8)宽方位角相干属性

$\mathrm{coherence}=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_s(i,j)}{3({\mathrm{\λ}}_1+{\mathrm{\λ}}_2+{\mathrm{\λ}}_3)}---(2-2-15)$

通过以上几种主要的宽方位角纹理属性体的计算，地层下不同的地质结构可以通过不同的属性体的计算来预测。例如溶孔、溶洞等地质构造可以在宽方位角的纹理属性中的Dip属性体中得到体现，如图9。而在宽方位角的纹理属性中的Continuety属性体中则能清晰明显的看到断层等地质结构，如图13。

下面是基于宽方位角的三维地震数据体的纹理属性的提取在四川某工区的实际应用，并与该工区的传统叠后数据体的地震纹理属性体做了对比，从结果可以看到，与传统叠后数据的地震纹理属性相比，方位地震纹理属性对裂缝和断层的预测更加明显清晰。

图9是由地质专家所绘制的纹理属性中的Dip属性图与溶孔、溶洞分布预测叠合图，可以看到，我们所计算的基于宽方位角地震资料的Dip属性与专家所预测的溶孔、溶洞位置几乎完全吻合，断层也非常清晰明显。

图10是由地质专家所绘制的产气井与断层的关系示意图，其中的线条代表断层，是由基于宽方位角地震资料的纹理属性中的Dip属性得到的。

从上述实施例可知，本发明提供的基于方位地震数据的纹理属性方法，从宽方位角地震资料出发，以宽方位角三维地震资料的波场特征为基础，利用振幅随方位角变化的特征，充分利用宽方位角三维地震资料的充足信息量(如偏移距信息和方位角信息等)来预测尺度较小、数量更庞大的小型或者微断裂系统，从而得到更加稳定、鲁棒性更好的方位纹理属性体。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (5)

1.基于方位地震数据的纹理属性方法，其特征在于，包括：

S1：计算宽方位角地震数据的梯度；

S2：根据步骤S1得到的梯度，获得各个方位角的三个方向的三维梯度数据体，在得到三维梯度数据体的基础上，获得各个方位角的三维地震数据体在x、y和z方向上的相关性矩阵T；

$T=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right];$

S3：对相关性矩阵T进行特征值分解：Tv＝λv；

其中，λ表示对角矩阵，对角线上的值为特征值，分别为λ₁，λ₂，λ₃，它们是由大到小排列的，即λ_i≥λ_i+1，v表示特征向量矩阵，列向量表示所述特征值λ₁，λ₂，λ₃对应的特征向量，分别为v₁，v₂，v₃；

S4：计算基于宽方位角GST的方位纹理属性体；

所述宽方位角的纹理属性包括：

(1)宽方位角连续属性C_con：

$C_{\mathrm{con}}=\frac{{\mathrm{\λ}}_1}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_2+{\mathrm{\λ}}_3};$

宽方位角连续性属性C_con的取值范围为0到1；

(2)宽方位角平面属性C_plane：

$C_{\mathrm{plane}}=\frac{{\mathrm{\λ}}_1-{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_2};$

(3)宽方位角线条属性C_line：

$C_{\mathrm{line}}=\frac{{\mathrm{\λ}}_2-{\mathrm{\λ}}_3}{{\mathrm{\λ}}_2+{\mathrm{\λ}}_3};$

所述宽方位角线条属性C_line的取值范围为0到1；

(4)宽方位角断裂系统属性C_fault：

$C_{\mathrm{fault}}=\frac{2{\mathrm{\λ}}_2({\mathrm{\λ}}_2-{\mathrm{\λ}}_3)}{({\mathrm{\λ}}_1+{\mathrm{\λ}}_2)({\mathrm{\λ}}_2+{\mathrm{\λ}}_3)};$

所述宽方位角断裂系统属性C_fault的取值范围为0到1；

(5)宽方位角混沌纹理属性C_chaos：

$C_{\mathrm{chaos}}=\frac{2{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1+{\mathrm{\λ}}_3}-1;$

(6)宽方位角倾角属性dip：

$\mathrm{dip}=a\tan (\sqrt{v_1{\left(1\right)}^2+v_1{\left(2\right)}^2}/v_1\left(3\right));$

其中，atan(·)表示求反正切值，v₁(1)表示特征向量v₁的第一个值，v₁(2)表示特征向量v₁的第二个值，v₁(3)表示特征向量v₁的第二个值；

(7)宽方位角方位属性azimuth：

azimuth＝atan2(v₁(2),v₁(1))；

(8)宽方位角相干属性coherence：

$\mathrm{coherence}=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_s(i,j)}{3({\mathrm{\λ}}_1+{\mathrm{\λ}}_2+{\mathrm{\λ}}_3)}.$

2.根据权利要求1基于方位地震数据的纹理属性方法，其特征在于，所述S1：计算宽方位角地震数据的梯度，包括：对每个宽方位角的三维地震数据体，分别在x、y和z方向进行求导运算，获得三个方向新的数据体，具体运算公式为：

$G(x_j,y_j,z_j;\mathrm{\σ})=\exp [(x_j^2+y_j^2+z_j^2)/\left(2{\mathrm{\σ}}^2\right)];$

其中，x_j，y_j和z_j分别表示沿着x，y和z轴的以待求导的分析点中心，σ为方差即尺度因子；

每个角度的三个方向的梯度数据体的计算表达式如下：

$G_{(x,\mathrm{\θ})}=S_{\mathrm{\θ}}(x,y,z)*(\∂G/\∂x);$

$G_{(y,\mathrm{\θ})}=S_{\mathrm{\θ}}(x,y,z)*(\∂G/\∂y);$

$G_{(z,\mathrm{\θ})}=S_{\mathrm{\θ}}(x,y,z)*(\∂G/\∂z);$

其中，S_θ(x,y,z)表示各个方位角的三维原始振幅数据体，θ∈{0°，30°，60°，90°，120°,150°}，*为卷积运算符号，和分别为x、y和z方向的三维差分高斯滤波器，G_(x,θ)、G_(y,θ)和G_(z,θ)即是所需要求得的x，y和z方向的三维梯度数据体。

3.根据权利要求1基于方位地震数据的纹理属性方法，其特征在于，所述相关性矩阵T它的矩阵因子具体表达式为，

$\left\{\begin{array}{c}{T}_{11}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(x,\mathrm{\θ})}(x,y,t)\×G_{(x,\mathrm{\θ})}(x,y,t)\\ T_{12}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(x,\mathrm{\θ})}(x,y,z)\×G_{(y,\mathrm{\θ})}(x,y,t)\\ T_{13}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(x,\mathrm{\θ})}(x,y,z)\×G_{(z,\mathrm{\θ})}(x,y,z)\\ T_{22}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(y,\mathrm{\θ})}(x,y,z)\×G_{(y,\mathrm{\θ})}(x,y,z)\\ T_{23}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(y,\mathrm{\θ})}(x,y,z)\×G_{(z,\mathrm{\θ})}(x,y,z)\\ T_{33}=\underset{t=z-N_t}{\overset{t=z+N_t}{\mathrm{\Σ}}}{G}_{(z,\mathrm{\θ})}(x,y,z)\×G_{(z,\mathrm{\θ})}(x,y,z)\end{array}\right.;$

其中，θ∈{0°，30°，60°，90°，120°,150°}，G_(x,θ)(x,y,t)、G_(y,θ)(x,y,t)、G_(z,θ)(x,y,t)表示梯度数据体，T_ij分别代表各个方向上的相关系数。

4.根据权利要求4基于方位地震数据的纹理属性方法，其特征在于，所述T_ij通过在时间上在各个角度的三个方向梯度数据体内取窗口长度为N_t方式计算获得。

5.根据权利要求1基于方位地震数据的纹理属性方法，其特征在于，所述特征值λ₁，λ₂，λ₃的大小关系，反应的目标点周围的地质结构；

当λ₁＝λ₂＝λ₃＝0时，表示目标点周围地质结构为没有差异，振幅强度一致；

当λ₁>λ₂＝λ₃＝0时，表示目标点周围地质结构为平面线性结构，在一个方向上有变化；

当λ₁>λ₂>λ₃＝0时，表示目标点周围地质结构为线性结构，同时当λ₁＝λ₂>λ₃＝0时，表示该线性结构是各向同性的；

当λ₁>λ₂>λ₃>0时，表示目标点周围地质结构为非线性的，同时当λ₁＝λ₂＝λ₃>0时，表示在目标点周围的结构为各向同性的。