一种联合CEEMD与广义S变换的地震数据去噪方法与流程

本发明涉及地球物理技术领域，尤其涉及一种联合ceemd与广义s变换的地震数据去噪方法。

背景技术

深层油气资源已成为我国油气勘探突破发现与规模增储的重要领域，随着勘探对象地面地质条件的复杂性以及勘探深度的增加，所获得地震数据的信噪比较低,尤其深层信号在强噪音背景下甚至无法分辨有效信息，这对去噪技术提出了更高的要求。近年来国内外学者就地震数据中的噪音压制进行了大量研究，并提出了一系列处理技术与方法，如小波变换域随机噪音衰减、多项式拟合去噪技术、时空双曲滤波等，这些去噪方法可以有效压制随机噪音，但他们采用确定性的基函数，难以准确表达地震信号的非平稳特征，去噪后信号会存在一定假象。此外，针对深层信号地震信号具有噪声能量强、有效信号能量弱的特点，常规去噪方法难以获得较好去噪效果，因此需要开发出了一些新的去噪方法如，张岩等依据图像块分组稀疏表示思想，利用邻近地震数据块之间多个地震记录道的相似性，提出基于多道相似组稀疏表示的地震数据随机噪音压制方法；韩文功、张军华通过对厚层和薄层模型的幅值和频谱特征进行研究，提出使用奇异值分解和曲波变换用于强噪音下地震弱信号的检测与识别；张华等提出使用二维曲波变换对三维地震数据的时间切片进行多尺度、多方向分解，在曲波域采用阈值法进行随机噪音压制；张恒磊等提出基于反射波各向异性而随机噪音各向同性特征，通过构造各向异性函数进行加权叠加的数据保真处理实现地震数据保真去噪；宋维琪针对地面微地震检测资料中噪声特点研究了地面微地震资料τ－p变换去噪方法。上述处理技术方法原理不同，适用条件不同，对随机噪音处理的效果均不理想。考虑到地震信号的非平稳性和去噪方法对非平稳信号的适应性，基于广义s变换(gst)、互补集合经验模态分解(ceemd)的地震数据去噪技术是近些年发展起来的地震资料去噪方法，其中基于广义s变换的去噪方法，能够根据实际地震信号的时频分布特点，通过灵活调节窗函数使其随频率的变化趋势呈现多种非线性变化特征，以更好适应具体信号的分析与处理。而基于ceemd的去噪方法，具有时变、分频和高保真特性，可以有效压制地震资料中的高频噪音。考虑到ceemd去噪方法在去除高频噪音的同时会压制高频有效信息，王姣等和zhangy等提出对含噪音的经验模态分量(imf)进行小波阈值去噪处理后，再与其他不含噪声的imf分量进行信号重构，数值模拟显示具有较好的去噪效果，但小波阈值去噪效果受不同小波基函数和处理参数影响较大，选取最优参数困难。gacis等对比分析了小波阈值去噪算法和eemd去噪算法在不同信噪比条件下的去噪效果，证明了基于eemd的去噪算法具有更好的去噪效果并且相对于小波阈值去噪算法eemd是一种数据自适应的方法。针对基于gst时间-频率域自适应滤波去噪会将大部分高频噪声压制同时，也会将能量较小的有效弱信号去除的问题，wangb等提出一种振幅补偿的广义s变换用于地震数据衰减补偿，通过改善gst时频谱的能量分布对能量较小的有效弱信号进行增强处理。黄捍东等提出在广义s变换实现时引入一系列函数库和快速傅里叶变换，使运算简洁高效，并通过选取合适参数组合对时频谱进行能量重新分配重构可以获得高分辨率的地震信号；邓攻等提出使用s变换与谱分解方法相结合用于深层反射地震弱信号的识别与提取，实际深部地震勘探数据试算表明由s变换谱分解方法获得的剖面可以有效提高深部弱反射信号信噪比和分辨率并且能够清晰刻画剖面上的低频细节特征；sattarih提出自适应快速稀疏的s变换谱分解方法用于地震数据高分辨率分析。

由于ceemd舍弃高频分量的去噪方法，容易造成随机噪声压制不完全或损失部分高频有效信号的问题，基于gst时间-频率域自适应去噪方法在压制高频随机噪声的同时会将能量较小的有效弱信号去除容易造成有效信号损失的不足。

技术实现要素：

本发明的目的在于解决上述现有技术存在的缺陷，尝试将ceemd与gst时频域去噪相结合，提供一种种联合ceemd与广义s变换的地震数据去噪方法。

有益效果：

本发明通过互补方式解决上述这两个去噪方法存在的问题，对ceemd去噪时需要舍弃的含噪音高频经验模态分量(imf)进行gst时间-频率域自适应滤波去噪，保留这些分量中的高频有效信息，最后与不含噪音的其他imf一起重构原信号，以保留原信号中的低频有效弱信号信息。模型和实际资料试算结果证明，本文提出方法在有效压制随机噪音的同时能够有效保留原始信号中的高频与低频有效信息，是一种有效的噪声压制方法，并且处理后剖面的信噪比获得较大改善。

附图说明

图1(a)为合成地震记录模拟分析模拟地震数据图；

图1(b)为合成地震记录模拟分析添加高斯白噪声后地震数据图；

图1(c)为合成地震记录模拟分析合成地震记录单道ceemd分解图；

图2(a)为合成地震记录基于gst时频域自适应滤波去噪结果图；

图2(b)为合成地震记录基于ceemd舍弃imf和imf2后的重构地震数据；

图2(c)为合成地震记录使用本发明方法去噪结果；

图3(a)为不同方法基于gst时频域自适应滤波去除噪声

图3(b)为不同方法ceemd方法去除的噪声

图3(c)为使用本发明方法去噪结果；

图4(a)为对应图1a中第15道合成地震记录时频谱；

图4(b)为对应图1b中第15道合成地震记录gst时频谱；

图4(c)为对应图2a中第15道合成地震记录时频谱；

图4(d)为对应图2b中第15道合成地震记录时频谱；

图4(e)为对应图2c中第15道合成地震记录时频谱；

图5(a)为实际地震记录剖面；

图5(b)为ceemd分解舍弃高频分量的去噪方法处理后的剖面；

图5(c)为gst域自适应时频滤波去噪后剖面；

图5(d)为本发明方法去噪处理后的剖面；

图6(a)为ceemd分解的imf1剖面；

图6(b)为ceemd分解的imf2剖面；

图7(a)为图4a第200道地震记录gst时频谱；

图7(b)为图4b第200道地震记录gst时频谱；

图7(c)为图4c第200道地震记录时频谱；

图7(d)为图4d第200道地震记录gst时频谱。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

所述ceemd方法的基本原理：

为了克服emd分解中存在模态混叠效应，huangne等提出一种噪声辅助信号分析的emd方法，即集合经验模态分解(ensembleempiricalmodedecomposition,eemd)。作为emd的改进型算法，eemd有效解决了emd的模态混叠问题，同时也保留了emd处理非平稳信号的优势，并且eemd分解和emd分解一样具有二进制滤波特性。eemd利用高斯白噪声频谱均匀分布的统计特性，向原始信号中加入不同的白噪声，使得信号在不同尺度上具有连续性，但这一方法需要多次加入不同辅助白噪声，然后通过平均的方式消除引入噪声的影响，最终使分解过程具有抗噪性。yehjr等对eemd方法作了进一步改进，采用加入正负成对的辅助噪声，可以很好消除重构信号中的残余噪声，并且加入噪声集合次数少，计算效率高，这种方法称为互补集合经验模态分解(completeensembleempiricalmodedecomposition,ceemd)。

对于目标信号x(t)进行ceemd分解步骤如下：

1)首先定义:算子ej,对于给定信号，通过emd方法求解第j个模态分量函数(imf)；wⁱ(n)表示单位方差的零均值高斯白噪声n(0,1)；i＝1,…,i；εk系数表示在每一个分解阶段需要选择的信噪比。

2)计算imf1(t)：使用添加的不同白噪声通过emd进行重复分解n次，计算n次分解结果的总平均值定义为目标信号x(t)的imf1(t)，如下式所示

3)计算i＝1时的一阶残差r1(t)，公式为

r1(t)＝x(t)-imf1(t)(2)

4)通过emd实现r1(t)+ε1e1(wⁱ(t))，满足第一个imf1(t)条件时，定义总体平均值为imf2(t)，公式为

5)计算i＝2,…,n时的残差ri(t)，公式为

ri(t)＝ri-1(t)-imfi(t)(4)

6)计算ri(t)+εiei(wⁱ(t))的imfk分量，计算总体平均值获得目标信号imf(k+1)，公式如下

7)重复步骤(5)(6),直到残差不能再分解为止，获得最终残差r(t)为

其中k是imf(t)的总数。因此，重构目标信号x(t)可以被表述为

ceemd能自适应的将一个复杂信号分解为一系列本征模态函数分量(imf)，imf分量满足从高频到低频系列分布，随机高频噪音一般分布于第一或前几个高频imf分量。基于ceemd的去噪方法通常是舍弃前几个高频imf分量，将剩余的imf分量重构原始信号实现滤波去噪。

所述于广义s变换的时频域自适应滤波去噪基本原理：

s变换由r.g.stockwell结合短时傅里叶变换(stft)和小波变换(wt)的特点，提出的一种非平稳信号时频分析方法。若信号为h(t)，则它的s变换定义为

其中τ和t是时间，τ为窗函数的中心点，控制窗函数在时间轴上的位置，f是频率。通常窗函数选取下式所示的高斯函数：

由于s变换与傅里叶变换关系为

因此s反变换定义为

常规基于广义s变换的时频滤波去噪，通常是将信号用s变换映射到时频域内，根据有效信号与噪声存在的时频差异进行阈值去噪，此过程需要进行人机交互确定噪声压制频带范围，具有一定人为主观经验性的缺陷，并且在去噪的同时也会将与噪声具有相同时频分布的有效信号去除。本文基于s变换域数据沿时间轴求和可以获得原始信号的傅里叶谱，在进行反傅里叶变换即可获得重构的原始信号，在s反变换获得重构时间域信号过程中采用数据自适应的时变滤波函数去除信号中噪声成分，这一时变滤波函数为数据驱动的自适应时变滤波去噪，无需人工干预。我们在robert等提出的基于广义s变换时频滤波去噪的方法基础上，采用非对称双曲线窗函数代替广义s变换中的高斯窗函数，并对时频滤波函数进行了改进，这一改进的自适应时变滤波函数如下

式中|s(τ,fi)|为信号由广义s变换获得的时频域数据,fi为单道地震记录gst数据中的每一频率成分。非对称双曲线函数表达式为

式中是一个向后参数，用于补偿非对称双曲线函数时窗反方向的衰减速度；是一个向前参数，用于增强非对称双曲线函数时窗前进方向的变化速度；λhy是一个曲率参数，用于控制非对称时窗函数向对称时窗随频率变化的速度。将gf(t-τ)代替式(8)中的高斯窗函数，得到这一窗口可调且非对称的广义s变换方法为

时变滤波函数通过(12)式求取信号经广义s变换获得时频谱中每一时频数据的比例因子，比例因子的大小随时频谱中时频数据的大小而变化。由(12)式可以看出，低频部分的时频数据比例因子较大，高频部分的时频数据比例因子较小。将信号经gst获得的时频域数据体与这一时变滤波函数在时频域相乘，低频部分由于与较大比例因子相乘可以明显突出，高频部分由于与较小比例因子相乘而被压制，因此通过这一时变滤波函数进而可以实现高频噪音的去除,之后再进行s反变换获得去噪重构后的时间域信号如下式

所述联合ceemd与广义s变换去噪方法原理：

对于含噪地震信号，当其中包含较小幅度的有效信号时，基于gst的时间-频率域自适应滤波去噪会将大部分高频噪声压制同时，也会将这些小幅度的有效信号去除；而单纯采用ceemd直接舍弃高频imf分量的去噪方法时，会在去除高频噪音的同时损失高频有效信息或造成噪声压制不完全。因此本文提出将ceemd分解与gst时频域自适应滤波相结合，通过互补方式解决上述这些问题,对ceemd去噪时需要舍弃的含噪音高频经验模态分量(imf)进行gst时间-频率域自适应滤波去噪，保留这些分量中的高频有效信息，最后与不含噪音的其他imf一起重构原信号，以保留原信号中的低频有效弱信号信息。这一方法的计算步骤如下

1)先对原始信号x(t)进行ceemd分解，获得各imf分量；

2)采用自相关法^[11]分析各imf分量，确定出前几个分量中需要进行去噪处理的imf分量；

3)对确定出的含噪imf分量进行基于gst的时频域自适应滤波去噪处理；

4)将去噪处理后的imf分量以及其他剩余imf分量进行信号重构，获得去噪后的信号x'(t)

式中:x'(t)为去噪后的信号；c'(t)为采用gst时频域自适应滤波去噪后的imf分量，c(t)为其他imf分量，r(t)为ceemd分解获得的残余分量。

本文采用如下信噪比(snr)计算表达式衡量各种方法的去噪效果，值越高，表明信号中剩余噪音越小，计算公式为

其中，snr为信噪比，单位为db，x(i)为原始数据，x'(i)为去噪后的数据。

模型数据测试分析

模型数据如图1(a)所示，该数据中地震子波是主频为35hz的零相位雷克子波。图1(b)为图1(a)中加入高斯白噪声(snr≈-1.97)后的数据，图1(c)为图1(a)单道地震数据经ceemd分解获得的一系列固有模态函数(imfs)，对这些imfs分量采用自相关法分析后发现确定imf1和imf2分量中含噪音较多。图2(a)为基于gst的时间-频率域自适应滤波去噪后的地震数据(snr≈0.80)。图2(b)为ceemd分解舍弃imf1和imf2分量后重构的地震数据(snr≈1.6)。图2(c)为采用本文提出的方法对ceemd去噪时需要舍弃的imf1和imf2分量在gst时频域进行自适应滤波去噪之后，再与其他imf分量重构原始信号后的地震数据(snr≈2.1)。

对比图2(a)、(b)和(c)可以看出ceemd直接舍弃imf1和imf2分量后，虽然可以较好压制随机噪音，提高了剖面的信噪比，但同时也损失了一部分高频有效信息，造成浅层因缺少高频成分(如图2(b)与图1(a)红色矩形框对比所示)分辨率较低，部分同相轴不连续；而基于gst的时间-频率域自适应滤波去噪，也可以较好压制随机噪音，但由于该方法的去噪特性，造成能量较弱的有效信号损失较严重(如图2(a)中200～400ms及700～800ms之间红色箭头所示)；而采用本文方法去噪处理后的剖面信噪比较高，浅层分辨率没有太大损失(如图2(c)红色矩形框所示)，较好的保持了原始信号中高频与低频有效信息，同相轴也相对更加连续，并且有效识别与提取出了在强噪音背景下的幅值较弱的有效信号(如图2(c)中200～400ms及700～800ms之间)。

图3为不同方法去除的噪声对比分析，从图3(a-c)对比可以看出，基于gst的时频域自适应滤波去噪方法，在去除大部分噪声的同时有效信号损失严重，去除的噪声数据中浅层和深层都包含部分同相轴信息如图3(a)中红色箭头所示；采用ceemd直接舍弃imf1和imf2分量的去噪方法，去除的噪声中只有在浅层损失部分有效信息(如图3(b)中红色箭头所示)；而采用本文提出的联合ceemd与广义s变换的去噪方法，去除的噪声中基本看不出有效信息损失(如图3(c)所示)。

图4分别对应图1与图2中地震数据的第15道地震记录及其gst时频谱，对比也可以看出基于gst时频域自适应滤波去噪可以去除大部分随机噪音，但也会压制能量较弱的低频有效信号，如图4(c)所示的时频谱与图4(a)原始信号的时频谱对比，可以明显看出0.1～0.2s、0.4～0.5s及0.7～0.8s之间缺失部分低频有效信号(如图4(c)中箭头所示)；含噪信号采用单纯ceemd分解去噪也可以去除大部分高频噪音成分，但去噪后的数据浅层缺失部分高频成分(如图4(d)中椭圆形框)所示；而采用本文方法处理后随机噪音能够被有效去除，并且去噪后地震记录中有效信号的高频部分与低频部分都能够有效保留如图4(e)所示。

综合对比分析可见，ceemd直接舍弃高频imf去噪方法容易导致高频有效信息损失；基于gst时频域自适应滤波去噪方法存在能量较弱低频有效信号损失的问题；本文提出的方法在有效压制随机噪音的同时，能够有效保护有效信号，是一种相对有效的去噪方法。

图5(a)为某探区的实际叠后地震剖面，从实际地震剖面中可以看出该地震数据中深层部分(如图5(a)中2～6s之间)存在较多随机噪音，地震剖面深层部分的信噪比、分辨率较低，有效信号与噪音难以区分。从图5(a)所示原始叠后地震资料中随机抽取几道数据进行ceemd分解及自相关分析，发现原始信号中随机噪音主要分布在imf1分量剖面(图6(a))和imf2分量剖面(图6(b))中所示，但也可以看出两个剖面也包含较多高频有效信息，故采用ceemd直接舍弃imf1和imf2分量去噪时会损失部分高频有效信号；而采用本文方法去噪时需要对imf1和imf2分量进行gst时频域自适应滤波去噪处理，然后与不含噪音的其他imf分量一起重构原始信号。

对实际地震资料进行不同去噪方法应用效果对比如图5(b-d)所示，其中图5(b)为ceemd直接舍弃imf1和imf2分量后的去噪剖面(snr≈4.385)；图5(c)为基于gst时频域自适应去噪处理后的剖面(snr≈6.018)；图5(d)为本文方法去噪处理后的剖面(snr≈10.696)。通过对比分析可以看出，采用单纯ceemd舍弃高频分量方法随机噪音可以被有效压制，深层有效信号的信噪比获得较大改善，尤其在图5(a)中被噪音堙没的有效弱信号的同相轴较为清晰(如图5(b)中蓝色椭圆标注部分)，但是也可以看到该方法容易造成地震数据中的部分高频有效信号缺失，造成去噪处理后剖面浅层的分辨率降低(见图5(b)中红色矩形框所示部分)，从处理后的单道时频谱中(图7(b))也可以看出浅层(0.5s～2s之间)缺失高频有效信号。图5(c)为基于gst时频域自适应去噪方法处理后的剖面中，可以看出压制噪音效果较好，但该方法在压制随机噪音的同时压制了能量较弱的低频有效信号,去噪处理后的剖面上能量较弱有效信号的同相轴连续性较差(如图5(c)中蓝色椭圆框所示)，弱能量的有效信号信息损失严重，从处理后的单道时频谱也可以看出(如图7(c)中红色箭头所示)缺失部分低频有效信号。如图4(d)采用本文方法处理后的剖面，可看出信噪比较高，高频有效信号没有太大损失浅层分辨率较高(见图5(d)中红色矩形框所示部分)，而且本文方法在处理深层地震数据时能够既能够有效压制随机噪音又能一定程度上保持原信号中的高频有效信号和能量较弱的低频有效信号成分(见图5(d)中蓝色椭圆形框所示部分)，去噪后剖面上有效信息较丰富，去噪效果相对较为理想，同相轴也相对更加清晰，并且从去噪后的单道地震记录的时频谱中(如图7(d)所示)也可以看出本文方法处理后的有效信号在时频谱中的时频分辨率相对较高，浅层和深层的有效信号的能量也相对于原始信号更为均衡。综合对比分析可见，本文提出的算法具有较好的噪声压制效果，并且具有较好的强噪声背景下有效弱信号的识别与提取能力。

去噪技术是地震数据处理中的关键技术之一，考虑到地震信号的非平稳性和去噪方法对非平稳信号的适应性，针对单纯ceemd分解舍弃高频分量的去噪方法和基于gst时频域自适应滤波存在的不足，本文将ceemd分解与gst时频域去噪相结合，提出了联合ceemd与广义s变换的去噪方法。该方法首先将信号经ceemd分解获得一系列imf分量，再通过自先关对比分析判别出含噪音较多的高频imf分量，然后对ceemd去噪时要舍弃的这些含噪音较多高频imf分量进行gst时频域自适应滤波去噪，保留这些分量中的高频有效信息，同时还避免了但当噪声能量与有效信号能量相接近时基于gst时频域滤波去噪方法会压制能量较弱有效信号的问题，达到了即能够有效压制随机噪音又能保持原信号中高频成分和低频弱信号信息不受损失的目的。模拟数据和实际资料测试分析结果表明，联合ceemd与gst时频域自适应去噪方法的去噪效果，优于单纯ceemd分解舍弃高频分量的去噪方法和基于广义s变换的时频域自适应滤波去噪方法，去噪效果证明了本文方法的有效性和优越性。

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