基于HHT方法的爆炸應變波時頻分析

楊仁樹，高祥濤，車玉龍，解北京

(1.中國礦業(yè)大學(北京)力學與建筑工程學院，北京 100083；2中國礦業(yè)大學(北京)深部巖土力學與地下工程國家重點實驗室，北京 100083；3中國礦業(yè)大學(北京)資源與安全工程學院，北京 100083)

固體介質中爆炸應變波信號為瞬變的非周期、非平穩(wěn)隨機超動態(tài)信號。在爆源近中區(qū)，固體介質中爆炸應變波動實質為介質內質點振動總體反映。該質點振動傳播的為體波，不同于通常的爆破地震波[1]，但其信號特征均屬非平穩(wěn)信號范疇。對固體介質中爆炸波頻譜分析早期主要用傅里葉變換[2-4]。因理論發(fā)展的限制，F(xiàn)ourier分析方法基于線性及平穩(wěn)假定。隨信號分析理論與技術的發(fā)展，數(shù)據(jù)信號分析出現(xiàn)諸多新方法。如短時快速傅里葉變換、Cohen類、小波分析及Wagner-Ville分布等均能用于線性非平穩(wěn)數(shù)據(jù)[5-7]分析，但其本質均未擺脫傅里葉分析的局限，存在基函數(shù)選擇及頻譜擴散問題[8]。此外，各種非線性時序分析方法被設計成非線性平穩(wěn)確定性系統(tǒng)，但實際所得大部分信號均具有非線性非平穩(wěn)特征，分析較困難。

HHT方法廣泛用于地震工程[9]、機械振動信號分析[10-11]、結構損傷分析[12]等領域，能較好處理工程中非線性非平穩(wěn)信號。其中對爆破振動信號HHT法的應用研究[1]取得較好效果。但用HHT方法分析研究固體介質近中區(qū)爆炸應變波信號未見報道。本文據(jù)爆炸波信號短時、瞬變、非周期、非線性、非平穩(wěn)特點，用HHT方法分析爆炸應變波信號時頻特征，并與STFT方法對比以證明HHT方法的優(yōu)越性。為研究固體介質中爆炸應力波傳播機理、工程結構的安全防護提供可靠依據(jù)。

1 HHT分析法[13]

1.1 EMD原理與算法

HHT信號分析法為全新分析技術，由EMD方法與Hilbert變換兩部分組成，核心為EMD。其據(jù)信號本身的時間尺度特性，將信號分解為含有不同時間尺度且滿足一定條件的一組IMF(Intrinsic Mode Function)分量。EMD算法為：① 找出信號x(t)上所有極值點，用三次樣條插值函數(shù)對所有極大值點進行插值，擬合出原始信號x(t)上包絡線xmax(t)。同理可得下包絡線xmin(t)。上、下兩條包絡線含所有信號數(shù)據(jù)。按順序連接上、下包絡線均值即得一條均值線m1(t)，將原數(shù)據(jù)序列減去m1(t)可得一去掉低頻的新數(shù)據(jù)序列h1(t)，用篩選過程終止準則判別不同信號h1k(t)是否IMF分量。② 符合終止準則條件的h1k(t)即為可表示信號數(shù)據(jù)中最高頻成分的第一IMFc1(t)，再用x(t)減去c1(t)獲得一新數(shù)據(jù)序列r1(t)。重復步驟①，獲得一系列cn(t)及最后一個不可分解序列rn(t)，rn(t)代表x(t)的均值或趨勢項，此時可終止模態(tài)分解過程。原始信號x(t)可由n階IMF分量及殘差rn(t)構成：

(1)

1.2 Hilbert變換與Hilbert譜

將EMD分解的每個IMF分量進行Hilbert變換，即可得每個IMF分量的瞬時頻譜，綜合所有IMF分量瞬時頻譜可得Hilbert譜，每個IMF分量進行Hilbert變換后，將x(t)表示成Hilbert譜形式：

(2)

將式(2)三維空間中信號幅度表達為時間與瞬時頻率的函數(shù)，即Hilbert譜。若振幅的平方對時間積分，可得Hilbert能量譜：

(3)

Hilbert能量譜可提供每個頻率的能量計算式，表達每個頻率在整個時間長度內累積的能量。對實驗所測爆炸應變信號進行Matlab程序語言設計，可得HHT分析結果。

2 巖石中爆炸應變測試實驗

實驗中采用在水泥砂漿試件中預埋應變磚方法測試爆源近中區(qū)爆炸應變信號。水泥砂漿試塊為圓柱體，尺寸Φ400 mm×450 mm；應變磚尺寸180 mm× 20 mm×10 mm。圓柱體試件與應變磚均由PC32.5普通硅酸鹽水泥、黃砂(細骨料)制成，水泥砂漿配比(質量比)為水泥：黃砂：水=1∶2∶0.5。試件養(yǎng)護28 d后測試其物理力學參數(shù)為密度2.3 g/cm3，抗壓強度56.5 MPa，彈性模量24.87 GPa，泊松比0.25，縱波速度4 445 m/s，橫波速度2 080 m/s。每塊模型預埋1個應變磚，埋深225 mm，理論上認為固體介質中距爆源8～150倍裝藥半徑范圍為衰減后應力波及爆生氣體膨脹作用的破裂區(qū)，該區(qū)在爆破過程中占重要地位。因此在該區(qū)域內應變磚的測點布置為距炮孔4 cm、8 cm、13 cm、18 cm，每個距離處分別設置徑向、橫向應變片。應變片為120 Ω箔式環(huán)氧酚醛類電阻應變片，敏感柵尺寸3×2(mm)，靈敏系數(shù)2.08%。實驗藥包為圓柱形，裝藥為二硝基重氮酚(C6H2(NO2)2N2O)，裝藥密度ρ0=0.82 g/cm3，爆速3 991 m/s，采用橋絲貫穿裝藥軸心瞬時起爆。炮孔用石英砂與AB膠混合后堵塞密實。對巖土中爆炸應力波測試，選動態(tài)信號數(shù)據(jù)采集儀采樣頻率為40 MHz，前置放大器頻率響應2.5 MHz。測試裝置示意圖見圖1。

圖1 測試裝置示意圖

3 計算結果及分析

3.1 實驗結果

相同藥量與起爆方式重復三組混凝土介質爆炸近中區(qū)測試實驗，獲得距爆源不同距離處爆炸應變波信號，限于篇幅僅以其中一組實驗中同一測點處兩個應變波形為例見圖2，其它測點波形相似。圖2中1#、2#爆炸應變波信號分別為距爆源中心4 cm處應變片所測徑向應變、橫向應變。由圖2可見，爆炸應變波在固體近中距離處同一測點同時產生徑向、橫向應變，徑向主要承受壓應變，橫向主要承受拉應變，徑向應變絕對值(8×104μs)略大于橫向應變絕對值(7.5×104μs)。

圖2 所測爆炸應變信號

3.2 爆炸波信號時頻特征分析

用HHT方法定量分析爆炸應變波信號的時頻分布特征。對圖2中1#、2#爆炸應變波信號分別進行EMD分解，并對分解的IMF分量進行快速傅里葉變換(FFT)，分析各分量頻域特征，結果見圖3。由圖3可見，EMD分解順序按頻率由高到低進行。圖3(a)、(c)中1#、2#爆炸應變波信號經EMD分解成7個IMF分量(從高頻到低頻)及1個殘余分量R(采集信號的整體變化趨勢)；圖3(b)中，1#徑向應變信號EMD分解出的7個IMF分量C1～C7頻率范圍主要集中在2.5～17.5 MHz、0～6 MHz、0～3 MHz、0～1.5 MHz、0～1 MHz、0～0.5 MHz、0～0.5 MHz；各IMF分量C1～ C7的平均頻率為7.25 MHz、3.1 MHz、1.14 MHz、0.627 MHz、0.314 MHz、0.196 MHz、0.157 MHz；圖3(d)中2#橫向應變信號EMD分解出的7個IMF分量C1～C7頻率范圍主要集中在2.5～8 MHz、1～5 MHz、0～3 MHz、0～1.5 MHz、0～1 MHz、0～0.5 MHz、0～0.5 MHz；各IMF分量C1～C7平均頻率為6.98 MHz、3.45 MHz、1.33 MHz、0.392 MHz、0.314 MHz、0.186 MHz、0.118 MHz。其它4組測點應變信號EMD分析結果基本相似，各測點信號IMF分量頻帶及相應能量分布見表1。由表1看出，或徑向應變或橫向應變，距爆源越近高頻IMF分量頻帶所占能量比例越大。如1#徑向IMF_C1分量(7.25 MHz)，其能量占總能量的53.6%；隨爆炸應力波向遠處傳播，其所占能量比例大的IMF分量頻帶逐漸向低頻率轉移，如7#徑向IMF_C5分量(0.275 MHz)其能量占總能量的56.2%。可見，炸藥在巖石介質中爆炸的頻率分布豐富，不同頻率波動對介質質點做功不同。每個固有模態(tài)分量(IMF)代表1個頻率段的作用結果。而爆源近中區(qū)爆炸應變波的主頻率在幾百千赫茲到幾兆赫茲范圍內。

圖3 爆炸應變波信號EMD分解及頻域圖

表1 各IMF分量所占原始信號能量百分比及頻率分布

3.3 爆炸波信號HHT法與STFT對比分析

HHT方法與STFT方法與非平穩(wěn)信號的優(yōu)越性對比已較清楚，但對瞬變非周期非平穩(wěn)隨機超動態(tài)爆炸應變波信號傳播細節(jié)特征，利用HHT法與STFT法對比研究巖石中爆炸應變波時頻特征研究未見報道。將實測爆炸應變波信號同時用HHT法、STFT法的分析結果見圖4。由圖4可見，與短時傅里葉變換時頻能量譜相比，希爾伯特能量譜能清晰表達能量隨時間頻率的分布細節(jié)，信號能量主要集中于小能量譜段，能準確刻畫信號特征，而短時傅里葉變換的能量譜范圍較寬，且頻率分辨率較低，未體現(xiàn)信號本身的非平穩(wěn)特性，只具有相對分布意義。此因為爆炸波傳播具有多頻率作用特性，短時傅里葉變換(STFT)所建頻譜圖為最簡單、最直觀的時頻分布。但在分析非平穩(wěn)信號時，受Heisenberg不確定性原理限制，時頻分辨率不能自適應改變，時間分辨率及頻率分辨率需要折中，且STFT使用固定的短時窗函數(shù)，較能降低頻率分辨率，不能顯示爆炸沖擊信號的能量細節(jié)。而爆炸應變信號的Hilbert能量譜可清晰反映隨時間頻率變化的具體分布，利用HHT方法可明顯表現(xiàn)出爆炸應變波傳播細節(jié)，從而印證爆炸應變波的多頻特征。

圖4 爆炸應變波信號HHT變換與STFT變換比較

3.4 爆炸波信號多頻率作用特性探討

為探討爆炸波信號多頻率作用特性，分別給出空氣、巖石介質中爆炸波的傳播數(shù)值模擬結果及有機玻璃中不耦合裝藥爆炸應力波傳播的動光彈實驗結果。圖5為爆炸波在不同介質中的傳播特征。圖5(a)為運用Euler-FCT格式模擬爆炸波傳播。由圖5(a)看出，爆炸波由不同層級壓力作用于周圍空氣介質，產生不同幅度波動。由于爆炸波動過程的不定常性，爆炸波與周圍介質作用在極短時間內往復振蕩[14]。此振蕩隨時間、距離衰減，介質內部受沖擊壓縮時出現(xiàn)電動勢、電磁輻射，必具有多頻率相互作用特點。圖5(b)為巖石中耦合裝藥爆炸應力波傳播，巖土介質運用分段的Druck-Prager本構模型，模擬結果清晰表明，爆炸波動能量隨時間、距離耗散，加之巖石介質的不均勻性，爆炸波會發(fā)生彌散現(xiàn)象，波動頻率更豐富。圖5(c)為有機玻璃中不耦合裝藥爆炸應力波動的多頻段特性[15]，實驗結果清晰表明，有機玻璃圍繞爆源周圍產生諸多疏密相間的爆炸應力波動條紋，其本質為因爆炸應力波動多頻段特性作用結果。

因此，固體介質中爆炸應力波信號由各種不同頻率波動共同作用，該復合疊加的波動作用于介質，若超過介質的固有頻率，則介質將發(fā)生拉伸或剪切破壞?？捎脭?shù)值模擬及文獻中數(shù)據(jù)佐證爆炸波動的多頻率特性。

圖5 不同介質中炸藥爆炸波傳播

4 結 論

爆炸應變波信號為典型的瞬變非平穩(wěn)隨機信號，且爆炸應力波頻譜構成復雜、影響因素眾多。HHT變換能較好用于具有非線性、非平穩(wěn)特性的爆炸應變波信號分析中，且能提取時程曲線的主要特征信息，結論如下：

(1) 通過爆炸應變波信號的EMD分解及時頻分析知，爆源近中區(qū)爆炸應變波的主頻率在幾百千赫茲到幾兆赫茲范圍內。

(2) 通過爆炸波信號HHT法與STFT時頻能對比分析知，Hilbert能量譜能較好、清晰描述爆炸應變波信號能量在時頻譜的分布，可獲得爆炸應變波信號的頻譜及功率譜隨時間變化的細節(jié)特征。

(3) 通過分析爆炸波在空氣、巖石介質中傳播數(shù)值模擬、爆炸應變波在有機玻璃介質傳播動光彈實驗結果知，爆炸波信號具有多頻率作用特性，即介質中爆炸波信號波形為由不同頻率波動共同作用結果。

[1]張義平，李夕兵，左宇軍.爆破震動信號的HHT分析與應用[M]. 北京：冶金工業(yè)出版社，2008.

[2]宋熙太，任建民.爆炸波分析及其波譜[J].爆炸與沖擊，1982,4:33-42.

SONG Xi-tai, REN Jian-min. Analysis of explosion waves and wave spectrum[J]. Explosion and Shock Waves, 1982, 4: 33-42.

[3]李廷芥，洪憫萱.瞬態(tài)波形的頻譜分析及在巖石動力學中的應用[J]. 爆炸與沖擊, 1985,5(4):10-16.

LI Ting-jie, HONG Min-xuan. Spectral analysis of transient waveforms and its applications in rock dynamics[J]. Explosion and Shock Waves, 1985,5(4):10-16.

[4]洪憫萱,李廷芥.爆炸波的頻譜分析[J].巖土力學,1987,8(1): 51-57.

HONG Min-xuan, LI Ting-jie. Spectral analysis of blast waveforms[J]. Rock and Soil Mechanics, 1987,8(1)：51-57.

[5]趙建平, 陳壽如.基于短時Fourier變換的爆炸波時頻能分布及分離[J]. 采礦技術,2011, 11(5): 84-86.

ZHAO Jian-ping, CHEN Shou-ru. Blast wave time-frequency energy distribution based on STFT[J]. Mining Technology, 2011, 11(5): 84-86.

[6]毛婧，韓焱.小波分析在爆炸沖擊波信號處理中的應用[J].彈箭與制導學報,2007,27(4):137-144.

MAO Jing, HAN Yan. Wavelet analysis in explosive shock wave signal processing[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2007,27(4):137-144.

[7]趙建平,林杭.基于時頻能分析技術的爆炸波能量分布及分離[J].巖石力學與工程學報,2012, 31(S1): 3278-3285.

ZHAO Jian-ping, LIN Hang. Energy distribution and separation of blast wave based on time-frequency energy analysis technology[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(S1): 3278-3285.

[8]趙楊，武岳，曹曙陽，等. 利用HHT方法對非平穩(wěn)風力的時頻分析[J]. 振動與沖擊, 2011, 30(2): 5-14.

ZHAO Yang, WU Yue, CAO Shu-yang,et al. Time-frequency analysis of a non-stationary wind pressure with HHT method[J]. Journal of Vibration and Shock, 2011, 30(2): 5-14.

[9]胡燦陽，陳清軍. 基于HHT法的地震地面運動局部譜密度估計[J]. 振動與沖擊，2007,26(10):126-131.

HU Can-yang, CHEN Qing-jun. Estimation of local spectral density of earthquake ground motion based on HHT theory[J]. Journal of Vibration and Shock, 2007,26(10):126-131.

[10]張帥，楊勇，韓清凱，等. 基于HHT的轉子系統(tǒng)定點碰摩實驗研究[J]. 振動與沖擊，2010,29(7)：121-125.

ZHANG Shuai, YANG Yong, HAN Qing-kai,et al. Experimental study on a rotor system with rub-impact at fixed limiter based on HHT[J]. Journal of Vibration and Shock, 2010,29(7)：121-125.

[11]唐先廣，郭瑜，丁彥春. 基于獨立分量分析與希爾伯特-黃變換的軸承故障特征提取[J]. 振動與沖擊,2011, 30(10):45-49.

TANG Xian-guang, GUO Yu, DING Yan-chun. Rolling element bearing fault feature extraction based on HHT and independent compoment analysis[J]. Journal of Vibration and Shock, 2011,30(10):45-49.

[12]黃天立.結構系統(tǒng)和損傷識別的若干方法研究[D].上海：同濟大學，2007.

[13]Huang N E, Shen S S P. The Hilbert-Huang transform and its applications[M]. World Scientific Publishing Co. Pte. Ltd., Toh Tuck Link, Singapore,2005.

[14]奧爾連科JI II,著. 孫承緯,譯.爆炸物理學[M].北京：科學出版社,2011.

[15]曹躍，趙泉，孔安.不耦合裝藥爆炸擴腔機理研究[J].爆破,2005,22(3):21-25.

CAO Yue, ZHAO Quan, KONG An. The springing mechanism of spherical charges in an eccentric decouple pattern[J]. Blasting, 2005,22(3):21-25.