基于SSEC-EWT的地震资料噪声压制算法

2021-10-23李亚南

石油地球物理勘探 2021年5期

孟娟高琴李亚南②

(①防灾科技学院电子科学与控制工程学院，河北三河 065201；②河北省地震灾害监测仪器与监测技术重点实验室，河北三河 065201)

0 引言

地震勘探是油气勘探开发的重要手段，且依赖于高质量的处理成果数据。受人类活动、环境、天气、仪器等多方面因素影响，实际采集的地震数据中往往存在各种噪声干扰，这无疑会影响地震资料处理和解释的准确性。为了提高地震勘探精度，使地震资料能更真实地反映地下地质情况，压制噪声干扰在地震数据处理中就显得极为重要了。

地震资料噪声包括规则噪声和随机噪声两大类。面波是其中一种能量较强且广泛存在的规则干扰，其频率低、衰减慢，大大降低了地震资料信噪比；而随机噪声则显得杂乱无规律，频率分布宽泛，空间上无确定的视速度。因此，可根据面波和随机噪声各自的特点探寻相应的适用去噪方法。

在面波压制领域，S变换[1]、广义S变换[2-3]、Curvelet变换[4-5]、小波变换[6-7]和f-x域滤波[8]等方法广泛使用，利用面波与有效波在变换域上的特征差异实施分离。这些方法有一定效果但也有局限性，且会不同程度地损害有效信号。

在随机噪声去除方面，经验模态分解(Empirical mode decomposition，EMD)[9]、小波阈值去噪[10]、深度学习[11]等方法被引入，但这些方法同样有各自优缺点。如EMD存在模态混叠、伪模态等问题，难以兼顾增强去噪效果与降低有效信号损伤；小波阈值去噪效果依赖于小波基的选取；深度学习需大量训练数据，计算量大。这些算法一般只针对单一类型噪声，难以同时处理面波和随机噪声这两种最普遍存在的干扰。

作为一种自适应的信号处理方法，EWT(经验小波变换)[12]结合小波变换和EMD技术优势，通过分割信号的傅里叶谱，并在每个分割区间构造正交小波基，将信号分解为若干分量，从而准确提取信号中固有模态分量，适用于非线性非平稳信号。

EWT推出后被广泛应用于噪声压制及故障诊断等方面[13-15]，在地震数据处理与分析领域，也倍受青睐。Liu等[16]基于EWT对多通道地震资料进行时频分析，分解出不同频率分量，并根据计算边界提取地震信号分量，分辨率明显高于传统连续小波变换。基于油气藏导致地震波衰减增大的原理，Hui等[17]利用EWT分析地震波衰减时频特性，圏定衰减异常位置，以厘清含气层厚度变化，实现油气检测。在地震数据去噪方面，覃发兵等[18]利用EWT将地震信号分解为若干分量，并计算各分量主频，剔除阈值外分量后重构实现随机噪声压制。Chen等[19]利用EWT将地震信号分解为若干IMF，再利用聚类算法对含噪分量做阈值化处理。

以上应用表明利用EWT可有效分析地震信号，但传统EWT根据信号频谱幅值进行频带划分和信号分解，易受噪声影响，因此有必要根据地震信号特点改进EWT频谱分割方法，使之更适用于地震信号。Liu等[20]基于尺度空间表示提取信号傅里叶谱的慢变分量，再根据尺度参数获取地震信号中包含的频率分量和边界，自适应分割频谱，但尺度参数通常依据经验获取，使频谱分割的准确性受限。孟娟等[21]对地震信号连续小波变换后计算各频率点能量，得到小波谱能量曲线(WSEC)，以该曲线局部极大值点的频率及ε邻域法确定频谱分割边界进行改进EWT，实现依地震信号频率和能量的自适应分解，取得良好效果，但连续小波变换中基小波函数的选取影响变换结果精度，且通过尺度参数调整小波函数形状，获得的并不是真实的时频谱。

为同时分离地震资料中的面波和随机噪声，针对传统EWT根据频谱极值分割频带抗干扰性差的不足，本文提出基于S谱能量曲线的改进EWT地震资料噪声压制算法，实现面波和随机噪声同步压制。通过对地震信号S变换后计算各频点能量得到S谱能量曲线(SSEC)，根据曲线极大值点分布进行频带分割，将原信号自适应分解为若干按频率和能量分布的IMF，并根据SSEC确定面波IMF；为保护面波IMF中的低频有效信号，对其进行带通滤波，以实现精准的面波剔除；最后根据有效波频率阈值和各IMF主频去除随机噪声，完成面波和随机噪声的压制，有效提高地震资料信噪比。

1 EWT基本原理

EWT通过计算信号的傅里叶频谱，并基于频谱幅值的极值点进行频谱分割，在每个分割区间内构建正交的小波滤波器组，提取不同频带内的调幅、调频分量，实现信号分解。其原理详述如下：

(1)获取原信号f(t)的傅里叶频谱F(ω)，ω∈[0，π]。

(2)获取F(ω)的M个按频率升序排列的局部极大值点ω1，…，ωM，将[0，π]频带划分成N(N≤M)个子频带，设以ωn为中心的分量子频带为[Ωn-1,Ωn]。获取边界Ω的常用方法有两种：一是根据频谱的M个极大值点，将相邻两个极大值中点定义为边界，Ωn-1=(ωn-1+ωn)/2，n=2，…,N，Ω0=0,ΩN=π；或先找到信号频谱的极大值点，再寻找两相邻极大值间的极小值点，以该极小值点为边界。这样将原信号频谱划分为N个连续、不重叠的子频带。

(1)

(2)

(3)

(4)

则得到的低频分量和高频分量分别为

(5)

(6)

须注意的是，传统EWT分解信号需预先设定经验模态分解阶数N，该参数选取会影响信号分解精度，过大或过小都不能较好表征原信号。

2 基于S谱能量曲线的改进EWT

2.1 S谱能量曲线SSEC

作为短时傅里叶变换和以Morlet小波为小波基的连续小波变换的发展，S变换采用带有频率变量的高斯窗函数截取信号，从而实现信号的局部分析。S变换免去了窗函数选择，且其窗函数可随频率需求自适应地变化，具有多分辨率的特点。与连续小波变换相比，S变换分辨率依赖于频率而不是尺度，因此更直观；且每一频率的绝对相位不变，与傅里叶谱直接关联，能精确描述信号各时刻的频谱，较好适应非平稳信号频率不规律变化的特点，因此在地震数据处理领域应用较广。其原理如下：

根据短时傅里叶变换，定义高斯窗函数为

(7)

式中σ为尺度因子，可控制窗函数宽度。它与频率直接联系，从而改变时频分辨率

(8)

代入短时傅里叶变换公式，即有

(9)

可得到其S变换

(10)

对信号进行S变换得到信号时频分布的S谱后，根据式(11)对频点i定义S谱能量，实现信号时频域到频率—能量域的转换。表示为

(11)

式中：fi为S变换后频点i的频率；S(fi,τ)为该频点对应的S变换值，即S谱。

对某信号按式(11)获取所有频点能量值，从而得到依频率分布的SSEC，该曲线能清晰描绘不同频率处信号的能量分布及不同频率成分的能量强弱。

2.2 基于SSEC的改进EWT

对含面波的地震信号进行S变换并获取SSEC，SSEC能直接反映不同频率成分的能量大小。由于面波的频率低且能量强，因此通过SSEC能区分面波与有效地震波。不同于传统EWT根据傅里叶频谱的局部极大值分割频谱，本文根据SSEC的极大值点自适应地分割频谱，即以SSEC极大值点对应的频率为初始边界，通过ε邻域法搜索重新确定边界，以进行改进EWT，根据信号不同频率分量的能量大小自适应分解。算法如下：

(1)对地震信号s(t)进行傅里叶变换和S变换，得到傅里叶频谱F(ω),ω∈[0,π]和S谱S(f,τ)。

(2)根据式(11)获取各频率点能量，绘制SSEC。

改进EWT可将原地震信号分解为一系列具有紧支撑频谱的IMF分量，对各IMF解调可得到其瞬时频率、幅度、相位等属性，从而获取原信号特征。

3 基于SSEC-EWT的噪声压制方法

与地震信号相比，面波频率范围低，频带一般为4～18Hz[2]，且能量较强。本文改进EWT是根据S变换后各频点能量大小进行频谱分割，各IMF是根据原信号不同频率成分的能量大小分解的，从而可将低频、高能量面波成分从原信号中分离出来。

而随机噪声在时域内广泛存在，且一般为高频率，可根据有效地震信号频率范围设定阈值，剔除阈值范围外的噪声IMF，实现随机噪声压制。基于此思路，本文提出基于SSEC-EWT的面波和随机噪声压制方法，其流程如图1所示。

图1 基于SSEC-EWT的面波和随机噪声压制方法流程

3.1 面波压制

根据有效波与面波在频率和能量上的不同，本文利用改进EWT进行地震信号分解，提取面波IMF后进行面波压制，流程如下：

(1)计算地震信号SSEC，以能量曲线极值点及ε邻域法分割频谱，进行改进的EWT，自适应得到依频率和能量分解的各模态分量IMF1，…，IMFN。

(2)对各IMF进行希尔伯特变换，得到瞬时频率、瞬时振幅，并据下式计算各IMF主频

(12)

式中：fm为IMF主频；f为瞬时频率；s(f)为瞬时振幅谱。

(3)结合SSEC，找到频率在4～18Hz范围内，且能量强于其他邻近子频带的子频段位置，确定面波所在的IMF。

(4)根据SSEC中面波频率峰值，对该IMF进行带通滤波。滤波器截止频率为面波IMF频谱峰值的6dB带宽频率，即以频谱功率下降为功率峰值25%、幅度下降为峰值50%时的频率为截止频率，对面波IMF进行带通滤波，以最大限度地保护面波IMF中的低频有效波不被破坏。

(5)从原记录剔除滤波后面波IMF，实现面波压制。

3.2 随机噪声压制

随机噪声一般以高频率广泛存在于地震资料中，常规地震勘探中检波器接收的有效地震波频率范围一般在150Hz以内[22]。因此，根据地震波与随机噪声的频率范围差异，可实现随机噪声衰减。本文基于改进EWT的随机噪声压制具体算法如下：

(1)设定有效地震波最高频率阈值为150Hz；

(2)计算完成面波压制后的其余IMF主频，剔除阈值范围外IMF，进行随机噪声压制；

(3)对剔除面波后IMF和随机噪声IMF 信号进行重构，得到噪声压制后的地震信号。

至此，完成面波和随机噪声的同时压制。对二维资料中所有道进行以上相同操作，即可得到最终去噪记录。

4 实验与分析

4.1 SSEC-EWT地震信号分解

先对模拟地震记录进行SSEC-EWT分解，并与传统的EWT及文献[21]中的小波谱能量曲线EWT(WSEC-EWT)进行地震信号分解对比。利用Ricker子波模拟单道地震信号，设0.5s处的面波是主频为10Hz的Ricker子波，幅值为5，能量较强，0.7、0.8、0.9s处为有效波，频率分别为20、35、50Hz，幅值分别为2、2、3，其合成波及SNR=3dB的含噪合成记录如图2所示。

对该含噪记录(图2b)进行S变换，根据S谱对每个频点按式(11)计算能量得到SSEC(图3a)。可见SSEC和WSEC(图3b)都能较好体现原信号不同频率处的能量大小。但因CWT是通过尺度参数改变小波函数形状，变换结果并不是真正的时频谱，从而导致信号特征精度受影响，而S变换反映的是信号最真实的频率构成，得到的是真实时频谱。因此，相比于WSEC，本文SSEC能更真实地体现原信号低频和高频部分各频率成分的能量分布。

图2 合成单道地震波(a)及含噪地震波(b)

图3 图2b数据的SSEC(a)和WSEC(b)

基于SSEC极大值点及ε邻域法进行子频带划分结果如图4a所示，可见较好地将面波频谱与邻近频率有效波频谱分割开；WSEC-EWT分割结果(图4b)中面波与有效波频谱有少许混叠；而传统EWT划分的子带(图4c)更不能很好地反映面波与有效波的差异，故不能直接应用于面波压制。相比WSEC，本文SSEC极大值点能更直接、准确地表示有效信号和面波的频率和能量大小，因此本文频谱分割方法能更显著体现地震信号中面波与有效波的能量、频率的不同，具有更好的频谱分割性能。

图4 SSEC-EWT(a)、WSEC-EWT(b)、传统EWT(c)三种方法的频谱分割结果

本文方法得到的IMF如图5a所示，可见SSEC-EWT和WSEC-EWT(图5b)都能较好地根据信号的频率和能量进行模态分解。图5a中IMF1主要为0.5s处的面波分量(含有微量有效波)，IMF2主要为0.7、0.8、0.9s处的有效波分量，而IMF4主要为高频随机噪声，有利于实现面波与随机噪声的有效分离；图5b中IMF1主要为面波，但其中包含部分有效波；而传统EWT分解结果(图5c)中，有效波与面波混叠不清，难以实现面波提取。因此，本文SSEC能完整描述原信号各频率点的能量特征，其曲线极大值点能更直接、更真实地表示有效信号和面波的频率和能量值，从而更好地进行频谱分割和信号分解，实现面波与有效波的精确分离。

图5 SSEC-EWT(a)、WSEC-EWT(b)、传统EWT(c)三种分解方法所得的IMF

为防止面波IMF中有效波被同时压制，根据SSEC面波频率峰值，对频段做带通滤波，保留面波频率附近的有效波。图6中滤波前的面波IMF中含有少量有效波，滤波后大部分有效波被保留，二者与理想纯面波的标准差分别为0.0225、0.0209，方差分别为0.1310、0.1154，可见能最大限度减少对有效波的损伤。从经面波和随机噪声提取后重构所得的最终去噪结果(图7)可见，面波压制较彻底，随机噪声也得到有效抑制，去噪后的SNR=8.5dB。

图6 滤波前、后的面波IMF

图7 理想合成波(a)及本文方法去噪结果(b)

4.2 SSEC-EWT抗噪性测试

为测试SSEC-EWT的抗噪性，定义

(13)

首先测试在SNR较低(图8)情况下SSEC-EWT分解的准确性。选取有效波频率分别为16、35、48Hz，面波频率为6Hz，添加随机噪声后SNR=-8.8dB，有效波淹没于噪声中难以分辨。

图8 模拟合成波(a)及其加噪信号(SNR=-8.8dB)(b)

从本文SSEC-EWT分解结果(图9a)可见，IMF1是面波IMF，其中含有很少量有效波，IMF2、IMF3主要为有效波，IMF4、IMF5主要为高频噪声。显然SSEC-EWT根据原信号的频率成分和能量大小进行信号分解，在SNR较低时仍具有较好的信号分解能力，具有较好的抗噪性。面波和随机噪声压制后所得结果(图9b)中，SNR=9.3dB，面波压制较彻底，随机噪声也得到有效压制，且有效波得到较好保护。可见在较低SNR下，本文算法仍具有较好性能。

图9 SNR较低时SSEC-EWT分解所得IMF(a)及去噪结果(b)

再针对不同SNR(-15～15dB)模拟信号，用本文算法进行处理(图10)并做对比，可见经本文算法处理后，SNR明显提高(幅度高达17.9dB)，尤其是在低SNR时，本文算法也能较好压制面波和随机噪声，具有较好抗噪性能。

图10 相同SNR条件下本文去噪效果

进一步测试有效波与面波的频率有重叠时，改进EWT的分解效果。图11a中面波主频设为6Hz(图11a上)，最大峰值为5，有效波主频为7、15、20Hz，最大峰值均为3，并添加随机噪声，面波与有效波的频率几近重叠(图11a下)。按式(11)得到的SSEC如图11b所示，该曲线较清晰地体现了不同频率地震波成分的能量大小，有助于区分有效波与面波。从基于该SSEC曲线进行频谱分割后所得IMF(图11c)可见，IMF1主要为面波，IMF4为随机噪声，由于面波与有效波频率非常接近，故面波IMF中尚含有少量的有效波。

图11 面波和有效波几近重叠的模拟地震波(a)及其SSEC(b)与分解所得IMF(c)

对IMF1进行带通滤波(图12a)后，除了极少部分有效波被损伤，绝大部分有效波得到有效保护。从最终去噪结果(图12b下)也可见，面波得到充分压制，有效波得到尽可能保护，显著提高了资料的SNR。

图12 面波IMF滤波前、后对比(a)及最终去噪结果(b)

4.3 合成二维地震道集

对同时含有面波和随机噪声的二维地震道集进行噪声压制测试。根据层状速度模型和波场模拟得到图13a所示地震道集。面波速度为200m/s，直达波速度为2000m/s，合成地震记录采样间隔为1ms，共计10道、1000个采样点，面波主频为8Hz，反射波频率为40Hz。添加信噪比为15dB的随机噪声(图13b)，面波与有效反射波叠加，导致同相轴模糊。

从应用本文算法压制面波和随机噪声后结果道集(图13c)可见，面波和随机噪声同时得到有效压制，有效地震信号更光滑、更平稳，且有效反射波与面波混叠的同相轴变得清晰。观察最终被去除的噪声残差(图13d)可知，其上基本为面波和随机噪声，仅损失极少量的有效波，且面波压制彻底。

图13 含面波(a)、且含随机噪声(b)地震道集、本文方法去噪结果(c)及去噪残差(d)

对经SSEC-EWT分解后所得面波IMF记录(图14a)进行滤波处理，可见滤波后记录(图14b)上大部分有效波被保留，仅极少部分有用信息被滤除。

图14 本文方法滤波前(a)、后(b)的面波IMF记录

将WSEC-EWT与本文SSEC-EWT提取面波的结果(图15)进行对比，可见二者都能较好地提取面波，但相比WSEC-EWT中面波IMF会包含“少量”有效波，而本文方法提取面波则只含有“微量”有效波，面波提取更精准。这是因为相比连续小波变换，S变换的时频分析精确度更高，而基于频率和能量进行频谱分割后能更好分解信号。

图15 WSEC-EWT(a)与本文算法(b)提取面波对比

为对比随机噪声压制效果，对图13b地震道集应用文献[18]算法去噪(图16a)，可见虽然随机噪声被压制，但面波未被压制；而本文算法去噪后(图13c)不仅能去除面波，还能有效压制随机噪声，而且能最大限度保护有效波。从文献[18]算法去噪残差(图16a)可知，部分有效波被损害；而本文算法去噪残差(图13d)中基本不含有效波。

图16 文献[18]去噪结果(a)及随机噪声残差(b)

不同于传统EWT按傅里叶频谱极值点分割频谱，本文算法是根据频点能量进行频谱分割，故可使面波、有效波、随机噪声根据各自能量大小自适应分割，进而能同时分离面波和随机噪声。

4.4 实际地震资料处理

图17所示实际地震数据共有95道，每道采样点数为500，采样间隔为2ms，图上可见典型的扫帚状面波，降低了地震资料的SNR。

图17 实际地震资料

观察从该记录抽取第17道信号进行S变换后根据SSEC所得子频带分割(图18a)情况，可见低频有效波和面波的频谱被分割为不同子频带。在根据频谱分割得到的IMF(图18b)中，IMF2为面波，IMF1为低频有效波，IMF3～IMF4为有效波，IMF5为随机噪声。对该记录中每道都进行改进EWT，并根据各IMF的频率和能量确定面波IMF，再做带通滤波，分离出记录中的面波。