CN110515128B - 基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法 - Google Patents

Abstract

基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法属地震勘探环境噪声的去噪技术领域，本发明包括：计算地震勘探数据的Kendall空间秩相关系数；利用相关噪声与同相轴空间持续性的不同，获得不同的Kendall空间秩相关系数；基于Kendall空间秩相关系数的斜坡保持复扩散噪声压制。本发明能有效解决斜坡保持复扩散方法去噪结果中残留大量相关噪声的问题，此外，利用空间Kendall秩相关系数构建的自适应扩散系数，可对背景噪声和有效信号进行不同程度的扩散，进一步改善地震勘探数据的去噪质量；模拟数据和真实沙漠地震勘探数据的实验结果均表明：本发明对沙漠地震勘探噪声具有明显的抑制作用，能有效去除相关噪声。

Description

基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法

技术领域

本发明属于地震勘探环境噪声的去噪技术领域，具体涉及基于地震勘探环境噪声中空间秩相关系数的复扩散去噪方法。

背景技术

随机噪声特性是地震数据处理和分析的关键因素。由于沙漠环境复杂多变，沙漠随机噪声呈现出特定的性质，如非高斯分布，非平稳和非线性。伴随着高频噪声的存在，出现了更多的低频噪声。此外，地表层低速区地震波的衰减是沙漠地区地震勘探数据分辨率低的重要因素。表面介质的不均匀性，沙漠地区表层沙丘的多样性，起伏的地形，使得相邻物理点之间的介质结构信息差异很大，这增加了信号衰减的复杂性。因此，沙漠地区的地震信号通常呈现低频和非平稳，具有与随机噪声相似的时间相关性和频率带宽。由于这些原因，传统的噪声抑制方法不再适用，阻碍了地质信息的采集以及去噪算法的研究，这对接下来的地震资料的解释带来了许多难题。

目前，已经提出了许多地震噪声抑制方法，例如f-x，小波阈值处理，时频峰值滤波(TFPF)，经验模态分解(EMD)，非局部均值(NLM)滤波等。这些方法可以对高频率噪声产生很好的抑制作用，例如白高斯噪声。而近年来，抑制沙漠低频噪声的研究主要是基于卷积神经网络(DnCNN)等方法，但这种方法在训练过程中经常需要花费大量时间。并且，在信号采集过程中，由于有效信号和噪声的频率相对接近，噪声与噪声、噪声和有效信号之间的相似性增加，各种变换域中的有效信号和随机噪声之间的差异被削弱。基于这些特性，噪声抑制效果减弱。虽然斜坡保持复扩散可以对沙漠噪声产生一定的抑制作用，但是这种方法在去噪的过程中，由于其利用的是梯度的不同改变不同位置处的扩散强度，而沙漠噪声中的相关噪声又具有与同向轴相似的性质，使得在去噪后会留下具有相关特性的噪声。

上述分析表明，沙漠噪声与同向轴具有某些相似的性质，会限制目前已有的利用相关性进行去噪的噪声压制技术，导致去噪后留下大量的相关噪声。因此，需要找到一种能够区分信号与相关噪声的方法，在对同向轴进行保留的同时，对相关噪声进行有效的抑制。

发明内容

本发明的目的在于提供一种能有效抑制噪声的基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法。

本发明的基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法，包括下列步骤：

1.1计算地震勘探数据的Kendall空间秩相关系数：

利用检波器对二维地震勘探数据进行采集，假设在采集到的二维地震勘探数据D中，有N条地震道，每道有M个样本点，取一个分析窗X，以当前样本点d_ij为原点，沿时间方向上取m个样本点作为窗的长度，沿着地震道方向取n条地震道作为窗的宽度，则窗X内数据为m×n大小的一个窗。以当前样本点水平方向上相邻的样本点d_i(j+1)为原点，以相同方式取与窗X相邻，大小为m×n的另一个窗作为窗Y的数据。

将两个窗X和Y中对应位置上的元素组成一组数据对(X_p,Y_p)，利用以上两个窗内的数据计算Kendall空间秩相关系数，为了方便计算将矩阵中的元素代入到公式中：

式中：

其中：sgn(·)是符号函数，τ_ij表示第i行第j列元素对应的Kendall空间秩相关系数，当滑动窗遍历二维矩阵d的全部数据后，得到相关系数矩阵τ，Kendall空间秩相关系数矩阵中的每个相关系数τ_ij表示以d_ij为原点的相邻窗X和Y的Kendall空间秩相关程度。

相关系数矩阵中的每一个元素，代表这个位置周围邻域的数据在空间上的相关程度的强弱，当相关系数为正，说明这部分地震数据具有空间正相关的性质，其空间变化趋势一致；当相关系数为负，说明这部分地震数据具有空间负相关的性质，其空间变化趋势不一致；当空间相关系数为0，说明数据之间呈现相互独立的状态。

1.2空间相关噪声与信号的区分

在计算Kendall空间秩相关系数的时候，根据单位周期内信号的时间长度选择窗X和窗Y的长度m，使得窗内的数据可以表征该段时间内数据的变化趋势；窗的宽度n分别取整数2≤n₁,n₂…,n_N≤10。

在利用不同宽度的窗对二维地震勘探数据计算Kendall空间秩相关系数的时候，固定其中一个窗X,另一个窗Y在自下而上一定区域内滑动。若在某个角度上的Kendall空间秩相关系数的值较大，且几乎不随着窗的宽度n的改变而改变，则说明该区域的数据为同相轴；若Kendall空间秩相关系数的值在某个角度上随着窗宽n的改变变化较大，则说明该区域的数据为相关噪声；若Kendall空间秩相关系数的值在各个角度上都比较小，而且不随窗宽n的变化而改变，说明这部分区域的数据为随机噪声。

根据Kendall空间秩相关系数的取值，选择窗的宽度n_i，使得同相轴与相关噪声计算出来的Kendall空间秩相关系数具有较大的差异。这时利用长度为m宽度为n_i的窗计算出Kendall空间秩相关系数。绘制Kendall空间秩相关系数的概率密度分布曲线，找到概率分布密度曲线中的波谷，波谷对应的Kendall空间秩相关系数就是区分同向轴和相关噪声的阈值。由于随机噪声的Kendall空间秩相关系数较小，所以在利用阈值过滤相关噪声的同时，也会有效的将随机噪声过滤掉。

1.3基于Kendall空间秩相关系数的斜坡保持去噪算法的噪声压制：

利用Kendall空间秩相关系数改进斜坡保持去噪算法中的扩散系数C(Im(I))的数学表达式为：

其中：I是地震勘探数据；Im(·)表示虚部因子；θ为角度参数；K＝e^-aλ(τ,Th)是一个与阈值函数λ(τ,Th)有关的函数，阈值函数的数学表达式为：

其中：a为实参数；τ为Kendall空间秩相关系数；Th是用来区分信号和噪声的阈值；

将改进后的扩散系数带入到斜坡保持去噪算法中，得到最终的去噪结果。

本发明通过对信号和噪声区域进行不同程度的扩散来抑制噪声，利用Kendall相关系数改变扩散强度，使得基于Kendall相关系数的复扩散方法可以有效的抑制具有相关性的背景噪声。

本发明利用噪声与同向轴空间连续性的不同，求得不同的Kendall空间秩相关系数。利用Kendall空间秩相关系数对同向轴与噪声加以区分。在利用斜坡保持复扩散方法进行噪声压制的过程中，利用相关系数作为区分信号与噪声的阈值，对信号与噪声部分进行不同程度的扩散去噪，使得去噪结果中，噪声部分被有效抑制，增强了对同向轴的保持能力。

本发明所解决的技术问题：斜坡保持复扩散可以有效的对部分沙漠噪声进行抑制，但是部分具有相关特性的噪声会在去噪结果中保留下来，导致沙漠噪声中的部分噪声不能被有效去除。针对这一问题，本发明提出利用Kendall空间秩相关系数对信号以及噪声部分加以区分，并构造不同大小的扩散系数，使得信号和噪声能够进行不同程度的扩散。这样，这保留同向轴的同时，对相关噪声也会产生有效的抑制作用，使得去噪结果中的相关噪声大大减少。

本发明的有益效果在于：研究斜坡保持复扩散在有效的对部分沙漠噪声进行抑制的同时，保留了部分具有相关特性的噪声，导致沙漠噪声中的部分噪声不能被有效去除等问题。针对这一问题，本发明提出利用Kendall空间秩相关系数对信号以及噪声部分加以区分，并构造不同大小的扩散系数，使得信号和噪声能够进行不同程度的扩散。这样，在保留同向轴的同时，对相关噪声也会产生有效的抑制作用，使得去噪结果中的相关噪声大大减少。本发明利用了Kendall空间秩相关系数具有衡量空间数据依赖程度的优点，对斜坡保持复扩散方法的扩散系数进行改进，以达到对相关噪声进行抑制，对同向轴进行保留的目的。

附图说明

图1为基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法的流程图

图2为模拟的纯净地震勘探记录图

图3为加沙漠噪声的模拟地震勘探记录图

图4为计算得到的Kendall空间秩相关系数矩阵

图5为利用不同宽度的窗计算同相轴位置处的Kendall空间秩相关系数图

图6为利用不同宽度的窗计算噪声位置处的Kendall空间秩相关系数图

图7为含噪数据的Kendall空间秩相关系数的概率密度分布图

图8为利用斜坡保持复扩散进行滤波的结果示意图

图9为基于Kendall空间秩相关系数的斜坡复扩散滤波结果示意图

图10为第41道数据利用斜坡保持复扩散滤波对比图

图11为第41道利用改进后的斜坡保持复扩散滤波结果对比图

图12为纯净信号的f-k谱

图13为含噪信号的f-k谱

图14为利用斜坡保持复扩散滤波后的f-k谱

图15为利用改进斜坡保持复扩散滤波后的f-k谱

图16为真实沙漠地震勘探噪声数据

图17为利用改进后的斜坡保持复扩散滤波后的结果图

图18为利用斜坡保持复扩散滤波后的结果图

具体实施方式

本发明的基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法，包括下列步骤：

1.1计算地震勘探数据的Kendall空间秩相关系数：

利用检波器对二维地震勘探数据进行采集，假设在采集到的二维地震勘探数据D中，有N条地震道，每道有M个样本点，取一个分析窗X，以当前样本点d_ij为原点，沿时间方向上取m个样本点作为窗的长度，沿着地震道方向取n条地震道作为窗的宽度，则窗X内数据为m×n大小的一个窗。以当前样本点水平方向上相邻的样本点d_i(j+1)为原点，以相同方式取与窗X相邻，大小为m×n的另一个窗作为窗Y的数据。

将两个窗X和Y中对应位置上的元素组成一组数据对(X_p,Y_p)，利用以上两个窗内的数据计算Kendall空间秩相关系数，为了方便计算将矩阵中的元素代入到公式中：

式中：

其中：sgn(·)是符号函数，τ_ij表示第i行第j列元素对应的Kendall空间秩相关系数，当滑动窗遍历二维矩阵d的全部数据后，得到相关系数矩阵τ，Kendall空间秩相关系数矩阵中的每个相关系数τ_ij表示以d_ij为原点的相邻窗X和Y的Kendall空间秩相关程度。

相关系数矩阵中的每一个元素，代表这个位置周围邻域的数据在空间上的相关程度的强弱，当相关系数为正，说明这部分地震数据具有空间正相关的性质，其空间变化趋势一致；当相关系数为负，说明这部分地震数据具有空间负相关的性质，其空间变化趋势不一致；当空间相关系数为0，说明数据之间呈现相互独立的状态。

1.2空间相关噪声与信号的区分

在计算Kendall空间秩相关系数的时候，根据单位周期内信号的时间长度选择窗X和窗Y的长度m，使得窗内的数据可以表征该段时间内数据的变化趋势；窗的宽度n分别取整数2≤n₁,n₂…,n_N≤10。

在利用不同宽度的窗对二维地震勘探数据计算Kendall空间秩相关系数的时候，固定其中一个窗X,另一个窗Y在自下而上一定区域内滑动。若在某个角度上的Kendall空间秩相关系数的值较大，且几乎不随着窗的宽度n的改变而改变，则说明该区域的数据为同相轴；若Kendall空间秩相关系数的值在某个角度上随着窗宽n的改变变化较大，则说明该区域的数据为相关噪声；若Kendall空间秩相关系数的值在各个角度上都比较小，而且不随窗宽n的变化而改变，说明这部分区域的数据为随机噪声。

根据Kendall空间秩相关系数的取值，选择窗的宽度n_i，使得同相轴与相关噪声计算出来的Kendall空间秩相关系数具有较大的差异。这时利用长度为m宽度为n_i的窗计算出Kendall空间秩相关系数。绘制Kendall空间秩相关系数的概率密度分布曲线，找到概率分布密度曲线中的波谷，波谷对应的Kendall空间秩相关系数就是区分同向轴和相关噪声的阈值。由于随机噪声的Kendall空间秩相关系数较小，所以在利用阈值过滤相关噪声的同时，也会有效的将随机噪声过滤掉。

1.3基于Kendall空间秩相关系数的斜坡保持去噪算法的噪声压制：

利用Kendall空间秩相关系数改进斜坡保持去噪算法中的扩散系数C(Im(I))的数学表达式为：

其中：I是地震勘探数据；Im(·)表示虚部因子；θ为角度参数；K＝e^-aλ(τ,Th)是一个与阈值函数λ(τ,Th)有关的函数，阈值函数的数学表达式为：

其中：a为实参数；τ为Kendall空间秩相关系数；Th是用来区分信号和噪声的阈值；

将改进后的扩散系数带入到斜坡保持去噪算法中，得到最终的去噪结果。

实施例

将本发明应用于模拟沙漠地区的含噪数据噪声压制：

图2是模拟出来的纯净的地震勘探数据，道数为88道，每道长度为2s。其中，同向轴的主频是20hz，起始到达时间为0.4s、0.8s和1.4s，速度分别为3500m/s、4000m/s和2600m/s，

图3是在图2中加入了模拟的沙漠噪声数据后得到的含噪数据。

结果分析：图4是利用长度为30个数据点，宽度为7条地震道的窗计算出来的Kendall空间秩相关系数。从相关系数矩阵中可以清晰的看出，同向轴位置处的相关系数，明显大于噪声位置处的相关系数。

图5是对信号位置处进行Kendall空间秩相关系数计算的结果，可以看出，当窗的宽度在5、7、9时，同相轴位置处的相关系数基本不随着穿宽度的改变而变化。

图6是对噪声位置处进行Kendall空间秩相关系数计算的结果，可以看出，当窗的宽度在5、7、9时，噪声位置处的相关系数随着穿宽度的增加而变小，宽度为7时可以与信号处的Kendall空间秩相关系数有明显差异。

图7是含噪数据经过计算Kendall空间秩相关系数后，统计的概率密度分布，可以看出在0.3位置处左右产生一个波谷，所以我们将0.3作为一个阈值来区分信号和噪声。

图8是利用斜坡保持复扩散对模拟沙漠噪声记录进行扩散去噪的结果，从去噪结果中可以看出，在大量噪声被有效压制的同时，留下了大量的相关噪声。

图9是利用基于Kendall空间相关系数的斜坡保持复扩散滤波方法对模拟含噪数据进行滤波的结果，可以看到有大量的相关噪声被有效的抑制，证明了本方法的有效性。

在给出整体去噪结果图的基础上，我们给出了两种方法的单道对比图。选取图8和图9的第41道数据，如图10和图11所示，从单道对比实验中也可以清晰的看出，当加入了Kendall空间相关系数后，会有更多的相关噪声被抑制。同时，我们给出了纯净信号与含噪数据的f-k谱(图12和图13)，以及去噪后(图8和图9)的f-k谱(图14和图15)，从f-k谱中也可以看出，当加入Kendall空间相关系数后，噪声压制效果会得到有效的增强。

将本发明应用到实际沙漠数据的去噪实验中。实际沙漠地震勘探数据如图16所示，图17为基于Kendall空间相关系数的斜坡保持复扩散去噪结果，图18为斜坡保持复扩散去噪结果。

结果分析：从实验结果图17和图18的对比中可以看出，原始的斜坡保持复扩散滤波方法对实际沙漠数据中的噪声具有一定的抑制作用，但是本发明所提出的方法获得的同向轴更加的清晰，对噪声的抑制作用效果更好。可以证明本发明在抑制沙漠地震勘探噪声上可以取得更好的效果。

Claims (1)

1.一种基于地震勘探环境噪声空间秩相关系数的复扩散去噪方法，其特征在于包括下列步骤：

1.1计算地震勘探数据的Kendall空间秩相关系数：

利用检波器对二维地震勘探数据进行采集，假设在采集到的二维地震勘探数据D中，有N条地震道，每道有M个样本点，取一个分析窗X，以当前样本点d_ij为原点，沿时间方向上取m个样本点作为窗的长度，沿着地震道方向取n条地震道作为窗的宽度，则窗X内数据为m×n大小的一个窗；以当前样本点水平方向上相邻的样本点d_i(j+1)为原点，以相同方式取与窗X相邻，大小为m×n的另一个窗作为窗Y的数据；

将两个窗X和Y中对应位置上的元素组成一组数据对(X_p,Y_p)，利用以上两个窗内的数据计算Kendall空间秩相关系数，为了方便计算将矩阵中的元素代入到公式中：

式中：

其中：sgn(·)是符号函数，τ_ij表示第i行第j列元素对应的Kendall空间秩相关系数，当滑动窗遍历二维矩阵d的全部数据后，得到相关系数矩阵τ，Kendall空间秩相关系数矩阵中的每个相关系数τ_ij表示以d_ij为原点的相邻窗X和Y的Kendall空间秩相关程度；

相关系数矩阵中的每一个元素，代表这个位置周围邻域的数据在空间上的相关程度的强弱，当相关系数为正，说明这部分地震数据具有空间正相关的性质，其空间变化趋势一致；当相关系数为负，说明这部分地震数据具有空间负相关的性质，其空间变化趋势不一致；当空间相关系数为0，说明数据之间呈现相互独立的状态；

1.2空间相关噪声与信号的区分

在计算Kendall空间秩相关系数的时候，根据单位周期内信号的时间长度选择窗X和窗Y的长度m，使得窗内的数据可以表征该段时间内数据的变化趋势；窗的宽度n分别取整数2≤n₁,n₂…,n_N≤10；

在利用不同宽度的窗对二维地震勘探数据计算Kendall空间秩相关系数的时候，固定其中一个窗X,另一个窗Y在自下而上一定区域内滑动；若在某个角度上的Kendall空间秩相关系数的值较大，且几乎不随着窗的宽度n的改变而改变，则说明该区域的数据为同相轴；若Kendall空间秩相关系数的值在某个角度上随着窗宽n的改变变化较大，则说明该区域的数据为相关噪声；若Kendall空间秩相关系数的值在各个角度上都比较小，而且不随窗宽n的变化而改变，说明这部分区域的数据为随机噪声；

根据Kendall空间秩相关系数的取值，选择窗的宽度n_i，使得同相轴与相关噪声计算出来的Kendall空间秩相关系数具有较大的差异；这时利用长度为m宽度为n_i的窗计算出Kendall空间秩相关系数；绘制Kendall空间秩相关系数的概率密度分布曲线，找到概率分布密度曲线中的波谷，波谷对应的Kendall空间秩相关系数就是区分同向轴和相关噪声的阈值；由于随机噪声的Kendall空间秩相关系数较小，所以在利用阈值过滤相关噪声的同时，也会有效的将随机噪声过滤掉；

1.3基于Kendall空间秩相关系数的斜坡保持去噪算法的噪声压制：

利用Kendall空间秩相关系数改进斜坡保持去噪算法中的扩散系数C(Im(I))的数学表达式为：

其中：I是地震勘探数据；Im(·)表示虚部因子；θ为角度参数；K＝e^-aλ(τ,Th)是一个与阈值函数λ(τ,Th)有关的函数，阈值函数的数学表达式为：

其中：a为实参数；τ为Kendall空间秩相关系数；Th是用来区分信号和噪声的阈值；

将改进后的扩散系数带入到斜坡保持去噪算法中，得到最终的去噪结果。

Images