基于小波和CEEMDAN的地震信号噪声压制

2021-05-28方根显

物探化探计算技术 2021年3期

殷明, 方根显

(东华理工大学地球物理与测控技术学院，南昌 330013)

0 前言

浅层地震勘探由于其精度高，在城市地下空间勘探中得到了广泛地应用，长期以来为工程建设、城市规划等提供了重要的支撑信息[1]，因此在城市地质调查中被作为常用的探测技术。由于城市特殊的地理环境的影响，地震信号不仅受各种各样的人为因素干扰，而且还有天然的环境因素干扰，因此浅层地震勘探的探测效果在复杂的环境中受干扰影响大。对于非平稳非线性的地震信号， huang等[2]提出了EMD(经验模态分解)克服了傅里叶变换的局限性，但使用EMD容易出现模态混淆[3]问题，在使用后期对其进行了改进算法-集合经验模态分解 (EEMD, ensemble empirical mode decomposition)[4]，但由此引入的白噪声并不能完全剔除。所以为了解决这一问题，CEEMDAN(自适应噪声完备集合经验模态分解)便被提出，这种方法的提出有效地克服了EEMD的不足。

针对地震信号中的噪声干扰不易压制的问题，使用小波-CEEMDAN相结合的方式，对其信号进行滤波。首先对信号进行小波分解，对其高频系数中存在的噪声信号，使用CEEMDAN进行分解，随后对分解的IMF分量求其自相关函数并判断噪声主要存在的分量并进行剔除，在对其与小波低频系数进行重构并求出信噪比和均方差，判断滤波后的效果，最后通过实际地震信号去噪对比效果,验证本文小波-CEEMDAN噪声压制方法的有效性。

1 浅层地震信号干扰及噪声压制方法

1.1 噪声干扰

在城市地震勘探中噪声干扰是不确定的，有规则干扰波与不规则干扰波。规则干扰波可以用常规的方法去压制噪声，如带通滤波、二维滤波、反褶积等[5]，但这些方法对于不规则干扰却不再适用，所以一般采用常规手段和小波变换相结合的方式进行噪声压制。

1)50 Hz工业电干扰。由于地震测线的上方存在高压输电线路时，会导致检波器感应到50 Hz的电压，出现电干扰[6-8]，这种50 Hz的电干扰，可以用单频干扰消除[8]。

2)地下建筑产生干扰次波。地下建筑物和松散层之间存在明显的物性差异，当地震波传播到这些界面时，会导致干扰波的形成(如折射，反射等)，所以通过改变激发和观察位置能明显减少其影响。

3)高、低频背景干扰。当地震波在松散介质中激发时，低频干扰(10 Hz～30 Hz)会随震源振动而形成。震源在硬介质中激发时，会导致80 Hz到200 Hz的高频干扰。减弱或消除高、低频干扰波对有效波的干扰，可以采用带通滤波截至高低频带的方法。

4)虚反射。它是由于地震波向上传播，当遇到地面后或者低速带底部时又转为向下反射传播，最后从下方反射界面再次反射到地面的现象[9]。由于该干扰波频带较宽, 可以利用反褶积来减弱或消除对有效波的干扰。

5)声波。它常常可以在地震记录中找到，它的频率较高(大于100 Hz)，特点是波速稳定(330 m/s ～340 m/s)。在浅层地震勘探中，特别是地面震源的声波干扰较为严重[10]，对此可以利用二维速度滤波方法消除。

6)面波。主要是地表传播的瑞雷面波，它在浅层地震勘探中是一种主要的干扰波。面波的主要特点是视速度小，频率低(10 Hz～30 Hz)、衰减慢、能量强，对此可通过提高低截止频率来衰减干扰。

7)微动。是由人为因素所导致的一系列不规则的振动或者大自然当中的一些噪声。这类型的噪声频带很宽(1 Hz～200 Hz)，并且难以区分，对此可以通过小波变换与CEEMDAN进行噪声压制。

1.2 小波变换

对于时频信号分析，小波变换是一种较好的方法，它具有信噪分离、提取弱信号特征、多分辨率的特性[11]，可以避免因傅氏变换出现的频泄现象，并且在时域和频域分析都有较好的效果，所以对于地震波这种非稳定信号的分析处理中可以使用小波变换。

对于任何的地震信号f(t),小波变换都可以将其表示为小波函数和地震信号f(t)内积:

(1)

在信号处理中常采用Mallat算法实现对信号的小波变换和信号重构，即将a尺度空间的剩余系数dj,k，经过滤波器系数h0(n)、h1(n)加权求和就得到j+1长度空间的剩余系数dj+1,k和小波系数Cj+1,k。可以表示为：

dj+1,k=∑mh0(m-2k)dj,m

(2)

Cj+1,k=∑mh1(m-2k)Cj,m

(3)

重构公式为：

Cj,n=∑mCj+1,kh1(m-2k)+

∑mdj+1,kh0(m-2k)

(4)

式中:h0、h1为滤波器系数;dj+1,k、Cj+1,k为剩余系数和小波系数。

从重构信号可见，选取小波系数和滤波手段，是压制噪声的关键。选用什么样的小波基函数和处理等手段是决定小波降噪效果的关键。

1.3 CEEMDAN原理

EMD(经验模态分解)是一种分析非平稳非线性的自适应方法[12]，但是容易导致模态混淆问题，所以Huang[13]提出了EEMD(集合经验模态分解)，通过添加白噪声来让EMD导致的模态混淆问题减少，然而在EEMD重构信号的同时会存在掺杂的噪声。在EEMD的基础上Torres[14]提出了CEEMDAN(自适应噪声完备集合经验模态分解)，CEEMDAN其实是在EEMD的基础上进行补充的一种信号处理方法。其主要在EEMD的算法上加入了有限次的自适应白噪声，在分解的每一个阶段都添加自适应白噪声，再通过计算获取每个IMF在比较少的次数下实现重构地震信号并且达到几乎没有误差，保证了分解的完整性同时也减少了计算的时间。

CEEMDAN分解如下：

1)通过对原始信号S(t)在添加不同的白噪声ωi(t)，原信号S(t)变为x(t)+ε0ωi(t)，ε为噪声系数。用EMD获取原始信号S(t)分解的第一个IMF分量公式表示为：

E1[x(t)+ε0ωi(t) ]

(5)

2)第一个计算残差信号r1为：

r1(t)=s(t)-IMF1[t]

(6)

将r1(t)和经过EMD分解的E1[ωi(t)]噪声分量相加，然后在对其结果进行EMD分解，即可以得到第二个IMF分量即:

(7)

式中：r1(t)为第一个残差分量；E1[ωi(t)]为信号通过EMD分解得到的第一个模态函数。

3)当IMFk[t],K=2、3、…、k时，计算K的残差便重复以上过程。

4)直到残差不能在进行分解便得到K个模态函数分量，其中K为完全自适应无需人为设置，最后的结果R(t)为式(8)。

(8)

5)重构信号x(t):

(9)

从上述步骤可以看出，CEEMDAN的分解能准确地重构目标信号，克服了EMD分解导致的模态混淆，有效地获得了分离谱，计算效率也较大地提高。

2 小波变换-CEEMDAN滤波

2.1 高频IMF分量的定位和选取

当使用传统的EMD类方法用于去噪时，根据以往的知识，即随机噪声信号通常分布在高频分量中，则直接丢弃最前的1个到2个高频IMF分量。再用剩下的IMF分量重构获得去噪后的信号。但是依靠经验选择几个高频IMF分量会都受主观因素的影响，高频分量舍弃过多或不足都会导致去噪效果的不佳。

故先通过CEEMDAN自适应分解得到最优的分解层数，然后在计算高频细节系数与分解的各个IMF分量之间的相关性系数的方法来确定高频IMF分量的准确位置，相关性系数大的IMF分量为地震信号中包含的噪声干扰将其剔除，其他的分量则为有效信号将其保留，从而实现噪声的准确定位。

2.2 地震信号噪声压制处理流程

小波变换降噪方式种类太多，大都采用软阀值、硬阀子。对此阀子的选取对于信号的降噪表现有很大影响，如果选择不好，要么消噪太多，要么就是不够，所以选取一个合适的阀子便显得格外重要。为了摆脱这种问题笔者将小波分解和CEEMDAN相结合对地震信号进行降噪。整体降噪实现步骤如下(图1)：

1)对地震原始信号S(t)进行小波分解，将分解后的每一层高频细节系数进行CEEMDAN分解。

2)使用CEEMDAN对高频细节系数分解出的IMF分量，从频率高低进行排列。并计算出各个IMF分量与经过小波分解后的高频系数的相关性系数。

3)对求出的相关性系数,分析系数之间的线性关系。

4)通过相关性分析，选择斜率较大的数据进行剔除，将保留下来的IMF分量进行重构得到处理后的小波高频分量。

5)对剩下的小波高频分量重复步骤2)、步骤3)、步骤4)。

6)重构原始信号。将所有处理后的小波高频细节系数与低频近似系数进行重构，得到去噪信号。

图1 小波—CEEMDAN降噪模型Fig.1 Wavelet-CEEMDAN noise reduction model

图2 地震原始信号Fig.2 The original seismic signal

图3 小波分解的高、低频部分Fig.3 High and low frequency part of wavelet decomposition(a)第一层高频细节系数;(b)第二层高频细节系数;(c)第三层高频细节系数;(d)第三层低频近似系数

3 实际应用

为了体现本方法的可行性，笔者选用了地震实测数据，数据来源于2018年南昌市多要素地质调查项目。图2是原始地震单道信号，可以看出在单道信号中夹杂者较多噪声干扰，不能明显地分辨有效反射波的准确位置。选择bior2.4小波基[15]对其进行小波分解。图3是经过小波分解后得到的地震低频近似信号和各层的高频细节信号。经过小波分解后的高频细节部分大都存在噪声，所以在对前三层高频信号进行CEEMDAN分解，图4是第一层高频细节信号经过CEEMDAN后自适应分解得到的12个IMF分量。通过对分解得到的12个IMF分量求出各阶自相关函数和相关性系数，准确定位高频IMF分量。

图5为各阶IMF自相关函数。通过图5的观察发现IMF1的自相关函数图像有用信号完全被噪声淹没。大量的噪声存在于IMF1分量中，而IMF2-MIF4分量中主要为有效信号，干扰噪声占比很小。表1为各阶IMF与高频细节系数的相关性系数。从表1可以看出IMF1的相关性系数高达0.986 4，此为极高相关分量。而之后分量的相关性系数大都比IMF1低2个数量级为极低相关分量。反应出大量的噪声信号存在于IMF1中。图6为相关性系数之间的线性关系。通过图6可以更加直观地反应出各相关性系数之间的关系，IMF1为要剔除的高频分量。在对剩下的IMF余量进行相加在与小波低频信号进行重构，得到处理后信号见图7。从图7上看，小波-CEEMDAN滤波取得了很好效果，对原始信号的干扰噪声有了有效的压制，突出了有效反射波，原始信号的主要变化特征也得到了很好的保留。

图4 CEEMDAN分解结果Fig.4 CEEMDAN decomposition results

图5 各阶IMF自相关函数Fig.5 Autocorrelation function of each IMF

表1 各阶IMF与高频细节系数的相关性系数

图6 相关性系数之间的线性关系Fig.6 Linear relationship between correlation coefficients

图7 小波-CEEMDAN滤波后信号Fig.7 Wavelet-CEEMDAN filtered signal

为了体现小波-CEEMDAN滤波的效果，笔者对信号进行信噪比(SNR)和均方差(RMSE)的量化比较评价：信噪比数值越大则信号滤波效果越好，均方差数值越小滤波效果越好。

SNR=10*lg{∑x(t)2/∑[x(t)-x(t)′]2}

(10)

(11)

式中:x(t)、x(t)′分别为原始信号和去噪后信号。

由式(10)、式(11)计算的信噪比和均方差结果见表2。

表2 信噪比和均方差结果

图8 使用小波-CEEMDAN滤波地震剖面图Fig.8 Seismic profile using wavelet- CEEMDAN filtering

图9 未使用小波-CEEMDAN滤波地震剖面图Fig.9 Seismic profile without wavelet- CEEMDAN filtering

从表2可以看出，使用小波-CEEMDAN相结合的滤波方法，比单独使用CEEMDAN进行滤波效果更好，大大提高了地震信号的信噪比，也降低了信号的均方差数值。因此可以说明地震信号通过小波-CEEMDAN滤波，在有效地去除信号的噪声信号的同时，也极大地保留了原始信号的主要特征。使用小波-CEEMDAN压制噪声效果要优于单独使用CEEMDAN方法。所以我们对整个剖面的地震信号进行小波-CEEMDAN滤波，滤波前后结果如图8、图9所示，可以看出，地震信号使用波-CEEMDAN滤波后，T0波组连贯性更好，更加平滑，横向分辨率也大大增加。T1波组反射波信号得到加强，反射界面更加突出。而未使用波-CEEMDAN滤波的地震信号T0波组地震信号横向分辨率较弱，连续性较弱。T1波地震反射信号较弱，连续性不好。从图8和图9的对比中可以看出，地震信号经过小波-CEEMDAN滤波后信号能量得到加强，消除了同相轴周边的噪声干扰，信噪比也大大提高，横向分辨率更好。

4 结论

1)对小波高频细节系数使用CEEMDAN分解，再通过IMF分量算取各自信号的自相关函数，能很好地对高频系数进行噪声去除，使用这种方法重构的小波信号能在提高数据拟合度的同时对信号的噪声进行压制。

2)基于小波-CEEMDAN处理的非平稳非线性地震信号比单独使用CEEMDAN处理的地震信号效果更好，使用小波-CEEMDAN滤波相结合的方式既可以避免小波去噪多种方式的选择和去噪阈值函数的选择，还可以最大限度的保留原始信号的特征特点，避免信号的失真，使处理后的信号更加真实，为地震信号噪声压制提供一种方法。

3)对于在强干扰地区的实测城市浅层地震记录经过小波-CEEMDAN滤波处理后可以看出，其对噪声压制取得了比较满意的效果，并且准确地反映了有效反射波的位置，数据的信噪比也有明显地提升，横向分辨率提高为以后的解释研究工作也提供了便利。