唐 杰,溫 雷,李 聰,戚瑞軒

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

由于微地震事件信號通常能量較弱,加之傳播過程中能量被地層吸收衰減,造成微地震事件淹沒在噪聲中,對后續(xù)的微地震事件精確定位產(chǎn)生了不利影響,因此提高微地震資料的信噪比是微地震數(shù)據(jù)處理和解釋的首要任務(wù)[1-2]。隨著定位精度和震源機(jī)制全矩張量反演需求的增加,對去噪技術(shù)的要求也逐步提高[3]。傳統(tǒng)壓制隨機(jī)噪聲的方法可分為空間域及變換域兩類,前者主要包括均值濾波、中值濾波及各向異性擴(kuò)散濾波等,后者主要包括傅里葉變換域?yàn)V波方法及基于小波變換、曲波變換等的閾值去噪方法等[4-5]。地面微地震資料具有強(qiáng)噪聲、弱有效信號的特點(diǎn),常規(guī)的去噪方法往往難以獲得較好的去噪效果[6]。針對經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)方法的不穩(wěn)定問題和模態(tài)混疊現(xiàn)象[7-9],總體經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)利用高斯白噪聲頻譜均勻分布的統(tǒng)計(jì)特性,在原始信號中加入不同的白噪聲,使得信號在不同尺度上具有連續(xù)性,但是該方法計(jì)算效率不高[10-11];完備總體經(jīng)驗(yàn)?zāi)B(tài)分解方法(CEEMD)通過加入正、負(fù)成對的輔助噪聲,能夠有效消除重構(gòu)信號中的殘余輔助噪聲,而且計(jì)算效率較EEMD有所提高[12],但是信號重構(gòu)精度會有所降低;而TORRES等[13]提出的改進(jìn)的完備總體經(jīng)驗(yàn)?zāi)B(tài)分解方法(ICEEMD)能夠精確重構(gòu)原始信號,有效減少虛假模態(tài)和模態(tài)中的噪聲,同時(shí)也能降低計(jì)算成本。

本文研究了基于ICEEMD的有效信號提取方法,將非平穩(wěn)信號分解為一系列相對平穩(wěn)的固有模態(tài)函數(shù),定義了一種自適應(yīng)間隔閾值,用來去除固有模態(tài)中的噪聲成分,通過重構(gòu)數(shù)據(jù)得到去噪后的結(jié)果,同時(shí)利用分解后的信號與Hilbert變換進(jìn)行初至檢測。最后用數(shù)值模擬數(shù)據(jù)及實(shí)際微地震資料驗(yàn)證了方法的有效性。

1 改進(jìn)的完備總體經(jīng)驗(yàn)?zāi)B(tài)分解閾值去噪

1.1 改進(jìn)的完備總體經(jīng)驗(yàn)?zāi)B(tài)分解(ICEEMD)

ICEEMD方法將待處理微地震信號看作初始數(shù)據(jù),在分解的每一階段添加一個(gè)特定白噪聲,計(jì)算一個(gè)固有模態(tài)函數(shù)(IMF)分量,然后得到一個(gè)殘差,進(jìn)而得到每個(gè)IMF分量,ICEEMD能自適應(yīng)地將一個(gè)復(fù)雜信號分解為一系列IMF分量,且該分量滿足從高頻到低頻的順序分布。ICEEMD克服了EMD原有的模態(tài)混疊等缺點(diǎn),保持了EMD的完備性,能夠進(jìn)行原始信號的精確重構(gòu),且有更好的收斂性[14-15]。

ICEEMD具體步驟如下。

2) 計(jì)算一級殘差r1=x-IMF1。

5) 計(jì)算k級殘差rk=rk-1-IMFk。

6) 重復(fù)步驟4)和步驟5)直至殘差不能被分解。

上述各式中,x表示輸入的原始信號;rk表示k級殘差;算子Ej(·)表示對給定信號通過EMD求得第j個(gè)模態(tài)分量;wi為單位方差均值為0的高斯白噪聲,i=1,2,…,m,m表示噪聲組總數(shù);xi=x+wi為加入不同噪聲后的信號;εk表示噪聲標(biāo)準(zhǔn)差,在每一個(gè)分解階段表示一個(gè)常數(shù),允許在每個(gè)階段選擇信噪比。

圖1分別給出了原始信號時(shí)域波形(全文需要分解的原始信號不做特殊說明均為各圖中的c0)及經(jīng)過EMD、EEMD、ICEEMD分解后的IMF分量。通過對比可以看出,EMD在分解信號時(shí)出現(xiàn)模態(tài)混疊,不能將各個(gè)模態(tài)分量有效分解出來,而EEMD與ICEEMD分解效果相近,均可以將不同模態(tài)分量有效地分解出來。

1.2 時(shí)頻域高斯噪聲分量檢測與去除

對主要包含高斯噪聲的IMF分量可以通過高階統(tǒng)計(jì)量(HOS)和Kurtosis準(zhǔn)則進(jìn)行檢測和去除,留下混合噪聲和信號的分量。包含N個(gè)數(shù)據(jù)的第i個(gè)IMF的Kurtosis峰度(Ki)可通過下式計(jì)算[16]:

(1)

式中:σIMFi和μIMFi分別是IMFi的標(biāo)準(zhǔn)差和平均值,i=1,2,…,K。區(qū)分高斯分布和非高斯分布的HOS準(zhǔn)則為:

(2)

式中:α是置信水平,RAVIER等[17]給出了優(yōu)化值α=90%。Ki若滿足(2)式即為高斯分布,需要在信號中濾除。該過程類似自動帶通濾波,去除地震信號頻帶相對較低和較高的頻率成分,通過此步驟可以去除高能量相干噪聲。

圖1 原始信號時(shí)域波形及分解后的IMF分量a EMD分解; b EEMD分解; c ICEEMD分解

1.3 自適應(yīng)間隔閾值去噪

閾值去噪方法的關(guān)鍵是閾值和閾值函數(shù)的選取。傳統(tǒng)的閾值函數(shù)主要有硬閾值函數(shù)和軟閾值函數(shù),若選取的閾值過大,雖然去噪徹底,但同時(shí)造成了有效信息的丟失;若選取的閾值過小,則導(dǎo)致去噪不徹底[18-19]。本文采用了自適應(yīng)間隔閾值:

(3)

式中:h(zi)是輸入信號相鄰零點(diǎn)之間的采樣間隔;h(ri)是這個(gè)間隔中的局部極值;i=1,2,…,I;0<γ<λ,0≤α≤1,γ是截?cái)嘀?截?cái)嘀狄韵聰?shù)值設(shè)為0;α決定了閾值函數(shù)的形狀。該閾值函數(shù)類似于硬閾值,但門檻值是平穩(wěn)過渡的,因此它可以視為硬閾值和軟閾值的組合。如果該高斯噪聲的方差為σ2,則對于J個(gè)采樣點(diǎn)的優(yōu)化閾值為[20]:

(4)

圖2a為含噪聲合成信號時(shí)域波形及分解后的IMF分量;圖2b為去噪后信號時(shí)域波形及分解后的IMF分量;圖3a為去噪前的時(shí)頻分布;圖3b為去噪后的時(shí)頻分布。合成記錄中含有大量噪聲,信噪比很低,對比去噪前后波形可看出,噪聲被基本去除。從分解后的模態(tài)分量可以看出,噪聲分量得到去除;從時(shí)頻分析結(jié)果也可以看出合成記錄中的噪聲得到有效去除,信噪比大大提高。

圖2 含噪聲合成信號時(shí)域波形及ICEEMD分解后的IMF分量(a)及去噪后信號時(shí)域波形及ICEEMD分解后的IMF分量(b)

圖4a為含噪聲微地震信號時(shí)域波形及分解后的IMF分量;圖4b為去噪后信號時(shí)域波形及分解后的IMF分量;圖5a為去噪前的時(shí)頻分布;圖5b為去噪后的時(shí)頻分布。所選微地震信號信噪比較低,無法有效識別初至。對比去噪前后波形可見,噪聲被基本去除,在分解后的各個(gè)分量中,噪聲分別得到有效去除;從時(shí)頻分析結(jié)果也可以看出去噪效果良好,噪聲基本得到去除,信噪比得到極大提高。

圖3 含噪聲合成信號去噪效果分析a 去噪前時(shí)頻分布; b 去噪后時(shí)頻分布

圖4 含噪聲微地震信號時(shí)域波形及分解后的IMF分量(a)及去噪后時(shí)域波形及分解后的IMF分量(b)

圖5 含噪聲微地震信號去噪效果分析a 去噪前時(shí)頻分布; b 去噪后時(shí)頻分布

2 基于Hilbert變換的信號檢測

微地震記錄的特點(diǎn)是頻率高、信噪比低,因此微地震事件的自動識別和初至到時(shí)拾取對實(shí)現(xiàn)海量微地震數(shù)據(jù)的自動處理具有重要意義。微地震事件的自動識別和初至到時(shí)拾取主要根據(jù)地震記錄中地震事件到達(dá)前后質(zhì)點(diǎn)振動性質(zhì)的差別構(gòu)建特征函數(shù)來實(shí)現(xiàn)[22-23]。劉玉海等[24]利用相鄰道有效信號相關(guān)性好、與隨機(jī)噪聲相關(guān)性差的特性來識別有效信號,壓制隨機(jī)噪聲。譚玉陽等[25]針對低信噪比事件容易被遺漏的情況,提出利用多道相似系數(shù)來檢測微地震事件,在實(shí)際資料的應(yīng)用中具有較高的準(zhǔn)確率。本文結(jié)合ICEEMD與Hilbert變換進(jìn)行信號初至檢測,利用ICEEMD將信號分解為一系列IMF分量,然后通過Hilbert變換將信號變換到頻率域,噪聲與信號可以更好地分開,進(jìn)而可以更加有效地進(jìn)行微地震有效信號檢測。

對分解出的分量,利用Hilbert變換獲得每個(gè)點(diǎn)的振幅(Hi):

(5)

式中:n=1,2,…,N,N表示分解出的IMF分量個(gè)數(shù);AMP表示取振幅;i=1,2,…,I,I表示信號大小。

然后計(jì)算持續(xù)能量比(ER1):

(6)

式中:L是τ測試點(diǎn)附近的能量求和窗口的長度。定義一個(gè)關(guān)于高能量地震信號的參數(shù)(ER2):

(7)

在獲得了ER2(τ)后,設(shè)定一個(gè)閾值,通過給定的閾值找到ER2的局部極大值,大于該閾值則說明存在微地震初至信息。閾值被設(shè)置為最大峰值的分?jǐn)?shù),定義了檢測靈敏度。

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

圖6為選取的某地面微地震資料和利用傳統(tǒng)的濾波方法和本文方法分別處理后的結(jié)果以及所對應(yīng)的時(shí)頻分布。其中圖6a是實(shí)際地震數(shù)據(jù),其信噪比比較低,無法確定準(zhǔn)確到達(dá)時(shí)間;圖6b和圖6c分別顯示了利用帶通濾波與本文方法去噪后的結(jié)果。帶通濾波可以去除某一頻帶外的噪聲,與有效信號頻率接近的噪聲無法完全去除;而本文方法去除了大部分噪聲,提高了信號的信噪比,使所選事件的信噪比更高。圖6d到圖6f給出了利用短時(shí)傅里葉變換將圖6a 到圖6c對應(yīng)信號變換到時(shí)頻域的結(jié)果。對比時(shí)頻分析結(jié)果可見,帶通濾波僅去除了特定頻率范圍之外的噪聲,而本文方法去除了大部分噪聲,信噪比得到了提高。圖7為去噪后的微地震數(shù)據(jù)以及對應(yīng)的特征函數(shù),由圖7可以看出,去噪后的有效信號(圖7a) 與其特征函數(shù)(圖7b)的峰值呈一一對應(yīng)關(guān)系,可以有效地檢測微地震信號初至。

圖6 對實(shí)際微地震資料利用帶通濾波與本文方法去噪效果對比a 實(shí)際微地震數(shù)據(jù); b 帶通濾波去噪結(jié)果; c 本文方法去噪結(jié)果; d,e,f 為對應(yīng)a,b,c的時(shí)頻分布

圖7 去噪后微地震數(shù)據(jù)(a)及對應(yīng)的特征函數(shù)(b)

4 結(jié)論

本文研究了基于ICEEMD的微地震噪聲壓制與初至檢測方法,理論模型和實(shí)際數(shù)據(jù)的應(yīng)用取得了很好的效果,得出的主要結(jié)論有:

1) ICEEMD方法可以將復(fù)雜信號有效地分解為一系列固有模態(tài)函數(shù),克服了EMD方法存在的模態(tài)混疊問題。閾值去噪中閾值的大小會影響去噪效果與有效信號的保護(hù),針對傳統(tǒng)閾值大小不好確定的問題定義了自適應(yīng)間隔閾值,可以看做是硬閾值與軟閾值的組合,與ICEEMD結(jié)合的去噪方法可以有效地壓制噪聲,數(shù)值模擬資料與實(shí)際微地震資料的應(yīng)用結(jié)果表明該方法在微地震信號噪聲壓制方面效果明顯,提高了數(shù)據(jù)信噪比,有較好的應(yīng)用價(jià)值;

2) ICEEMD與Hilbert變換相結(jié)合的信號檢測方法可以顯著提高微地震事件的探測精度,利于后續(xù)震源定位及反演。