王 磊
內蒙古自治區地震局 內蒙古 呼和浩特 010010
基于傅里葉變換的數字信號處理是信號分析領域中的重要手段之一,也是地震學中分析監測數據的經典方法之一。傅里葉變換基于信號整體進行處理,將記錄為時間序列的信號變換至頻率域從而顯示信號整體所擁有的頻率特征。當輸入信號的頻率成分不隨時間改變即平穩信號時,利用此方法進行頻譜分析可以獲得正確的信號頻率參數,但當輸入信號是頻率隨時間變化的非平穩信號時,傅里葉變換并不能有效確定不同頻率的信號所出現的時間。
為了進一步分析非平穩信號中頻率隨時間的變化規律,學者們嘗試在時間與頻率聯合域內進行信號分析,簡稱時頻分析。目前的時頻分析手段有短時傅里葉變換、小波變換、Winger-Ville分布、Radon變換、Gabor變換等等多種分析方法,這些時頻分析方法都是以傅里葉變換為基礎,通過對輸入信號加窗、改變基函數等處理方式將輸入信號分解。由于傅里葉變換在處理非平穩信號時因信號易出現虛假信號和假頻現象,因此上述方法在處理非平穩信號時也會出現同樣的問題。
1998年,Huang首次提出了基于固有模態函數(Intrinsic Mode Function,簡稱IMF)利用經驗模態分解(Empirical Mode Decomposition,簡稱EMD)和希爾伯特譜分析(Hilbert Spectrum Analysis)進行信號處理分析的方法——希爾伯特—黃變換(Hilbert-Huang Transform,簡稱HHT),為非平穩信號的分析開辟了一條新的道路。隨后在Huang又對HHT的算法和譜分析作了詳細說明。HHT沒有參照傅里葉變換的理論基礎采用固定基函數對輸入信號進行分解,而是給出了一種新的信號分量(IMF)的定義,之后利用經驗模態分解方法將輸入信號分解為若干個滿足定義的IMF分量,再對每個分量進行Hilbert變換得到信號的時頻信息。
本文先簡要說明希爾伯特黃變換的理論,隨后介紹其在地震監測資料處理中的應用。
自希爾伯特—黃變換提出以來,有許多學者對其理論和應用都做了深入研究。由于Huang在提出HHT時是從實際應用的角度出發而非先建立完善的理論體系后推廣到實際應用,因此其具體理論仍處于完善階段。目前希爾伯特—黃變換的數據處理分析主要可分為兩部分:經驗模態分解(EMD)和希爾伯特譜分析。
希爾伯特黃變換的創新性之一在于提出了固有模態函數這一概念,具體定義為:①整個信號的零點數和極值點數相等或最多相差1個;②對于信號上任意一點,由極大值確定的上包絡線和由極小值確定的下包絡線的均值為0,即整個信號關于時間軸對稱。在進行時頻分析之前,需要對輸入信號作經驗模態分解EMD,將原時間序列x(t)分解為若干IMF分量與殘余信號R(t)的和(如式2-1),經驗模態分解的具體方法如下:

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

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

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

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

之后構造信號:

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

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

進而將原始信號x(t)轉化為有限個解析信號的和,即:

由式2-9可知,HHT利用IMF分量將原始信號分解為有限個解析信號之和,并且每個IMF信號的瞬時頻率及瞬時幅值信息都得以體現。而傅里葉變換一方面將信號分為無數解析信號之和,其理論基礎就決定了在實際應用中并不能獲得輸入信號的全部信息;另一方面,對于每個分量僅能獲得對整體區域的貢獻,而不能清晰的了解不同時間段內不同頻率信號的貢獻量。盡管目前短時傅里葉變換及小波變換等基于傅里葉理論的方法同樣能取得與HHT一樣的非平穩信號處理效果,但是從理論基礎上來看HHT仍有很大的發展和提升空間。
希爾伯特—黃變換以其較強自適應性和高效等優點,受到國內外學者廣泛關注。自Huang提出這一方法后,諸多學者也在針對方法存在的不足不斷完善和改進。鐘佑明詳細分析了HHT在理論基礎、包絡線擬合、篩選準則、端點效應等方面的問題。由于固有模態函數IMF在提出時完全基于實際應用,缺少理論方面的定義及論證,鐘佑明、秦樹人于2006年提出了Hilbert變換的局部乘積定理并從理論層面和物理含義層面分析論證,最終結合定理對IMF定義、EMD、瞬時頻率分析等部分作了解釋。雖然這一解釋在一定程度上填補了HHT的部分理論的空白,但HHT的理論體系仍需進一步完善和發展。
EMD的第一步便是依據信號極值確定輸入數據的上下包絡線,因此包絡線的擬合方法將直接影響到后續數據分析結果。EMD的包絡線擬合中實際上是在兩極值點間插值,要在避免引入額外分量的同時保證曲線高階光滑,因此Wu 和Huang綜合考慮后選擇使用光滑性較好的三次樣條插值方法擬合曲線。有不少學者提出了其他算法,包括埃爾米特( Hermite) 插值法,多項式擬合插值法、高次樣條插值法、多項式擬合法、分段冪函數插值法等,不少方法都取得了較為理想的結果。
作為獲得IMF的方法,EMD的篩選方法決定了各IMF的信號質量及計算效率。為了解決EMD存在的模態混疊問題,Wu和Huang提出引入白噪聲用來緩解模態混疊問題的集合經驗模態分析(ensemble empirical mode decomposition,簡稱EEMD),但這一方法引入了大量的額外計算步驟,增加了時間成本。Yeh 等在 2010 年提出了互補的集合經驗模態分解算法(complementary ensemble empirical mode decomposition,簡稱CEEMD)。為了減少EEMD算法中邊界數據對整體數據的“污染”,有學者選擇在邊界的第一個極值點處增加正弦信號,該方法稱為 Slope-EEMD,簡稱S-EEMD,即帶斜度的集合經驗模態分解。Huang于2013年提出了導數優化的 EMD 算法( Derivativeoptimized EMD,DEMD),改用埃爾米特多項式來獲得上下包絡的一階差分作為參數來進行優化。
由于信號的兩個端點僅有一側的信息,無法確定是否為極值點,因此在作EMD時需要在盡量保持曲線變化特征的同時對數據進行合理延拓,以此來獲取完整的包絡線。目前主要的延拓方法有鏡像閉合延拓法,極點對稱延拓法,自回歸模型法,正交多項式擬合法及神經網絡法等。
目前在地震的監測預測研究工作中主要使用兩種數據:一種是由數字地震計記錄的地震動位移或位移的加速度,最終轉換為地震波形數據統一存儲;另一類是由大地磁場、大地電場、地殼形變、地下水位、地下水溫、重力等多種間接監測手段的數據組成,統稱為定點地球物理監測數據。這些監測數據中大部分因區域環境、儀器、人工等多種原因,存在不規則干擾及長期趨勢變化,屬于非平穩數據。因此有不少學者使用HHT對這些數據進行分析處理時,取得了一定的成果。
天然地震波形數據主要用于地震發生的時間、位置、震級的判定,其中一項工作就是從波形中分辨出來自不同反射層的不同傳播類型的波,簡稱震相識別。馮紅武基于HHT的較高時頻分辨能力在震相自動識別方法和技術上做了一定研究。王玥琪使用HHT對陜、甘、寧、蒙交界地區寬頻帶數字地震儀記錄到的不同類型的ML2.0 級以上地震數據進行分析,結果顯示其能較好地描述和分辨爆破事件、天然地震事件、塌陷事件,以及臺站背景噪聲等優勢頻率分布情況。賈瑞生等人則研究了如何自動拾取低信噪比波形數據的P波震相初至時間。
HHT在定點地球物理監測數據的處理中同樣取得了不錯的效果,安張輝等利用HHT對2013年蘆山Ms7.0 地震前甘肅省平涼臺、四川省成都臺和云南省元謀臺地電場觀測資料進了分析研究,發現地震前1~2個月地電場時頻圖出現能量增強現象, 認為HHT分析增加了觀測到震前地電場異常變化的可能性。張國清等通過HHT對寶昌、大同地電阻率數據的處理實驗證明了HHT在地電阻率數據的降噪和異常信息提取中均能取得較好的效果。孫和平等采用EEMD方法提取了全球5個超導重力儀觀測數據的重力極潮,消除了儀器漂移誤差,并且通過與前人成果對比分析認為EEMD得到的結果正確切更具有實際物理意義。
1)希爾伯特黃變換在理論上有不同于傳統時頻分析方法的優勢,在實際應用中也體現出較好的干擾識別和信息提取能力。同時諸多學者對于其理論體系的完善和計算方法的改進也使得希爾伯特黃變換的數據處理效果和計算效率都有顯著的提升。從整體上看,希爾伯特黃變換有較好的發展前景
2)雖然目前有很多改進的算法和判別準則,但希爾伯特黃變換的理論體系尚不健全,對固有模態函數、瞬時頻率的定義等問題仍需進一步討論,因此該方法無論從理論還是實際應用方面都有很大的改進和提升的空間。
3)目前大多數地震監測數據都是帶有各種干擾信息的非平穩信號,希爾伯特黃變換正是針對這類信號的時頻分析方法。隨著理論與實際應用研究,一定可以將此方法進一步的應用地震監測預測工作中,獲得更多的應用成果。