CN103675902A - 一种最优方向边缘监测方法 - Google Patents

Abstract

本发明提供了一种最优方向边缘监测方法，属于石油地球物理勘探领域。所述方法利用结构导向平滑滤波器结合结构导向边缘检测滤波提取三维地震数据中的地质体边界信息。与第一代相干算法相比，本发明具有自适应寻找方位信息的优点，另外，本发明使用的扩散滤波方程具有压制噪音的功能，对资料的适应程度更高。

Description

一种最优方向边缘监测方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种最优方向边缘监测方法。

背景技术

检测地震数据中地质体的边界信息是地球物理油藏检测描述的一项重要工作，从上世纪70年代开始，相干技术在油藏描述中得到了广泛应用，三维相干技术经历了三代算法的发展，第一代算法以互相关为基础，其优点是计算比较简单，缺点是仅适用于高质量的地震资料；第二代算法是以协方差矩阵为基础，通过沿协方差矩阵中各个检测倾角/方位计算多道相干，其优点是抗噪能力强，缺点是分辨率较低；第三代算法是计算协方差矩阵中特征值，其优点是抗噪能力强，其缺点是不适应陡倾角地区。

基于扩散方程的噪声衰减方法是上世纪90年代初兴起的一种数字图像增强技术，由于其特有的边缘保持性能和结构特征的倾向去噪，很快在多个领域得到广泛应用。Fhemers和Hocker(2003)首先将扩散滤波技术引入地震资料的处理解释中，提出构造约束、边缘保护的各向异性扩散滤波技术，用于增强地震数据的结构特征使其易于解释。国内学者孙夕平、杨长春也对非线性各向异性扩散滤波技术做了大量研究。

但是，目前边缘检测算法都采用倾角扫描来约束边缘检测，不能适应陡倾角。而扩散方程在方位信息约束下可以提取该方位的特征信息，并且能够较好地适应陡倾角。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种最优方向边缘监测方法，将基于扩散方程的结构平滑滤波器与相似系数计算公式相结合，即结构导向边缘检测滤波，检测三维地震数据中地质体边缘信息(例如河道、断层、断裂)，利用结构张量的方位信息及各向异性扩散的方向滤波特性，自适应地寻找边界信息，同时具有压制背景噪音的特点。

本发明是通过以下技术方案实现的：

一种最优方向边缘监测方法，利用结构导向平滑滤波器结合结构导向边缘检测滤波提取三维地震数据中的地质体边界信息。

所述方法包括以下步骤：

(1)求取三维地震数据中的每个样点的梯度 $\mathrm{grad}=\left[\begin{array}{ccc}\frac{\mathrm{df}}{\mathrm{dx}}& \frac{\mathrm{df}}{\mathrm{dy}}& \frac{\mathrm{df}}{\mathrm{dt}}\end{array}\right],$ 利用奇异值分解计算所述三维地震数据的结构张量T(x)＝grad^Tgrad＝λ_uuu^T+λ_vvv^T+λ_www^T；其中，f表示原始三维地震数据，x，y，t分别表示所述三维地震数据的三个方向，λ_u、λ_v和λ_w分别表示结构张量T的奇异值，T是T(x)中的一个元素，u、v和w分别为特征向量；规定特征值降序排列，λ_u≥λ_v≥λ_w≥0，其中特征向量u代表梯度最大的方向，垂直于地质体平面特征；特征向量w代表梯度最小的方向，平行于线性地质体；特征向量v和w都平行于地质体平面；

(2)利用结构导向平滑滤波器计算结构导向边缘检测滤波的分子：

所述结构导向平滑滤波器为(2)式：

$g\left(x\right)-\mathrm{\α}\▿\·D\left(x\right)\▿g\left(x\right)=f\left(x\right)---\left(2\right)$

其中f(x)为输入的三维地震数据，g(x)为平滑输出结果；常数变量α为时间步长，α＝0时，f(x)＝g(x)，D(x)与T(x)具有相同特征向量，不同特征值；

所述结构导向边缘检测滤波为(14)式：

$s_{\mathrm{vw},u}=\frac{{<{<f>}_{\mathrm{vw}}^2>}_u}{{<{<f^2>}_{\mathrm{vw}}>}_u}---\left(14\right)$

(14)式中，s_vw，u为计算得到的边缘系数；

(3)利用结构导向平滑滤波器计算结构导向边缘检测滤波的分母；

(4)计算步骤(2)得到的分子与步骤(3)得到的分母的比值，该比值即为地震数据中的地质体边界信息，也就是边缘检测结果。

所述步骤(2)是这样实现的：

首先对三维地震数据利用(2)式在

时进行计算得到中间结果，然后对中间结果进行平方，再次对平方后的中间结果利用(2)式在时进行计算得到输出结果，该输出结果即为结构导向边缘检测滤波的分子。

所述步骤(3)是这样实现的：

首先对三维地震数据进行平方，然后对平方后的三维地震数据利用(2)式在

时进行计算得到中间结果，再对中间结果进行平方，再次对平方后的中间结果利用(2)式在

时进行计算得到输出结果，该输出结果即为结构导向边缘检测滤波的分母。

与现有技术相比，本发明的有益效果是：

(1)与第一代相干算法相比，本发明具有自适应寻找方位信息的优点；

(2)本发明使用的扩散滤波方程具有压制噪音的功能，对资料的适应程度更高。

附图说明

图1-1是现有共轭梯度算法的步骤图。

图1-2是本发明对共轭梯度算法的改进示意图。

图2是利用第一代相干算法处理的结果。

图3是利用本发明方法处理的结果。

图4本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明方法用到的原理如下：

结构导向平滑滤波器及边缘检测滤波器需要地震数据中地质体的方位信息，利用结构张量的奇异值分解就可以获得地质体的各向异性方位信息。

对于三维地震数据，结构张量T(方向梯度的外积)为3×3正定-半正定对称矩阵，其奇异值分解为T＝λ_uuu^T+λ_vvv^T+λ_www^T

式中λ_u、λ_v和λ_w分别表示结构张量的奇异值，u、v和w分别为特征向量。规定特征值降序排列，λ_u≥λ_v≥λ_w≥0，其中特征向量u代表梯度最大的方向，垂直于地质体平面特征。特征向量w代表梯度最小的方向，平行于线性地质体。特征向量v和w都平行于地质体平面。

(1)结构导向平滑滤波器

利用结构张量场的特征向量可设计各向异性平滑滤波器，这种平滑滤波器可以沿着面特征对地震数据进行平滑，避免了常规加权平滑算子破坏了地质体的边缘特征。

Weick和Fehmers设计了结构导向平滑滤波器。设输入图像为f(x)，输入结构张量场为T(x)，则各向异性扩散方程的解g(x)即为平滑结果：

$\frac{\&PartialD;g(x;\mathrm{\&tau;})}{\&PartialD;\mathrm{\&tau;}}=\&dtri;\&CenterDot;D\left(x\right)\&dtri;g(x;\mathrm{\&tau;})---\left(1\right)$

解的初始值g(x；τ＝0)＝f(x)。对τ＞0，g(x；τ)为输入f(x)的一平滑解。D(x)与T(x)具有相同的特征向量。如果要沿特征向量v(x)进行平滑，则需令λ_u(x)＝0和λ_v(x)＝1。Dave(2009)提出了一种新颖的数值解法：

$g\left(x\right)-\mathrm{\&alpha;}\&dtri;\&CenterDot;D\left(x\right)\&dtri;g\left(x\right)=f\left(x\right)---\left(2\right)$

其中f(x)为输入地震数据，g(x)为平滑输出结果。常数变量α相当一时间步长，α＝0时，f(x)＝g(x)。与常规扩散方程相比，新的方程并不依赖时间τ，并且能自适应地寻找最优平滑方向。

本发明中(2)式采用中心差分格式，具体格式如下：

$\stackrel{\&OverBar;}{x}\left(a11\stackrel{\&OverBar;}{x}u\right)=\frac{a11_{\frac{1}{2}i,j,k}{u}_{1i,j,k}{u}_{i,j,k}{a11}_{\frac{1}{2}i,j,k}{u}_{i,j,k}{u}_{1i,j,k}}{{\mathrm{dx}}^2}---\left(3\right)$

$\stackrel{\&OverBar;}{y}\left(a22\stackrel{\&OverBar;}{y}u\right)=\frac{{a22}_{i,\frac{1}{2}j,k}{u}_{i,1j,k}{u}_{i,j,k}{a22}_{i,\frac{1}{2}j,k}{u}_{i,j,k}{u}_{i,1j,k}}{{\mathrm{dy}}^2}---\left(4\right)$

$\stackrel{\&OverBar;}{z}a33\stackrel{\&OverBar;}{z}u=\frac{{a33}_{i,j,\frac{1}{2}k}{u}_{i,j,1k}{u}_{i,j,k}{a33}_{i,j,\frac{1}{2}k}{u}_{i,j,k}{u}_{i,j,1k}}{{\mathrm{dz}}^2}---\left(5\right)$

$\stackrel{\&OverBar;}{x}a12\stackrel{\&OverBar;}{y}u=\frac{{a12}_{i1,j,k}{u}_{i1,j1,k}{u}_{i1,j1,k}{a12}_{i1,j,k}{u}_{i1,j1,k}{u}_{i1,j1,k}}{4\mathrm{dxdy}}---\left(6\right)$

$\stackrel{\&OverBar;}{x}a13\stackrel{\&OverBar;}{z}u=\frac{a13_{i1,j,k}{u}_{i1,j,k1}{u}_{i1,j,k1}{a13}_{i1,j,k}{u}_{i1,j,k1}{u}_{i1,j1,k1}}{4\mathrm{dxdz}}---\left(7\right)$

$\stackrel{\&OverBar;}{y}a21\stackrel{\&OverBar;}{x}u=\frac{a{21}_{i,j1,k}{u}_{i1,j1,k}{u}_{i1,j,k}{a21}_{i,j1,k}{u}_{i1,j1,k}{u}_{i1,j1,k}}{4\mathrm{dxdy}}---\left(8\right)$

$\stackrel{\&OverBar;}{y}a23\stackrel{\&OverBar;}{z}u=\frac{{a23}_{i,j1,k}{u}_{i,j1,k1}{u}_{i,j,k1}{a23}_{i,j1,k}{u}_{i,j1,k1}{u}_{i,j1,k1}}{4\mathrm{dydz}}---\left(9\right)$

$\stackrel{\&OverBar;}{z}a31\stackrel{\&OverBar;}{x}u=\frac{{a31}_{i,j,k1}{u}_{i1,j,k1}{u}_{i1,j,k1}{a31}_{i,j,k1}{u}_{i1,j,k1}{u}_{i1,j,k1}}{4\mathrm{dxdz}}---\left(10\right)$

$\stackrel{\&OverBar;}{z}a32\stackrel{\&OverBar;}{y}u=\frac{{a32}_{i,j,k1}{u}_{i,j1,k1}{u}_{i,j1,k1}{a32}_{i,j,k1}{u}_{i,j1,k1}{u}_{i,j1,k1}}{4\mathrm{dxdz}}---\left(11\right)$

其中a_mn为3×3矩阵D(x)的元素，i，j，k分别为三维地震数据每个采样点x、y、t方向的采样序号，dx，dy，dz分别为采样间隔(均等于1)。

(2)各向异性扩散方程的数值实现

求解线性平滑方程(2)需要求解一大型稀疏方程。从数值求解的观点来看，新方程等价于一个单一的大时间步长，并且在数值上是无条件稳定的。求解大型稀疏矩阵有许多方法，直接求解算法如LU分解，ILU分解，SVD，迭代算法如LSQR，共轭梯度(CG)、预条件共轭梯度(preconditioned-CG)、伪牛顿算法、多重网格算法等。直接求解算法具有计算效率高的优点，但是所需存储内存较大。在整个数值算法实现过程中，遇到的最大困难就内存。为减少计算过程对内存的需求，因此本发明对共轭梯度算法进行了改进，图1-1是原始的共轭梯度算法，图1-2是本发明修改后的共轭梯度算法(该算法是用来求解结构导向滤波计算公式(14)的)。对于共轭梯度算法而言，初始解与真解越接近，收敛速度越快，而本发明的平滑真解与原始地震剖面存在相似性，因此采用原始地震剖面数据作为初始解，采用梯度的模作为收敛条件，如果两次迭代过程中收敛条件比值小于一定值就退出迭代，这样就避免了迭代次数的设置，大大减少了计算工作量。实际资料表明只需迭代2～4次满足迭代条件。。

(3)结构导向边缘检测滤波

加权相似系数计算公式：

$s\left[i\right]=\frac{{\left(\underset{j}{\mathrm{\&Sigma;}}h\right[i-j\left]f\right[j\left]\right)}^2}{\underset{j}{\mathrm{\&Sigma;}}h[i-j]f^2\left[j\right]}---\left(12\right)$

式中h[i]为加权系数，它相当于一个平滑滤波器。

利用加权相似系数相应的概念，定义结构导向相干系数计算公式为

$s_{\mathrm{2,1}}=\frac{{\&lang;{\&lang;f\&rang;}_2^2\&rang;}_1}{{\&lang;{\&lang;f^2\&rang;}_2\&rang;}_1}---\left(13\right)$

式中<·>₁表示沿第一坐标轴平滑，<·>₂表示沿第二坐标轴平滑。<·>₁有助于稳定计算，而<·>₂对微小值进行加权叠加。从结构导向意义上讲，两次平滑分别沿两个正交方向对地震数据进行平滑。对于三维地震数据而言，如果沿特征向量V和W所在平面进行平滑，则定义的相似计算公式为

$s_{\mathrm{vw},u}=\frac{{<{<f>}_{\mathrm{vw}}^2>}_u}{{<{<f^2>}_{\mathrm{vw}}>}_u}---\left(14\right)$

该公式定义的是一种平面相干属性体。

本发明方法的具体实现步骤如下：

(1)求取每个样点的梯度 $\mathrm{grad}=\left[\begin{array}{ccc}\frac{\mathrm{df}}{\mathrm{dx}}& \frac{\mathrm{df}}{\mathrm{dy}}& \frac{\mathrm{df}}{\mathrm{dt}}\end{array}\right]$

f表示原始三维地震数据，x，y，t分别表示所述三维地震数据的三个方向，利用奇异值分解计算结构张量T＝grad^Tgrad＝λ_uuu^T+λ_vvv^T+λ_www^T，该式中的T表示共轭转置；

其中，λ_u、λ_v和λ_w分别表示结构张量的奇异值，u、v和w分别为特征向量，表示在每个样点处利用原始三维地震数据的结构张量求取特征向量，这三个向量可以描述相应的方向；

规定特征值降序排列，λ_u≥λ_v≥λ_w≥0，其中特征向量u代表梯度最大的方向，垂直于地质体平面特征。特征向量w代表梯度最小的方向，平行于线性地质体。特征向量v和w都平行于地质体平面；

(2)利用结构导向平滑滤波器(即公式(2))计算结构导向边缘检测滤波(即公式(14))的分子，详见附图4，即首先对原始三维地震数据利用平滑方程(2)在时计算中间结果，再对中间结果平方后，再次利用平滑方程(2)在

时计算输出结果1，即(14)式的分子；

(3)利用结构导向平滑滤波器(即公式(2))计算结构导向边缘检测滤波(即公式(14))的分母，详见附图4，即首先对原始三维地震数据平方后，利用平滑方程(2)在

时计算中间结果，再对中间结果平方后，再次利用平滑方程(2)在

时计算输出结果2，即(14)式的分母；

计算分母时是将原始数据平方后进行计算的，计算分子时是直接用原始数据进行计算，分子和分母的计算就差在这第一步，因为原始数据平方前和平方后进行计算，输出是具有不同特征的数据。

(4)计算步骤(2)得到的分子与步骤(3)得到的分母的比值，即得到公式(14)的结果，该比值即为边缘检测结果。

具体实施时，在Cygwin下利用fortran语言对本发明方法进行了实现。采用某区地震资料对本方法进行了验证。附图2是用第一代相干算法计算得到的结果，附图3是用本发明方法计算得到的结果。从附图3上可以看出，本发明算法在断裂及溶洞刻画上优于第一代相干算法。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (4)

1.一种最优方向边缘监测方法，其特征在于：所述方法利用结构导向平滑滤波器结合结构导向边缘检测滤波提取三维地震数据中的地质体边界信息。

2.根据权利要求1所述的最优方向边缘监测方法，其特征在于：所述方法包括以下步骤：

(1)求取三维地震数据中的每个样点的梯度 $\mathrm{grad}=\left[\begin{array}{ccc}\frac{\mathrm{df}}{\mathrm{dx}}& \frac{\mathrm{df}}{\mathrm{dy}}& \frac{\mathrm{df}}{\mathrm{dt}}\end{array}\right],$ 利用奇异值分解计算所述三维地震数据的结构张量T(x)＝grad^Tgrad＝λ_uuu^T+λ_vvv^T+λ_www^T；

其中，f表示原始三维地震数据，x，y，t分别表示所述三维地震数据的三个方向，λ_u、λ_v和λ_w分别表示结构张量T的奇异值，u、v和w分别为特征向量；规定特征值降序排列，λ_u≥λ_v≥λ_w≥0，其中特征向量u代表梯度最大的方向，即垂直于地质体平面特征；特征向量w代表梯度最小的方向，平行于线性地质体；

特征向量v和w都平行于地质体平面；

(2)利用结构导向平滑滤波器计算结构导向边缘检测滤波的分子：

所述结构导向平滑滤波器为(2)式：

$g\left(x\right)-\mathrm{\&alpha;}\&dtri;\&CenterDot;D\left(x\right)\&dtri;g\left(x\right)=f\left(x\right)---\left(2\right)$

其中f(x)为输入的三维地震数据，g(x)为平滑输出结果；常数变量α为时间步长，α＝0时，f(x)＝g(x)，D(x)与T(x)具有相同特征向量，不同特征值；

所述结构导向边缘检测滤波为(14)式：

$s_{\mathrm{vw},u}=\frac{{<{<f>}_{\mathrm{vw}}^2>}_u}{{<{<f^2>}_{\mathrm{vw}}>}_u}---\left(14\right)$

(14)式中，s_vw，u为计算得到的边缘系数；

(3)利用结构导向平滑滤波器计算结构导向边缘检测滤波的分母；

(4)计算步骤(2)得到的分子与步骤(3)得到的分母的比值，该比值即为地震数据中的地质体边界信息。

3.根据权利要求2所述的最优方向边缘监测方法，其特征在于：所述步骤(2)是这样实现的：

首先对三维地震数据利用(2)式在

时进行计算得到中间结果，然后对中间结果进行平方，再次对平方后的中间结果利用(2)式在

时进行计算得到输出结果，该输出结果即为结构导向边缘检测滤波的分子。

4.根据权利要求2所述的最优方向边缘监测方法，其特征在于：所述步骤(3)是这样实现的：

首先对三维地震数据进行平方，然后对平方后的三维地震数据利用(2)式在

时进行计算得到中间结果，再对中间结果进行平方，再次对平方后的中间结果利用(2)式在

时进行计算得到输出结果，该输出结果即为结构导向边缘检测滤波的分母。

Images