基于小波和CEEMDAN的地震信號(hào)噪聲壓制

殷 明, 方根顯

(東華理工大學(xué) 地球物理與測控技術(shù)學(xué)院，南昌 330013)

0 前言

淺層地震勘探由于其精度高，在城市地下空間勘探中得到了廣泛地應(yīng)用，長期以來為工程建設(shè)、城市規(guī)劃等提供了重要的支撐信息[1]，因此在城市地質(zhì)調(diào)查中被作為常用的探測技術(shù)。由于城市特殊的地理環(huán)境的影響，地震信號(hào)不僅受各種各樣的人為因素干擾，而且還有天然的環(huán)境因素干擾，因此淺層地震勘探的探測效果在復(fù)雜的環(huán)境中受干擾影響大。對(duì)于非平穩(wěn)非線性的地震信號(hào)， huang等[2]提出了EMD(經(jīng)驗(yàn)?zāi)B(tài)分解)克服了傅里葉變換的局限性，但使用EMD容易出現(xiàn)模態(tài)混淆[3]問題，在使用后期對(duì)其進(jìn)行了改進(jìn)算法-集合經(jīng)驗(yàn)?zāi)B(tài)分解 (EEMD, ensemble empirical mode decomposition)[4]，但由此引入的白噪聲并不能完全剔除。所以為了解決這一問題，CEEMDAN(自適應(yīng)噪聲完備集合經(jīng)驗(yàn)?zāi)B(tài)分解)便被提出，這種方法的提出有效地克服了EEMD的不足。

針對(duì)地震信號(hào)中的噪聲干擾不易壓制的問題，使用小波-CEEMDAN相結(jié)合的方式，對(duì)其信號(hào)進(jìn)行濾波。首先對(duì)信號(hào)進(jìn)行小波分解，對(duì)其高頻系數(shù)中存在的噪聲信號(hào)，使用CEEMDAN進(jìn)行分解，隨后對(duì)分解的IMF分量求其自相關(guān)函數(shù)并判斷噪聲主要存在的分量并進(jìn)行剔除，在對(duì)其與小波低頻系數(shù)進(jìn)行重構(gòu)并求出信噪比和均方差，判斷濾波后的效果，最后通過實(shí)際地震信號(hào)去噪對(duì)比效果,驗(yàn)證本文小波-CEEMDAN噪聲壓制方法的有效性。

1 淺層地震信號(hào)干擾及噪聲壓制方法

1.1 噪聲干擾

在城市地震勘探中噪聲干擾是不確定的，有規(guī)則干擾波與不規(guī)則干擾波。規(guī)則干擾波可以用常規(guī)的方法去壓制噪聲，如帶通濾波、二維濾波、反褶積等[5]，但這些方法對(duì)于不規(guī)則干擾卻不再適用，所以一般采用常規(guī)手段和小波變換相結(jié)合的方式進(jìn)行噪聲壓制。

1)50 Hz工業(yè)電干擾。由于地震測線的上方存在高壓輸電線路時(shí)，會(huì)導(dǎo)致檢波器感應(yīng)到50 Hz的電壓，出現(xiàn)電干擾[6-8]，這種50 Hz的電干擾，可以用單頻干擾消除[8]。

2)地下建筑產(chǎn)生干擾次波。地下建筑物和松散層之間存在明顯的物性差異，當(dāng)?shù)卣鸩▊鞑サ竭@些界面時(shí)，會(huì)導(dǎo)致干擾波的形成(如折射，反射等)，所以通過改變激發(fā)和觀察位置能明顯減少其影響。

3)高、低頻背景干擾。當(dāng)?shù)卣鸩ㄔ谒缮⒔橘|(zhì)中激發(fā)時(shí)，低頻干擾(10 Hz～30 Hz)會(huì)隨震源振動(dòng)而形成。震源在硬介質(zhì)中激發(fā)時(shí)，會(huì)導(dǎo)致80 Hz到200 Hz的高頻干擾。減弱或消除高、低頻干擾波對(duì)有效波的干擾，可以采用帶通濾波截至高低頻帶的方法。

4)虛反射。它是由于地震波向上傳播，當(dāng)遇到地面后或者低速帶底部時(shí)又轉(zhuǎn)為向下反射傳播，最后從下方反射界面再次反射到地面的現(xiàn)象[9]。由于該干擾波頻帶較寬, 可以利用反褶積來減弱或消除對(duì)有效波的干擾。

5)聲波。它常常可以在地震記錄中找到，它的頻率較高(大于100 Hz)，特點(diǎn)是波速穩(wěn)定(330 m/s ～340 m/s)。在淺層地震勘探中，特別是地面震源的聲波干擾較為嚴(yán)重[10]，對(duì)此可以利用二維速度濾波方法消除。

6)面波。主要是地表傳播的瑞雷面波，它在淺層地震勘探中是一種主要的干擾波。面波的主要特點(diǎn)是視速度小，頻率低(10 Hz～30 Hz)、衰減慢、能量強(qiáng)，對(duì)此可通過提高低截止頻率來衰減干擾。

7)微動(dòng)。是由人為因素所導(dǎo)致的一系列不規(guī)則的振動(dòng)或者大自然當(dāng)中的一些噪聲。這類型的噪聲頻帶很寬(1 Hz～200 Hz)，并且難以區(qū)分，對(duì)此可以通過小波變換與CEEMDAN進(jìn)行噪聲壓制。

1.2 小波變換

對(duì)于時(shí)頻信號(hào)分析，小波變換是一種較好的方法，它具有信噪分離、提取弱信號(hào)特征、多分辨率的特性[11]，可以避免因傅氏變換出現(xiàn)的頻泄現(xiàn)象，并且在時(shí)域和頻域分析都有較好的效果，所以對(duì)于地震波這種非穩(wěn)定信號(hào)的分析處理中可以使用小波變換。

對(duì)于任何的地震信號(hào)f(t),小波變換都可以將其表示為小波函數(shù)和地震信號(hào)f(t)內(nèi)積:

(1)

在信號(hào)處理中常采用Mallat算法實(shí)現(xiàn)對(duì)信號(hào)的小波變換和信號(hào)重構(gòu)，即將a尺度空間的剩余系數(shù)dj,k，經(jīng)過濾波器系數(shù)h0(n)、h1(n)加權(quán)求和就得到j(luò)+1長度空間的剩余系數(shù)dj+1,k和小波系數(shù)Cj+1,k。可以表示為：

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

(2)

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

(3)

重構(gòu)公式為：

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

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

(4)

式中:h0、h1為濾波器系數(shù);dj+1,k、Cj+1,k為剩余系數(shù)和小波系數(shù)。

從重構(gòu)信號(hào)可見，選取小波系數(shù)和濾波手段，是壓制噪聲的關(guān)鍵。選用什么樣的小波基函數(shù)和處理等手段是決定小波降噪效果的關(guān)鍵。

1.3 CEEMDAN原理

EMD(經(jīng)驗(yàn)?zāi)B(tài)分解)是一種分析非平穩(wěn)非線性的自適應(yīng)方法[12]，但是容易導(dǎo)致模態(tài)混淆問題，所以Huang[13]提出了EEMD(集合經(jīng)驗(yàn)?zāi)B(tài)分解)，通過添加白噪聲來讓EMD導(dǎo)致的模態(tài)混淆問題減少，然而在EEMD重構(gòu)信號(hào)的同時(shí)會(huì)存在摻雜的噪聲。在EEMD的基礎(chǔ)上Torres[14]提出了CEEMDAN(自適應(yīng)噪聲完備集合經(jīng)驗(yàn)?zāi)B(tài)分解)，CEEMDAN其實(shí)是在EEMD的基礎(chǔ)上進(jìn)行補(bǔ)充的一種信號(hào)處理方法。其主要在EEMD的算法上加入了有限次的自適應(yīng)白噪聲，在分解的每一個(gè)階段都添加自適應(yīng)白噪聲，再通過計(jì)算獲取每個(gè)IMF在比較少的次數(shù)下實(shí)現(xiàn)重構(gòu)地震信號(hào)并且達(dá)到幾乎沒有誤差，保證了分解的完整性同時(shí)也減少了計(jì)算的時(shí)間。

CEEMDAN分解如下：

1)通過對(duì)原始信號(hào)S(t)在添加不同的白噪聲ωi(t)，原信號(hào)S(t)變?yōu)閤(t)+ε0ωi(t)，ε為噪聲系數(shù)。用EMD獲取原始信號(hào)S(t)分解的第一個(gè)IMF分量公式表示為：

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

(5)

2)第一個(gè)計(jì)算殘差信號(hào)r1為：

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

(6)

將r1(t)和經(jīng)過EMD分解的E1[ωi(t)]噪聲分量相加，然后在對(duì)其結(jié)果進(jìn)行EMD分解，即可以得到第二個(gè)IMF分量即:

(7)

式中：r1(t)為第一個(gè)殘差分量；E1[ωi(t)]為信號(hào)通過EMD分解得到的第一個(gè)模態(tài)函數(shù)。

3)當(dāng)IMFk[t],K=2、3、…、k時(shí)，計(jì)算K的殘差便重復(fù)以上過程。

4)直到殘差不能在進(jìn)行分解便得到K個(gè)模態(tài)函數(shù)分量，其中K為完全自適應(yīng)無需人為設(shè)置，最后的結(jié)果R(t)為式(8)。

(8)

5)重構(gòu)信號(hào)x(t):

(9)

從上述步驟可以看出，CEEMDAN的分解能準(zhǔn)確地重構(gòu)目標(biāo)信號(hào)，克服了EMD分解導(dǎo)致的模態(tài)混淆，有效地獲得了分離譜，計(jì)算效率也較大地提高。

2 小波變換-CEEMDAN濾波

2.1 高頻IMF分量的定位和選取

當(dāng)使用傳統(tǒng)的EMD類方法用于去噪時(shí)，根據(jù)以往的知識(shí)，即隨機(jī)噪聲信號(hào)通常分布在高頻分量中，則直接丟棄最前的1個(gè)到2個(gè)高頻IMF分量。 再用剩下的IMF分量重構(gòu)獲得去噪后的信號(hào)。但是依靠經(jīng)驗(yàn)選擇幾個(gè)高頻IMF分量會(huì)都受主觀因素的影響，高頻分量舍棄過多或不足都會(huì)導(dǎo)致去噪效果的不佳。

故先通過CEEMDAN自適應(yīng)分解得到最優(yōu)的分解層數(shù)，然后在計(jì)算高頻細(xì)節(jié)系數(shù)與分解的各個(gè)IMF分量之間的相關(guān)性系數(shù)的方法來確定高頻IMF分量的準(zhǔn)確位置，相關(guān)性系數(shù)大的IMF分量為地震信號(hào)中包含的噪聲干擾將其剔除，其他的分量則為有效信號(hào)將其保留，從而實(shí)現(xiàn)噪聲的準(zhǔn)確定位。

相關(guān)系數(shù)計(jì)算方法：

設(shè)Si(n)為第i個(gè)IMF分量Ui(n)與X(n)的相關(guān)系數(shù)，則Si(n)可以表達(dá)為：

(10)

式中:Ui(n)、X(n)為IMF分量Ui(t)和高頻細(xì)節(jié)信號(hào)X(t)的離散化表示，n=1、2、3、…、m。

在大多數(shù)情況下求得的相關(guān)系數(shù)的取值范圍應(yīng)在0～1之間，但通過實(shí)驗(yàn)發(fā)現(xiàn)，個(gè)別的IMF分量在與高頻細(xì)節(jié)信號(hào)求得的相關(guān)系數(shù)會(huì)出現(xiàn)小于“0”的情況，由于欲通過各階IMF的相關(guān)系數(shù)來反應(yīng)噪聲信號(hào)在信號(hào)中的變化趨勢，所以當(dāng)相關(guān)系數(shù)為負(fù)值時(shí)，對(duì)其取絕對(duì)值。

在剔除IMF分量的選取上，第1階IMF分量到第i階IMF分量Ui(n)和小波高頻細(xì)節(jié)系數(shù)X(n)之間的相關(guān)系數(shù)Si(n),總體來說是呈現(xiàn)單調(diào)遞減的函數(shù)關(guān)系，且它們各點(diǎn)之間的斜率穩(wěn)定且大體不會(huì)發(fā)生較大的變化。直到第i+1階IMF分量時(shí)斜率出現(xiàn)折點(diǎn)且驟然變小之后趨于平滑。這說明在第1階到i階，IMF分量中主要為噪聲信號(hào)。從第i+1階開始有效信號(hào)為主要信號(hào)。此時(shí)剔除前i個(gè)IMF分量，從第i+1階開始重構(gòu)信號(hào)。

2.2 地震信號(hào)噪聲壓制處理流程

小波變換降噪方式種類太多，大都采用軟閥值、硬閥子。對(duì)此閥子的選取對(duì)于信號(hào)的降噪表現(xiàn)有很大影響，如果選擇不好，要么消噪太多，要么就是不夠，所以選取一個(gè)合適的閥子便顯得格外重要。為了擺脫這種問題筆者將小波分解和CEEMDAN相結(jié)合對(duì)地震信號(hào)進(jìn)行降噪。整體降噪實(shí)現(xiàn)步驟如下(圖1)：

1)對(duì)地震原始信號(hào)S(t)進(jìn)行小波分解，將分解后的每一層高頻細(xì)節(jié)系數(shù)進(jìn)行CEEMDAN分解。

2)使用CEEMDAN對(duì)高頻細(xì)節(jié)系數(shù)分解出的IMF分量，從頻率高低進(jìn)行排列。并計(jì)算出各個(gè)IMF分量與經(jīng)過小波分解后的高頻系數(shù)的相關(guān)性系數(shù)。

3)對(duì)求出的相關(guān)性系數(shù),分析系數(shù)之間的線性關(guān)系。

4)通過相關(guān)性分析，選擇斜率較大的數(shù)據(jù)進(jìn)行剔除，將保留下來的IMF分量進(jìn)行重構(gòu)得到處理后的小波高頻分量。

5)對(duì)剩下的小波高頻分量重復(fù)步驟2)、步驟3)、步驟4)。

6)重構(gòu)原始信號(hào)。將所有處理后的小波高頻細(xì)節(jié)系數(shù)與低頻近似系數(shù)進(jìn)行重構(gòu)，得到去噪信號(hào)。

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

圖2 地震原始信號(hào)Fig.2 The original seismic signal

圖3 小波分解的高、低頻部分Fig.3 High and low frequency part of wavelet decomposition(a)第一層高頻細(xì)節(jié)系數(shù);(b)第二層高頻細(xì)節(jié)系數(shù);(c)第三層高頻細(xì)節(jié)系數(shù);(d)第三層低頻近似系數(shù)

3 實(shí)際應(yīng)用

為了體現(xiàn)本方法的可行性，筆者選用了地震實(shí)測數(shù)據(jù)，數(shù)據(jù)來源于2018年南昌市多要素地質(zhì)調(diào)查項(xiàng)目。圖2是原始地震單道信號(hào)，可以看出在單道信號(hào)中夾雜者較多噪聲干擾，不能明顯地分辨有效反射波的準(zhǔn)確位置。選擇bior2.4小波基[15]對(duì)其進(jìn)行小波分解。圖3是經(jīng)過小波分解后得到的地震低頻近似信號(hào)和各層的高頻細(xì)節(jié)信號(hào)。經(jīng)過小波分解后的高頻細(xì)節(jié)部分大都存在噪聲，所以在對(duì)前三層高頻信號(hào)進(jìn)行CEEMDAN分解，圖4是第一層高頻細(xì)節(jié)信號(hào)經(jīng)過CEEMDAN后自適應(yīng)分解得到的12個(gè)IMF分量。通過對(duì)分解得到的12個(gè)IMF分量求出各階自相關(guān)函數(shù)和相關(guān)性系數(shù)，準(zhǔn)確定位高頻IMF分量。

圖5為各階IMF自相關(guān)函數(shù)。通過圖5的觀察發(fā)現(xiàn)IMF1的自相關(guān)函數(shù)圖像有用信號(hào)完全被噪聲淹沒。大量的噪聲存在于IMF1分量中，而IMF2-MIF4分量中主要為有效信號(hào)，干擾噪聲占比很小。表1為各階IMF與高頻細(xì)節(jié)系數(shù)的相關(guān)性系數(shù)。從表1可以看出IMF1的相關(guān)性系數(shù)高達(dá)0.986 4，此為極高相關(guān)分量。而之后分量的相關(guān)性系數(shù)大都比IMF1低2個(gè)數(shù)量級(jí)為極低相關(guān)分量。反應(yīng)出大量的噪聲信號(hào)存在于IMF1中。圖6為相關(guān)性系數(shù)之間的線性關(guān)系。通過圖6可以更加直觀地反應(yīng)出各相關(guān)性系數(shù)之間的關(guān)系，IMF1為要剔除的高頻分量。在對(duì)剩下的IMF余量進(jìn)行相加在與小波低頻信號(hào)進(jìn)行重構(gòu)，得到處理后信號(hào)見圖7。從圖7上看，小波-CEEMDAN濾波取得了很好效果，對(duì)原始信號(hào)的干擾噪聲有了有效的壓制，突出了有效反射波，原始信號(hào)的主要變化特征也得到了很好的保留。

圖4 CEEMDAN分解結(jié)果Fig.4 CEEMDAN decomposition results

圖5 各階IMF自相關(guān)函數(shù)Fig.5 Autocorrelation function of each IMF

表1 各階IMF與高頻細(xì)節(jié)系數(shù)的相關(guān)性系數(shù)

圖6 相關(guān)性系數(shù)之間的線性關(guān)系Fig.6 Linear relationship between correlation coefficients

圖7 小波-CEEMDAN濾波后信號(hào)Fig.7 Wavelet-CEEMDAN filtered signal

為了體現(xiàn)小波-CEEMDAN濾波的效果，筆者對(duì)信號(hào)進(jìn)行信噪比(SNR)和均方差(RMSE)的量化比較評(píng)價(jià)：信噪比數(shù)值越大則信號(hào)濾波效果越好，均方差數(shù)值越小濾波效果越好。

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

(10)

(11)

式中:x(t)、x(t)′分別為原始信號(hào)和去噪后信號(hào)。

由式(10)、式(11)計(jì)算的信噪比和均方差結(jié)果見表2。

表2 信噪比和均方差結(jié)果

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

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

從表2可以看出，使用小波-CEEMDAN相結(jié)合的濾波方法，比單獨(dú)使用CEEMDAN進(jìn)行濾波效果更好，大大提高了地震信號(hào)的信噪比，也降低了信號(hào)的均方差數(shù)值。因此可以說明地震信號(hào)通過小波-CEEMDAN濾波，在有效地去除信號(hào)的噪聲信號(hào)的同時(shí)，也極大地保留了原始信號(hào)的主要特征。使用小波-CEEMDAN壓制噪聲效果要優(yōu)于單獨(dú)使用CEEMDAN方法。所以我們對(duì)整個(gè)剖面的地震信號(hào)進(jìn)行小波-CEEMDAN濾波，濾波前后結(jié)果如圖8、圖9所示，可以看出，地震信號(hào)使用波-CEEMDAN濾波后，T0波組連貫性更好，更加平滑，橫向分辨率也大大增加。T1波組反射波信號(hào)得到加強(qiáng)，反射界面更加突出。而未使用波-CEEMDAN濾波的地震信號(hào)T0波組地震信號(hào)橫向分辨率較弱，連續(xù)性較弱。T1波地震反射信號(hào)較弱，連續(xù)性不好。從圖8和圖9的對(duì)比中可以看出，地震信號(hào)經(jīng)過小波-CEEMDAN濾波后信號(hào)能量得到加強(qiáng)，消除了同相軸周邊的噪聲干擾，信噪比也大大提高，橫向分辨率更好。

4 結(jié)論

1)對(duì)小波高頻細(xì)節(jié)系數(shù)使用CEEMDAN分解，再通過IMF分量算取各自信號(hào)的自相關(guān)函數(shù)，能很好地對(duì)高頻系數(shù)進(jìn)行噪聲去除，使用這種方法重構(gòu)的小波信號(hào)能在提高數(shù)據(jù)擬合度的同時(shí)對(duì)信號(hào)的噪聲進(jìn)行壓制。

2)基于小波-CEEMDAN處理的非平穩(wěn)非線性地震信號(hào)比單獨(dú)使用CEEMDAN處理的地震信號(hào)效果更好，使用小波-CEEMDAN濾波相結(jié)合的方式既可以避免小波去噪多種方式的選擇和去噪閾值函數(shù)的選擇，還可以最大限度的保留原始信號(hào)的特征特點(diǎn)，避免信號(hào)的失真，使處理后的信號(hào)更加真實(shí)，為地震信號(hào)噪聲壓制提供一種方法。

3)對(duì)于在強(qiáng)干擾地區(qū)的實(shí)測城市淺層地震記錄經(jīng)過小波-CEEMDAN濾波處理后可以看出，其對(duì)噪聲壓制取得了比較滿意的效果，并且準(zhǔn)確地反映了有效反射波的位置，數(shù)據(jù)的信噪比也有明顯地提升，橫向分辨率提高為以后的解釋研究工作也提供了便利。