CN109100788B - 地震数据非局部均值去噪方法 - Google Patents

Abstract

本发明提供的地震数据非局部均值去噪方法，获取待检测地震数据中的待检测数据点；计算待检测数据点与邻域数据点的振幅差，计算待检测数据点的SUSAN面积；利用SUSAN面积计算待检测数据点的边缘响应；重复上述步骤，计算所有数据点的边缘响应，得到待检测地震数据的边缘响应；划分待检测地震数据中的优化区域和常规区域；采用非局部均值去噪方法对常规区域进行去噪处理；对优化区域进行优化去噪。能够利用SUSAN边缘检测，提取出数据结构的边缘信息，进而从数据中划分出优化区域。优化区域内，独立进行数据的相似性计算，能够提升优化区域内数据之间的加权效果，弥补了常规非局部均值去噪方法在数据的结构边缘区域即优化区域去噪不彻底的不足。

Description

地震数据非局部均值去噪方法

技术领域

本发明属于数据处理技术领域，具体涉及地震数据非局部均值去噪方法。

背景技术

噪声的压制是地震数据处理流程中的重要组成部分。为能够最大限度的降低噪声影响，提高地震数据数据处理精度，首先需要明确噪声的类型。地震信号中的噪声通常可划分为相干噪声和非相干噪声。

地震数据的噪声压制方法通常是根据数据中有效信号和噪声在一定条件下特性的不同实现信号分离，进而达到噪声压制的目的。例如带通滤波，f-k和k_x-k_y滤波等方法是通过对地震数据进行傅里叶变换后，根据有效信号和噪声信号在频率-波数域或波数域中能量分布形态的不同实现信噪分离。

非局部均值去噪算法最初被应用于图像处理领域，该算法可看做是一种随机噪声的衰减滤波器，结合图像的自相似性质，引入权重因子，根据图像中构造相似性的强弱进行加权计算，压制随机噪声。由于该方法对任意数值点去噪时，均需与数据内所有点进行相似性计算，会产生巨大的计算量。因此，诸多学者对该方法的优化策略主要集中于降低算法的计算耗时，Mahmoudi等在实现非局部均值去噪算法的过程中，提出选取目标点附近区域进行相似性计算的策略，有效降低了计算量；Sheng等根据非局部均值去噪算法高度并行化的特点，提出基于GPU加速的非局部均值去噪算法，提升了算法的计算效率；张权等结合SUSAN(Smallest Univalue Segment Assimilating Nucleus，最小同值分割吸收核)边缘检测结果，优化非局部均值去噪算法参数，明显改善了去噪效果。目前非局部均值去噪算法在医学影像、雷达数据和音频信号等方面也均取得了较为成功的应用。

但是由于地震数据中的有效信号富含连续的同相轴，相似性信息丰富，而目前常规的非局部均值去噪方法在选取相似性数据时，仅是对全局数据或目标点周围数据的一次遍历选取，并未充分发挥数据的自相似性结构，噪声压制效果不佳。

发明内容

针对现有技术中的缺陷，本发明提供一种地震数据非局部均值去噪方法，充分考虑了地震数据中相似数据的比重，提高了去噪效果。

一种地震数据非局部均值去噪方法，包括以下步骤：

S1：输入待检测地震数据，并获取待检测地震数据中的待检测数据点；

S2：计算待检测数据点与邻域数据点的振幅差，通过设置振幅阈值，计算待检测数据点的SUSAN面积；

S3：利用待检测数据点的SUSAN面积计算待检测数据点的边缘响应；

S4：重复步骤S1～S3，计算待检测地震数据中所有数据点的边缘响应，得到待检测地震数据的边缘响应；

S5：划分待检测地震数据中的优化区域和常规区域；

S6：采用非局部均值去噪方法对常规区域进行去噪处理；

S7：计算优化区域内数据点的相似度，进行优化去噪。

进一步地，所述步骤S2具体包括：

设置伪圆形模板；

利用伪圆形模板计算待检测数据点与邻域数据点的振幅差A；

A＝A(r₀)-A(r)；

其中，r₀、r分别为待检测数据点和邻域数据点的坐标，A()为振幅值；

确定待检测数据点与邻域数据点的相似度：

其中，t₁为预设的振幅阈值，c(r₀,r)＝1表示待检测数据点与邻域数据点相似，c(r₀,r)＝0表示待检测数据点与邻域数据点非相似；

计算待检测数据点的SUSAN面积n(r₀)：

进一步地，所述步骤S3具体包括：

定义SUSAN面积阈值g＝t₂n_max；

其中，t₂为比例系数，n_max为37；

计算待检测数据点的边缘响应R(r₀)：

进一步地，所述步骤S5具体包括：

对待检测地震数据的边缘响应进行平滑处理；

根据平滑后的边缘响应数据中的最大值确定阈值比例系数t₃，对待检测地震数据进行区域划分：

其中，S(r₀)为1的坐标集为优化区域，S(r₀)为0的坐标集为常规区域，R为所有待检测点的边缘响应的集合。

进一步地，所述步骤S7具体包括：

定义在有界区域I内的含噪声数据体d(i)，其中i∈I，假设含噪声数据体由有效数据s(i)和噪声数据n(i)叠加构成：

d(i)＝s(i)+n(i)；

对任意位置i的数据点进行去噪，去噪后的数据

可表示为：

式中，ω(i,j)为分别以目标点位置i和其他任意位置j为中心的局部数据体之间的相似性系数，且满足0≤ω(i,j)≤1，

Z(i)为规则化因子，

h为常数参数，D²(i,j)为高斯加权欧几里得距离的平方项，即

式中，表示高斯加权欧几里得距离的平方，a为高斯函数的标准差；{N_k}_k∈I为邻域系统，L为邻域系统内坐标的集合，l为邻域系统内任一点的坐标，G_a代表标准差为a的高斯函数；

设置窗口，对窗口内的数据进行相似性系数的计算，并采用非局部均值去噪方法去噪。

由上述技术方案可知，本发明提供的地震数据非局部均值去噪方法，利用SUSAN边缘检测，提取出数据结构的边缘信息，进而从数据中划分出优化区域。优化区域内，独立进行数据的相似性计算，能够提升优化区域内数据之间的加权效果，理论上等同于重新定义了优化区域中邻域N的形状和选取范围，弥补了常规非局部均值去噪方法在数据的结构边缘区域即优化区域去噪不彻底的不足。在数据结构变化明显的区域，本发明提出方法充分考虑了地震数据中相似数据的比重，能够获得更好的去噪效果，具有广泛的研究价值和实用前景。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍。在所有附图中，类似的元件或部分一般由类似的附图标记标识。附图中，各元件或部分并不一定按照实际的比例绘制。

图1为地震数据非局部均值去噪方法的流程图。

图2为伪圆形模板的示意图。

图3为伪圆形模板检测数据示意图。

图4为SUSAN边缘检测过程示意图。

图5为SUSAN边缘检测结果。其中a为原始数据；b为噪声数据；c为含噪声数据；d为SUSAN边缘响应数据；e为平滑后的边缘相应数据；f为优化区域分布。

图6为高斯加权欧几里得距离计算示意图。

图7为非局部均值去噪计算计算示意图。

图8为含噪声速度模型的去噪结果对比。a为非局部均值去噪法的去噪结果；b为优化非局部均值去噪法的去噪结果；c为非局部均值去噪法预测的噪声；d为优化非局部均值去噪法预测的噪声。

具体实施方式

下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只作为示例，而不能以此来限制本发明的保护范围。需要注意的是，除非另有说明，本申请使用的技术术语或者科学术语应当为本发明所属领域技术人员所理解的通常意义。

实施例一：

本申请提供的地震数据非局部均值去噪方法的流程如图1所示，主要包括两部分，数据输入后，左侧流程为：SUSAN边缘信息检测和优化去噪区提取；右侧流程为：根据提取的优化去噪区进行的非局部均值去噪的优化计算。

1、SUSAN边缘检测。

SUSAN边缘检测算法，根据待检测数据点与邻域数据点的振幅差，通过设置振幅阈值，获得待检测数据点的SUSAN面积，进而确定待检测数据点是否处于数据的结构边缘，完成边缘检测计算。

输入待检测地震数据，并获取待检测地震数据中的待检测数据点。

由于计算机数据均以矩形网格方式排列，因此设定由37个数据点构成的伪圆形模板，如图2所示。图2中，中间方块为待检测数据点，其他方块为邻域数据点，计算待检测数据点与邻域数据点的振幅差A(r0)-A(r)。r0、r分别为待检测数据点和邻域数据点的坐标，A()为振幅值。根据下式确定待检测数据点与邻域数据点的相似度：

c(r₀,r)＝1表示待检测数据点与邻域数据点相似，c(r₀,r)＝0表示待检测数据点与邻域数据点非相似。式(1)中，t₁为预设的振幅阈值，然后根据函数：

计算SUSAN面积n(r₀)。定义SUSAN面积阈值g＝t₂n_max，t₂为比例系数，通常为0.75，n_max为37，即可计算待检测数据点的边缘响应R(r₀)：

对待检测地震数据中所有数据点重复上述过程，即得到数据体的边缘响应。

为便于理解，假设伪圆形模板检测数据如图3所示，假设图中中间方形区域和外围区域的数据振幅差超过振幅差阈值t₁，分别对a-e五个点进行边缘检测，对数据套用如图3所示的伪圆形模板，图4为SUSAN边缘检测过程示意图，图4的十字即为待检测数据点位置，对应图2的中间方块。经由伪圆形模板选取的a-e点的邻域数据如图4中的第一行圆形所示，根据(2)式计算出各点的SUSAN面积，如图4中的第二行圆形所示，其中灰色区域标识的为SUSAN面积。根据(3)式判断a、b和e三点为边缘点，c、d为非边缘点。

2、优化去噪区提取。

结合SUSAN边缘检测算法提取出数据的边缘响应信息，从数据中划分出优化区域，对优化区域内的数据进行独立去噪，可充分利用优化区域内数据的相似性信息，改善去噪效果。

为便于理解，图5给出一个SUSAN边缘检测的例子。图5a为原始数据，加以如图5b所示的信噪比20dB的高斯噪声后，所得数据如图5c所示。对于含噪声数据，首先进行SUSAN边缘检测，获取数据的边缘响应信息。需要注意的是，如果相邻数据点的噪声振幅差过大，在进行SUSAN边缘检测时，会被误检测为边缘点。根据如图5c所示数据，所获得的边缘响应数据中，除明显的边缘信息外，仍存在较多孤立点。如图5d所示，这些噪点并非边缘检测算法的期望输出。因此，需对SUSAN边缘检测结果进行平滑。这样做一方面弱化孤立点的振幅值。另一方面，扩展检测出的有效边缘响应的面积，平滑后数据如图5e所示。根据平滑后的边缘响应数据中的最大值MAX(R)，设定阈值比例系数t₃，对数据进行区域划分：

设定S(r₀)为1的坐标集为优化区域，对应图5f中线条区域，S(r₀)为0的坐标集为常规区域，对应图5f中其余区域。

3、优化的地震数据非局部均值去噪原理。

定义在有界区域I内的含噪声数据体d(i)，其中i∈I，假设含噪声数据体由有效数据s(i)和噪声数据n(i)叠加构成：

d(i)＝s(i)+n(i) (5)

对任意位置i的数据点进行去噪，去噪后的数据

可表示为：

式中，ω(i,j)为分别以目标点位置i和其他任意位置j为中心的局部数据体之间的相似性系数，且满足0≤ω(i,j)≤1，

ω(i,j)的计算公式如下：

式中，Z(i)为规则化因子：

以确保

h为常数参数，用于控制式(7)、(8)中指数函数衰减趋势，亦可理解h为通过影响相似性系数ω(i,j)的变化趋势，进而控制滤波强度的参数。D²(i,j)为高斯加权欧几里得距离的平方项，即，

式中，

表示高斯加权欧几里得距离的平方，a为高斯函数的标准差。为能够量化表示以i，j为中心的局部数据体的相似性，定义邻域系统{N_k}_k∈I，通过对比N_i、N_j区域内的数据即可获得相似性信息。而对于一个给定的数据点，它的邻域范围的选取通常为正方形或正方体，这主要取决于计算数据的维度和算法需要。式(9)中L为邻域内坐标的集合，l为邻域内任意一点的坐标。G_a代表标准差为a的高斯函数，以二维数据体为例，高斯函数G_a可表示为：

式中，x₀，y₀为高斯函数的中心，坐标x，y对应式(9)中的l项。若假设N_i、N_j为7*7的正方形，则(9)式的计算可表示如图6所示。

在应用式(9)计算某点与数据内其他所有点的相似性时，会引入巨大的计算量，因此可在目标点位置i附近设置窗口，仅对窗口内部的数据进行相似性系数的计算和去噪，从而降低计算量。如图7所示，对黑色实心圆处的数据进行去噪计算，黑色虚线框为设置的窗口，假设N_i、N_j均为3*3的正方形，如图7中黑色实线框和灰色实线框所示，根据式(7)、(9)依次计算窗口内以各点为中心的邻域数据N_j与以目标点为中心的邻域数据N_i的相似性系数；然后应用(6)式，逐点进行加权计算，完成去噪。

结合本专利提出的优化方法，对优化区域内的数据去噪时，公式(9)中j的选取范围为优化区域内的所有数据。在常规区域内的数据，j的选取范围仍为以i为中心的窗口内的数据。

由此可看出，本申请方法与常规非局部均值去噪方法的不同之处在于对数据的区域划分和公式(9)中参数j的选取范围上。该方法基于数据的结构边缘信息，对数据进行区域划分，并对优化区域内数据进行独立去噪计算。这样做的原因在于，优化区域内通常对应地震数据中，振幅差异明显或同相轴连续性较强的区域，而该区域内的数据具有更强的自相似性，将j的选取限定在优化区域内更能充分利用数据中的相似信息，对非相干噪声压制效果更好。

为了验证本申请方法的去噪效果，对含噪声速度模型进行去噪测试，如图5c所示。所用参数如表1所示。图8a、b分别为非局部均值去噪法和本申请提出的优化非局部均值去噪法去噪结果，通过对比可看出，本方法对盐丘边界处的去噪效果更好，在相同色标刻度下，去噪结果更接近原始数据。

图8c、d分别为非局部均值去噪法和本申请方法去除的噪声，非局部均值去噪法去除的噪声在盐丘边界处存在较明显的痕迹，且振幅值为0，即去噪不足。本申请方法去除的噪声无明显痕迹，去除的噪声呈随机分布，更接近理论值。

表1速度模型去噪测试参数表

该方法利用SUSAN边缘检测，提取出数据结构的边缘信息，进而从数据中划分出优化区域。优化区域内，独立进行数据的相似性计算，能够提升优化区域内数据之间的加权效果，理论上等同于重新定义了优化区域中邻域N的形状和选取范围，弥补了常规非局部均值去噪方法在数据的结构边缘区域即优化区域去噪不彻底的不足。在数据结构变化明显的区域，本发明提出方法充分考虑了地震数据中相似数据的比重，能够获得更好的去噪效果，具有广泛的研究价值和实用前景。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。

Claims (4)

1.一种地震数据非局部均值去噪方法，其特征在于，包括以下步骤：

S1：输入待检测地震数据，并获取待检测地震数据中的待检测数据点；

S2：计算待检测数据点与邻域数据点的振幅差，通过设置振幅阈值，计算待检测数据点的SUSAN面积；

S3：利用待检测数据点的SUSAN面积计算待检测数据点的边缘响应；

S4：重复步骤S1～S3，计算待检测地震数据中所有数据点的边缘响应，得到待检测地震数据的边缘响应；

S5：划分待检测地震数据中的优化区域和常规区域；所述步骤S5具体包括：

对待检测地震数据的边缘响应进行平滑处理；

根据平滑后的边缘响应数据中的最大值确定阈值比例系数t₃，对待检测地震数据进行区域划分：

其中，S(r₀)为1的坐标集为优化区域，S(r₀)为0的坐标集为常规区域，R为所有待检测点的边缘响应的集合；

S6：采用非局部均值去噪方法对常规区域进行去噪处理；

S7：计算优化区域内数据点的相似度，进行优化去噪。

2.根据权利要求1所述地震数据非局部均值去噪方法，其特征在于，所述步骤S2具体包括：

设置伪圆形模板；

利用伪圆形模板计算待检测数据点与邻域数据点的振幅差A；

A＝A(r₀)-A(r)；

其中，r₀、r分别为待检测数据点和邻域数据点的坐标，A()为振幅值；

确定待检测数据点与邻域数据点的相似度：

其中，t₁为预设的振幅阈值，c(r₀,r)＝1表示待检测数据点与邻域数据点相似，c(r₀,r)＝0表示待检测数据点与邻域数据点非相似；

计算待检测数据点的SUSAN面积n(r₀)：

3.根据权利要求2所述地震数据非局部均值去噪方法，其特征在于，所述步骤S3具体包括：

定义SUSAN面积阈值g＝t₂n_max；

其中，t₂为比例系数，n_max为37；

计算待检测数据点的边缘响应R(r₀)：

4.根据权利要求1所述地震数据非局部均值去噪方法，其特征在于，所述步骤S7具体包括：

定义在有界区域I内的含噪声数据体d(i)，其中i∈I，假设含噪声数据体由有效数据s(i)和噪声数据n(i)叠加构成：

d(i)＝s(i)+n(i)；

对任意位置_i的数据点进行去噪，去噪后的数据可表示为：

式中，ω(i,j)为分别以目标点位置_i和其他任意位置j为中心的局部数据体之间的相似性系数，且满足0≤ω(i,j)≤1，

Z(i)为规则化因子，

h为常数参数，D²(i,j)为高斯加权欧几里得距离的平方项，即式中，

表示高斯加权欧几里得距离的平方，a为高斯函数的标准差；{N_k}_k∈I为邻域系统，L为邻域系统内坐标的集合，l为邻域系统内任一点的坐标，G_a代表标准差为a的高斯函数；

设置窗口，对窗口内的数据进行相似性系数的计算，并采用非局部均值去噪方法去噪。

Images