楊光亮

(1)中國地震局地震研究所,武漢 430071 2)地殼運(yùn)動實(shí)驗(yàn)室,武漢 430071)

一種去除地震背景噪音的新方法＊

楊光亮1,2)

(1)中國地震局地震研究所,武漢 430071 2)地殼運(yùn)動實(shí)驗(yàn)室,武漢 430071)

結(jié)合希爾伯特-黃變換方法,根據(jù)信號背景噪音時間延續(xù)性假設(shè),提出了一種自適應(yīng)地震信號去噪算法。利用該算法實(shí)現(xiàn)了臺站地震監(jiān)測信號的去噪分析。分析表明：1)該去噪算法能根據(jù)背景信號特性有效去除信號的低頻干擾;2)該算法可自適應(yīng)完成信號分解去噪,計(jì)算效率高,無時間分辨率和頻率分辨率問題;3)該算法處理高頻干擾存在缺陷,有待進(jìn)一步完善。

模態(tài)分解;固有模態(tài)函數(shù);相關(guān)分析;特性時間延續(xù);去噪

1 引言

天然地震監(jiān)測臺網(wǎng)接收到的有用信號淹沒在大量的干擾噪聲中,有些地震由于震級較小、信號信噪比低,將導(dǎo)致無法識別。微小地震反映殼內(nèi)部的微破裂,如果能精確監(jiān)測到微破裂的發(fā)展過程,對地震預(yù)警將起到重要作用。另外,由于大量的干擾噪音也影響對地震 P波初至的識別。因此,對地震監(jiān)測弱信號的提取、去噪的研究具有重要實(shí)際意義。

地震監(jiān)測信號是非平穩(wěn)非各態(tài)歷經(jīng)的隨機(jī)信號。目前廣泛用于地震信號去噪的方法,如傅立葉頻譜分析、功率譜分析、小波變換等方法,都是將地震信號近似為平穩(wěn)隨機(jī)信號進(jìn)行處理。非平穩(wěn)信號的主要特征是其頻率是時間的函數(shù)。這種時變頻率雖然可以通過已有的時頻聯(lián)合分析方法,如短時傅立葉變換 (STFT)、W igner-Ville分布、小波變換(WT)等進(jìn)行分析,但這些時頻分析法的最終理論依據(jù)都是傅立葉變換,在本質(zhì)上還無法擺脫 Fourier變換的局限性,而且受 Heisenberg不確定原理的限制。

希爾伯特-黃變換 (HHT)是一種新的數(shù)據(jù)分析方法[1],它從根本上突破了傅立葉變換理論的限制,首次建立了一種基于瞬時頻率的信號分析方法,它可以比較精確地描繪出非平穩(wěn)信號中頻率隨時間的變化規(guī)律,更容易避免假頻和多余信號分量等冗余問題的出現(xiàn)。它依據(jù)數(shù)據(jù)本身的時間尺度特征將信號分解為有限個固有模態(tài)函數(shù) (I MF,然后對各模態(tài)分量進(jìn)行變換從而得到信號能量在時間尺度上的分布規(guī)律,實(shí)現(xiàn)信號特征的量化提取。相對于傳統(tǒng)的數(shù)據(jù)分析方法,HHT有完全自適應(yīng)性,能處理非線性非平穩(wěn)數(shù)據(jù),不受 Heisenberg測不準(zhǔn)原理制約等優(yōu)點(diǎn),在客觀性和分辨率方面都具有明顯的優(yōu)越性,能提取到更多、更接近實(shí)際的信號特性[2-5]。

2 HHT去噪模型及計(jì)算算法

地震信號具有短時、突變等特點(diǎn),是一種典型的非平穩(wěn)隨機(jī)信號。對于地震監(jiān)測臺站接收到的信號,可以認(rèn)為是背景噪音與有用信號的混合。假設(shè)背景噪音在一定時間范圍內(nèi)是平穩(wěn)的,并假設(shè)其信號特征在時間上有一定的延續(xù)性,據(jù)此可以將背景噪音作為信號基底 (圖 1),通過對基底和分析區(qū)信號的特征分析,可以尋找到兩者的內(nèi)在關(guān)聯(lián)。

圖1 信號基底與分析區(qū)示意圖Fig.1 Schematic diagram of the backround signal and the analyzed area

HHT算法依據(jù)數(shù)據(jù)本身的時間尺度特征可以將信號分解為有限個固有模態(tài)函數(shù) (I MF),經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)的依據(jù)為：

1)整個信號中,零點(diǎn)數(shù)與極點(diǎn)數(shù)相等或至多相差 1。這一限制條件近似于傳統(tǒng)平穩(wěn)高斯過程中關(guān)于窄帶的定義。

根據(jù)上述設(shè)想,我們設(shè)計(jì)了如下地震監(jiān)測信號去噪算法,算法流程如圖 2所示。該算法需要解決3個關(guān)鍵問題：1)EMD分解時,包絡(luò)插值算法的選擇;2)EMD分解終止條件的選擇;3)相關(guān)閾值的設(shè)置。這里忽略了邊界效應(yīng)的討論。Huang提出HHT時,采用三次樣條插值,但在實(shí)際使用過程中容易出現(xiàn)插值“過沖”和“欠沖”。這里我們采用線性插值算法計(jì)算包絡(luò)線,并采用 Rilling[6]提出的終止條件終止模態(tài)分解。

圖2 信號去噪流程Fig.2 Flowchart of denoising

相關(guān)閾值與信號長度密切相關(guān)。理論上,當(dāng)分析信號和基底較長時,相互間的信號特征延續(xù)性較差,當(dāng)信號較短時,則延續(xù)性較好。因此,相關(guān)閾值的選擇應(yīng)根據(jù)信號長度及當(dāng)?shù)刭Y料特性,經(jīng)多次試驗(yàn)來獲得。

3 地震臺站監(jiān)測信號震例分析

圖3與圖 4中的信號是上海某地震監(jiān)測臺一段連續(xù)垂直向地震信號,按上述算法分成基底區(qū)和分析區(qū)兩部分,作為分析數(shù)據(jù),信號采樣頻率 100 Hz,信號總長度 600 s,分段后基底區(qū)與分析區(qū)均為 300 s。

圖3 基底原始信號Fig.3 Original background signal

圖4 待去噪的原始信號Fig.4 Orignal signal to be denoised

通過 EMD分解獲得了基底和分析信號的模態(tài)(圖 5)。基底分解后得到 13個基本模態(tài)和一個余量,即趨勢項(xiàng)。分析信號分解得到 11個模態(tài)和一個余量。原始信號經(jīng)過模態(tài)分解后,得到從高頻到低頻的一系列模態(tài)分量 imf1,imf2,i mf3,imf4,…,這里并不是說 imf1的頻率一定比 imf2的高,而是 imf1中的某個局部的頻率比 imf2中相同局部的頻率要高,在相同時間的局部,imf1頻率一定大于 imf2。至于分解過程造成的誤差,主要是包絡(luò)方式的選取、邊界效應(yīng)的處理和濾波停止條件的設(shè)計(jì),會不斷累積到下一層分解中,并不一定是最后一個余量 (趨勢項(xiàng))。

根據(jù) EMD模態(tài)分解結(jié)果,各模態(tài)頻率成分由高到低排列,如果我們能從分析信號中識別噪音干擾,就可以據(jù)此對信號進(jìn)行去噪。只需要將干擾信號模態(tài)直接剔出,然后將其他模態(tài)信號累加,就可得到去噪后的結(jié)果。在這種意義上說,模態(tài)分解相當(dāng)于一個濾波器。

圖5 基底信號與分析信號模態(tài)分解結(jié)果Fig.5 Mode decomposition results of background signal and the analysis signal

圖6 基底與分析區(qū)信號相關(guān)系數(shù)Fig.6 Correlation coefficients between signals of background and analyzed area

圖7 噪音頻譜Fig.7 Noise spectrum

通過對分解出來的各固有模態(tài)作相關(guān)分析,計(jì)算得到相關(guān)系數(shù)關(guān)系如圖 6所示。相關(guān)系數(shù)有 4個大于 0.3的局部峰值,其中 2個集中在 0.3以上,另2個集中在 0.6以上,而其他峰值相對較小,因此,設(shè)置相關(guān)閾值分別為 0.3和 0.6時,可以得到 2組噪音模態(tài),圖 7是兩組噪音模態(tài)的頻譜圖。相關(guān)閾值為 0.6時,剔除分析信號中最后 2個模態(tài),其他模態(tài)通過模態(tài)累加就可以復(fù)合出去噪后信號,當(dāng)相關(guān)閾值設(shè)置為 0.3時,剔除最后 4個模態(tài),其他模態(tài)累加,可得到另一組復(fù)合后的分析信號(圖 8)。

4 討論與結(jié)論

該去噪算法是在假設(shè)背景噪音在一定時間范圍內(nèi)是平穩(wěn)的,并假設(shè)其信號特征在時間上有一定的延續(xù)性而作出的。一般認(rèn)為時間序列越長,其信號特性延續(xù)性越差。為此,我們另外取了連續(xù)的兩段信號利用上述算法進(jìn)行了計(jì)算,相關(guān)系數(shù)關(guān)系如圖8所示,x軸和 y軸分別為信號分解后的模態(tài)序數(shù)。從圖 8中可看到相關(guān)峰值明顯增多,相關(guān)平面不再平坦,其中相關(guān)系數(shù)大于 0.3的有 5個。該段分析信號去噪后時域和頻域結(jié)果如圖 9。通過實(shí)際資料的處理,表明該算法利用背景基底特征,能有效去除觀測信號中的部分噪音干擾。

從以上處理結(jié)果看,去除的干擾信息都集中在低頻端,這可能是低頻干擾信號更穩(wěn)定,持續(xù)時間更長,更易被識別。由此可知,該去噪算法對高頻端的處理還有待進(jìn)一步探索。

通過以上實(shí)例的分析,初步得出結(jié)論：1)該去噪算法能根據(jù)背景信號特性有效去除信號的低頻干擾;2)該算法去噪可自適應(yīng)完成,計(jì)算效率高;3)該算法對高頻干擾存在缺陷,有待進(jìn)一步完善。

圖8 去噪后基底與分析區(qū)信號相關(guān)系數(shù)Fig.8 Correlation coefficients between signals of background and analyzed area after denoising

圖9 10s信號去噪結(jié)果及頻譜圖Fig.9 Spectrum and denoising results of 10 s signal

1 HuangN E.The empiricalmode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Procedures of the Royal Cociety of London, 1998,A454(3)：903-995.

2 Qiuhui Chen and Norden Huang.A B-spline approach for empirical mode decompositions[J].Advances in ComputationalMathematics,2006,24：171-195.

3 Norden E Huang Zhengshen and Steven R Long.A new view of nonlinearwaterwaves：the Hilbert spectrum[J].Annual Review of FluidMechanics,1999,31：417-457.

4 張立,等.基于 HHT提取重力固體潮的地震前兆信息[J].地震學(xué)報(bào),2007,29(2)：222-226.(ZhangLi,et al. HHT-based extraction of gravity tide of earthquake precursor information[J].Acta Seismologica Sinica,2007,29(2)：222-226)

5 楊志華,等.一種基于 HHT的信號周期性分析方法及應(yīng)用[J].中山大學(xué)學(xué)報(bào) (自然科學(xué)版),2005,44(2)：14-18.(Yang Zhihua,et al.Method and applications of periodic signal based on HHT[J].Journal of Sun Yat-sen University (science),2005,44(2)：14-18)

6 Gabriel Rilling and Paulo Gon?alvés.Empirical mode decomposition as a filter bank[J].IEEE Signal ProcessingLetters,2004,11(2)：112-114.

SEISM IC SIGNAL DENO ISING ALGORITHM BASED ON HHT AND ITS APPL ICATI ON

Yang Guangliang1,2)

(1)Institute of Seism ology,CEA,W uhan 430071 2)CrustalM ovem ent Laboratory,W uhan 430071)

On the basis of the Hilbert-Huang Transform(HHT)method and according to the assumption of the background noise signal continuity in time domain,an adaptive algorithm of seis mic signal denoising is presented. The algorithm,as used to analyze the seis mic wave,has certain positive results.Through the analysis following conclusions can be drawn.1)According to the background noise signal characteristics,the algorithm can effectively remove the low-frequency noise signal from the original signal.2)By the use of the adaptive algorithm we can decompose the signal and denoising,and avoid the ti me resolution and frequency resolution problem.3)The algorithm has the flaw for dealingwith high-frequency interferenc,and so needs to be improved.

EMD(EmpiricalMode Decomposition);I MF(Intrinsic Mode Function);correlation analysis;characteristics of the time extension;denoising

1671-5942(2011)04-0090-04

2011-03-17

中國地震局地震研究所重點(diǎn)基金(IS200916004)

楊光亮,男,博士,主要從事重力動態(tài)變化與重力殼幔結(jié)構(gòu)反演研究工作.E-mail：vforyang@gmail.com

P207

A