用于地震数据处理的迭代倾角导向中值滤波器

摘要

提供了用于对在地震数据中示出有冲突倾角的地震数据进行随机噪声衰减的迭代倾角导向中值滤波器。所述数据的处理窗口或样本内的若干支配倾角通过傅里叶径向变换在频率-波数域中被识别。之后沿着所述支配倾角施加中值滤波器以去除噪声，滤波之后的剩余信号被保留用于进一步的中值滤波器迭代。重复迭代以在所述剩余数据中沿着主要支配倾角施加中值滤波器。在随后的迭代中继续处理直到所有选定的倾角被处理完。每次迭代的剩余信号之后被累加用于最终输出。

说明书

技术领域

本发明涉及用于地下地层结构分析的地震数据处理，特别是涉 及用于这种处理的迭代倾角导向中值滤波器。

背景技术

在地震数据处理中，噪声衰减对于产生相关区域中的地震勘测 结果的更加精确的表征非常重要。地震数据的叠前道集中的信噪比的 优化处理可产生更佳的后续处理、成像和解释。

在努力抑制噪声方面已采用若干的信号处理技术。这些技术被 分类成三个主要组：频-空(f-x)域预测滤波、奇异值分解(SVD) 法和中值滤波。

在若干年前，在f-x域中采取通过预测反褶积来进行随机噪声 衰减。F-x反褶积是基于假设在每个单一频率上的空间信号是由有限 数量的复谐波所构成的。在噪声存在的情况下，自回归模型适合于预 测谐波的迭加。F-x反褶积对衰减随机噪声有效并且能够处理被称为 冲突倾角(conflicting dip)的事件。冲突倾角指的是地震数据示 出可能会出现多于一个可能的倾角的情形。然而，f-x反褶积已知会 在存在极强噪声时使信号电平显著失真。由此建议采取投影滤波器来 取代预测滤波器。

奇异值分解(SVD)是优化地震道集中的横向相关事件的另一工 具。在采用特征值分解以提取相关事件前，奇异值分解形成数据协方 差矩阵。SVC能通过只累加表示横向相关信号的最大奇异值的贡献来 有效地抑制随机噪声。扩展的SVD应用已用于f-x域中的地震数据。

中值滤波同样已经在石油和天然气工业中被广泛接受。在 Bednar,J.B.于1983年发表在Geophysics第48期第1598-1610页上 的文章“Applications of median filtering to deconvolution, pulse estimation,and statistical editing of seismic data” 中讨论了中值滤波在地震勘探中的一些应用。在Duncan,G.和 G.Beresford于1995年发表在Geophysics第60期第1157-1168页 上的文章“Some analyses of2D median f-k filters”中介绍了 2D中值f-k滤波器，其使用f—k滤波器的截短的脉冲响应系数作为 加权系数用于加权的中值处理。在Zhang,R.和T.J.Ulrych于2003 年发表在SEG扩充文摘第1949-1952页的文章“Multiple  suppression based on the migration operator and a hyperbolic  median filter”中探讨了使用双曲线中值滤波器来进行多重抑制， 而在Liu,C.、Y.Liu、B.Yang、D.Wang和J.Sun于2006年发表在 Geophysics第71期V105-V110上的文章“A2D multistage median  filter to reduce random seismic noise”中采用2D多级中值滤波 器来抑制陆地地震数据中的随机噪声。在Liu,Y.、Y.Luo和Y.Wang 于2009年发表在SEG扩充文摘29第3342-3346页的文章“Vector  median filter and its Applications in Geophysics”中提出在地 质物理中应用矢量中值滤波器(VMF)。在Huo,S.、Luo,Y.和 P.G.Kelamis于2009年发表在SEG扩充文摘28第31-35页的文章 “Simultaneous sources separation via multi-directional  vector-median filter”中讨论了将矢量中值滤波器扩展为多方向矢 量中值滤波器(MD-VMF)以分隔混合的现场地震数据。

中值滤波器处理假设已经预先对相关事件进行平整化，而 MD-VMF滤波假设在操作窗口中存在单一倾角。因此，目前已知的， 中值滤波器或多方向矢量中值滤波器方法都不能够处理具有冲突倾 角的地震道集。

当地震数据示出冲突倾角时使用被称为频率-波数(或F-K)滤 波的技术。经验显示频率-波数滤波器不是边缘保持滤波器，因此在 滤波之后具有边缘效应。F-K滤波器不能令人满意地执行尖峰或脉冲 噪声的消除。如上所记载的，中值滤波器，作为边缘保持滤波器，一 般适用于衰减尖峰噪声，但是只对平整的事件有用。

发明内容

简而言之，本发明提供了一种新的和改进的用于衰减其中时间- 空间处理结果指示有冲突倾角的地震数据中的随机噪声的处理地震 数据的计算机实施的方法。根据本发明的计算机实施的方法，将地震 数据的时间-空间处理结果在计算机中聚集以形成该地震数据的一系 列重叠的时间-空间样本窗口。将单个地震数据的样本窗口的数据转 换为频率-波数域数据。随后，将通过转化样本窗口数据所获得的频 率-波数域数据转换为傅里叶径向域数据。确定样本窗口的傅里叶径 向域数据的峰值图，对选定的傅里叶径向域数据的峰值进行转换，以 指示出样本窗口数据中的支配倾角。沿着指示出的支配倾角中的一个 所选定的支配倾角将中值滤波器施加到地震数据的样本窗口中的数 据，从而对该样本窗口中的数据进行噪声衰减。对中值滤波后的数据 的剩余信号进行存储。

本发明还提供了一种新的和改进的用于衰减其中时间-空间处 理结果指示有冲突倾角的地震数据中的随机噪声的地震数据的数据 处理系统。该数据处理系统包括存储地震数据的时间-空间处理结果 的数据存储器，以及处理器。该数据处理系统的处理器将地震数据的 时间-空间处理结果在计算机中聚集以形成该地震数据的一系列重叠 的时间-空间样本窗口，将单个地震数据的样本窗口的数据转换为频 率-波数域数据。随后，处理器将通过转化样本窗口数据所获得的频 率-波数域数据转换为傅里叶径向域数据，并确定样本窗口的傅里叶 径向域数据的峰值图。随后，处理器根据样本窗口数据中指示出的支 配倾角对选定的傅里叶径向域数据的峰值进行转换。处理器沿着所指 示出的支配倾角的一个选定的支配倾角将中值滤波器施加到样本窗 口的数据，从而对该样本窗口中的数据进行噪声衰减。处理器然后对 样本窗口的中值滤波后的数据的剩余信号进行存储。

本发明还提供了一种新的和改进的在计算机可读介质中存储计 算机可操作指令的数据存储装置，所述计算机可操作指令使数据处理 系统对其中时间-空间处理结果指示有冲突倾角的地震数据进行处理 以衰减该地震数据中的随机噪声。存储在数据存储装置中的指令引发 数据处理系统将地震数据的时间-空间处理结果在计算机中聚集以形 成该地震数据的一系列重叠的时间-空间样本窗口。该指令还引发数 据处理系统将单个地震数据的样本窗口的数据转换为频率-波数域数 据，并将通过转化样本窗口数据所获得的频率-波数域数据转换为傅 里叶径向域数据。该指令还引发数据处理系统确定样本窗口的傅里叶 径向域数据的峰值图，并根据样本窗口数据中指示出的支配倾角对选 定的傅里叶径向域数据的峰值进行转换。该指令引发数据处理系统沿 着所指示出的支配倾角中的一个选定的支配倾角将中值滤波器施加 到地震数据的样本窗口中的数据，从而对样本窗口中的数据进行噪声 衰减。随后，该指令引发数据处理系统对样本窗口的中值滤波后数据 的剩余信号进行存储。

附图说明

图1是根据本发明的在数据处理系统中执行的用于迭代倾角导 向中值滤波器的数据处理步骤组的功能框图。

图2是根据本发明的用于迭代倾角导向中值滤波器的数据处理 系统的示意框图。

图3A、图3B和图3C是在不同域中的样本平面波的变换示意图。

图4A、图4B和图4C是处理后的地震数据的绘图组。

图5A、图5B、图5C、图5D和图5E是处理后的地震数据的另一 绘图组。

图6A、图6B和图6C是处理后的地震数据的又一绘图组。

图7A、图7B、图7C、图7D和图7E是处理后的地震数据的追加 绘图组。

具体实施方式

本发明提供了在用于相关地下结构分析的地震数据中衰减随机 噪声的迭代倾角导向滤波器。如将要阐述的，在对数据执行傅里叶径 向变换之后通过在数据窗口中搜寻多个支配倾角来处理数据。之后以 降序序列施加中值滤波。通过这种方法，可以消除具有冲突倾角的数 据中的随机噪声。

由一组数据处理步骤所组成的流程图F(图1)示出了在计算机 程序软件中实现的本发明的逻辑结构。流程图F是高级逻辑流程图， 其示出了根据本发明的通过根据本发明的迭代倾角导向中值滤波来 处理地震数据的方法。应当理解的是，该流程图示出了根据本发明运 行的计算机程序码元的结构。本发明在其必要实施例中通过计算机组 件来实行，所述计算机组件以这样一种形式使用程序代码指令，即， 指示数字数据处理系统D(图2)执行对应于流程图F中所示步骤的 处理步骤序列。图1的流程图F包括针对用于随机噪声衰减的在数据 处理系统D中迭代倾角导向中值滤波地震数据的计算机实施处理数 据的优选步骤序列。

流程图F是高级逻辑流程图，其示出了根据本发明的处理方法。 本发明的方法在数据处理系统D的计算机10(图2)中执行并且可以 使用图1的存储在存储器12中且由计算机10的系统处理器14执行 的计算机程序步骤实施。如将要阐述的，流程图F示出了采用数据处 理系统D来迭代倾角导向中值滤波地震数据以衰减地震数据中的不 需要的随机噪声及用于其他目的的优选实施例。

在流程图F的步骤20(图1)期间，数据处理系统D的存储器 中的输入地震数据被分隔并聚集到若干个时间-空间重叠的小数据样 本窗口中。应当理解的是，用于根据本发明进行处理的输入地震数据 可以是叠前道集或叠后剖面。输入数据被聚集到的作为开始的第一个 时间-空间窗口以及相继的其他数据窗口中，处理如图1所示的后续 序列继续进行。

在步骤22期间，将被处理的地震数据的选定时间-空间样本窗 口从时间-空间域变换为频率-波数(f-k)域。在步骤24期间，对样 本窗口数据执行步骤22所产生的f-k域数据被变换为傅里叶径向域， 优选地使用公知的环绕径向道变换。

在步骤26期间，通过累加从步骤24产生的傅里叶径向变换数 据中所示出的每个径向形成峰值图。在步骤28中，评估在步骤26 期间所形成的傅里叶径向域峰值图中的选定峰值以示出用于样本窗 口的原始t-x域数据中的支配倾角。按照峰值图中示出的它们各自的 峰值来评估倾角的支配度。

在步骤30期间，之后沿着在步骤28期间示出的支配倾角的第 一选定支配倾角的角度施加中值滤波器以衰减来自信号的噪声。以在 峰值图中示出的不同倾角的各自的峰值的降序序列以不同的示出倾 角迭代地执行中值滤波步骤30。由于通常存在多个支配倾角，本发 明使用以示出峰值的降序序列按照倾角进行滤波的迭代方法。对于中 值滤波的每个处理迭代，作为步骤30的一部分所施加的中值滤波器 被沿着在剩余倾角中具有最大示出峰值的剩余倾角的对应角度施加。 处理针对每个选定倾角继续进行。在步骤30期间，在每个迭代之后， 将前面预估的信号数据从输入数据中减去，剩余或残留的数据被用作 在步骤30期间的后续中值滤波的样本窗口数据的输入数据。残留数 据同样存储在计算机存储器中。当沿着由峰值图示出的每个选定支配 倾角的中值滤波结束时，剩余的噪声衰减后的信号被累加以形成用于 样本数据窗口的最终输出值。结果是当前被处理的样本窗口的噪声衰 减后的信号的倾角导向中值滤波最终输出。针对样本窗口的在中值滤 波步骤30期间的迭代数量取决于数据的复杂性(即，存在于样本窗 口的数据中的示出冲突倾角的数量)。

之后在步骤32期间，将步骤30产生的针对当前被处理的样本 窗口的中值滤波后的数据存储到数据处理系统D的存储器中。在步骤 34期间，确认输入数据中的每个样本窗口是否都已经被处理。如果 否，则在步骤36期间增加数据窗口序列或计数到新的相关数据窗口 并且处理回归到步骤22以继续以上述方式进行处理。针对连续样本 窗口的处理数据的序列以这样的方式继续直到相关的输入时间-空间 数据的每个数据样本窗口被处理。

当在步骤34期间确定该处理完成时，步骤38使得每个时间-空 间数据样本窗口的信号输出被从数据处理系统D的存储器中聚合。步 骤38通过显示器26产生倾角导向中值滤波后的噪声衰减的地震数据 的输出地震显示用于分析、解释和评估。

基于下面的分析，可进一步理解图1中所示出及上面所描述的 处理序列。在t-x域中平面波d(x，t)可以表达为

$d(x,t)=w\left(t\right)\otimes \delta (t-t_{\theta }-\mathrm{px}),---\left(1\right)$

其中小波w(t)与具有倾角p和截断t₀的线性事件卷积。图3A是数据 样本时间窗口中的多个平面波在t-x域中的表示。如可见，平面波 组40和42具有冲突倾角。

变换方程(1)到f-x域得到：

$D(x,\omega )=W\left(\omega \right)e^{-\mathrm{t\omega px}}{e}^{-\mathrm{tw}{t}_0}---\left(2\right)$

其中ω表示时间频率。由方程(2)所表示的图3A的平面波之后在 步骤22期间被转移到f-k域，得到：

$D(k_x,\omega )=W\left(\omega \right)e^{-\mathrm{tw}{t}_0}\delta (k_x+\mathrm{p\omega })---\left(3\right)$

其中k_x表示波数。图3B是基于在图3A中所示的f-x域中的平面波 在频率-波数域或f-k域中的平面波表示。

从f-k笛卡尔坐标系移动到r-θ极坐标系，得到：

$D(r,\theta )=W\left(r\cos \theta \right)e^{-t{t}_{0}r\cos \theta }\delta (\mathrm{rstn\theta }+\mathrm{pr}\cos \theta )---\left(4\right)$

其中ω＝rcosθ、k_x＝rsinθ。在方程式(4)中孤立θ，在步骤24期 间得到傅里叶径向变换：

$D(t,\theta )=W\left(r\cos \theta \right)e^{-t{t}_{0}r\cos \theta }\delta (\theta -\frac{\pi }{2}-{\tan }^{-1}p)---\left(5\right)$

图3C是基于图3A和图3B中所示的平面波的在傅里叶径向域或 r-θ域中的平面波的表示。

如在图3C中所见，在傅里叶径向域中，数据中的原始t-x事件 转移到θ＝tan^-1p处的线条，其不根据r变化。

之后在步骤26期间通过累加每个示出θ的能量来产生支配倾角 图。

$M\left(\theta \right)={\Sigma }_{n=1}^{N_r}\left|\right|D(r_N,\theta )\left|\right|---\left(6\right)$

在函数M(θ)中的峰值是在t-x域中具有支配能量的倾角或斜坡 的指示器。当峰值大于其两个相邻数值时选择该峰值。可以识别多个 峰值，但是只需要若干个最大值。峰值的数量可以根据在t-x窗口中 的冲突倾角的数量来决定。

在步骤28期间在傅里叶径向域中识别了M(θ)的峰值之后，在 t-x域中的对应倾角被确定。这些倾角之后被用作在步骤30期间在 t-x域引导中值滤波来预估信号数据。

通常，中值滤波的标量中值可以通过排序来确定。给定标量组 {x_i|i＝1,…,N}，该组可以以升序或降序来排序，排序后位于中间的 倾角数值可以被预估或设置为中值。

滤波器还可以定义为从中值到组中其他成员的累积距离的最小 值。可以如下描述：

$D\left(x_j\right)=\underset{i=1}{\overset{N}{\Sigma }}{\left|\right|x_j-x_i\left|\right|}_l,x_j\in \left\{x_i\right|i=1,...,N\}---\left(7\right)$

方程式(7)定义了从倾角数值组的一个成员到该组中剩余成员 的距离之和，其中j＝1,…,N而l表示范数排序。中值标量x_m则可以从 满足等式(8)的组中选取。

$x_m=\underset{x_j}{\arg \min }D\left(x_j\right)---\left(8\right)$

如果使用L1范数(l＝1)，则使用两种方法所定义的中值数值 是相同的。然而，如果非-L1范数被采用，则基于排序或最小距离确 定的中值数值可能不同。例如，对于由{1,2,3,4,1000}所组成的 一组数字，根据L1-范数最小距离或排序确定的中值数值是3，而基 于L-2范数最小距离定义的中值数值是4。应当认识到的是，由于范 数排序根据具体应用可以调整，在实践中最小距离定义可以比基于排 序的定义更加灵活(Astola等人于1990年4月发表在IEEE会报第 78卷第678-689页的文章“Vector Median Filters”)。

如图2所示，数据处理系统D包括具有处理器12和与处理器12 耦接以存储操作指令、控制信息以及数据库记录于其中的存储器14 的计算机10。如果需要，计算机10可以是具有节点的多核处理器， 例如来自HP、Intel公司或先进微电子器件(AMD)的多核处理器。 计算机20还可以是任意传统类型的具有适当处理能力的大型计算 机，例如可从N.Y,Armonk，国际商用机器(IBM)或其他来源获得 的那些大型计算机。如下记载，同样也可以使用其他数字处理器。

应当注意到的是，可以使用其他数字处理器，例如膝上型计算 机、笔记本型计算机形式的个人计算机，或其他合适的程序化的或可 编程的数字数据处理设备。

处理器10的形式一般是具有用户界面44和用于显示根据本发 明的地震数据测量的处理的输出数据或记录的输出显示46的计算 机。输出显示46包括诸如打印机和输出显示屏幕之类的能够以图表、 数据表、图形图像、数据绘图等形式提供打印输出信息或可视化显示 作为输出记录或图像的组件。

计算机10的用户界面44还包括合适的用户输入装置或输入/输 出控制单元48用于提供用户访问控制或访问信息和数据库记录以及 操作计算机10。数据处理系统D还包括存储在存储器(可以是内部存 储器12，或者如在相关数据库服务器56中的54处所示的外部的、 网络或非网络的存储器)中的数据库50。

数据处理系统D包括存储在计算机10的存储器12中的程序代 码60。如前面所述，根据本发明的程序代码60以计算机可操作指令 的形式存在，该计算机可操作指令引发数据处理器14形成用于分析、 解释和评估的倾角导向中值滤波后的地震数据。

应当注意的是程序代码60的形式可以是提供控制数据处理系统 D的功能并指示其操作的特定有序操作集的微码、程序、例行程序或 符号化计算机可操作语言。程序代码60的指令可以存储在计算机10 的存储器12中，或者存储在计算机磁盘、磁带、传统硬盘、电只读 存储器、光存储装置中，或者存储在其他适当的具有非易失性计算机 可使用介质存储于其上的数据存储装置上。程序代码60还可以作为 计算机可读介质包含在数据存储装置上，例如所示的服务器56。

考虑如图4A所示的由多个冲突倾角所组成的噪声合成炮集。观 察该数据，可以识别信号的五个不同的倾角。然而，如果在数据上的 滑坡数据被定义为25道和50时间样本，则在每个窗口只找到两个冲 突倾角。因此，根据本发明的具有5点中值滤波的两个迭代被用于处 理。图4B描述了根据本发明的滤波结果，其中噪声被很好地衰减。 图4C示出了被消除的能量。可以观察到不存在被错误抑制的相关事 件(信号)。

图5A示出了具有间断点的噪声合成道集。图5B和图5D分别是 根据本发明的倾角导向中值滤波器和传统的f-x预测反褶积的滤波 结果。虽然可以看出两种方法都将随机噪声衰减到了不同程度，但是 图5B的倾角导向中值滤波器数据比图5D的数据更好地保留了间断 点。类似结论可以通过比较图5C和图5E的数据显示而得出。图5C 是图5A和图5B之间的数据内容的差异图，而图5E是图5A和图5D 之间的数据内容的差异图。

图6A是实际陆地地震数据的样本炮集。图6B展示了根据本发 明的3-迭代倾角导向中值滤波器的滤波结果，图6C是衰减的能量的 绘图。在图6B中可看到信噪比被很好地优化，而在图6C中可观察到 很少的信号泄露。

图7A到图7D是进一步示出根据本发明的倾角导向中值滤波器 和传统f-x预测反褶积之间的比较。图7A是极端噪声炮集。图7B 是通过根据本发明的倾角导向中值滤波器对图7A中所示的数据进行 滤波之后剩余信号的显示。图7D是通过现有技术f-x预测反褶积对 图7A的数据处理之后的剩余信号的显示。两个结果都是每个滤波器 的一系列实验参数中最佳选择。可以观察到倾角导向中值滤波比f-x 预测反褶积保留了更多的信号。图7C是图7A和图7B之间的数据内 容的差异图，而图7E是图7A和图7D之间差异图。

对于陆地地震数据处理而言，获得更好的信噪比一直是个挑战。 图6A和图7A示出了低S/N比地震道集的示例。尽管传统中值滤波器 可以消除如图6A和图7A中所示的那样的随机噪声，但是其对倾角事 件，特别是具有冲突倾角的事件执行不佳。本发明解决了这些问题并 优化了数据中的信噪比，同时很好地保留了信号的有意义的信息内 容。

如上所述，可以看到本发明提供了频率-波形滤波器和采用随后 的中值滤波器的地震数据的傅里叶径向变换。因此，本发明避免了在 地震数据中传统地使横向事件平整的中值滤波效应。通过在累加期间 沿着示出的支配倾角在数据中执行迭代，本发明能够处理数据中存在 的冲突倾角，而目前已知的其他去噪声方法则不能处理。

本发明已经充分的描述，使得本领域的普通技术人员可以复现 并获得本发明于此所提及的结果。尽管如此，本领域的任何技术人员 根据本文中的本发明可能进行未在此描述的修改，将这些修改应用到 确定的处理方法，或使用其结果，要求在下列权利要求中的请求主题； 这类修改应该由本发明的范围覆盖。

应当注意和认识到的是，可以在不脱离如所附权利要求所阐述 的本发明的精神或范围的前提下对本发明的上述细节进行改进和修 改。