CN103364833A - 一种高精度地层倾角估计方法 - Google Patents

Abstract

本发明公开了一种高精度地层倾角估计方法，首先计算梯度向量，然后向量场翻转，在进行向量场翻转、计算向量集合距离、计算向量滤波加权、最后进行向量加权滤波处理。该方法能够保持在同相轴单一方向延伸区域的地层倾角估计的空间一致性，提高倾角估计方法在断层、角度不整合以及褶皱形态的地层结构处的地层倾角估计精确度；该技术方案的算法内容易于实现，计算效率高；同时可以灵活调整分析窗口尺寸，满足不同倾角估计效果的需求，对噪声干扰严重的地震数据，可采用大的分析窗提高倾角估计结果的空间一致性，而对高信噪比的地震数据，可采用小的分析窗提高断层和角度不整合区域的倾角估计结果的精确度。

Description

一种高精度地层倾角估计方法

技术领域

本发明属于地震勘探技术领域，涉及一种地层倾角估计方法，尤其是一种利用加权向量方向滤波器高精度估计地层倾角和方位角。

背景技术

地震同相轴的倾角和方位角是继时间构造和层位拾取成图之后三维地震数据解释的重要产物，地层倾角和方位角属性体不仅自身包含有重要的地震地层学和地理学信息，可以直接用于勘探工区的解释，而且作为一种中间处理结果，可以用于后续的地震数据处理解释环节。三维地震数据体中每个采样位置的倾角和方位角对定义了一个局部的地层平面，从而允许沿着这个局部平面分析地层不连续性，也可利用局部地层连续性度量指导数据滤波。Barnes利用估计的倾角和方位角属性发展了一系列计算机自动生成的地震纹理属性，包括对地层发散、汇聚、平行以及混沌形态度量。Barnes通过光源照明分析模型将倾角和方位角属性联系起来给出一种阴影浮雕属性，该方法与利用二维颜色映射表表示倾角和方位角属性类似，以时间切片的方式显示，可突出细小断层和一些沉积特征，而这些特征往往与地层微小挤压或者反射波形的细微变化有关。由于倾角和方位属于地震数据的一阶导数属性，而曲率属于二阶导数属性，Al-Dossary和Marfurt给出直接基于倾角和方位角属性的体积曲率属性提取方法。因而，如果能够准确地估算三维地震数据体的倾角和方位角属性，不仅可以提高倾角和方位角自身作为解释属性的有效性，而且可以提高边缘检测和构造导向滤波的性能以及各种相关属性的计算精度。

Milkereit采用倾斜叠加方法计算二维地震剖面内同相轴的视倾角；Fomel与Schleicher等利用平面波分解方法获取局部地层倾角信息；Barnes利用二维复地震道分析理论计算瞬时倾角和方位角属性，随后，Barnes通过对瞬时地震属性进行局部加权平均，去除脉冲干扰、快变成分等杂乱变化，提高了瞬时地震倾角和方位角的计算效果；Marfurt等通过利用一系列的试验倾角截取地震数据道构造分析窗计算地震波形的相干值，将取得最大相干值对应的试验倾角记为局部地层倾角。Marfurt结合Kuwahara局部邻域多窗遍历的思想给出一种更为可靠的地层倾角估计方法，该方法避免了常规的中心分析窗横跨断层和角度不整合时造成不准确的倾角估计。然而，这种方法随着邻域分窗的增加，其计算复杂度也成倍增加；另外，通过对主测线和联络测线进行试验倾角扫描搜索真实地层倾角的方式，当倾角扫描间隔过于密集时也会加剧其计算复杂度，并且倾角扫描间隔较大时会遗漏一些细小结构。

利用有限差分方法计算反射波阵面的梯度向量进而估计地层倾角和方位角是一种简单而又有效的途径。Randen等利用梯度结构张量法对三维地震数据的梯度向量场进行平滑滤波获得局部地层倾角和方位角估计结果。Luo等提出基于加权结构张量法的梯度向量场平滑方法，该方法提高了倾角和方位角估计结果的一致性。Wang等给出一种向量翻转迭代滤波的方法来平滑梯度向量场。上述所有梯度向量场的平滑滤波方法都采用中心分析窗思想，在断层和角度不整合等地层边缘区域，倾角估计结果是两侧地层的平均取向，不能准确地反应当前分析点所处地层的取向。

梯度结构张量（GST）方法被广泛应用于自然图像处理分析中，用于确定信号结构和光流场方向估计、进行纹理特征分析以及提取物体特征点等。Randen等首先利用梯度结构张量方法分析三维叠后地震数据，估计局部波阵面的倾角和方位角。考虑图1a中的平面波波阵面，其真实方向向量为v_n，与局部波阵面相垂直。在理想情况下，波阵面的梯度向量与其真实方向向量v_n平行。然而，在实际情况下，一方面利用有限差分方法计算梯度向量不可避免会产生数值误差；另外，数据往往存在各种噪声干扰，使得计算的梯度向量也和真实方向向量之间存在偏差。根据图1b，假设某个分析邻域中心点x处的真实方向向量为v_n(x)，而在该局部区域内计算的梯度向量为g(y)，那么梯度向量和真实方向向量之间的偏差可表示为

e(x,y)=||e(x,y)||₂=||g(y)-(g(y)^Tv_n(x))v_n(x)||₂       (1)

在式(1)中，e(x,y)为误差向量，(g(y)^Tv_n(x))v_n(x)为梯度向量g(y)在v_n(x)方向的投影向量，易知e(x,y)度量梯度向量g(y)和方向向量v_n(x)之间的垂直距离。分析单一梯度向量与真实方向向量之间的偏差往往存在一定的随机性。因而，需要以点x为中心的高斯邻域内衡量两种方向向量之间的整体偏差程度：

ξ(x,y)=∫_Ωe²(x,y)G_ρ(y-x)dy        (2)

将式（2）进行简单的数学推导，整体误差函数ξ(x,y)可转化为

ξ(x,y)=∫_Ωg(y)^Tg(y)G_ρ(y-x)dy-∫_Ωv_n(x)^T(g(y)g(y)^T)v_n(x)G_ρ(y-x)dy.

(3)

根据式（3），如果已经求得以点x为中心的高斯邻域G_ρ(y-x)内各点的梯度向量，显然通过最小化整体误差函数ξ(x,y)可以获得真实方向向量v_n(x)的某种估计。由于右边项的第一项∫_Ωg(y)^Tg(y)G_ρ(y-x)dy与v_n(x)无关，因而最小化整体误差函数ξ(x,y)即等价于最大化右边项的第二项：

$\underset{v_n\left(x\right)}{\max }{\∫}_{\mathrm{\Ω}}{v}_n{\left(x\right)}^T\left(g\left(y\right)g{\left(y\right)}^T\right)v_n\left(x\right)G_{\mathrm{\ρ}}(y-x)\mathrm{dy},s.t.v_n{\left(x\right)}^T{v}_n\left(x\right)=1---\left(4\right)$

在式（4）中，由于目标函数的积分变量为y，与v_n(x)无关，可将目标函数重写为 $v_n^T\left({\∫}_{\mathrm{\Ω}}\left({\mathrm{gg}}^T\right)G_{\mathrm{\ρ}}\mathrm{dy}\right)v_n=v_n^T{J}_{\mathrm{\ρ}}{v}_n$

其中，J_ρ即为具有尺度参数ρ的梯度结构张量，

J_ρ=∫_Ω(gg^T)G_ρdy=Gρ*(gg^T)            (5)

此时，利用拉格朗日乘子法可将式(4)中带有约束的最优化问题化为一个无约束的最优化问题：

$\mathrm{\λ}(1-v_n^T{v}_n)+v_n^T{J}_{\mathrm{\ρ}}{v}_n.---\left(6\right)$

式(6)是关于v_n的二次凸问题，对其求导并令其导数等于零易得

J_ρv_n=λv_n.         (7)

因此，当利用梯度向量计算点x处的梯度结构张量时，该点反射波阵面方向向量的最佳估计对应于梯度结构张量的一个特征向量，并且根据式(4)最大化

J_ρv_n的优化目标可知，估计的方向向量是梯度结构张量的最大特征值对应的特征向量。

以上现有技术具有以下缺点：

(1)当遇到复杂地层结构或者低信噪比的地震数据时，梯度结构张量的三个特征值往往在幅值上比较接近，没有一个明显偏大的特征值，不好确定局部地层的方向估计；

(2)当遇到断层、角度不整合以及褶皱形态的地层结构时，由于在局部分析邻域内没有明确的地层延伸方向，梯度结构张量方法不能准确地估计分析点所处地层结构的方向。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种高精度地层倾角估计方法，其能够保持倾角估计方法在同相轴单一方向延伸区域的地层倾角估计的空间一致性，提高倾角估计方法在断层、角度不整合以及褶皱形态的地层结构处的地层倾角估计精确度，且易于实现，计算效率高。

本发明的目的是通过以下技术方案来解决的：

该种高精度地层倾角估计方法，包括以下步骤：

1）计算梯度向量

首先将低阶差分算子进行参数化表示，然后应用参数化的梯度算子计算合成平面波数据的梯度向量，得到带有参数的梯度向量表示，最后将这种参数化的梯度向量与平面波数据梯度向量的解析解形式进行比较，从而确定该差分算子中各项的未知参数；

2）向量场翻转

将梯度向量表示的方向信息转化为方位信息，即统一将梯度向量调整为半空间指向；

3）计算向量集合距离

针对翻转以后的向量场，计算滤波窗口W内任意向量

与所有其它向量之间的角度距离

$\mathrm{\α}({\stackrel{\‾}{V}}_i,{\stackrel{\‾}{V}}_j)={\cos }^{-1}\left(\frac{{\stackrel{\‾}{V}}_i^T{\stackrel{\‾}{V}}_j}{\left|{\stackrel{\‾}{V}}_i\right|\left|{\stackrel{\‾}{V}}_j\right|}\right),---\left(14\right)$

然后在滤波窗口W内，向量

对应的集合角度距离A_j可表示为该向量与所有其它向量之间的角度距离的平均：

$A_j=\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\α}({\stackrel{\‾}{V}}_i,{\stackrel{\‾}{V}}_j)\right)/N.---\left(15\right)$

4）计算向量滤波加权

构造向量滤波加权使其反比于被加权向量和滤波窗口W内其它相邻向量的集合距离，当分析向量取自滤波窗口W内取向最为集中向量子集，因其具有较小的集合距离而被赋予较大的加权值；当分析向量

受到噪声干扰表现为异常取向，因其具有较大的集合距离而被赋予较小的加权值；当在断层和角度不整合等地层不连续区域，滤波窗口W内同时包含两种或者多种延展方向的反射同相轴，那么需要使得占主导地位的同相轴方向向量被赋予较大的加权值，而其它处于次要地位的同相轴方向向量被赋予较小的加权值；根据以上要求，计算分析向量

的滤波加权μ_j为

${\mathrm{\μ}}_j=\frac{R^{\mathrm{\λ}-1}{(\mathrm{\π}-A_j)}^{\mathrm{\λ}}}{R^{\mathrm{\λ}-1}{\left(\mathrm{\π}{-A}_j\right)}^{\mathrm{\λ}}+{(1-R)}^{\mathrm{\λ}-1}{\left(A_j\right)}^{\mathrm{\λ}}}---\left(16\right)$

其中，滤波加权μ_j为关于集合角度距离A_j的单调递减的S形映射函数，R为S形映射函数的过渡点位置，且R∈(0,1)；λ控制S形映射函数过渡区域的性态，且λ≥1；

5）向量加权滤波处理

对于三维地震数据，滤波窗口W为以目标点为中心、尺寸为N=N_x×N_y×N_z的立方体，滤波窗口W内翻转的梯度向量为

对于二维地震数据，滤波窗口W为以目标点为中心、尺寸为N=N_x×N_z或N_y×N_z的矩形窗口，滤波窗口W内翻转的梯度向量为

在上述两种情况下，向量加权滤波处理具有如下统一的加权滤波形式：

$\stackrel{\widehat{\‾}}{V}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_j{\stackrel{\‾}{V}}_j=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{\μ}}_j/\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_j){\stackrel{\‾}{V}}_j---\left(13\right)$

式中，

为当前滤波窗口W内目标位置的滤波输出；w_j为向量

的归一化加权值；μ_j为向量

的加权值。

进一步的，在以上步骤1）中，依据下式计算地震数据的梯度向量：

（1）对于二维数据U，利用强化各向同性梯度算子计算的梯度向量的各个分量为：

$\left\{\begin{array}{c}{V}_x(i,j)=(U_{i+1,j}-U_{i-1,j})+0.25(U_{i+1,j-1}-U_{i-1,j-1}+U_{i+1,j+1}-U_{i-1,j+1})\\ V_z(i,j)=(U_{i,j+1}-U_{i,j-1})+0.25(U_{i-1,j+1}-U_{i-1,j-1}+U_{i+1,j+1}-U_{i+1,j-1})\end{array}\right..---\left(9\right)$

式中，V_x(i,j)为二维地震数据U在位置x_i,j处沿x方向的梯度分量；V_z(i,j)为二维地震数据U在位置x_i,j处沿z方向的梯度分量；

（2）对于三维数据，应用强化各向同性梯度算子计算的梯度向量的各个分量为：

V_x(i,j,k)=(U_i+1,j,k-U_i-1,j,k)+0.245*(U_i+1,j-1,k-U_i-1,j-1,k+U_i+1,j+1,k

-U_i-1,j+1,k)+0.245*(U_i+1,j,k-1-U_i-1,j,k-1+U_i+1,j,k+1-U_i-1,j,k+1)

+0.085*(U_i+1,j-1,k-1-U_i-1,j-1,k-1+U_i+1,j+1,k-1-U_i-1,j+1,k-1)

+0.085*(U_i+1,j-1,k+1-U_i-1,j-1,k+1+U_i+1,j+1,k+1-U_i-1,j+1,k+1),

(10)

V_y(i,j,k)=(U_i,j+1,k-U_i,j-1,k)+0.245*(U_i-1,j+1,k-U_i-1,j-1,k+U_i+1,j+1,k

-U_i+1,j-1,k)+0.245*(U_i,j+1,k-1-U_i,j-1,k-1+U_i,j+1,k+1-U_i,j-1,k+1)

+0.085*(U_i-1,j+1,k-1-U_i-1,j-1,k-1+U_i+1,j+1,k-1-U_i+1,j-1,k-1)

+0.085*(U_i-1,j+1,k+1-U_i-1,j-1,k+1+U_i+1,j+1,k+1-U_i+1,j-1,k+1),

(11)

V_z(i,j,k)=(U_i,j,k+1-U_i,j,k-1)+0.245*(U_i-1,j,k+1-U_i-1,j,k-1+U_i+1,j,k+1

-U_i+1,j,k-1)+0.245*(U_i,j-1,k+1-U_i,j-1,k-1+U_i,j+1,k+1-U_i,j+1,k-₁)

+0.085*(U_i-1,j-1,k+1-U_i-1,j-1,k-1+U_i+1,j-1,k+1-U_i+1,j-1,k-1)

+0.085*(U_i-1,j+1,k+1-U_i-1,j+1,k-1+U_i+1,j+1,k+1-U_i+1,j+1,k-₁).

(12)

。

进一步的，上述步骤2）具体为：首先选择具有零倾角的水平地层的方向向量作为参考向量V_R，该向量垂直指向下半空间；然后逐个比较梯度向量和参考向量V_R之间的夹角，如果它们之间的夹角大于90°，那么就将该梯度向量的各个分量均乘以-1进行180°的方向翻转，使得翻转处理后的梯度向量统一指向下半空间。

本发明具有以下有益效果：

本发明能够保持在同相轴单一方向延伸区域的地层倾角估计的空间一致性，提高倾角估计方法在断层、角度不整合以及褶皱形态的地层结构处的地层倾角估计精确度；该技术方案的算法内容易于实现，计算效率高；同时可以灵活调整分析窗口尺寸，满足不同倾角估计效果的需求，对噪声干扰严重的地震数据，可采用大的分析窗提高倾角估计结果的空间一致性，而对高信噪比的地震数据，可采用小的分析窗提高断层和角度不整合区域的倾角估计结果的精确度。

附图说明

图1为现有技术中利用梯度结构张量估计地层方向原理图；

图2为本发明梯度度向量场翻转与均值平滑处理示意图，

其中（a）是根据参考向量V_R翻转梯度向量；（b）是对翻转梯度向量进行平滑得到

图3为加权向量方向滤波器的S形加权映射函数图；

图4为加权向量方向滤波器在断层区域的滤波操作示意图；

其中:(a)含有断层的合成模型;(b)对(a)中断层区域局部放大；

图5为二维合成模型的地层倾角估计图；

其中：（a）原始模型；（b）信噪比为11dB的加噪模型；（c）GST方法的倾角估计结果；（d）BCDF方法的倾角估计结构；（e）WVDF方法的倾角估计结果；

图6为图5b-e中标记的黑色矩形框区域对应的放大显示图；

图7为某油田三维地震资料的倾角属性分析图。

具体实施方式

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

梯度向量计算与向量场翻转

有多种有限差分方法可以计算三维地震波场的梯度，如中心差分、Sobel和Prewitt差分算子等，对于三维地震数据U中任意位置x_i,j,k的梯度向量V=[V_x,V_y,V_t]可以利用中心差分法求得

$\left\{\begin{array}{c}{V}_x(i,j,k)=(U_{i+1,j,k}-U_{i-1,j,k})\\ V_y(i,j,k)=(U_{i,j+1,k}-U_{i,j-1,k})\\ V_t(i,j,k)=(U_{i,j,k+1}-U_{i,j,k-1})\end{array}\right..---\left(8\right)$

然而，这些常用的差分方法都是低阶算子，由于没有考虑到信号结构的频率成分和方向特性，往往会导致梯度向量的计算不精确。本发明采用强化各向同性梯度算子来计算梯度向量。强化各向同性梯度算子首先将常规的低阶差分算子进行参数化表示，然后应用参数化的梯度算子计算合成平面波数据的梯度向量，得到带有参数的梯度向量表示，最后将这种参数化的梯度向量与平面波数据梯度向量的解析解形式进行比较，从而确定该差分算子中各项的未知参数，使得由于信号的频率成分和方向特性造成的数值偏差最小化。对于二维数据U，利用强化各向同性梯度算子计算的梯度向量的各个分量为

$\left\{\begin{array}{c}{V}_x(i,j)=(U_{i+1,j}-U_{i-1,j})+0.25*(U_{i+1,j-1}-U_{i-1,j-1}+U_{i+1,j+1}-U_{i-1,j+1})\\ V_t(i,j)=(U_{i,j+1}-U_{i,j-1})+0.25*(U_{i-1,j+1}-U_{i-1,j-1}+U_{i+1,j+1}-U_{i+1,j-1})\end{array}\right..---\left(9\right)$

式中：V_x(i,j)为二维地震数据U在位置x_i,j处沿x方向的梯度分量；

V_z(i,j)为二维地震数据U在位置x_i,j处沿z方向的梯度分量；

而在三维情况下，应用强化各向同性梯度算子计算的梯度向量的各个分量为

V_x(i,j,k)=(U_i+1,j,k-U_i-1,j,k)+0.245*(U_i+1,j-1,k-U_i-1,j-1,k+U_i+1,j+1,k

-U_i-1,j+1,k)+0.245*(U_i+1,j,k-1-U_i-1,j,k-1+U_i+1,j,k+1-U_i-1,j,k+1)     (10)

+0.085*(U_i+1,j-1,k-1-U_i-1,j-1,k-1+U_i+1,j+1,k-1-U_i-1,j+1,k-1)

+0.085*(U_i+1,j-1,k+1-U_i-1,j-1,k+1+U_i+1,j+1,k+1-U_i-1,j+1,k+1),

V_y(i,j,k)=(U_i,j+1,k-U_i,j-1,k)+0.245*(U_i-1,j+1,k-U_i-1,j-1,k+U_i+1,j+1,k

-U_i+1,j-1,k)+0.245*(U_i,j+1,k-1-U_i,j-1,k-1+U_i,j+1,k+1-U_i,j-1,k+1)    (11)

+0.085*(U_i-1,j+1,k-1-U_i-1,j-1,k-1+U_i+1,j+1,k-1-U_i+1,j-1,k-1)

+0.085*(U_i-1,j+1,k+1-U_i-1,j-1,k+1+U_i+1,j+1,k+1-U_i+1,j-1,k+1),

V_z(i,j,k)=(U_i,j,k+1-U_i,j,k-1)+0.245*(U_i-1,j,k+1-U_i-1,j,k-1+U_i+1,j,k+1

-U_i+1,j,k-1)+0.245*(U_i,j-1,k+1-U_i,j-1,k-1+U_i,j+1,k+1-U_i,j+1,k-1)     (12)

+0.085*(U_i-1,j-1,k+1-U_i-1,j-1,k-1+U_i+1,j-1,k+1-U_i+1,j-1,k-₁)

+0.085*(U_i-1,j+1,k+1-U_i-1,j+1,k-1+U_i+1,j+1,k+1-U_i+1,j+1,k-1).

基于有限差分方法的求导操作会放大地震资料中的高频噪声，但是又不能对梯度向量场进行简单的均值平滑来压制干扰噪声，这是因为波阵面有两个完全相反的方向，例如水平地层的倾角既可以为0°，也可以为180°，如果采用简单的均值平滑滤波的话，滤波结果很有可能是90°的错误倾角估计。因此，对梯度向量场进行平滑滤波前，需要将梯度向量表示的方向信息转化为方位信息，即统一将梯度向量调整为半空间指向，这就是梯度向量场的翻转处理。根据图2中的梯度向量场翻转的示意图，首先选定一个参考向量V_R，该技术方案选择具有零倾角的水平地层的方向向量作为参考向量，该向量垂直指向下半空间；然后逐个比较梯度向量和参考向量之间的夹角，如果它们之间的夹角大于90°，那么就将该梯度向量的各个分量均乘以-1进行180°的方向翻转，使得翻转处理后的梯度向量统一指向下半空间。

Wang等提出对翻转梯度向量场进行迭代均值平滑滤波处理的方法，迭代过程中对参考向量V_R进行不断修正。分析图2b中对翻转后的梯度向量场进行均值平滑处理时，如果待处理向量取向非常接近时，采用均值平滑滤波方法得到的平均向量将具有向量集合中近似中间位置的指向，即取得预期的滤波效果。然而，如果待处理向量中大部分的取向较为接近，而小部分向量的指向与多数向量取向的偏离较大，超过90°甚至接近180°夹角时，采用均值平滑方法就会得到错误的方向向量估计结果。例如大部分向量指向接近0°，而小部分向量的指向靠近180°，这种情形会在大倾角地层情况时遇到，采用均值平滑滤波获得地层方向估计结果必然会大幅度偏离0°，从而会错误地反映为小角度倾斜地层结构。因此，有必要为翻转后的梯度向量场设计更为合理的滤波方法，提高地层方向估计的准确度；另外，由于基于有限差分的地层方向估计方法都采用中心窗分析方式，致使在断层和角度不整合等地层不连续性区域给出错误的方向估计结果，因而有必要提高现有方法对地层不连续性区域的方向估计效果，从而更有利于对这些区域的分析刻画。

加权向量方向滤波方法

向量滤波方法在彩色图像处理领域得到广泛的研究和应用，这类方法包括向量中值滤波器（Vector Median Filter，VMF），基本向量方向滤波器（Basic VectorDirectional Filter，BVDF）以及两者衍生的各种向量滤波方法等。其中，基本向量方向滤波器仅与待分析向量的方向信息有关，等价于将待分析向量归一化后的向量中值滤波结果。因而，应用基本向量方向滤波器对翻转梯度向量场进行滤波，将非常有利于剔除分析窗内异常取向的向量，提高地层方向估计的准确性；另外，对于断层和角度不整合等地层不连续区域，基本向量方向滤波器可以较好地保持地层方向变化的边界，降低这些区域估计结果的不确定性。然而，当地震数据受到噪声干扰或者对于复杂地层结构的区域，分析窗内的梯度向量往往都会偏离真实的地层结构方向，基本向量方向滤波器选择分析窗内的中值方向向量作为滤波输出，如果该输出向量仍然偏离真实的地层方向，就会导致估计的地层方向结果出现抖动变化的现象，缺乏空间一致性。为了解决这一问题，本发明提出利用向量方向距离作为权值的加权向量方向滤波器（Weighted Vector Directional Filter，WVDF）来平滑翻转的梯度向量场。

对于三维数据体中以目标点为中心、尺寸为的立方体分析窗N_x×N_y×N_z，将分析窗内翻转的梯度向量记为其中N=N_x×N_y×N_z。同样地，用于二维地震剖面分析的矩形分析窗，若其尺寸为N_x×N_z或者N_y×N_z，即有N=N_x×N_z或N_y×N_z。对于分析窗W内翻转的梯度向量，WVDF方法具有如下的加权滤波形式：

$\stackrel{\widehat{\‾}}{V}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_j{\stackrel{\‾}{V}}_j=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{\μ}}_j/\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_j){\stackrel{\‾}{V}}_j---\left(13\right)$

其中，

为WVDF在当前分析窗内目标位置的滤波输出；w_j为向量

的归一化加权值；μ_j为向量

的预加权值。

根据式(13)，WVDF的滤波输出为分析窗内所有向量的加权平均，而分析窗内任意向量

的加权系数μ_j与该向量在分析窗内与其他向量的某种集合距离有关。为了获得良好的地层方向估计效果，通常期望滤波器的滤波输出主要来自分析窗内取向最为集中的向量子集的贡献，即取向最为集中的向量子集都具有较大的加权值，而其它向量具有较小的加权值。因而，这就要求自适应加权值本身作为一种置信度量，反映分析窗内的向量是否属于取向最为集中向量子集。本发明构造自适应加权值使其反比于被加权向量和分析窗内其它相邻向量的集合距离，当分析向量取自滤波窗内取向最为集中向量子集，因其具有较小的集合距离而被赋予较大的加权值；而当分析向量受到噪声干扰等表现为异常取向，会因其具有较大的集合距离而被赋予较小的加权值，使其对滤波结果的影响较小；另外，在断层和角度不整合等地层不连续区域，如果分析窗内同时包含两种或者多种延展方向的反射同相轴，那么需要使得占主导地位的同相轴方向向量被赋予较大的加权值，而其它处于次要地位的同相轴方向向量被赋予较小的加权值，从而使得滤波输出结果主要反映占主导地位的同相轴方向，降低这些区域内地层方向估计的不确定性。

针对地层方向估计性能的要求，在WVDF的设计中，本发明利用向量角度距离作为向量之间的距离度量，考虑如下两个向量之间的夹角角度

$\mathrm{\α}({\stackrel{\‾}{V}}_i,{\stackrel{\‾}{V}}_j)={\cos }^{-1}\left(\frac{{\stackrel{\‾}{V}}_i^T{\stackrel{\‾}{V}}_j}{\left|{\stackrel{\‾}{V}}_i\right|\left|{\stackrel{\‾}{V}}_j\right|}\right),---\left(14\right)$

那么在滤波窗W内，向量

对应的集合角度距离可表示为

$A_j=\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\α}({\stackrel{\‾}{V}}_i,{\stackrel{\‾}{V}}_j)\right)/N.---\left(15\right)$

在式(15)中，向量

的集合角度距离A_j在区域[0,π]内取值，度量向量

相对滤波窗W内所有其它向量的整体偏离程度。当A_j>π/2时表示该向量为异常向量，偏离滤波窗内的绝大多数向量；而当A_j→0时表示该窗内的向量取向较为集中，并且该向量也属于分析窗内集中取向的向量子集。由于这种向量距离度量方法完全符合地层方向滤波的要求，因此可将其用于分析地层梯度向量。

在确定滤波窗内每个向量的集合距离度量之后，就需要建立这种集合距离度量和该向量的加权系数之间的映射关系。根据本发明对WVDF的设计目标，一种可行的映射关系要将最小集合角度距离转化为最大的向量加权值，而将最大集合角度距离转化为最小的向量加权值。因此，根据这种映射关系的要求，本发明为WVDF构造一种单调递减的S形映射函数，从而向量

的加权系数μ_j计算为

${\mathrm{\μ}}_j=\frac{R^{\mathrm{\λ}-1}{(\mathrm{\π}-A_j)}^{\mathrm{\λ}}}{R^{\mathrm{\λ}-1}{\left(\mathrm{\π}{-A}_j\right)}^{\mathrm{\λ}}+{(1-R)}^{\mathrm{\λ}-1}{\left(A_j\right)}^{\mathrm{\λ}}}---\left(16\right)$

其中，R为S形映射函数的过渡点位置，且R∈(0,1)；λ控制S形映射函数过渡区域的性态，且λ≥1。

图3中为选择三组不同参数R和λ的S形加权映射函数曲线，从图中可以看出，当集合角度距离A_j→0时，向量加权值μ_j→1；而当A_j→π时，μ_j→0。这种加权映射函数在加权向量方向滤波操作中起着模糊分类器的作用。如图4所示，根据向量的集合角度距离A_j，加权映射函数将滤波窗内的向量分为取向最为集中的向量子集和异常取向的向量子集，对取向最为集中的向量子集中的向量分配较大的加权值，使其对滤波输出起着主要作用；而对异常取向的向量子集中的向量分配较小的加权值，使其几乎不影响地层方向的估计结果。这里，加权映射函数的参数R起着子集分割阈值界限的作用，而参数λ控制着两类向量子集之间过渡的平缓程度。从图3中可见，大的λ取值会导致更为陡峭的过渡区域，而小的λ取值使得两个向量子集之间平缓过渡；大的R取值会使加权映射函数对取向最为集中的向量子集的选择更为严格，而小的R取值会使更多的向量包括到取向最为集中的向量子集中。因此，可以灵活调整两个参数的取值，使得WVDF适用于不同条件地震数据的处理分析。当地震数据的信噪比偏高时，可以选取大的R和大的λ，使得加权映射函数具有更强的选择性，从而提高断层和角度不整合等地层不连续区域的地层方向估计性能；而对信噪比偏低的地震数据，需要选择较小的R和较小的λ，使得取向最为集中的向量子集中包含更多的向量，从而保证了WVDF足够的滤波性能，可以提高地层方向估计的空间一致性。

影响WVDF的滤波性能的另一个因素是数据分析窗的尺寸大小。根据不同的解释和处理目标，N_x、N_y和N_z可在3点到21点之间取值。一般地，对于大的数据分析窗可以得到具有高度空间一致性的地层倾角和方位角估计结果，适用于识别大尺度的地层结构；相反地，对于小的数据分析窗可以得到具有高分辨率的地层倾角和方位角属性，有利于一些细微结构的检测。如果地震数据经过精细的保边缘滤波处理，由于不需要重点考虑对梯度向量场的噪声衰减处理，可以优选小的数据分析窗对地震数据进行倾角和方位角属性分析，获得高分辨率的地层方向估计结果。

分析加权向量方向滤波过程，首先需要计算分析窗内每个向量的集合角度距离，即遍历分析窗内的每个向量，并计算当前向量和历经向量之间的内积，进一步得到它们之间的夹角以及当前向量的集合角度距离，最后将各个向量的集合角度距离转变为加权值得到加权滤波结果。可以看出，加权向量方向滤波操作具有较高的计算复杂度。图中，红色箭头指示根据垂直向下的参考向量翻转后的梯度向量，（b）中的绿色框表示以目标点C为中心的分析窗，蓝色表示大的加权值，透明部分表示非常小的加权值。集合角度距离的计算步骤。不失一般性，考虑一个立方体的分析窗，其尺寸为N_x=N_y=N_z=M，那么对每个分析窗内计算向量的集合角度距离所需的内积运算次数为

n_vip=M³(M³-1)/2～O(M⁶).            (17)

从式(17)可以看出，总的内积运算次数n_vip接近M⁶的量级，将是一个相对较大的数目。然而，考虑到相邻两点的分析窗具有非常高的重叠度，当前点对应分析窗内一部分向量的角度距离可由其相邻位置的分析窗内已经得到的结果获得，从而可以降低总体的计算量。容易知道两个相邻分析点对应的分析窗有M²(M-1)个内积运算完全相同，只有M²个内积运算不同。如果考虑相邻分析点之间运算的重叠度，那么每个分析窗内所必须的内积运算次数变为

n_vip-r=M²(M-1)M²+M²(M²-1)/2～O(M⁵).         (18)

另外，还可以利用地震数据呈层状结构分布的特点对WVDF的滤波过程进行一定程度的简化，从而进一步降低其计算复杂度。本发明对WVDF的滤波过程作出如下的简化：首先对每个立方体分析窗沿着垂直方向进行简单的均值平滑操作，得到一个尺寸为N_x×N_y的矩形向量场；然后应用WVDF处理得到的矩形向量场，获得当前分析点的地层方向估计结果。由于地震数据具有呈层状分布的特点，而且地层的不连续性往往表现为大角度结构，因而这种简化操作在多数情况是可行的，并且几乎不影响地震数据横向不连续性的刻画和细节特征的分析。

数值仿真结果

合成模型资料

为了验证WVDF估计地层方向的有效性，本发明利用平面波信号得到二维的含有断层的合成模型。该模型采用的平面波信号定义为

u(x,z)=sin(2πf(-x·sinθ+z·cosθ)),       (19)

其中，f为平面波信号的频率，即合成地震同相轴的频率；θ为平面波信号的波前方向，即合成地震同相轴的倾角。

如图5a所示，该合成模型的左侧为具有地层倾角θ=0°的水平地层，而右侧为具有地层倾角θ=15°的倾斜地层，两种不同倾角的地层由一处弯曲的断层分割开。分别向合成模型加入8dB、11dB和14dB的加性高斯白噪声得到含噪的合成模型。图5b为11dB的含噪合成数据剖面，从图中可以看出，这种噪声强度虽然不会影响对各种地层结构的分析和识别，但是经过差分求导运算后，噪声干扰会被放大，从而会严重影响对地层方向的估计效果。

为WVDF方法选择参数R=0.1和λ=4，对应的S形加权映射函数如图3中的实线表示，然后将WVDF用于处理含噪合成模型的翻转梯度向量场。为了说明本发明方法的有效性，进一步将梯度结构张量方法（GST）和基于基本向量方向滤波器（BVDF）的翻转向量场滤波方法用于估计该含噪合成模型的地层方向，并选择这两种方法的滤波分析窗尺寸与WVDF相同，均为9×9的分析窗。为了对比三种方法对地层倾角估计的性能，本发明分别从得到倾角结果的视觉效果以及它们与真实地层倾角之间的均方根误差进行分析；另外，为了更为细致地对比各种处理效果，本发明将合成模型剖面分为断层区域和非断层区域分别对待。考虑到分析窗的尺寸大小，本发明简单地将两种不同倾角地层分界线两侧各取半个分析窗大小的区域，即每边各取4个样点，联合构成断层区域，而将其他区域作为非断层区域。这样可以得到每种方法的估计结果在断层和非断层区域的均方根误差，从而很容易对比每种方法在非断层区域的滤波性能以及在断层区域对地层倾角突变的刻画能力。针对三种不同信噪比水平的合成数据剖面，三种方法得到的倾角估计结果的均方根误差如表1所示，并且将信噪比为11dB的合成数据剖面的倾角估计结果显示在图5c-e。比较表1中均方根误差结果可以看出，WVDF方法和BVDF方法在断层区域的倾角估计性能均要好于GST方法，而WVDF方法和GST方法在非断层区域比BVDF方法获得更好的滤波效果，这种结论可以从图5c-e的结果得到验证。为了从视觉效果上进行更加细致地对比，将图5b-e中黑色框选择的断层区域在图6中放大显示。对比图6b-d可以清晰地看出，GST方法在断层区域由于将断层作为一种大倾角结构而造成错误的地层倾角估计结果，而WVDF方法和BVDF方法可以较为精确地定位不同倾向地层的倾角突变边界；对于非断层区域，WVDF方法和GST方法可以得到较好的空间一致性的地层倾角估计，而BVDF方法由于滤波性能的限制，估计的倾角结果中存在较为严重的块状现象，可解释性较差。因而可见，本发明提出的WVDF方法不仅兼顾了BVDF方法对地层倾角变化边界精确定位的优势，而且还保持了GST方法优异的平滑滤波性能。

表1三种不同倾角估计方法处理含噪合成模型对应的均方根误差（RMSE）

实际地震资料

将提出的WVDF方法用于估计某油田三维地震资料的地层方向。该地震数据的时间采样频率为1kHz，道间距和线间距都是20.0m。在该三维区块中，一个复杂的断层网络横截几条河道砂体沉积结构，因而要求地层方向估计方法具有精确的倾角变化边缘和细节结构保持能力，从而有助于提高后继地震资料处理和分析的精度。与合成模型的处理类似，为WVDF方法的S形加权映射函数选择参数R=0.1和λ=4，并选取滤波器的分析窗尺寸为9×9×9。另外，为了进一步说明本发明方法的有效性，还利用GST方法、瞬时倾角估计方法和基于相似体度量的倾角扫描方法估计该三维地震资料的地层方向，其中，GST方法的分析窗尺寸也取为9×9×9，而倾角扫描方法的分析窗的空间尺寸为9×9，垂直方向尺寸近似取为地震波的波长长度，为19样点。图7a是该三维纯波地震资料过测线AA′的垂直剖面，图7b-e为四种方法估计的地层倾角结果。对比各种方法估计的倾角剖面可以看出，瞬时倾角估计方法对噪声干扰非常敏感，得到的地层倾角的空间一致性较差，并且在大倾角地层处由于发生频率混叠而给出错误的倾角估计结果；基于相似体度量的倾角扫描方法会造成条带状的假象，并且由于倾角扫描步长和扫描范围的限制，会在地层倾角发生突变的断层区域给出错误的倾角估计结果；GST方法和WVDF方法由于它们优异的滤波性能受噪声干扰的影响小，都能够得到具有较高空间一致性的地层倾角估计结果，但是分析两种方法在地层倾角发生突变的断层区域的估计效果，可以看出GST方法在断层区域得到的是断层两侧的地层方向平均取向，而本发明提出的WVDF方法可以较好地区别断层两侧的不同地层取向，能够较为精确地定位地层方向发生突变的边界。因此，将本发明的WVDF方法得到的地层倾角和方位角估计结果进一步用于后继的处理解释环节，可以得到更为精细的地层不连续性刻画效果。

Claims (3)

1.一种高精度地层倾角估计方法，其特征在于，包括以下步骤：

1)计算梯度向量

首先将低阶差分算子进行参数化表示，然后应用参数化的梯度算子计算合成平面波数据的梯度向量，得到带有参数的梯度向量表示，最后将这种参数化的梯度向量与平面波数据梯度向量的解析解形式进行比较，从而确定该差分算子中各项的未知参数；

2)向量场翻转

将梯度向量表示的方向信息转化为方位信息，即统一将梯度向量调整为半空间指向；

3)计算向量集合距离

针对翻转以后的向量场，计算滤波窗口W内任意向量

与所有其它向量之间的角度距离

$\mathrm{\α}({\stackrel{\‾}{V}}_i,{\stackrel{\‾}{V}}_j)={\cos }^{-1}\left(\frac{{\stackrel{\‾}{V}}_i^T{\stackrel{\‾}{V}}_j}{{|\stackrel{\‾}{V}}_i\left|\right|{\stackrel{\‾}{V}}_j|}\right),---\left(14\right)$

然后在滤波窗口W内，向量

对应的集合角度距离A_j可表示为该向量与所有其它向量之间的角度距离

的平均：

$A_j=\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\α}({\stackrel{\‾}{V}}_i,{\stackrel{\‾}{V}}_j)\right)/N.---\left(15\right)$

4)计算向量滤波加权

构造向量滤波加权使其反比于被加权向量和滤波窗口W内其它相邻向量的集合距离，当分析向量

取自滤波窗口W内取向最为集中向量子集，因其具有较小的集合距离而被赋予较大的加权值；当分析向量

受到噪声干扰表现为异常取向，因其具有较大的集合距离而被赋予较小的加权值；当在断层和角度不整合等地层不连续区域，滤波窗口W内同时包含两种或者多种延展方向的反射同相轴，那么需要使得占主导地位的同相轴方向向量被赋予较大的加权值，而其它处于次要地位的同相轴方向向量被赋予较小的加权值；根据以上要求，计算分析向量

的滤波加权μ_j为

${\mathrm{\μ}}_j=\frac{R^{\mathrm{\λ}-1}{(\mathrm{\π}-A_j)}^{\mathrm{\λ}}}{R^{\mathrm{\λ}-1}{(\mathrm{\π}-A_j)}^{\mathrm{\λ}}+{(1-R)}^{\mathrm{\λ}-1}{\left(A_j\right)}^{\mathrm{\λ}}}---\left(16\right)$ 其中，滤波加权μ_j为关于集合角度距离A_j的单调递减的S形映射函数，R为S形映射函数的过渡点位置，且R∈(0，1)；λ控制S形映射函数过渡区域的性态，且λ≥1；

5)向量加权滤波处理

对于三维地震数据，滤波窗口W为以目标点为中心、尺寸为N＝N_x×N_y×N_z的立方体，滤波窗口W内翻转的梯度向量为

j＝1，2，…，N；对于二维地震数据，滤波窗口W为以目标点为中心、尺寸为N＝N_x×N_z或N_y×N_z的矩形窗口，滤波窗口W内翻转的梯度向量为

j＝1，2，…，N；在上述两种情况下，向量加权滤波处理具有如下统一的加权滤波形式：

$\stackrel{\widehat{\‾}}{V}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_j{\stackrel{\‾}{V}}_j=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{\μ}}_j/\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_j){\stackrel{\‾}{V}}_j---\left(13\right)$

式中，

为当前滤波窗口W内目标位置的滤波输出；w_j为向量

的归一化加权值；μ_j为向量

的加权值。

2.根据权利要求1所述的高精度地层倾角估计方法，其特征在于，在步骤1)中，依据下式计算地震数据的梯度向量：

(1)对于二维数据U，利用强化各向同性梯度算子计算的梯度向量的各个分量为：

$\left\{\begin{array}{c}{V}_x(i,j)=(U_{i+1,j}-U_{i-1.j})+0.25(U_{i+1,j-1}-U_{i-1,j-1}+U_{i+1,j+1}-U_{i-1,j+1})\\ V_z(i,j)=(U_{i,j+1}-U_{i,j-1})+0.25(U_{i-1,j+1}-U_{i-1,j-1}+U_{i+1,j+1}-U_{i+1,j-1})\end{array}\right..---\left(9\right)$

式中，V_x(i，j)为二维地震数据U在位置x_i，j处沿x方向的梯度分量；V_z(i，j)为二维地震数据U在位置x_i，j处沿z方向的梯度分量；

(2)对于三维数据，应用强化各向同性梯度算子计算的梯度向量的各个分量为：

V_x(i，j，k)＝(U_i+1，j，k-U_i-1，j，k)+0.245*(U_{i+1，j-1，k}-U_{i-1，j-1，k}+U_{i+1，j+1，k}-U_{i-1，j+1，k})+0.245*(U_{i+1，j，k-1}-U_{i-1，j，k-1}+U_{i+1，j，k+1}-U_{i-1，j，k+1)}

                                                                 (10)+0.085*(U_{i+1，j-1，k-1}-U_{i-1，j-1，k-1}+U_{i+1，j+1，k-1}-U_{i-1，j+1，k-1})+0.085*(U_{i+1，j-1，k+1}-U_{i-1，j-1，k+1}+U_{i+1，j+1，k+1}-U_{i-1，j+1，k+1})，

V_y(i，j，k)＝(U_i，j+1，k-U_i，j-1，k)+0.245*(U_{i-1，j+1，k}-U_{i-1，j-1，k}+U_{i+1，j+1，k}-U_{i+1，j-1，k})+0.245*(U_{i，j+1，k-1}-U_{i，j-1，k-1}+U_{i，j+1，k+1}-U_{i，j-1，k+1})

                                                                      (11)+0.085*(U_{i-1，j+1，k-1}-U_{i-1，j-1，k-1}+U_{i+1，j+1，k-1}-U_{i+1，j-1，k-1})+0.085*(U_{i-1，j+1，k+1}-U_{i-1，j-1，k+1}+U_{i+1，j+1，k+1}-U_{i+1，j-1，k+1})，

V_z(i，j，k)＝(U_i，j，k+1-U_i，j，k-1)+0.245*(U_{i-1，j，k+1}-U_{i-1，j，k-1}+U_{i+1，j，k+1}-U_{i+1，j，k-1})+0.245*(U_{i，j-1，k+1}-U_{i，j-1，k-1}+U_{i，j+1，k+1}-U_{i，j+1，k-1})

                                                                    (12)+0.085*(U_{i-1，j-1，k+1}-U_{i-1，j-1，k-1}+U_{i+1，j-1，k+1}-U_{i+1，j-1，k-1})+0.085*(U_{i-1，j+1，k+1}-U_{i-1，j+1，k+1}+U_{i+1，j+1，k+1}-U_{i+1，j+1，k-1}).

。

3.根据权利要求1所述的高精度地层倾角估计方法，其特征在于，步骤2)具体为：首先选择具有零倾角的水平地层的方向向量作为参考向量V_R，该向量垂直指向下半空间；然后逐个比较梯度向量和参考向量V_R之间的夹角，如果它们之间的夹角大于90°，那么就将该梯度向量的各个分量均乘以-1进行180°的方向翻转，使得翻转处理后的梯度向量统一指向下半空间。

Images