基于CEEMDAN·MPE-NHT的爆破地震波信號時頻分析*

孫 苗,吳 靜,吳 立,楊鈞凱,覃亞男

(1.湖北國土資源職業(yè)學(xué)院 環(huán)境與工程學(xué)院,武漢 430090;2.中國地質(zhì)大學(xué)(武漢) a.巖土鉆掘與防護教育部工程研究中心;b.工程學(xué)院,武漢 430074;3.湖北工程學(xué)院 土木工程學(xué)院,孝感 432000)

希爾伯特-黃變換(Hilbert-Huang Transform,HHT)由經(jīng)驗?zāi)B(tài)分解(Ensemble Empirical Mode,EMD)[1]和希爾伯特變換(Hilbert Transform,HT)[2]組成,其中EMD以其能識別數(shù)據(jù)內(nèi)在屬性并根據(jù)數(shù)據(jù)本身特征進行分解在信號分析領(lǐng)域得到了廣泛應(yīng)用[3]。對EMD分解結(jié)果進行Hilbert變換,更是實現(xiàn)了將時域信號轉(zhuǎn)變?yōu)轭l域信號,幫助爆破工程技術(shù)人員對爆破地震波信號傳播特征、內(nèi)在屬性進行識別,對爆破地震波危害控制起到了一定的指導(dǎo)作用[4-6]。

但EMD對含有噪聲的爆破地震波信號相對敏感,而實際爆破地震波監(jiān)測又無法避免噪聲信號的混入[7],這便造成了實際爆破地震波信號EMD得到的固有模態(tài)函數(shù)(Intrinsic Mode Function,IMF)[8,9]存在嚴重模態(tài)混淆[10-12],而Hilbert變換處理這類IMF分量會產(chǎn)生錯誤的時頻分析結(jié)果,這樣的結(jié)果對爆破地震波危害效應(yīng)分析意義不大,有時候甚至?xí)蔀樽璧K爆破工程技術(shù)人員作出判斷的干擾項。因此如何對HHT處理含噪爆破地震波信號遇到導(dǎo)致分析精度受損的問題進行處理,是目前爆破地震波危害效應(yīng)分析亟待解決的問題。

鑒于此,對影響EMD-Hilbert時頻分析精度的因素進行逐一改進。對EMD進行改進,使其可降低對噪聲的敏感性。具體實現(xiàn)通過兩步,首先在EMD中添加自適應(yīng)白噪聲,得到自適應(yīng)補充集合經(jīng)驗?zāi)B(tài)分解(Complete Ensemble Empirical Mode Decomposition with Adaptive Noise,CEEMDAN)[13-15],用于抵抗監(jiān)測中混入噪聲的低頻部分;再者引入多尺度排列熵(Multiscale Permutation Entropy,MPE)[16-18]到CEEMDAN中,充分發(fā)揮MPE隨機性檢測能力,抑制高頻噪聲對分解精度的影響。以上兩步得到的IMF可認為是降噪處理后的IMF,對CEEMDAN·MPE得到的IMF進行歸一化Hilbert變換(Normalized Hilbert transform,NHT)[19,20]使得傳統(tǒng)的Hilbert變換不受Bedrosian定理的約束,降低負值瞬時頻率出現(xiàn)的概率。通過上述三步可實現(xiàn)HHT時頻分析精度提升。

最后進行含噪仿真爆破振動信號和實測含噪水下鉆孔爆破地震波信號CEEMDAN·MPE-NHT時頻分析算法和HHT算法對比研究,以驗證CEEMDAN·MPE-NHT算法能有效提高HHT時頻精度,得到反映真實爆破振動屬性的時頻特征參數(shù),對爆破振動特征識別和制定科學(xué)的防護措施具有一定的指導(dǎo)作用。

1 CEEMDAN·MPE-NHT時頻分析算法

CEEMDAN·MPE-NHT時頻分析算法可分兩步構(gòu)建,第一步建立CEEMDAN·MPE算法,其后對CEEMDAN·MPE算法得到的IMF進行NHT。

CEEMDAN·MPE算法的核心是通過剔除噪聲來抑制EMD模態(tài)混淆。它將MPE的隨機性檢測功能與CEEMDAN的自適應(yīng)性相結(jié)合,可有效地減少噪聲的干擾。有關(guān)CEEMDAN和MPE的詳細介紹,請參閱文獻[13-15]和[16-18]。在這里,僅通過流程圖1分析CEEMDAN·MPE的操作過程。

圖1 CEEMDAN·MPE算法流程圖Fig. 1 CEEMDAN·MPE algorithm flow chart

Hilbert變換得到具有實際物理意義的瞬時頻率通常需滿足十分苛刻的要求,模態(tài)混淆的IMF分量往往不能滿足要求。鑒于此Huang等提出NHT[20],以此來提高瞬時頻率的計算精度。

NHT的原始操作是對IMF的調(diào)頻分量進行Hilbert變換,IMF來源是EMD或集合經(jīng)驗?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)。不難發(fā)現(xiàn),本文得到的IMF物理意義較EMD和EEMD得到的IMF更加清晰,真實性更高。將文獻[20]中的IMF替換成CEEMDAN·MPE得到的IMF,便是本文的實現(xiàn)途徑。進一步分析,得到CEEMDAN·MPE-NHT時頻分析算法運算流程圖,見圖2。

圖2 CEEMDAN·MPE-NHT算法流程圖Fig. 2 CEEMDAN·MPE-NHT algorithm flow chart

2 仿真信號CEEMDAN·MPE-NHT時頻分析

為突出CEEMDAN·MPE-NHT時頻分析算法相比HHT可有效提高含噪爆破地震波信號時頻特征參數(shù)的提取精度,特進行CEEMDAN·MPE-NHT和HHT的含噪仿真振動信號時頻分析對比研究。

含噪仿真振動信號S(t)=x1(t)+x2(t),其中x1(t)=wgn(1,N,0.1),即功率為0.1的噪聲信號,如圖3所示;x2(t)=sin(2×pi×100×t),即頻率為100 Hz的正弦信號,如圖4所示;仿真信號如圖5所示。采樣點數(shù)N=512,采樣時間t=1/N∶1/N∶1。

圖3 功率為0.1的噪聲信號圖Fig. 3 Noise signal with power of 0.1

圖4 頻率為100 Hz的正弦信號圖Fig. 4 Sine signal with frequency of 100 Hz

為突出CEEMDAN·MPE對EMD模態(tài)混淆的抑制作用,分別采用EMD和CEEMDAN·MPE對S(t)進行模態(tài)分解。圖6為EMD分解結(jié)果,不難發(fā)現(xiàn)IMF1是難以除去的噪聲信號;高頻模態(tài)混淆嚴重,如IMF2在0.53 s,0.64 s及0.85 s附近均出現(xiàn)了高頻向中頻發(fā)散的趨勢;中頻模態(tài)混淆有所緩解,如IMF3在0.41 s及右端點附近存在向低頻發(fā)展的趨勢;低頻相對穩(wěn)定,如IMF5和IMF6。圖7為CEEMDAN·MPE分解結(jié)果,分解得到的IMF從高頻向低頻依次排列,未見明顯模態(tài)混淆現(xiàn)象,模態(tài)顯示相比EMD結(jié)果清晰且穩(wěn)定。

圖7 仿真信號經(jīng)CEEMDAN·MPE得到的IMFFig. 7 IMFs of simulated signal by CEEMDAN·MPE

為突出NHT對瞬時頻率具備清晰物理意義的貢獻,進行Hilbert變換和NHT對比分析。對EMD得到的IMF進行Hilbert變換,得到的時頻譜見圖8。對CEEMDAN·MPE得到的IMF進行NHT,得到的時頻譜見圖9。

圖8 EMD-HT得到的仿真信號時頻譜圖Fig. 8 Time-frequency spectrum of simulated signal obtained by EMD-HT

圖9 CEEMDAN·MPE-NHT得到的仿真信號時頻譜圖Fig. 9 Time-frequency spectrum of simulated signal obtained by CEEMDAN-INHT

圖8為EMD-HT時頻譜,模態(tài)混淆分量經(jīng)Hilbert變換得到的時頻譜分辨率不高,出現(xiàn)了100 Hz以上難以識別的虛假分量。圖9為CEEMDAN·MPE-NHT時頻譜,該算法通過CEEMDAN·MPE來抑制S(t)中的x1(t);CEEMDAN·MPE得到的IMF經(jīng)NHT處理后得到的頻率分布在0～100 Hz,和x2(t)具有對應(yīng)性,時頻參數(shù)可清晰識別。側(cè)面說明了CEEMDAN·MPE得到的IMF經(jīng)過NHT能夠得到具有實際物理意義的時頻信息,即CEEMDAN·MPE-NHT時頻分析算法不僅可有效抑制噪聲引起的EMD模態(tài)混淆,同時得到時頻分辨率雙高的信號頻譜圖。

3 炸礁爆破地震波信號CEEMDAN·MPE-NHT時頻分析

以三峽-葛洲壩兩壩間蓮沱河段航道整治炸礁工程為研究對象,選擇下岸溪～丁頭鎮(zhèn)區(qū)間LT5炸礁區(qū)為研究對象。該區(qū)總的炸礁工程量2561.0 m3,清碴工程量2484.2 m3。LT5區(qū)周邊環(huán)境復(fù)雜,炸礁區(qū)距離最近的民房70 m、距離公路(陡紙線)81 m,還有滑坡、崩塌和泥石流等不良地質(zhì)情況,對工程影響較大。

選取炸礁爆破施工期間處于較危險地帶的民房作為研究對象,對民房進行現(xiàn)場監(jiān)測得到一系列監(jiān)測數(shù)據(jù)。觀察監(jiān)測數(shù)據(jù)發(fā)現(xiàn)水下鉆孔爆破地震波信號三個方向振動速度呈現(xiàn)出有規(guī)律的一致性,水平徑向峰值振動速度最大,水平切向次之,而垂直方向振動速度最小?？紤]文章篇幅,僅對水下鉆孔爆破地震波監(jiān)測信號的水平徑向峰值振動速度進行CEEMDAN·MPE-NHT時頻分析。有代表性水平徑向地震波監(jiān)測信號如圖10所示。

圖10 水平徑向地震波監(jiān)測信號波形圖Fig. 10 Waveform of monitored radial seismic signal

對圖10信號進行CEEMDAN·MPE-NHT時頻分析。先對圖10信號進行CEEMDAN·MPE,得到如圖11所示的IMF分量。可發(fā)現(xiàn),IMF從高頻到低頻排列,每個IMF都攜帶了爆破地震波信號的特定時頻信息。其中IMF3攜帶圖10所示信號的主要能量,其次是IMF4和IMF5。

圖11 水平徑向地震波信號CEEMDAN·MPE結(jié)果Fig. 11 The results of CEEMDAN·MPE of the radial seismic signal

進一步分析,對圖11中的IMF1～IMF6執(zhí)行NHT,獲得每個IMF的時頻能量特征圖,其中圖12～圖17為對應(yīng)IMF的邊際譜;圖18～圖23為對應(yīng)IMF的時頻譜;這里未對IMF7和R進行變換,是因為結(jié)合圖11不難發(fā)現(xiàn)IMF7和R占有的能量可忽略不計。觀察圖12～圖23可發(fā)現(xiàn),CEEMDAN·MPE-NHT獲得的時頻譜在時域和頻域都具有高分辨率,這點和仿真信號時頻分析得到的結(jié)果一致。IMF1的頻率范圍為100～250 Hz,持續(xù)時間為0.24～0.38 s;IMF2的頻率范圍為50～100 Hz,持續(xù)時間為0.23～0.42 s;IMF3的頻率范圍為20～50 Hz,持續(xù)時間為0.22～0.47 s;IMF4的頻率范圍為10～30 Hz,持續(xù)時間為0.28～1.17 s;IMF5的頻率范圍為5～10 Hz,持續(xù)時間為0.00～1.20 s;IMF6頻率最低,約為0～5 Hz,持續(xù)時間為0.00～1.20 s。

圖12 IMF1邊際譜Fig. 12 Marginal spectrum of IMF1

圖13 IMF2邊際譜Fig. 13 Marginal spectrum of IMF2

圖14 IMF3邊際譜Fig. 14 Marginal spectrum of IMF3

圖15 IMF4邊際譜Fig. 15 Marginal spectrum of IMF4

圖16 IMF5邊際譜Fig. 16 Marginal spectrum of IMF5

圖17 IMF6邊際譜Fig. 17 Marginal spectrum of IMF6

圖18 IMF1時頻譜Fig. 18 Time-frequency spectrum of IMF1

圖19 IMF2時頻譜Fig. 19 Time-frequency spectrum of IMF2

圖20 IMF3時頻譜Fig. 20 Time-frequency spectrum of IMF3

圖21 IMF4時頻譜Fig. 21 Time-frequency spectrum of IMF4

圖22 IMF5時頻譜Fig. 22 Time-frequency spectrum of IMF5

圖23 IMF6時頻譜Fig. 23 Time-frequency spectrum of IMF6

通過圖12～圖23的分析,可發(fā)現(xiàn)信號在低頻停留時間大于高頻,主要能量集中在50 Hz以下。

進一步分析,獲得圖10中信號的三維時頻能量譜,如圖24所示。從圖24可看出,本次水下炸礁工程監(jiān)測得到的水下鉆孔爆破地震波信號的主頻為26.83 Hz,次主頻為19.32 Hz。

圖24 信號時間-頻率-能量三維圖Fig. 24 Time-frequency-energy spectrum of the signal

根據(jù)結(jié)構(gòu)抗震知識,當(dāng)爆破地震波的頻率與房屋的固有頻率相同時,結(jié)構(gòu)的振幅將達到最大,即發(fā)生共振危險。物體的受迫振動是根據(jù)特定規(guī)律進行的,即由于形狀和結(jié)構(gòu)的不同,它們具有不同的固有頻率。為了探索建筑結(jié)構(gòu)的固有頻率,有必要分析結(jié)構(gòu)的振動特性,如振動特征頻率與頻率對應(yīng)的振型,即模態(tài)分析。通過PKPM結(jié)構(gòu)軟件的SETWE模塊,對距離爆區(qū)最近的2層房屋進行了模態(tài)分析。房屋的第一個6階振型圖如圖25所示,每個振型對應(yīng)的固有頻率如表1所示。

表1 民房前6階陣型對應(yīng)的自振頻率(單位:Hz)Table 1 Natural frequency corresponding to the first 6-order array of civil houses(unit:Hz)

圖25 民房前6階陣型圖Fig. 25 The first six formations of civilian houses

從表1中可看出,受保護房屋的第六階陣型對應(yīng)的自振頻率為27.08 Hz,本次爆破的主頻為26.83 Hz;被保護房屋的第五階陣型對應(yīng)的自振頻率為20.17 Hz,本次爆破的次主頻為19.32 Hz。根據(jù)結(jié)構(gòu)共振的知識,不難發(fā)現(xiàn)此次爆破產(chǎn)生的地震波可能會引起該房屋的共振。

進一步定量分析不同IMF分量對民房每一階陣型的影響程度,IMF對不同的建筑結(jié)構(gòu)具有不同的放大倍數(shù),根據(jù)式(1),可計算IMF分量對民房每一階陣型的放大倍數(shù),計算結(jié)果見表2,β為單個IMF主頻和房屋各階陣型固有頻率之比,λ為阻尼比,一般建筑都λ取0.05?？紤]到本次爆破地震波能量主要在IMF3中,其次在IMF4及IMF5中,所以僅對IMF3～IMF5對房屋不同陣型的放大性進行討論。

表2 IMF分量對不同陣型下民房的放大系數(shù)DTable 2 The IMF component corresponds to the magnification factor D under different formations of civil houses

(1)

分析表2可發(fā)現(xiàn),房屋固有頻率和信號主頻越接近即單個IMF主頻和房屋各階陣型固有頻率之比β趨近于1時,D越大。D值大小反應(yīng)在一定的外界爆破信號激勵下,結(jié)構(gòu)對應(yīng)質(zhì)點產(chǎn)生的質(zhì)點振動速度放大效應(yīng),D值越大對應(yīng)的爆破作用放大效應(yīng)越強,對結(jié)構(gòu)的損傷越大,結(jié)構(gòu)發(fā)生破壞的概率越大。

不難發(fā)現(xiàn)IMF3可使民房第6階陣型房屋質(zhì)點速度產(chǎn)生9.9倍的放大效應(yīng)。當(dāng)IMF主頻和陣型自振頻率相等時,即β=1,此時D=10,結(jié)構(gòu)振動幅度達到最大,這便解釋了為什么民房處測得的速度峰值僅0.49 cm/s,但在民房處依舊產(chǎn)生了大量的裂縫,IMF分量對民房每一階陣型的放大倍數(shù)是不一樣的,實際結(jié)構(gòu)在共振的作用下振幅呈現(xiàn)出放大的特征,可據(jù)此解釋房屋在監(jiān)測低振速環(huán)境下發(fā)生開裂的原因。因此,實際工作中安全評估不能以單一速度峰值作為判別標(biāo)準(zhǔn)。

繼續(xù)觀察表2還可發(fā)現(xiàn),對于民房結(jié)構(gòu)而言由于結(jié)構(gòu)前六階陣型相對較小,而高頻信號對此種結(jié)構(gòu)的影響相對50 Hz以下的低頻影響相對小。針對此現(xiàn)象建議施工采用降低單段藥量法、微差起爆、優(yōu)化裝藥結(jié)構(gòu)等提高爆破地震波信號頻率的方法。

綜上分析可發(fā)現(xiàn),基于CEEMDAN·MPE-NHT的爆破地震波信號時頻分析方法,不僅有助于抑制含噪監(jiān)測信號引起的EMD模態(tài)混淆,同時CEEMDAN·MPE得到的IMF經(jīng)NHT處理后,得到的時頻譜在分辨率和精度上都得到了提升,分析結(jié)果有助于爆破振動特征的識別和爆破振動危害控制。

4 結(jié)論

(1)為提高HHT含噪爆破地震波信號時頻分析精度,得到反映真實爆破振動屬性的時頻能量參數(shù)。對HHT進行了改進,得到CEEMDAN·MPE-NHT時頻分析算法。該算法通過控制噪聲來改善EMD模態(tài)混淆,對IMF調(diào)頻分量進行Hilbert變換,降低負值瞬時頻率出現(xiàn)的概率,從算法原理上實現(xiàn)提升信號時頻分析精度。

(2) 通過HHT和CEEMDAN·MPE-NHT算法仿真爆破振動信號時頻分析對比研究可發(fā)現(xiàn),CEEMDAN·MPE-NHT算法將CEEMDAN的自適應(yīng)性和MPE隨機性檢測能力相結(jié)合,降低EMD對噪聲的敏感性。不僅可有效抑制噪聲引起的EMD模態(tài)混淆,同時結(jié)合NHT可得到時頻分辨率雙高的信號頻譜圖。

(3) 根據(jù)結(jié)構(gòu)共振放大系數(shù)D的原理,發(fā)現(xiàn)單個IMF分量主頻和房屋不同陣型下固有頻率越接近時,對應(yīng)質(zhì)點振動速度放大效應(yīng)越明顯。當(dāng)出現(xiàn)極限情況,即當(dāng)β=1,對應(yīng)質(zhì)點振動速度增幅可高達10倍?？蓳?jù)此解釋房屋在低振速情況下,出現(xiàn)大量的爆破振動裂縫的情形。