殷 明, 方根顯
(東華理工大學(xué) 地球物理與測(cè)控技術(shù)學(xué)院,南昌 330013)
淺層地震勘探由于其精度高,在城市地下空間勘探中得到了廣泛地應(yīng)用,長(zhǎng)期以來(lái)為工程建設(shè)、城市規(guī)劃等提供了重要的支撐信息[1],因此在城市地質(zhì)調(diào)查中被作為常用的探測(cè)技術(shù)。由于城市特殊的地理環(huán)境的影響,地震信號(hào)不僅受各種各樣的人為因素干擾,而且還有天然的環(huán)境因素干擾,因此淺層地震勘探的探測(cè)效果在復(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]問(wèn)題,在使用后期對(duì)其進(jìn)行了改進(jìn)算法-集合經(jīng)驗(yàn)?zāi)B(tài)分解 (EEMD, ensemble empirical mode decomposition)[4],但由此引入的白噪聲并不能完全剔除。所以為了解決這一問(wèn)題,CEEMDAN(自適應(yīng)噪聲完備集合經(jīng)驗(yàn)?zāi)B(tài)分解)便被提出,這種方法的提出有效地克服了EEMD的不足。
針對(duì)地震信號(hào)中的噪聲干擾不易壓制的問(wèn)題,使用小波-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)并求出信噪比和均方差,判斷濾波后的效果,最后通過(guò)實(shí)際地震信號(hào)去噪對(duì)比效果,驗(yàn)證本文小波-CEEMDAN噪聲壓制方法的有效性。
在城市地震勘探中噪聲干擾是不確定的,有規(guī)則干擾波與不規(guī)則干擾波。規(guī)則干擾波可以用常規(guī)的方法去壓制噪聲,如帶通濾波、二維濾波、反褶積等[5],但這些方法對(duì)于不規(guī)則干擾卻不再適用,所以一般采用常規(guī)手段和小波變換相結(jié)合的方式進(jìn)行噪聲壓制。
1)50 Hz工業(yè)電干擾。由于地震測(cè)線的上方存在高壓輸電線路時(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)致干擾波的形成(如折射,反射等),所以通過(guò)改變激發(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]。由于該干擾波頻帶較寬, 可以利用反褶積來(lái)減弱或消除對(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ì)此可通過(guò)提高低截止頻率來(lái)衰減干擾。
7)微動(dòng)。是由人為因素所導(dǎo)致的一系列不規(guī)則的振動(dòng)或者大自然當(dāng)中的一些噪聲。這類(lèi)型的噪聲頻帶很寬(1 Hz~200 Hz),并且難以區(qū)分,對(duì)此可以通過(guò)小波變換與CEEMDAN進(jìn)行噪聲壓制。
對(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)過(guò)濾波器系數(shù)h0(n)、h1(n)加權(quán)求和就得到j(luò)+1長(zhǎng)度空間的剩余系數(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)可見(jiàn),選取小波系數(shù)和濾波手段,是壓制噪聲的關(guān)鍵。選用什么樣的小波基函數(shù)和處理等手段是決定小波降噪效果的關(guān)鍵。
EMD(經(jīng)驗(yàn)?zāi)B(tài)分解)是一種分析非平穩(wěn)非線性的自適應(yīng)方法[12],但是容易導(dǎo)致模態(tài)混淆問(wèn)題,所以Huang[13]提出了EEMD(集合經(jīng)驗(yàn)?zāi)B(tài)分解),通過(guò)添加白噪聲來(lái)讓EMD導(dǎo)致的模態(tài)混淆問(wèn)題減少,然而在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)白噪聲,再通過(guò)計(jì)算獲取每個(gè)IMF在比較少的次數(shù)下實(shí)現(xiàn)重構(gòu)地震信號(hào)并且達(dá)到幾乎沒(méi)有誤差,保證了分解的完整性同時(shí)也減少了計(jì)算的時(shí)間。
CEEMDAN分解如下:
1)通過(guò)對(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)過(guò)EMD分解的E1[ωi(t)]噪聲分量相加,然后在對(duì)其結(jié)果進(jìn)行EMD分解,即可以得到第二個(gè)IMF分量即:
(7)
式中:r1(t)為第一個(gè)殘差分量;E1[ωi(t)]為信號(hào)通過(guò)EMD分解得到的第一個(gè)模態(tài)函數(shù)。
3)當(dāng)IMFk[t],K=2、3、…、k時(shí),計(jì)算K的殘差便重復(fù)以上過(guò)程。
4)直到殘差不能在進(jìn)行分解便得到K個(gè)模態(tài)函數(shù)分量,其中K為完全自適應(yīng)無(wú)需人為設(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ì)算效率也較大地提高。
當(dāng)使用傳統(tǒng)的EMD類(lèi)方法用于去噪時(shí),根據(jù)以往的知識(shí),即隨機(jī)噪聲信號(hào)通常分布在高頻分量中,則直接丟棄最前的1個(gè)到2個(gè)高頻IMF分量。 再用剩下的IMF分量重構(gòu)獲得去噪后的信號(hào)。但是依靠經(jīng)驗(yàn)選擇幾個(gè)高頻IMF分量會(huì)都受主觀因素的影響,高頻分量舍棄過(guò)多或不足都會(huì)導(dǎo)致去噪效果的不佳。
故先通過(guò)CEEMDAN自適應(yīng)分解得到最優(yōu)的分解層數(shù),然后在計(jì)算高頻細(xì)節(jié)系數(shù)與分解的各個(gè)IMF分量之間的相關(guān)性系數(shù)的方法來(lái)確定高頻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之間,但通過(guò)實(shí)驗(yàn)發(fā)現(xiàn),個(gè)別的IMF分量在與高頻細(xì)節(jié)信號(hào)求得的相關(guān)系數(shù)會(huì)出現(xiàn)小于“0”的情況,由于欲通過(guò)各階IMF的相關(guān)系數(shù)來(lái)反應(yīng)噪聲信號(hào)在信號(hào)中的變化趨勢(shì),所以當(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),總體來(lái)說(shuō)是呈現(xiàn)單調(diào)遞減的函數(shù)關(guān)系,且它們各點(diǎn)之間的斜率穩(wěn)定且大體不會(huì)發(fā)生較大的變化。直到第i+1階IMF分量時(shí)斜率出現(xiàn)折點(diǎn)且驟然變小之后趨于平滑。這說(shuō)明在第1階到i階,IMF分量中主要為噪聲信號(hào)。從第i+1階開(kāi)始有效信號(hào)為主要信號(hào)。此時(shí)剔除前i個(gè)IMF分量,從第i+1階開(kāi)始重構(gòu)信號(hào)。
小波變換降噪方式種類(lèi)太多,大都采用軟閥值、硬閥子。對(duì)此閥子的選取對(duì)于信號(hào)的降噪表現(xiàn)有很大影響,如果選擇不好,要么消噪太多,要么就是不夠,所以選取一個(gè)合適的閥子便顯得格外重要。為了擺脫這種問(wèn)題筆者將小波分解和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)過(guò)小波分解后的高頻系數(shù)的相關(guān)性系數(shù)。
3)對(duì)求出的相關(guān)性系數(shù),分析系數(shù)之間的線性關(guān)系。
4)通過(guò)相關(guān)性分析,選擇斜率較大的數(shù)據(jù)進(jìn)行剔除,將保留下來(lái)的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ù)
為了體現(xiàn)本方法的可行性,筆者選用了地震實(shí)測(cè)數(shù)據(jù),數(shù)據(jù)來(lái)源于2018年南昌市多要素地質(zhì)調(diào)查項(xiàng)目。圖2是原始地震單道信號(hào),可以看出在單道信號(hào)中夾雜者較多噪聲干擾,不能明顯地分辨有效反射波的準(zhǔn)確位置。選擇bior2.4小波基[15]對(duì)其進(jìn)行小波分解。圖3是經(jīng)過(guò)小波分解后得到的地震低頻近似信號(hào)和各層的高頻細(xì)節(jié)信號(hào)。經(jīng)過(guò)小波分解后的高頻細(xì)節(jié)部分大都存在噪聲,所以在對(duì)前三層高頻信號(hào)進(jìn)行CEEMDAN分解,圖4是第一層高頻細(xì)節(jié)信號(hào)經(jīng)過(guò)CEEMDAN后自適應(yīng)分解得到的12個(gè)IMF分量。通過(guò)對(duì)分解得到的12個(gè)IMF分量求出各階自相關(guān)函數(shù)和相關(guān)性系數(shù),準(zhǔn)確定位高頻IMF分量。
圖5為各階IMF自相關(guān)函數(shù)。通過(guò)圖5的觀察發(fā)現(xiàn)IMF1的自相關(guān)函數(shù)圖像有用信號(hào)完全被噪聲淹沒(méi)。大量的噪聲存在于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)系。通過(guò)圖6可以更加直觀地反應(yīng)出各相關(guān)性系數(shù)之間的關(guān)系,IMF1為要剔除的高頻分量。在對(duì)剩下的IMF余量進(jìn)行相加在與小波低頻信號(hào)進(jìn)行重構(gòu),得到處理后信號(hào)見(jiàn)圖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é)果見(jiàn)表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ù)值。因此可以說(shuō)明地震信號(hào)通過(guò)小波-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)過(guò)小波-CEEMDAN濾波后信號(hào)能量得到加強(qiáng),消除了同相軸周邊的噪聲干擾,信噪比也大大提高,橫向分辨率更好。
1)對(duì)小波高頻細(xì)節(jié)系數(shù)使用CEEMDAN分解,再通過(guò)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í)測(cè)城市淺層地震記錄經(jīng)過(guò)小波-CEEMDAN濾波處理后可以看出,其對(duì)噪聲壓制取得了比較滿意的效果,并且準(zhǔn)確地反映了有效反射波的位置,數(shù)據(jù)的信噪比也有明顯地提升,橫向分辨率提高為以后的解釋研究工作也提供了便利。