基于改進(jìn)的希爾伯特-黃變換的地震強震動特征提取

郭亞茹,豐繼華,楊黎薇,鐘玉盛

(1.云南省地震局,昆明 650224;2.云南民族大學(xué),昆明 650031)

0 引言

近半個世紀(jì)以來,伴隨新出現(xiàn)時頻分析方法被不斷成功地應(yīng)用于地球物理、環(huán)境科學(xué)和地質(zhì)學(xué)等領(lǐng)域,現(xiàn)代時頻分析理論得到了長足發(fā)展。在信號頻譜研究和分析過程中,使用最多的是快速傅里葉變換(Fast Fourier Transform,FFT)。根據(jù)信號處理過程,可將以傅里葉理論為基礎(chǔ)的時頻分析方法歸納為三類[1]：第一類首先對信號加窗,其次進(jìn)行傅里葉變換,前提條件為假設(shè)信號局部平穩(wěn),典型代表為STFT和Gabor變換等[2-3];第二類是直接更換傅里葉變換的基函數(shù),典型代表為小波變換等[4-5];第三類是從信號中獲得雙線性函數(shù),然后執(zhí)行傅里葉變換,但存在交叉項的問題,如Winger-Ville分布等。

傅里葉分析的理論是將信號分解為正弦信號的加權(quán)和,當(dāng)信號僅由幾個有限的正弦信號組成時,使用傅里葉分析是可行的,但如果信號沒有規(guī)律,使用傅里葉分解時需要大量的虛假的正弦信號去拼湊。此外,對于非線性分布的信號,始終存在來自交叉項的干擾。再者,受到不確定性原理的限制,傳統(tǒng)以傅里葉理論為基礎(chǔ)的時頻分析方法在時間和頻率分辨率并不理想,所以需要一種能自適應(yīng)準(zhǔn)確反映頻率隨時間變化規(guī)律的新方法。1998年,美國華裔科學(xué)家HUANG博士將經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)算法引入到希爾伯特變換中,提出了希爾伯特-黃變換(Hilbert-Huang Transform,簡稱“HHT”)。該方法被認(rèn)為是近二個世紀(jì)應(yīng)用數(shù)學(xué)領(lǐng)域中最重大的發(fā)現(xiàn),一經(jīng)提出,便受到了研究者們的普遍關(guān)注。

HHT旨在服務(wù)于非平穩(wěn)和非線性信號分析,如在地震波分析、腦電波分析研究以及機械故障檢測等。在地震研究領(lǐng)域內(nèi),希爾伯特-黃變換的獨特性在于分析信號過程中不需要基函數(shù)[6-7],其基本原理是將原始信號自適應(yīng)分解為多個本征模態(tài)函數(shù)(IMF),再將IMF進(jìn)行Hilbert變換得出具有物理意義的物理參量。2001年HUANG[8]利用HHT對1999年集集地震中Tcu129地震臺站的加速度記錄進(jìn)行分析,與短期傅里葉變換和Morlet小波分析相比,得到的希爾伯特譜能準(zhǔn)確反映地震動能量的時頻分布;2003年ZHANG等[9]使用 EMD對地震波傳播進(jìn)行研究;2009年劉強等[10]通過HHT變換計算結(jié)構(gòu)地震反應(yīng)的能量;2012年WANG等[11]運用EMD方法對地震信號進(jìn)行時頻分析,證實了該算法的有效性;2016年梁岳[12]將改進(jìn)HHT運用于儲層預(yù)測性能中。

本文選取了部分云南漾濞Ms6.4地震中不同震中距的強震動臺站所獲取的三分量加速度記錄數(shù)據(jù),基于改進(jìn)后的希爾伯特黃算法對其進(jìn)行時頻分析,以不同角度客觀分析強震動的尺度和頻域特征。通過對實際地震動記錄進(jìn)行特征提取后,獲取的希爾伯特時頻譜信號可清晰反映出強震能量隨時頻變化分布特征的能量譜,同時實現(xiàn)了時頻局部化的目的。與梁宏等[13]利用模態(tài)混疊問題及集合經(jīng)驗?zāi)Ｊ椒纸夥椒?EEMD)分析九寨溝M7.0地震加速度記錄相比,本文改進(jìn)后的模態(tài)混疊問題及集合經(jīng)驗?zāi)Ｊ椒纸夥椒?彌補了其幅值系數(shù)不確定性的不足,其對地震強震動特征的提取過程可為類似破壞性淺源地震的強震動特征分析提供更有效的信息參數(shù),為現(xiàn)代信號應(yīng)用分析于強震動特征提供另一種新思路。

1 希爾伯特-黃變換的改進(jìn)

1.1 引用數(shù)據(jù)來源

根據(jù)中國地震臺網(wǎng)正式測定,北京時間2021年5月21日21時48分,云南省大理州漾濞縣發(fā)生Ms6.4級地震,震中位于東經(jīng)99.87°,北緯25.67°漾濞縣蒼山西鎮(zhèn)境內(nèi),震源深度為8 km[14]。地震共計造成了3人遇難,34人受傷,部分房屋破裂甚至倒塌,同時云南多地州有強烈震感。此次地震是云南繼Ms6.5魯?shù)榈卣鸷蚆s6.6景谷地震后近10 a以來發(fā)生的又一次破壞性淺源地震。震中位于滇西北地區(qū)維西-喬后斷裂帶附近,震源機制結(jié)果顯示：此次地震是一次走滑型破裂事件,與區(qū)域構(gòu)造特征一致。從地震序列的類型上可以看出：此次漾濞地震屬于“前震-主震-余震型”地震[15-18]。

云南地震局強震動臺網(wǎng)獲取強震動數(shù)據(jù)25組75條強震動數(shù)據(jù)(三分向記錄)。獲取強震動數(shù)據(jù)的25個強震動臺站,土層臺有23個,基巖臺有2個。按照國家強震動觀測常規(guī)數(shù)據(jù)的處理流程對25個臺站的記錄進(jìn)行數(shù)據(jù)處理[19-21]。為了闡明提取記錄加速度特性效果,本文根據(jù)峰值加速度以及震中等因素,分別選取最具代表性的4個臺站：震中距為8.0 km的漾濞臺(YBX)、震中距為31.1 km的月溪井臺(DLY)、震中距為66.8 km的太和臺(BTH)、震中距為129.6 km的施甸臺(SDX)。其中：太和臺(BTH)為基巖臺。圖1為該地震中獲取加速度記錄的強震動臺站位置分布,具體記錄情況詳見表1。

圖1 云南強震動臺站位置分布圖Fig. 1 Location distribution map of strong motion stations in Yunnan

表1 云南漾濞Ms6.4強震動記錄Table 1 Ms6.4 strong ground motion records of Yangbi,Yunnan

1.2 希爾伯特-黃變換簡要概述

希爾伯特-黃變換(HHT)主要包含兩個部分：經(jīng)驗?zāi)B(tài)分解(EMD)與希爾伯特變換(Hilbert transform)[22]。其原理為：先利用EMD方法分解原始信號得出本征模態(tài)函數(shù),再將本征模態(tài)函數(shù)進(jìn)行希爾伯特變換得出具有實際物理意義的參量。

1) 經(jīng)驗?zāi)B(tài)分解(EMD)基本原理

EMD是希爾伯特-黃變換的核心步驟,可將原始混合信號分解為一系列本征模態(tài)函數(shù)(IMF)組件。EMD不同于快速傅里葉變換(FFT),其以原有信號為基礎(chǔ)直接進(jìn)行分解,不需要基函數(shù),適用于任意形式的信號,能夠應(yīng)用于多種領(lǐng)域。EMD的分解過程一般基于以下條件進(jìn)行：

(a)用于分解的數(shù)據(jù)信號是由幾種不同的本征模態(tài)函數(shù)(IMF)構(gòu)成的,至少含有兩個以上的極值點。

(b)IMF同時具有線性與非線性的特點。

(c)如果數(shù)據(jù)信號沒有極值點只有拐點,可以利用二次或者多次求導(dǎo)取得極值,再進(jìn)行積分。

EMD分解法是先通過尋找原始信號的極大極小值,分別取得極大值包絡(luò)mmax與極小值包絡(luò)mmin,再求其包絡(luò)均值m(t)[22]。設(shè)原始信號為x(t),h(t)為準(zhǔn)模態(tài)函數(shù),在第一次分解中的第一次篩分為：

m1(t)=(mmax(t)+mmax(t))/2

(1)

h1(t)=x(t)-m1(t)

(2)

如果第一次篩分結(jié)果h1(t)不滿足IMF定義,便開始第二次篩分,此時h1(t)作為被篩分信號重復(fù)以上過程,直到滿足黃鍔博士等提出的停止準(zhǔn)則,即可跳出篩分循環(huán),完成第一次分解,得到第一階本征模態(tài)函數(shù)c1(t),采用的停止準(zhǔn)則為：

(3)

當(dāng)Sd滿足0.2≤Sd≤0.3,即可跳出第一次分解,獲得第一階本征模態(tài)函數(shù)c1(t)：

c1(t)=h1k(t)

(4)

此時令：r1(t)=x(t)-c1(t)作為第二次分解的原始信號進(jìn)行篩分循環(huán)：

r2(t)=r1(t)-c2(t),……,rn(t)=rn-1(t)-cn(t)

當(dāng)rn為一個常數(shù)或者一個單調(diào)遞增(遞減)的函數(shù)再或者只含有一個極點的函數(shù)時結(jié)束分解,此時rn為剩余量,

(5)

在IMF生成過程中,EMD在每個時間段提取頻率最高的分量。因此,每個組成部分都不能間歇,否則,將發(fā)生模態(tài)混疊現(xiàn)象。

2) 希爾伯特變換(Hilbert transform)基本概念

假設(shè)本征模態(tài)函數(shù)x(t)的希爾伯特變換為?i(t),那么：

(6)

式中：P表示柯西(Cauchy)主值。根據(jù)這個定義,ci(t)和?i(t)形成一個復(fù)共軛對,可以得到一個解析信號zi(t)：

zi(t)=ci(t)+j?i(t)=ai(t)·ejθi(t)

(7)

式中：

(8)

式中：ai(t)為反映ci(t)能量隨時間變化的瞬時幅度,θi(t)為ci(t)的瞬時相位。而相位求導(dǎo)即可獲得瞬時頻率,如公式(9)所示：

(9)

原始數(shù)據(jù)x(t)(不包含余量rn(t))可表示為：

(10)

希爾伯特振幅譜H(ω,t)可表示為：

(11)

由希爾伯特振幅譜H(ω,t)進(jìn)而獲得：邊際譜h(ω)和瞬時能量密度級IE：

(12)

(13)

3) 模態(tài)混疊問題及集合經(jīng)驗?zāi)Ｊ椒纸?EEMD)原理

HHT最明顯的問題是由過包絡(luò)、欠包絡(luò)、模式混合以及干擾信號頻域特征提取的噪聲導(dǎo)致的模態(tài)混疊。集合經(jīng)驗?zāi)Ｊ椒纸?Ensemble Empirical Mode Decomposition,EEMD)由HUANG等針對模態(tài)混疊現(xiàn)象提出的,其原理為：利用相位-幅度耦合方法來量化高頻帶和低頻帶之間的相互作用,分析獨立成份消除模態(tài)混合,改善每個本征模態(tài)分量的正交性[23]。

設(shè)x(t)為分析信號。EEMD的原理過程如下：

(a)信號xi(t)為：

xi(t)=x(t)+βω(i)(t)

(14)

式中：β為增加的白噪聲的方差,ω(i)(i=1,2,……)i是EMD傳導(dǎo)的數(shù)量,表示均值為0,方差為1的白噪聲。

(c)計算每個最終的IMF：

(15)

4) 改進(jìn)后的模態(tài)混疊問題及集合經(jīng)驗?zāi)Ｊ椒纸?EEMD)原理及方法

模態(tài)混疊問題及集合經(jīng)驗?zāi)Ｊ椒纸夥椒?EEMD)具有幅值系數(shù)不確定性,EEMD中EMD執(zhí)行的總次數(shù)也是添加噪聲的次數(shù)為N,添加白噪聲的幅值系數(shù)為ε。WU和HUANG曾在論文中討論：幅度的選擇會直接影響分解效果,如果增加的噪聲幅度太小,其可能不會引入EMD所依賴的極端變化,如果增加的噪聲幅度太大,則分解結(jié)果可能不是原始數(shù)據(jù),并且可能出現(xiàn)冗余IMF。根據(jù)算法提出的研究人員的建議與實驗經(jīng)驗,建議在實驗過程中選用N=100,ε取值范圍為0.01～0.5較為合適。但是ε取值的具體情況應(yīng)根據(jù)實驗中的具體情況而定,并沒有遵循某種特定的規(guī)則。

由此嘗試在兩個強震記錄的加速度信號中加入不同幅值系數(shù)(信噪比)的白噪聲,通過更改所加白噪聲的標(biāo)準(zhǔn)差觀察相關(guān)系數(shù)的變化,進(jìn)而得出最佳標(biāo)準(zhǔn)差下最能消除模態(tài)混疊現(xiàn)象的白噪聲的分布。

在本文中,經(jīng)過對比分析發(fā)現(xiàn)對于連續(xù)信號,正態(tài)分布信號,具有線性關(guān)系的信號,皮爾遜相關(guān)系數(shù)作為評判指標(biāo)最為合適,當(dāng)然用斯皮爾曼相關(guān)系數(shù)也可以,但是就效率而言沒有皮爾遜相關(guān)系數(shù)高。所以最終我們采用了皮爾遜相關(guān)系數(shù),即ρX,Y定義為：

(16)

式中：cov(X,Y)是X和Y的協(xié)方差,σX為X的標(biāo)準(zhǔn)差,σY為Y的標(biāo)準(zhǔn)差,μX、μY分別表示X和Y的數(shù)學(xué)期望,所以相關(guān)系數(shù)可以表示為：

(17)

由圖2加入噪聲后的相關(guān)系數(shù)曲線可以直觀得到對于漾濞臺(YBX)、月溪井臺(DLY)、太和臺(BTH)和施甸臺(SDX)各個分量的最佳系數(shù),選取結(jié)果即表2所列。為便于敘述將添加特定幅值系數(shù)白噪聲的EMD定義為AEEMD。為更為直觀比較EMD、EEMD與AEEMD分解效果,以太和臺(BTH)為例,對太和臺(BTH)EW向量進(jìn)行分解,分別計算每種方法分解得到的各個IMF的方差σ與貢獻(xiàn)率mσ,以此作為模態(tài)混疊現(xiàn)象消除程度的評判標(biāo)準(zhǔn)。

圖2 加入噪聲后的相關(guān)系數(shù)曲線Fig.2 Correlationcoefficientcurveafteraddingnoise表2 噪聲系數(shù)選取Table2 Noisefigureselection臺站代碼方向最佳系數(shù)YBXEW0.30NS0.50UD0.50DLYEW0.43NS0.45UD0.49BTHEW0.46NS0.46UD0.38SDXEW0.38NS0.50UD0.11

第k層本征模態(tài)函數(shù)ck(t)的方差σk為：

σk=E(ck(t)2)-(E(ck(t)))2

(18)

則第k層IMF的貢獻(xiàn)率為：

(19)

式中：σk與mσk越大說明所得IMF分量與原始信號最為相近,分解效果越好,反之效果則越差。分解效果如圖3、表3所示。

表3 三種分解方法中IMF的方差(BTHEW)Table 3 IMFs’ variance of three decomposition methods (BTHEW)

圖3 分解譜(BTHEW)Fig. 3 Decomposition spectrum(BTHEW)

圖3(a)為EMD分解結(jié)果,信號分為10層,圖中Data為原始信號,C1-C10為10層IMF,R10為分解10層之后留下來的殘余信息。從圖中可以看到IMF1到IMF10都包含原始信號,我們無法分辨哪個層具有最原始信號的信息。

圖3(b)為EEMD分解結(jié)果,添加根據(jù)算法提出的研究人員的建議與實驗經(jīng)驗：100次幅值系數(shù)為0.2的高斯白噪聲。

圖3(c)為AEEMD即調(diào)整幅值系數(shù)為0.46后分解結(jié)果。

表3-表4所列為三種分解方法中IMF的方差與貢獻(xiàn)率。對于EMD方法：表3-表4均顯示IMF2的方差最大與貢獻(xiàn)率最高,僅從理論可知IMF1應(yīng)與原始信號最為相近,因此最為合理的解釋為EMD方法出現(xiàn)了模態(tài)混合現(xiàn)象。對于EEMD方法：IMF1方差最大和貢獻(xiàn)率最高,應(yīng)該具有原始信號信息最多,圖3(b)中IMF1與原始信號最為相似恰巧印證這一點。AEEMD方法同樣顯示IMF1方差最大和貢獻(xiàn)率最高,并且AEEMD中IMF1的貢獻(xiàn)率高于EEMD分解方法,所以AEEMD相較于EEMD更為突出的改善了模態(tài)混疊現(xiàn)象。

表4 三種分解方法中IMF的貢獻(xiàn)率(BTHEW)Table 4 IMFs’ variance of three decomposition methods (BTHEW)

上述實驗表明：添加的白噪聲幅值系數(shù)與所分解數(shù)據(jù)有關(guān),只有添加適當(dāng)?shù)陌自肼?才可以最大程度改善模態(tài)混疊現(xiàn)象,提高分解的準(zhǔn)確性。

2 特征提取

Hilbert-Huang變換(HHT)類似于傅里葉變換,兩者均是將時域中的信號進(jìn)行解析,以便于頻域分析。在此,我們用希爾伯特-黃變換對強震動加速度進(jìn)行分析：一是彌補傅里葉變換在信號時頻特性提取中的不足;二是通過HHT量化提取信號能量的時頻分布特征。

2.1 時域分解特征

圖4表示經(jīng)過改進(jìn)后希爾伯特黃變換之后的漾濞臺(YBX)、月溪井臺(DLY)、太和臺(BTH)和施甸臺(SDX)各個分量的三維采樣點-IMF層數(shù)-幅值譜,圖形將4個臺站每分量的采樣點-IMF層數(shù)-幅值進(jìn)行了立體展現(xiàn),從圖中可以看出：在任意一層本征模態(tài)函數(shù)均存在不同的頻率成分,各頻率成分的幅值(即其的相對含量)也不盡相同,與每個IMF分量對應(yīng)的瞬時頻率之間往往沒有明顯的界限;區(qū)別于小波分析中頻帶必須嚴(yán)格劃分,并且隨著震中距的逐漸增大,每個本征模態(tài)函數(shù)的振幅波峰差距逐漸減小,而基巖臺加速度記錄相對更為平穩(wěn)。獲取的希爾伯特時域譜信號可清晰反映出強震能量隨時頻變化分布特征的能量譜,同時實現(xiàn)了時頻局部化的目的,為類似破壞性淺源地震的強震動特征分析提供更有效的信息參數(shù)。

圖4 三維采樣點-IMF層數(shù)-幅值譜Fig. 4 3D sample point-IMF layer number-amplitude spectrums

2.2 Hilbert邊際譜分析

圖5是經(jīng)過改進(jìn)希爾伯特黃變換后得到的漾濞臺(YBX)、月溪井臺(DLY)、太和臺(BTH)和施甸臺(SDX)各個分量的邊際譜與FFT譜歸一化后的對比圖。以YBXEW為例,圖5 YBXEW為漾濞臺東西分量處理結(jié)果,震中距為8.0 km。對比FFT 譜和HHT邊際譜可以直觀發(fā)現(xiàn)：FFT譜在低頻時會低估信號幅值,而在高頻時放大信號幅值。漾濞臺(YBX)其他兩個分量(北南向NS、垂直向UD)、月溪井臺(DLY)三分量(東西EW、南向NS、垂直向UD)、施甸臺(SDX)三分量(東西EW、南向NS和垂直向UD)同樣印實了這一點。

圖5 Hilbert邊際譜與傅里葉譜對比Fig. 5 Comparison of Hilbert Marginal Spectrum and Fourier Spectrum

對比4個站點Hilbert邊際譜與傅里葉譜可以發(fā)現(xiàn)：隨著震中距的增加,最大地震動的峰值加速度出現(xiàn)明顯的衰減,而加速度記錄的主頻分量也出現(xiàn)相應(yīng)增加。與梁宏等[13]利用EEMD分析九寨溝M7.0地震加速度記錄相比,抑制模態(tài)混疊現(xiàn)象效果更為突出,經(jīng)本文改進(jìn)過的希爾伯特黃變換所得到的Hilbert邊際譜與傅里葉譜對比圖,更為清晰刻畫出在震中距影響下FFT譜對信號低頻低估高頻放大現(xiàn)象的變化：隨著震中距不斷增加FFT譜與邊際譜差異性逐漸變小。

2.3 Hilbert變換與譜分析

圖6為經(jīng)過改進(jìn)希爾伯特黃變換后得到的漾濞臺(YBX)、月溪井臺(DLY)、太和臺(BTH)和施甸臺(SDX)三分量的譜分析結(jié)果二維顯示,其中太和臺(BTH)為基巖臺,另外三個為土層臺,時頻譜橫坐標(biāo)為頻率,縱坐標(biāo)為HHT的瞬時能量,圖像顯示了能量分布特征：當(dāng)震中距為50 km以內(nèi)時能量集中于20 Hz左右,隨著震中距的增大,最大加速度峰值出現(xiàn)急速衰減,能量的分布規(guī)律也顯現(xiàn)出特定形態(tài)。聯(lián)系時域分解特征與邊際譜特征分析,該方法在強震記錄數(shù)據(jù)處理時頻分析方面展現(xiàn)了特有的實效性。

圖6 HHT時頻譜二維顯示Fig. 6 2D HHT spectrums

3 結(jié)語

本文選取了云南漾濞Ms6.4地震中不同震中距的強震動臺站所獲取的三分量加速度記錄數(shù)據(jù),采用了基于改進(jìn)后的AEEMD算法對其進(jìn)行時頻分析。通過分析可知：

1) 該方法可在保證信號分解穩(wěn)定性的同時,抑制模態(tài)混疊現(xiàn)象。通過本文對強震動數(shù)據(jù)的分析擬合,不但成功將加速度數(shù)據(jù)分解為本征模態(tài)函數(shù),還獲取了希爾伯特時頻譜信號。

2) 本文獲得的邊際譜可更加準(zhǔn)確的反映信號的實際頻率成分,更適于處理非平穩(wěn)信號。

3) 本文獲得的能量譜可清晰反映出強震能量隨時頻變化分布特征,實現(xiàn)了時頻同時局部化的目的。

經(jīng)分析對比發(fā)現(xiàn)：傳統(tǒng)的FFT譜普遍存在著低頻時低估信號幅值并且高頻時又過度放大信號幅值的問題,本文所改進(jìn)的方法能在一定程度上修正這種誤差。實驗結(jié)果表明：基于改進(jìn)的希爾伯特-黃變換在地震強震動特征提取過程中具有顯著優(yōu)點,與前人利用EEMD分析地震加速度記錄相比,抑制模態(tài)混疊現(xiàn)象效果更為突出,其得到的邊際譜和時頻譜可在最大程度上保留原始數(shù)據(jù)中的主要特性,可為類似破壞性淺源地震的強震動特征分析提供更有效的信息參數(shù),更適于處理非平穩(wěn)信號,為現(xiàn)代信號應(yīng)用分析于強震動特征提供另一種新思路。