王 磊

內(nèi)蒙古自治區(qū)地震局 內(nèi)蒙古 呼和浩特 010010

1 引言

基于傅里葉變換的數(shù)字信號處理是信號分析領(lǐng)域中的重要手段之一，也是地震學(xué)中分析監(jiān)測數(shù)據(jù)的經(jīng)典方法之一。傅里葉變換基于信號整體進行處理，將記錄為時間序列的信號變換至頻率域從而顯示信號整體所擁有的頻率特征。當(dāng)輸入信號的頻率成分不隨時間改變即平穩(wěn)信號時，利用此方法進行頻譜分析可以獲得正確的信號頻率參數(shù)，但當(dāng)輸入信號是頻率隨時間變化的非平穩(wěn)信號時，傅里葉變換并不能有效確定不同頻率的信號所出現(xiàn)的時間。

為了進一步分析非平穩(wěn)信號中頻率隨時間的變化規(guī)律，學(xué)者們嘗試在時間與頻率聯(lián)合域內(nèi)進行信號分析，簡稱時頻分析。目前的時頻分析手段有短時傅里葉變換、小波變換、Winger-Ville分布、Radon變換、Gabor變換等等多種分析方法，這些時頻分析方法都是以傅里葉變換為基礎(chǔ)，通過對輸入信號加窗、改變基函數(shù)等處理方式將輸入信號分解。由于傅里葉變換在處理非平穩(wěn)信號時因信號易出現(xiàn)虛假信號和假頻現(xiàn)象，因此上述方法在處理非平穩(wěn)信號時也會出現(xiàn)同樣的問題。

1998年，Huang首次提出了基于固有模態(tài)函數(shù)（Intrinsic Mode Function，簡稱IMF）利用經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition，簡稱EMD）和希爾伯特譜分析（Hilbert Spectrum Analysis）進行信號處理分析的方法——希爾伯特—黃變換（Hilbert-Huang Transform，簡稱HHT），為非平穩(wěn)信號的分析開辟了一條新的道路。隨后在Huang又對HHT的算法和譜分析作了詳細說明。HHT沒有參照傅里葉變換的理論基礎(chǔ)采用固定基函數(shù)對輸入信號進行分解，而是給出了一種新的信號分量（IMF）的定義，之后利用經(jīng)驗?zāi)B(tài)分解方法將輸入信號分解為若干個滿足定義的IMF分量，再對每個分量進行Hilbert變換得到信號的時頻信息。

本文先簡要說明希爾伯特黃變換的理論，隨后介紹其在地震監(jiān)測資料處理中的應(yīng)用。

2 希爾伯特—黃變換（HHT）理論簡述

自希爾伯特—黃變換提出以來，有許多學(xué)者對其理論和應(yīng)用都做了深入研究。由于Huang在提出HHT時是從實際應(yīng)用的角度出發(fā)而非先建立完善的理論體系后推廣到實際應(yīng)用，因此其具體理論仍處于完善階段。目前希爾伯特—黃變換的數(shù)據(jù)處理分析主要可分為兩部分：經(jīng)驗?zāi)B(tài)分解（EMD）和希爾伯特譜分析。

2.1 經(jīng)驗?zāi)B(tài)分解

希爾伯特黃變換的創(chuàng)新性之一在于提出了固有模態(tài)函數(shù)這一概念，具體定義為：①整個信號的零點數(shù)和極值點數(shù)相等或最多相差1個；②對于信號上任意一點，由極大值確定的上包絡(luò)線和由極小值確定的下包絡(luò)線的均值為0，即整個信號關(guān)于時間軸對稱。在進行時頻分析之前，需要對輸入信號作經(jīng)驗?zāi)B(tài)分解EMD，將原時間序列x(t)分解為若干IMF分量與殘余信號R(t)的和（如式2-1），經(jīng)驗?zāi)B(tài)分解的具體方法如下：

①求取時間序列信號x(t)的極值，并分別將極大值和極小值擬合成兩條光滑曲線作為信號上包絡(luò)線U(t)和下包絡(luò)線D(t)，之后依據(jù)式2-2求取上、下包絡(luò)線均值線m(t)：

②用x(t)減去m(t)得h1(t)，此時依據(jù)定義判斷h1(t)是否為IMF，若h1(t)不屬于IMF，則再對h1(t)作EMD直至得到滿足條件的IMF1(t)作為第一階IMF，

③用x(t)減去IMF1(t)得x1(t)，若x1(t)滿足殘余信號的要求，則EMD結(jié)束，否則再以x1(t)作為輸入信號重復(fù)上述過程，直至獲得R(t)為止。

2.2 希爾伯特譜分析

在通過EMD將原始信號分解為各階IMF與殘余信號R(t)之后，由于殘余信號在定義上僅代表數(shù)據(jù)整體變化趨勢，因此僅需要對各階IMF分量做希爾伯特變換并分析。對于任意階固有模態(tài)函數(shù)IMFi，其Hilbert變換為：

之后構(gòu)造信號：

其中幅值ai(t)與相位φi(t)分別為：

定義信號的瞬時頻率fi(t)為：

進而將原始信號x(t)轉(zhuǎn)化為有限個解析信號的和，即：

由式2-9可知，HHT利用IMF分量將原始信號分解為有限個解析信號之和，并且每個IMF信號的瞬時頻率及瞬時幅值信息都得以體現(xiàn)。而傅里葉變換一方面將信號分為無數(shù)解析信號之和，其理論基礎(chǔ)就決定了在實際應(yīng)用中并不能獲得輸入信號的全部信息；另一方面，對于每個分量僅能獲得對整體區(qū)域的貢獻，而不能清晰的了解不同時間段內(nèi)不同頻率信號的貢獻量。盡管目前短時傅里葉變換及小波變換等基于傅里葉理論的方法同樣能取得與HHT一樣的非平穩(wěn)信號處理效果，但是從理論基礎(chǔ)上來看HHT仍有很大的發(fā)展和提升空間。

3 HHT的發(fā)展與完善

希爾伯特—黃變換以其較強自適應(yīng)性和高效等優(yōu)點，受到國內(nèi)外學(xué)者廣泛關(guān)注。自Huang提出這一方法后，諸多學(xué)者也在針對方法存在的不足不斷完善和改進。鐘佑明詳細分析了HHT在理論基礎(chǔ)、包絡(luò)線擬合、篩選準則、端點效應(yīng)等方面的問題。由于固有模態(tài)函數(shù)IMF在提出時完全基于實際應(yīng)用，缺少理論方面的定義及論證，鐘佑明、秦樹人于2006年提出了Hilbert變換的局部乘積定理并從理論層面和物理含義層面分析論證，最終結(jié)合定理對IMF定義、EMD、瞬時頻率分析等部分作了解釋。雖然這一解釋在一定程度上填補了HHT的部分理論的空白，但HHT的理論體系仍需進一步完善和發(fā)展。

EMD的第一步便是依據(jù)信號極值確定輸入數(shù)據(jù)的上下包絡(luò)線，因此包絡(luò)線的擬合方法將直接影響到后續(xù)數(shù)據(jù)分析結(jié)果。EMD的包絡(luò)線擬合中實際上是在兩極值點間插值，要在避免引入額外分量的同時保證曲線高階光滑，因此Wu 和Huang綜合考慮后選擇使用光滑性較好的三次樣條插值方法擬合曲線。有不少學(xué)者提出了其他算法，包括埃爾米特( Hermite) 插值法，多項式擬合插值法、高次樣條插值法、多項式擬合法、分段冪函數(shù)插值法等，不少方法都取得了較為理想的結(jié)果。

作為獲得IMF的方法，EMD的篩選方法決定了各IMF的信號質(zhì)量及計算效率。為了解決EMD存在的模態(tài)混疊問題，Wu和Huang提出引入白噪聲用來緩解模態(tài)混疊問題的集合經(jīng)驗?zāi)B(tài)分析（ensemble empirical mode decomposition，簡稱EEMD），但這一方法引入了大量的額外計算步驟，增加了時間成本。Yeh 等在 2010 年提出了互補的集合經(jīng)驗?zāi)B(tài)分解算法（complementary ensemble empirical mode decomposition，簡稱CEEMD）。為了減少EEMD算法中邊界數(shù)據(jù)對整體數(shù)據(jù)的“污染”，有學(xué)者選擇在邊界的第一個極值點處增加正弦信號，該方法稱為 Slope-EEMD，簡稱S-EEMD，即帶斜度的集合經(jīng)驗?zāi)B(tài)分解。Huang于2013年提出了導(dǎo)數(shù)優(yōu)化的 EMD 算法( Derivativeoptimized EMD，DEMD)，改用埃爾米特多項式來獲得上下包絡(luò)的一階差分作為參數(shù)來進行優(yōu)化。

由于信號的兩個端點僅有一側(cè)的信息，無法確定是否為極值點，因此在作EMD時需要在盡量保持曲線變化特征的同時對數(shù)據(jù)進行合理延拓，以此來獲取完整的包絡(luò)線。目前主要的延拓方法有鏡像閉合延拓法，極點對稱延拓法，自回歸模型法，正交多項式擬合法及神經(jīng)網(wǎng)絡(luò)法等。

4 HHT在地震監(jiān)測數(shù)據(jù)處理中的應(yīng)用

目前在地震的監(jiān)測預(yù)測研究工作中主要使用兩種數(shù)據(jù)：一種是由數(shù)字地震計記錄的地震動位移或位移的加速度，最終轉(zhuǎn)換為地震波形數(shù)據(jù)統(tǒng)一存儲；另一類是由大地磁場、大地電場、地殼形變、地下水位、地下水溫、重力等多種間接監(jiān)測手段的數(shù)據(jù)組成，統(tǒng)稱為定點地球物理監(jiān)測數(shù)據(jù)。這些監(jiān)測數(shù)據(jù)中大部分因區(qū)域環(huán)境、儀器、人工等多種原因，存在不規(guī)則干擾及長期趨勢變化，屬于非平穩(wěn)數(shù)據(jù)。因此有不少學(xué)者使用HHT對這些數(shù)據(jù)進行分析處理時，取得了一定的成果。

天然地震波形數(shù)據(jù)主要用于地震發(fā)生的時間、位置、震級的判定，其中一項工作就是從波形中分辨出來自不同反射層的不同傳播類型的波，簡稱震相識別。馮紅武基于HHT的較高時頻分辨能力在震相自動識別方法和技術(shù)上做了一定研究。王玥琪使用HHT對陜、甘、寧、蒙交界地區(qū)寬頻帶數(shù)字地震儀記錄到的不同類型的ML2.0 級以上地震數(shù)據(jù)進行分析，結(jié)果顯示其能較好地描述和分辨爆破事件、天然地震事件、塌陷事件，以及臺站背景噪聲等優(yōu)勢頻率分布情況。賈瑞生等人則研究了如何自動拾取低信噪比波形數(shù)據(jù)的P波震相初至?xí)r間。

HHT在定點地球物理監(jiān)測數(shù)據(jù)的處理中同樣取得了不錯的效果，安張輝等利用HHT對2013年蘆山Ms7.0 地震前甘肅省平?jīng)雠_、四川省成都臺和云南省元謀臺地電場觀測資料進了分析研究，發(fā)現(xiàn)地震前1～2個月地電場時頻圖出現(xiàn)能量增強現(xiàn)象， 認為HHT分析增加了觀測到震前地電場異常變化的可能性。張國清等通過HHT對寶昌、大同地電阻率數(shù)據(jù)的處理實驗證明了HHT在地電阻率數(shù)據(jù)的降噪和異常信息提取中均能取得較好的效果。孫和平等采用EEMD方法提取了全球5個超導(dǎo)重力儀觀測數(shù)據(jù)的重力極潮，消除了儀器漂移誤差，并且通過與前人成果對比分析認為EEMD得到的結(jié)果正確切更具有實際物理意義。

5 結(jié)論與展望

1）希爾伯特黃變換在理論上有不同于傳統(tǒng)時頻分析方法的優(yōu)勢，在實際應(yīng)用中也體現(xiàn)出較好的干擾識別和信息提取能力。同時諸多學(xué)者對于其理論體系的完善和計算方法的改進也使得希爾伯特黃變換的數(shù)據(jù)處理效果和計算效率都有顯著的提升。從整體上看，希爾伯特黃變換有較好的發(fā)展前景

2）雖然目前有很多改進的算法和判別準則，但希爾伯特黃變換的理論體系尚不健全，對固有模態(tài)函數(shù)、瞬時頻率的定義等問題仍需進一步討論，因此該方法無論從理論還是實際應(yīng)用方面都有很大的改進和提升的空間。

3）目前大多數(shù)地震監(jiān)測數(shù)據(jù)都是帶有各種干擾信息的非平穩(wěn)信號，希爾伯特黃變換正是針對這類信號的時頻分析方法。隨著理論與實際應(yīng)用研究，一定可以將此方法進一步的應(yīng)用地震監(jiān)測預(yù)測工作中，獲得更多的應(yīng)用成果。