基于EMD分解的爆破振動信號主成分提取與分析

張聲輝，高文學，劉江超，李小帥，胡 宇，張小軍

(北京工業大學城市建設學部，北京 100124)

由于爆炸能量的瞬間釋放，產生的爆炸應力波迅速沿著介質傳播，致使很大范圍內的周邊建(構)筑物都會受到爆破振動的影響，所以其帶來的隱患一直以來皆受到學者們的重點關注[1-2]。爆破振動監測可作為施工方案和參數設計優化的衡量手段，對工程的經濟安全生產有著重要的指導意義和反饋價值[3]。

爆破振動信號受到復雜環境條件和監測儀器自身的限制，具有瞬時非平穩性的特點，基本都會包含噪聲和基線偏移等雜亂的干擾信號成分，而它們為信號的關鍵分析帶來較多困難[4-7]?，F如今，各領域信號的處理方法和分析研究也漸為深入。由剛開始的傅立葉變換(FFT)，隨繼出現短時傅里葉變換(STFT)、小波變換(WT)、小波包變換(WPT)、希爾伯特-黃變換(HHT)等技術，這些方法皆被應用到爆破振動信號領域[8-15]。邱賢陽等[9]針對金銅礦邊坡爆破獲取了單段振速信號，利用EMD分解和Hilbert能量譜得到并探討了不同延時下兩段疊加信號的振動特性，分析了短延時爆破在不同條件下的信號疊加降振；曹曉立等[10]基于Hilbert-Huang 變換分析很好地反映了魯坨路路塹邊坡爆破信號的時頻特征及能量分布情況；張其虎等[11]利用HHT法與EMD分解研究了不同藥量對金沙礦業官房礦段地下淺孔爆破振動多參量的影響；錢守一等[12]和龔敏等[13-14]也將HHT瞬時能量對延時爆破的振動信號進行了精確識別與分析。上述學者也通過研究證明EMD分解相較于小波(包)分析有無需選定基函數的優點，因此，近年來通過HHT分析爆破振動信號較為廣泛。

本文結合希爾伯特-黃(HHT)分析方法自適應性好的優點，基于崇禮國家跳臺滑雪中心邊坡巖體爆破現場振動監測采集的信號，通過EMD分解算法去除高頻噪聲和低頻趨勢項來實現爆破振動信號主成分的重點提取，并對爆破振動信號主成分展開了分析研究。

1 爆破振動監測試驗

1.1 工程概況

崇禮國家跳臺滑雪中心是2022年冬奧會主要比賽建設場地，距離城區20 km左右。施工建設區域的主要巖性為不同分化程度的花崗巖，根據局部山體開挖設計要求，崇禮跳臺滑雪中心北側邊坡由上而下逐級采用臺階爆破開挖。臺階高度為10 m，邊坡坡率為1∶0.9和1∶0.5，高程1 692 m的臺階寬度為4 m，其余均為2 m；爆破開挖的邊坡設計如圖1所示。

圖1 爆破開挖邊坡橫斷面Fig.1 Cross section of slope excavated by blasting

此次爆破開挖采用梅花形布孔方式和2#巖石乳化炸藥，炸藥單耗q=0.4 kg/m3，密度1.1 g/cm3，炮孔直徑90 mm，孔深10 m，超深1 m，填塞3.0 m左右，孔距3.0 m，排距2.5 m。預裂孔間距1.0 m，線裝藥密度0.5 kg/m。采用逐排起爆網路，最大段藥量550 kg。

1.2 爆破振動監測信號

爆破振動信號用TC-4850N型測振儀進行采集，測點布置在爆區上部邊坡高程1 702 m臺階平臺上，且距離坡頂1 m的位置處，位置確定后，將傳感器用石膏固定粘在清理干凈的基巖上，以保證監測效果的良好性。采集到的爆破振動信號如圖2所示。

圖2 爆破振動監測信號Fig.2 Blasting vibration monitoring signal

2 HHT變換原理

2.1 經驗模態分解算法與去噪

HHT分析方法是由Norden E Huang等[16]提出，該方法可針對非平穩隨機信號實現對稱與平穩化的技術處理，在信號分析領域取得了關鍵突破。假設有一爆破振動信號，可利用HHT法將此信號通過EMD分解，即可得到若干表征信號的固有模態函數(IMF)，而其中所有IMF都是單分量幅值或者頻率調制信號，再對IMF分量通過Hilbert變換處理后獲取信號時頻和能量分布特征，進而也可局部分析爆破信號在時間與頻率中的不同能量分布規律。

為研究如爆破開挖產生的瞬態且非平穩振動信號，頻率應該假定為與時間相關的函數，對監測的爆破振動信號進行EMD分解的要求是一步步篩選而獲得IMF。IMF是要符合以下假設[17]:①信號數據序列中極值點與零點個數的差值為0或±1；②信號任一點通過局部極大值得到的上包絡線與極小值得到的下包絡線的均值是0。

EMD分解算法[17]是對監測振動信號X(t)求解其一切極值點，利用三次樣條插值函數曲線對極大、小值點各自求插值，擬合得到爆破振動信號X(t)的上、下包絡線Xmax(t)、Xmin(t)，再求它們均值線m1(t)為

m1(t)=[Xmax(t)+Xmin(t)]/2

(1)

h1(t)=X(t)-m1(t)

(2)

信號的不一樣也許使h1(t)是一個IMF分量，但如若不是，此時應將h1(t)當成原信號重復上述步驟k-1次，則有

hk(t)=hk-1(t)-mk(t)

(3)

根據IMF假設和篩選停止的準則[17]，直到可把hk(t)認為是第一階IMF，記為c1(t)，殘差r1(t)=X(t)-c1(t)，并記為第二次篩選的爆破振動信號，重復n次運算，當cn(t)或rn(t)小于設定誤差閾值，或殘差rn(t)變成單調函數時，此時分解便終止不再繼續。至此，原信號X(t)可由n階IMF分量和殘差rn(t)組成，也就是

(4)

信號經EMD分解得到由各頻段構成的信號ci(t)，去除噪聲高頻信號分量，可組合不同頻段的IMF分量實現提取我們所需要的爆破振動主成分信號。

2.2 Hilbert變換與譜

對上述EMD分解的爆破振動信號IMF分量的Hilbert變換：

(5)

式中：PV為柯西主值，構造如下解析信號為

z(t)=c(t)+jH[c(t)]=a(t)ejΦ(t)

(6)

在相位函數的基礎上定義瞬時頻率為

(7)

IMF分量經Hilbert變換之后，爆破振動信號Hilbert譜可寫成：

(8)

式(8)省去了殘余函數r，Re為取實部。若H(ω,t)對時間積分，即可得到Hilbert邊際譜：

(9)

邊際譜表示了爆破振動信號各頻率在全局上累加的幅值或能量，同時，可定義Hilbert瞬時能量為IE(t)；實際上，若振幅的平方對時間進行積分，則能獲取Hilbert能量譜，再在能量譜基礎上定義Hilbert邊際能量為E(ω)，兩者表達式如下：

(10)

(11)

這里，瞬時能量和邊際能量分別提供了信號能量隨時間和頻率的變化情況。

3 爆破振動信號主成分提取結果分析

3.1 爆破振動信號主成分提取

基于現場爆破振動監測采集的信號(見圖2)，采用上述EMD分解算法編制MATLAB程序進行運算，所得的各階IMF分量呈現出的模態(見圖3)，且具有明晰的物理意義，并計算了不同IMF分量能量分布特征如圖4所示。

圖3 EMD分解后的各階IMF分量Fig.3 IMF components of each order after EMD decomposition

圖4 IMF各階分量能量分布Fig.4 Energy distribution of each order IMF component

由圖3和4可知，監測的爆破振動信號經EMD分解得到c1～c9及余量共計10個IMF分量。一般情況下，這些IMF各階分量是按局部頻帶區間由高至低依次排列，研究結果充分說明了EMD分解具有可使信號反映出局部性強的特點；c1分量作為高頻噪聲，出現在整個信號采集的過程中，其波形震蕩密集并表現出高頻、低幅值、波長短等特點，是原始信號在復雜環境監測條件下受到噪聲干擾的成分，而殘余項曲線呈現了振動監測過程中產生的零點漂移與信號偏離基線的趨勢，歸為振動信號成分中的低頻趨勢項。c2～c6分量的振速幅值比較大，涵蓋了信號的絕大部分能量，該些分量會對周圍邊坡巖體與建(構)筑物造成一定強度的擾動，不可與其他分量等同視之，而應加大關注；c7～c9分量的爆破振動頻率衰減較快，波長也隨之拉大，直至降至成微弱的振速和能量。

3.2 爆破振動信號主成分分析

為對爆破振動信號主成分提取后的結果展開研究，在上述對信號經EMD分解得到的IMF分量的前提下，對c2～c9的IMF分量進行重構，得到去噪和消除殘余項的主成分信號(見圖5)，再通過Hilbert變換能求解爆破信號能量與時間或頻率之間的對應關系，得到Hilbert能量譜(見圖6)、瞬時能量(見圖7)和邊際能量(見圖8)。

圖5 爆破振動監測信號主成分提取后的對比情況Fig.5 Comparison of blasting vibration monitoring signal after extracting principal components

圖6 Hilbert能量譜Fig.6 Hilbert energy spectrum

圖7 瞬時能量Fig.7 Instantaneous energy

由圖5可知，在爆破振動監測信號尾振段有很明顯的高頻振蕩噪聲干擾成分，信號主成分提取后的峰值振速變化很小，有效地保留了信號的形態與且曲線光滑性提高，爆破振動信號主成分提取效果比較理想。

由圖6可知，巖質邊坡爆破振動信號的頻率較為豐富，其Hilbert能量譜中偏紅色點對應的頻率主要集中在60 Hz以下，說明高頻分成能量較少；爆破振動信號主成分提取前后的Hilbert能量譜在頻率為10 Hz以下的顏色有明顯的改變，研究發現提取了爆破振動信號主成分的能量降低了，可以表明是消除了低頻趨勢項的結果。

由圖7可以看出，瞬時能量隨監測時間的變化情況，爆破振動信號主成分提取前后的瞬時能量峰值均出現在0.155 s左右，分別為112.12、108.70 cm2·Hz/s2，誤差為3%左右，相對較小，研究結果與原始爆破振動信號峰值振速分布規律基本一致；同時信號的瞬時能量主要集中在0.1～0.3 s，結果和圖2的監測波形或圖6的Hilbert能量分布情況一樣，有良好的對應效果；爆破振動信號主成分提取前的瞬時能量在I，II，III，IV處的光滑性較差，經過去噪和消除趨勢項的處理之后，即信號主成分提取后的瞬時能量在I，II，III，IV處的光滑性有一定的改善，其曲線整體過渡性光滑度有所提高，信號主成分的提取效果良好。

由圖8可知，信號主成分提取前后邊際能量在頻段內有一些波動，但它們都主要集中在60 Hz以下，爆破振動信號主成分提取后的頻率分布稍有變窄，且在10 Hz以下的邊際能量顯著降低，與Hilbert能量譜的分析結果一致。

圖8 邊際能量Fig.8 Marginal energy

因此，在分析和評估邊坡巖體爆破開挖產生的振速和能量時，我們應充分考慮并對信號高頻噪聲與低頻趨勢項進行剔除，獲取爆破振動信號主成分，可提高信號的局部時頻特征，以準確評估建(構)筑物的爆破振動，并為保證爆破施工安全與邊坡穩定提供參考依據。

4 結論

1)EMD分解的各階IMF分量在局部頻帶區間從高到低依次排列，呈現出一定尺度范圍的模態，具有明晰的物理意義，去除高頻噪聲和低頻趨勢項可以實現爆破振動信號主成分的重點提取。

2)Hilbert能量譜能直觀體現信號能量與時頻之間的關系，信號能量主要在60 Hz以下的低頻段和0.1～0.3 s的時間段，爆破振動信號主成分提取前后的瞬時能量峰值分別為112.12、108.70 cm2·Hz/s2，相對誤差較小。

3)提取主成分后的爆破振動信號頻率分布稍有變窄，且在10 Hz以下的邊際能量顯著降低，而瞬時能量和振動速度曲線整體過渡性光滑度有所提高，信號主成分提取效果良好。因此，在分析邊坡巖體爆破開挖產生的振速和能量時，獲取爆破振動信號主成分，可提高信號的局部時頻特征，以準確評估建(構)筑物的爆破振動，并為保證爆破施工安全與邊坡穩定提供參考依據。