希爾伯特—黃變換及其在地震監(jiān)測(cè)數(shù)據(jù)處理中的應(yīng)用

王 磊

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

1 引言

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

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

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

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

2 希爾伯特—黃變換（HHT）理論簡(jiǎn)述

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

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

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

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

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

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

2.2 希爾伯特譜分析

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

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

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

定義信號(hào)的瞬時(shí)頻率fi(t)為：

進(jìn)而將原始信號(hào)x(t)轉(zhuǎn)化為有限個(gè)解析信號(hào)的和，即：

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

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

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

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

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

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

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

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

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

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

5 結(jié)論與展望

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

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

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