CN111694056A - 一种压制地震资料异常噪声的方法、存储介质及设备 - Google Patents

Abstract

本发明公开了一种压制地震资料异常噪声的方法、存储介质及设备，读取含有高振幅异常噪声的原始地震资料，根据地震资料中噪声的复杂程度选取混合高斯分布的高斯成分个数和迭代次数；针对原始地震资料中有效信号的特点选取固定稀疏变换基，利用线性化Bregman方法得到地震资料中有效信号的稀疏表示，根据稀疏表示与固定稀疏变换基的乘积和原始地震资料得到关于有效信号的估计和异常噪声的估计，将异常噪声建模为混合高斯分布，利用EM算法更新每个采样点上异常噪声的权重矩阵；直到原始地震资料中异常噪声压制完成或者达到迭代次数。本发明的压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法，可以有效压制异常噪声，提高地震资料信噪比。

Description

一种压制地震资料异常噪声的方法、存储介质及设备

技术领域

本发明属于地震勘探数据处理技术领域，具体涉及一种压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法、存储介质及设备。

背景技术

按照噪声在地震剖面上出现的特征，将其分为随机噪声和相关噪声。地震资料成像过程中，检波器接收到的反射信号含有大量的随机噪声。随机噪声的存在，直接影响了地震资料时间剖面同相轴的连续性，极大增加了地震资料解释的难度。为了满足目前地震资源勘探的要求，压制随机噪声，提高地震资料信噪比是地震数据处理的首要任务。

异常噪声属于随机噪声。陆上和海上地震资料中通常包含异常噪声。气流、电力线和其他人文噪声、记录和奇偶性误差、未修正的极性反转、孤立噪声或耦合不良的检波器、风、雨和其他影响都可能导致空间方向上的异常噪声。异常噪声比一般随机噪声幅值大，通常为非高斯分布，从而对高斯分布假设下的经典去噪方法产生了无限制的失真影响。

现有技术：

中值滤波方法。该方法基于噪声特性设定一个窗口长度控制阈值,再根据地震数据与阈值的关系选取合适的窗口进行中值滤波,以取得有效信号与噪声之间的最佳平衡效果。

现有技术的缺点：

1、中值滤波方法的滤波窗口的选取不能随噪声的不同而变化，会损失有效信号或噪声不能够被压制。

2、迭代初值难以选取，不同初值的选取会极大影响异常噪声压制效果。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法、存储介质及设备，将地震资料中的异常噪声建模为混合高斯分布以自适应地对异常点赋予较小的权重，从而降低异常噪声的影响，混合高斯分布中的各个参数以及权重矩阵可由EM算法得到；基于得到的权重矩阵，利用线性化Bregman算法求得变换域下的鲁棒稀疏表示，最终得到压制异常噪声以后的地震资料。

本发明采用以下技术方案：

一种压制地震资料异常噪声的方法，包括以下步骤：

S1、读取含有高振幅异常噪声的原始地震资料，根据地震资料中噪声的复杂程度选取混合高斯分布的高斯成分个数和迭代次数；

S2、针对原始地震资料中有效信号的特点选取固定稀疏变换基；

S3、利用线性化Bregman方法得到地震资料中有效信号的稀疏表示，根据稀疏表示x^*与固定稀疏变换基的乘积和原始地震资料得到关于有效信号的估计和异常噪声估计，将异常噪声建模为混合高斯分布，利用EM算法更新每个采样点上异常噪声的权重矩阵w；

S4、重复步骤S3直到原始地震资料中异常噪声压制完成，。

具体的，步骤S1中，将原始地震资料中二维地震剖面y表示为有效信号s与异常噪声e的线性叠加；利用混合高斯分布拟合异常噪声e，根据异常噪声e的复杂程度确定迭代次数。

具体的，步骤S2中，基于wrapping方法的离散Curvelet，其正变换步骤如下：

S201、对二维信号做二维FFT，得到信号的二维傅里叶矩阵如下：

其中，(n₁,n₂)表示数据中的采样点；

S202、对每个尺度和角度方向的频率矩阵做加窗处理如下：

其中，

为窗函数，j表示第j个尺度，l表示第l个角度；

S203、对步骤S202中得到的每个矩阵做环绕处理，标准遵循Candes准则得到如下：

其中，W表示对加窗后的频率矩阵做环绕处理；

S204、对步骤S203中得到的

做二维FFT处理获得Curvelet变换系数矩阵C^D(j,l,k)，其中，k表示位移。

进一步的，Curvelet反变换的实现步骤如下：

S205、将每个尺度和角度方向的C^D(j,l,k)序列做二维FFT变换得到如下傅里叶序列如下：

S206、将步骤S205中得到的每个尺度和方向对的傅里叶矩阵与对应窗的环绕矩阵

相乘得到数据如下：

S207、将步骤S206中的数据进行解环绕处理，将得到的数据相加得到原始信号的傅里叶矩阵如下：

S208、对步骤S207中得到的数据进行IFFT处理得到重构的原始信号。

具体的，步骤S3中，迭代过程具体为：

混合高斯分布中高斯成分个数K，总迭代次数I，曲波变换字典D；每步迭代j增加1，并计算x^j；当

小于预设的值或迭代次数j＝I时，迭代终止；输出稀疏表示系数x^*＝x^j和有效信号

进一步的，稀疏系数的目标函数如下：

其中，D∈R^m×n为固定字典变换基，x∈R^n×1为对应于字典D的稀疏系数，⊙为Hadamard乘积，λ为，y为原始地震资料中二维地震剖面。

进一步的，异常噪声e建模为混合高斯分布如下：

其中，

表示第k个以0为均值、

为方差的高斯分布，共K个高斯分布；π_k为噪声混合比例。

进一步的，混合高斯分布的权重矩阵w更新公式为：

其中，i＝1,2,...,m，k＝1,2,...K，m为含噪数据y中的元素个数，γ_ik为自定义变量。

一种存储一个或多个程序的计算机可读存储介质，所述一个或多个程序包括指令，所述指令当由计算设备执行时，使得所述计算设备执行所述的方法中的任一方法。

一种计算设备，包括：

一个或多个处理器、存储器及一个或多个程序，其中一个或多个程序存储在所述存储器中并被配置为所述一个或多个处理器执行，所述一个或多个程序包括用于执行所述的方法中的任一方法的指令。

与现有技术相比，本发明至少具有以下有益效果：

本发明一种压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法，将异常噪声建模为混合高斯分布，可以利用EM算法自适应地得到每个采样点中异常噪声的权重，从而对高幅值异常噪声赋予较小的权重，以降低高幅值噪声对整个目标函数的影响；利用地震资料中有效信号在某些变换域下具有稀疏表示这一特性，对有效信号具有高保真性，能够较大程度地保护有效信号；使用固定字典，计算效率较高。

进一步的，读取含有高振幅异常噪声的原始地震资料，根据地震资料中噪声的复杂程度选取混合高斯分布的高斯成分个数和迭代次数；选取混合高斯分布的高斯成分个数的目的是可以更加准确地估计高振幅异常噪声的实际分布，选取迭代次数的目的是可以在噪声压制干净的同时尽可能地节省处理时间。

进一步的，针对原始地震资料中有效信号的特点选取固定稀疏变换基；这样做的目的是使得地震资料有效信号有更加稀疏的表示，同时固定变换基可以使得处理速度更快。

进一步的，利用线性化Bregman方法得到地震资料中有效信号的稀疏表示，根据稀疏表示x^*与固定稀疏变换基的乘积和原始地震资料得到关于有效信号的估计和异常噪声估计，将异常噪声建模为混合高斯分布，利用EM算法更新每个采样点上异常噪声的权重矩阵w；这样做的目的是线性化Bregman方法可以快速有效求得地震资料有效信号的稀疏表示，将异常噪声建模为混合高斯分布可以更加精确地拟合异常噪声的实际分布，EM算法可以有效求解目标函数。

进一步的，将异常噪声建模为混合高斯分布可以更加精确地拟合异常噪声的实际分布，从而更加精准地压制异常噪声。

进一步的，更新每个采样点上异常噪声的权重矩阵w，可以给含噪数据中噪声能量大的元素赋予较小的权重，给噪声能量小的元素赋予较大的元素，从而降低噪声能量大的元素对整个处理结果的影响。

综上所述，本发明利用曲波变换固定变换基对有效信号进行稀疏表示，极大程度地保留有效信号，对有效信号的损伤极小，将地震资料中的异常噪声建模为混合高斯分布以自适应地对异常点赋予较小的权重，从而降低异常噪声的影响。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为Curvelet变换原子示意图；

图2为模拟地震数据图，其中，(a)为含有三根同相轴的模拟地震数据；(b)为在图a基础上加上高斯随机噪声和表示坏道的高振幅异常噪声的模拟地震数据；

图3本发明方法得到的有效信号示意图，其中，(a)为有效信号；(b)为异常噪声；

图4为含有高振幅异常噪声的叠前原始单炮地震记录；

图5为发明方法得到的第二种有效信号示意图，其中，(a)为有效信号；(b)为异常噪声；

图6为图4中原始单炮地震记录与图5a中有效信号的振幅谱；

图7为本发明流程图。

具体实施方式

本发明提供了一种压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法，将地震资料中的异常噪声建模为混合高斯分布以自适应地对异常点赋予较小的权重，从而降低异常噪声的影响，混合高斯分布中的各个参数以及权重矩阵可由EM算法得到；基于得到的权重矩阵，利用线性化Bregman算法求得变换域下的鲁棒稀疏表示，从而对异常噪声进行压制，得到最终的有效信号。

请参阅图7，本发明一种压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法，包括以下步骤：

S1、读取含有高振幅异常噪声的原始地震资料，根据地震资料中噪声的复杂程度选取混合高斯分布的高斯成分个数和迭代次数；

将原始地震资料中二维地震剖面y表示为有效信号s与异常噪声e的线性叠加；

y＝s+e

其中，y∈R^m×1为含有异常噪声的观测信号，s∈R^m×1为有效信号，e∈R^m×1为异常噪声；

考虑利用混合高斯分布来拟合异常噪声e，由混合高斯分布的特点可知，异常噪声e越复杂，用来拟合e的混合高斯分布的高斯成分个数越多。因此，应该根据异常噪声e的复杂程度来选取混合高斯分布的高斯成分个数。

类似地，根据异常噪声e的复杂程度确定迭代次数，e越复杂，选取的迭代次数应越多。

S2、针对原始地震资料中有效信号的特点选取固定稀疏变换基，如曲波变换，二维离散余弦变换等；

原始地震资料中的有效信号在某些固定变换基张成的空间下具有稀疏表示，选择曲波变换作为原始地震资料中有效信号的稀疏表示字典。

基于wrapping方法的离散Curvelet，其正变换步骤如下：

S201、对二维信号做二维FFT，得到信号的二维傅里叶矩阵

其中，(n₁,n₂)表示数据中的采样点；

S202、对每个尺度和角度方向的频率矩阵做加窗处理

其中，

为窗函数，j表示第j个尺度，l表示第l个角度；

S203、对步骤S202中得到的每个矩阵做环绕处理，标准遵循Candes准则得到

其中，W表示对加窗后的频率矩阵做环绕处理；

S204、对步骤S203中得到的

做二维FFT处理获得Curvelet变换系数矩阵C^D(j,l,k)，其中，k表示位移；

Curvelet反变换的实现步骤如下：

S205、将每个尺度和角度方向的C^D(j,l,k)序列做二维FFT变换得到如下傅里叶序列

S206、将步骤S205中得到的每个尺度和方向对的傅里叶矩阵与对应窗的环绕矩阵

相乘得到如下形式的新数据

S207、将步骤S206中的数据进行解环绕处理，将得到的数据相加得到原始信号的傅里叶矩阵

S208、对步骤S207中得到的数据进行IFFT处理得到重构的原始信号。

请参阅图1，为曲波(Curvelet)变换原子。图1中为具有不同尺度不同取向的时空域Curvelet原子。Curvelet变换原子适宜于具有各向异性的曲线状结构信号的多尺度和多方向分析，因此选取Curvelet变换原子作为稀疏固定变换字典，可以做到对地震记录中有效信号的稀疏表示。

S3、利用线性化Bregman方法得到地震资料中有效信号的稀疏表示，根据该稀疏表示与固定稀疏变换基的乘积和原始地震资料得到关于有效信号的估计和异常噪声的估计，将异常噪声建模为混合高斯分布，利用EM算法更新每个采样点上异常噪声的权重矩阵；

提出求解稀疏系数的目标函数如下：

其中，D∈R^m×n为固定字典变换基，x∈R^n×1为对应于字典D的稀疏系数。⊙为Hadamard乘积，代表了前后两项对应元素直接相乘的运算，w表示每个采样点上异常噪声的权重矩阵。

由于将异常噪声e建模为混合高斯分布，即：

其中，

表示以0为均值、

为方差的高斯分布，共K个高斯分布；π_k为噪声混合比例；

为第k个高斯成分噪声方差。

根据EM算法，可得到混合高斯分布的每个参数和权重矩阵w的更新公式：

其中，i＝1,2,...,m，k＝1,2,...K，γ_ik为自定义，p(B|A)表示在A条件下B的概率，D_i为固定字典变换基D∈R^m×n的第i列。

基于如上参数更新公式，可得到如下迭代过程。

具体步骤为：

初始化：混合高斯分布中高斯成分个数K，总迭代次数I，曲波变换字典D；

迭代：每步迭代j增加1，并计算：

其中，i＝1,2,...,m，k＝1,2,...K，D∈R^m×n为固定字典变换基，D_i为D∈R^m×n的第i列，x∈R^n×1为对应于字典D的稀疏系数，

表示l₂范数，||x||₁表示l₁范数。

终止条件：当

小于预设的值或迭代次数j＝I时，迭代终止；

输出：稀疏表示系数x^*＝x^j和有效信号

S4、重复步骤S3直到原始地震资料中异常噪声压制完成，涌浪噪声压制完成后，反射信号能量和连续性明显加强，信噪比会显著提高，从而使得后续的反演结果更加可信，也更加便于地震资料解释人员对地下信息的解释。

由于一次迭代并不能完全压制涌浪噪声，因此需要重复步骤S3直到原始地震资料中异常噪声压制完成。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

将基于本发明的压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法应用到含有高振幅异常噪声的模拟地震数据和含有高振幅涌浪噪声的叠前原始单炮地震记录上，来达到压制异常噪声，保护有效信号的目的。应用结果表明，本发明可以在保护有效信号的同时压制大部分异常噪声，为后续地震数据的处理提供高信噪比的地震资料。

请参阅图2，为含有高斯随机噪声以及表示坏道的高振幅异常噪声的模拟地震数据，利用压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法得到的有效信号和异常噪声如图3a和3b所示。从图3A中可以看出，有效信号里几乎没有涌浪噪声的残留，从图3b的异常噪声中也看不到有效信号的存在，另外，用3a中的有效信号与图2a中的真实无噪的有效信号对比可以看出二者基本相同，均说明本发明提出的方法可以对异常噪声进行有效压制。

图4为含有高振幅异常噪声的叠前原始单炮地震记录，利用压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法得到的有效信号和异常噪声如图5a与图5b所示。从图5a中可以看出，有效信号里几乎没有异常噪声的残留，从图5b的异常噪声中也几乎看不到有效信号的存在，说明本发明提出的方法可以对高振幅涌浪噪声进行有效压制。

请参阅图6，为了进一步说明本发明方法的有效性，使用图4原始单炮地震记录与图5a中有效信号的振幅谱进行验证。有效信号的振幅谱与原始单炮地震记录的振幅谱相较可以看出，本发明方法并没有损失有效信号的频率成分，说明了本发明方法的保频性。

以上的模拟数据与叠前原始单炮地震记录算例说明了，利用本发明的压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法，可以实现对高振幅异常噪声的有效压制。

综上所述，本发明一种压制地震资料异常噪声的混合高斯鲁棒稀疏表示方法，将地震资料中的异常噪声建模为混合高斯分布，从而可以更加精确地拟合异常噪声的实际分布，对异常噪声进行更加精准地压制；同时利用地震资料有效信号在curvelet变换基下具有稀疏性这一特点，可以求得有效信号的稀疏表示，极大程度地保留有效信号，对有效信号的损伤极小。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (10)

1.一种压制地震资料异常噪声的方法，其特征在于，包括以下步骤：

S1、读取含有高振幅异常噪声的原始地震资料，根据地震资料中噪声的复杂程度选取混合高斯分布的高斯成分个数和迭代次数；

S2、针对原始地震资料中有效信号的特点选取固定稀疏变换基；

S3、利用线性化Bregman方法得到地震资料中有效信号的稀疏表示，根据稀疏表示x^*与固定稀疏变换基的乘积和原始地震资料得到关于有效信号的估计和异常噪声估计，将异常噪声建模为混合高斯分布，利用EM算法更新每个采样点上异常噪声的权重矩阵w；

S4、重复步骤S3直到原始地震资料中异常噪声压制完成，。

2.根据权利要求1所述的方法，其特征在于，步骤S1中，将原始地震资料中二维地震剖面y表示为有效信号s与异常噪声e的线性叠加；利用混合高斯分布拟合异常噪声e，根据异常噪声e的复杂程度确定迭代次数。

3.根据权利要求1所述的方法，其特征在于，步骤S2中，基于wrapping方法的离散Curvelet，其正变换步骤如下：

S201、对二维信号做二维FFT，得到信号的二维傅里叶矩阵如下：

其中，(n₁,n₂)表示数据中的采样点；

S202、对每个尺度和角度方向的频率矩阵做加窗处理如下：

其中，

为窗函数，j表示第j个尺度，l表示第l个角度；

S203、对步骤S202中得到的每个矩阵做环绕处理，标准遵循Candes准则得到如下：

其中，W表示对加窗后的频率矩阵做环绕处理；

S204、对步骤S203中得到的

做二维FFT处理获得Curvelet变换系数矩阵C^D(j,l,k)，其中，k表示位移。

4.根据权利要求3所述的方法，其特征在于，Curvelet反变换的实现步骤如下：

S205、将每个尺度和角度方向的C^D(j,l,k)序列做二维FFT变换得到如下傅里叶序列如下：

S206、将步骤S205中得到的每个尺度和方向对的傅里叶矩阵与对应窗的环绕矩阵

相乘得到数据如下：

S207、将步骤S206中的数据进行解环绕处理，将得到的数据相加得到原始信号的傅里叶矩阵如下：

S208、对步骤S207中得到的数据进行IFFT处理得到重构的原始信号。

5.根据权利要求1所述的方法，其特征在于，步骤S3中，迭代过程具体为：

混合高斯分布中高斯成分个数K，总迭代次数I，曲波变换字典D；每步迭代j增加1，并计算x^j；当

小于预设的值或迭代次数j＝I时，迭代终止；输出稀疏表示系数x^*＝x^j和有效信号

6.根据权利要求5所述的方法，其特征在于，稀疏系数的目标函数如下：

其中，D∈R^m×n为固定字典变换基，x∈R^n×1为对应于字典D的稀疏系数，⊙为Hadamard乘积，λ为，y为原始地震资料中二维地震剖面。

7.根据权利要求5所述的方法，其特征在于，异常噪声e建模为混合高斯分布如下：

其中，

表示第k个以0为均值、

为方差的高斯分布，共K个高斯分布；π_k为噪声混合比例。

8.根据权利要求5所述的方法，其特征在于，混合高斯分布的权重矩阵w更新公式为：

其中，i＝1,2,...,m，k＝1,2,...K，m为含噪数据y中的元素个数，γ_ik为自定义变量。

9.一种存储一个或多个程序的计算机可读存储介质，其特征在于，所述一个或多个程序包括指令，所述指令当由计算设备执行时，使得所述计算设备执行根据权利要求1至8所述的方法中的任一方法。

10.一种计算设备，其特征在于，包括：

一个或多个处理器、存储器及一个或多个程序，其中一个或多个程序存储在所述存储器中并被配置为所述一个或多个处理器执行，所述一个或多个程序包括用于执行根据权利要求1至8所述的方法中的任一方法的指令。

Images