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

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

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

希爾伯特-黃變換(Hilbert-Huang Transform,HHT)由經驗模態分解(Ensemble Empirical Mode,EMD)[1]和希爾伯特變換(Hilbert Transform,HT)[2]組成,其中EMD以其能識別數據內在屬性并根據數據本身特征進行分解在信號分析領域得到了廣泛應用[3]。對EMD分解結果進行Hilbert變換,更是實現了將時域信號轉變為頻域信號,幫助爆破工程技術人員對爆破地震波信號傳播特征、內在屬性進行識別,對爆破地震波危害控制起到了一定的指導作用[4-6]。

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

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

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

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

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

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

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

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

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

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

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

為突出CEEMDAN·MPE-NHT時頻分析算法相比HHT可有效提高含噪爆破地震波信號時頻特征參數的提取精度,特進行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所示。采樣點數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模態混淆的抑制作用,分別采用EMD和CEEMDAN·MPE對S(t)進行模態分解。圖6為EMD分解結果,不難發現IMF1是難以除去的噪聲信號;高頻模態混淆嚴重,如IMF2在0.53 s,0.64 s及0.85 s附近均出現了高頻向中頻發散的趨勢;中頻模態混淆有所緩解,如IMF3在0.41 s及右端點附近存在向低頻發展的趨勢;低頻相對穩定,如IMF5和IMF6。圖7為CEEMDAN·MPE分解結果,分解得到的IMF從高頻向低頻依次排列,未見明顯模態混淆現象,模態顯示相比EMD結果清晰且穩定。

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

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

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

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

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

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

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

進一步分析,對圖11中的IMF1～IMF6執行NHT,獲得每個IMF的時頻能量特征圖,其中圖12～圖17為對應IMF的邊際譜;圖18～圖23為對應IMF的時頻譜;這里未對IMF7和R進行變換,是因為結合圖11不難發現IMF7和R占有的能量可忽略不計。觀察圖12～圖23可發現,CEEMDAN·MPE-NHT獲得的時頻譜在時域和頻域都具有高分辨率,這點和仿真信號時頻分析得到的結果一致。IMF1的頻率范圍為100～250 Hz,持續時間為0.24～0.38 s;IMF2的頻率范圍為50～100 Hz,持續時間為0.23～0.42 s;IMF3的頻率范圍為20～50 Hz,持續時間為0.22～0.47 s;IMF4的頻率范圍為10～30 Hz,持續時間為0.28～1.17 s;IMF5的頻率范圍為5～10 Hz,持續時間為0.00～1.20 s;IMF6頻率最低,約為0～5 Hz,持續時間為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的分析,可發現信號在低頻停留時間大于高頻,主要能量集中在50 Hz以下。

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

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

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

表1 民房前6階陣型對應的自振頻率(單位: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中可看出,受保護房屋的第六階陣型對應的自振頻率為27.08 Hz,本次爆破的主頻為26.83 Hz;被保護房屋的第五階陣型對應的自振頻率為20.17 Hz,本次爆破的次主頻為19.32 Hz。根據結構共振的知識,不難發現此次爆破產生的地震波可能會引起該房屋的共振。

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

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

(1)

分析表2可發現,房屋固有頻率和信號主頻越接近即單個IMF主頻和房屋各階陣型固有頻率之比β趨近于1時,D越大。D值大小反應在一定的外界爆破信號激勵下,結構對應質點產生的質點振動速度放大效應,D值越大對應的爆破作用放大效應越強,對結構的損傷越大,結構發生破壞的概率越大。

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

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

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

4 結論

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

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

(3) 根據結構共振放大系數D的原理,發現單個IMF分量主頻和房屋不同陣型下固有頻率越接近時,對應質點振動速度放大效應越明顯。當出現極限情況,即當β=1,對應質點振動速度增幅可高達10倍?？蓳私忉尫课菰诘驼袼偾闆r下,出現大量的爆破振動裂縫的情形。