CN104777513B - 地震数据梯度信息不连续性边界检测方法 - Google Patents

Abstract

本发明涉及一种地震数据处理方法，具体为地震数据梯度信息不连续性边界检测方法，输入三维地震数据，计算梯度向量体，建立二维滑动分析时窗，二维滑动分析时窗分为5个子窗口，计算滑动时窗内每个样点的梯度幅值；再计算子窗口的峰态值，比较每个子窗口的峰态值；将峰态值最大的那个子窗口内所有的梯度向量构成梯度张量矩阵，求取梯度结构张量矩阵的特征值，并将三个特征值进行排序，利用排序后的特征值大小构建梯度张量属性，将求取的梯度张量属性赋给位于二维滑动分析时窗内中心的分析点；返回计算下一个分析点的梯度张量属性值。该方法既考虑了时窗内位于不同区域位置数据，又考虑了时窗内数据本身的不稳定性，充分利用了各种信息。

Description

地震数据梯度信息不连续性边界检测方法

技术领域

本发明涉及一种地震数据处理方法，具体为地震数据梯度信息不连续性边界检测方法。

背景技术

地震数据中含有大量丰富的地层和岩性信息，这些信息很多都是储层中断层、河道、裂缝发育区以及其他一些不连续性结构特征的响应，但是利用常规的地震数据解释方法并不能够得到清晰直观的不连续性信息。

2002年，Cohen和Coifman提出了基于局部结构熵的相干性算法。该算法把由地震数据估计的局部结构熵作为相干测定，先构造一个分析数据体，并将它分为4个子数据体，并分别生成4个列向量，利用4个列向量的互相关生成一个4×4的相关矩阵，将该矩阵的归一化道作为局部熵估计。该算法从计算上讲是非常有效的，因为它避免了协方差矩阵中大量本征值的计算，不足之处是没有考虑地下构造的倾角对局部结构熵估计的影响。

2005年，Lu Wenkai和Li Yandong等提出了基于高阶累积量的相干估计方法(HOSC)。高阶累积量(HOS)成功用于信号处理，而相干分析的基本问题就是延迟估计，这正是HOS方法的特点。HOSC方法同时利用3个地震道来计算拥有零延迟相关的归一化4阶矩的二维切片，将二维切片上最大相关点作为相干估计。另外，Lu Wenkai提出了超级地震道(ST-supertrace)概念。先定义一个分析时窗，它包含分析点周围的许多地震道，然后重新排列多道为超道，最后在整个三维数据体中移动分析时窗，这样就把原始地震数据体转换为超级地震道数据体。ST数据体保持了原始数据体中所拥有的地层倾角信息，故在ST数据体中利用倾角扫描很容易获取最佳倾角估计。把ST与HOSC算法组合起来可以得到有效的相干估计算法ST-HOSC。该算法在检测不连续性信息方面比C1相干体和HOS相干体更好。

以上两种现有的检测不连续性边缘的算法，都是在一个整体时窗内所有的地震数据进行分析处理，并没有考虑整体时窗内数据本身的分布情况和结构特征，而且没有考虑时窗内不同位置区域的数据对中心分析点的贡献。

发明内容

针对上述技术问题，本发明提供一种新的不连续性边界的检测方法。该不连续性边界检测算法从地震数据梯度信息的角度出发，在构建滑动时窗时既考虑时窗内位于不同区域位置数据对中心分析点的贡献程度，又考虑时窗内数据本身的不稳定性。

具体的技术方案是：

地震数据梯度信息不连续性边界检测方法，包括以下步骤：

(1)输入三维地震数据，并利用式(1)～(3)计算地震数据在x，y，z三个方向上的梯度向量体；

上式中，f(x，y，z)为三维地震数据；Δx、Δy和Δz分别为地震数据在x，y，z三个方向上的面元大小；g_x、g_y和g_z分别为地震数据在x、y、z三个方向上的梯度向量体；

(2)建立二维滑动分析时窗，二维滑动分析时窗分为五个3×3的子窗口，分别位于二维滑动分析时窗的中心位置、左上位置、右上位置、左下位置和右下位置；

(3)利用式(4)计算滑动时窗内每个样点的梯度幅值；再利用式(5)计算二维滑动分析时窗内5个子窗口的峰态值：

上式中，g为计算的梯度幅值；g_i为子窗口内第i个样点的梯度幅值大小；为子窗口内所有梯度幅值的平均值；K和J分别为子窗口的峰态值和样点数；

(4)比较二维滑动分析时窗内每个子窗口的峰态值；

(5)将峰态值最大的那个子窗口内所有的梯度向量构成梯度张量矩阵，即式(6)；

上式中g_x、g_y和g_z分别为子窗口内的地震数据在x、y、z三个方向上的梯度向量；J为子窗口内的样点个数；

(6)利用QR分解法求取梯度结构张量矩阵的特征值，并将三个特征值进行排序，λ₁＞λ₂＞λ₃；

(7)利用排序后的特征值大小构建梯度张量属性，即式(7)；

上式中λ₁为最大特征值；N为特征值的个数；为所有特征值的平均值；

(8)将求取的梯度张量属性赋给位于二维滑动分析时窗内中心的分析点；返回步骤(3)，滑动二维滑动分析时窗并计算下一个分析点的梯度张量属性值；

(9)将整个三维地震数据依据以上过程计算完毕，得到最后的三维梯度张量属性体。

本发明提供的地震数据梯度信息不连续性边界检测方法，将一个二维滑动分析时窗分为五个子窗口进行分析处理，可以避免单一地对一个整体分析窗口进行处理，保证了时窗内位于不同区域位置数据对中心分析点的贡献，具有很好的边缘贡献效应。峰态值的大小反应了一组数的不连续性程度，取峰态值最大的那个子窗口内的数据来进行分析，可以最大程度地量化空间的不连续性度，这样可以更好地刻画地震数据中的不连续性边界。

本发明提供的地震数据梯度信息不连续性边界检测方法，构建滑动时窗时既考虑了时窗内位于不同区域位置数据对中心分析点的贡献程度，又考虑了时窗内数据本身的不稳定性，充分利用了时窗内地震数据的各种信息。

附图说明

图1是本发明实施例的x方向样本；

图2是本发明实施例的y方向样本；

图3是本发明实施例的z方向样本；

图4是本发明实施例的二维滑动时窗之一；

图5是本发明实施例的二维滑动时窗之二；

图6是本发明实施例的二维滑动时窗之三；

图7是本发明实施例的二维滑动时窗之四；

图8是本发明实施例的二二维滑动时窗之五。

具体实施方式

结合实施例说明本发明的具体实施过程。

实例处理的数据为x方向、y方向和z方向的三组样本，如图1、图2和图3。

一种基于最大峰态约束下梯度张量属性的不连续性边界检测方法，依次包括以下步骤：

(1)确定需要处理的三维地震数据体。利用三维地震数据在x方向、y方向和z方向的差值，分别计算出这三个方向的梯度向量体。

x方向：

y方向：

z方向：

(2)在每个梯度向量体中建立二维滑动时窗。二维滑动时窗分为5个小子窗口，如图4、图5、图6、图7和图8所示，分别计算每个子窗口中所有样点的梯度幅值以及子窗口的峰态值。

梯度幅值：

峰态值：

(3)比较二维滑动时窗中每个子窗口的峰态值大小。将峰态值最大的子窗口内所有的梯度向量值作为梯度张量矩阵的元素。

梯度张量矩阵：

(4)利用QR分解法求取梯度张量矩阵的三个特征值，并将特征值进行排序：

λ₁＝5＞λ₂＝-1＝λ₃＝-1

(5)排序后的特征值构成二维滑动时窗内中心点的梯度张量属性：

(6)沿着纵向方向滑动二维分析时窗，完成纵向上所有地震数据的梯度张量属性计算。沿着横向方向滑动二维分析时窗，完成横向上所有地震数据的梯度张量属性计算。

(7)完成一个二维地震数据的梯度张量属性计算后，接着对下一个二维地震数据采用相同的方法逐步二维滑动分析时窗完成梯度张量属性计算。

(8)完成整个三维地震数据体的梯度张量属性计算。

Claims (1)

1.地震数据梯度信息不连续性边界检测方法，其特征在于：包括以下步骤：

(1)输入三维地震数据，并利用式(1)～(3)计算地震数据在x，y，z三个方向上的梯度向量体；

$g_x=\frac{f(x_{i+1},y,z)-f(x_i,y,z)}{\mathrm{\Δ}x}---\left(1\right);$

$g_y=\frac{f(x,y_{j+1},z)-f(x,y_j,z)}{\mathrm{\Δ}y}---\left(2\right);$

$g_z=\frac{f(x,y,z_{k+1})-f(x,y,z_k)}{\mathrm{\Δ}z}---\left(3\right);$

上式中，f(x，y，z)为三维地震数据；Δx、Δy和Δz分别为地震数据在x，y，z三个方向上的面元大小；g_x、g_y和g_z分别为地震数据在x、y、z三个方向上的梯度向量体；

(2)建立二维滑动分析时窗，二维滑动分析时窗分为五个3×3的子窗口，分别位于二维滑动分析时窗的中心位置、左上位置、右上位置、左下位置和右下位置；

(3)利用式(4)计算二维滑动分析时窗内每个样点的梯度幅值；再利用式(5)计算二维滑动分析时窗内5个子窗口的峰态值：

$g=\sqrt{g_x^2+g_y^2+g_z^2}---\left(4\right);$

$K=\frac{1}{n}\underset{m=1}{\overset{n}{\Σ}}{(g_m-\stackrel{\‾}{g})}^4---\left(5\right);$

上式中，g为计算的梯度幅值；g_m·为子窗口内第m个样点的梯度幅值大小；为子窗口内所有梯度幅值的平均值；K和n分别为子窗口的峰态值和样点数；

(4)比较二维滑动分析时窗内每个子窗口的峰态值；

(5)将峰态值最大的那个子窗口内所有的梯度向量构成梯度张量矩阵，即式(6)；

$\stackrel{\‾}{T}=\stackrel{\‾}{g\·g^T}=\frac{1}{n}\left\{\begin{array}{ccc}\underset{m=1}{\overset{n}{\Σ}}{g}_x^2& \underset{m=1}{\overset{n}{\Σ}}{g}_x{g}_y& \underset{m=1}{\overset{n}{\Σ}}{g}_x{g}_z\\ \underset{m=1}{\overset{n}{\Σ}}{g}_y{g}_x& \underset{m=1}{\overset{n}{\Σ}}{g}_y^2& \underset{m=1}{\overset{n}{\Σ}}{g}_y{g}_z\\ \underset{m=1}{\overset{n}{\Σ}}{g}_z{g}_x& \underset{m=1}{\overset{n}{\Σ}}{g}_z{g}_y& \underset{m=1}{\overset{n}{\Σ}}{g}_z^2\end{array}\right\}---\left(6\right)$

上式中g_x、g_y和g_z分别为地震数据在x、y、z三个方向上的梯度向量；n为子窗口内的样点个数；

(6)利用QR分解法求取梯度结构张量矩阵的特征值，并将三个特征值进行排序，λ₁＞λ₂＞λ₃；

(7)利用排序后的特征值大小构建梯度张量属性，即式(7)；

$A_G=\frac{{\mathrm{\λ}}_1}{\frac{1}{q}\underset{p=1}{\overset{q}{\Σ}}{({\mathrm{\λ}}_i-\stackrel{\‾}{\mathrm{\λ}})}^3}---\left(7\right)$

上式中λ₁为最大特征值；q为特征值的个数；为所有特征值的平均值；

(8)将求取的梯度张量属性赋给位于二维滑动分析时窗内中心的分析点；返回步骤(3)，滑动二维滑动分析时窗并计算下一个分析点的梯度张量属性值；

(9)将整个三维地震数据依据以上过程计算完毕，得到最后的三维梯度张量属性体。