HHT在爆炸振動信號處理中的應用

蘇秀紅

（中國工程物理研究院總體工程研究所，四川綿陽621900）

HHT在爆炸振動信號處理中的應用

蘇秀紅

（中國工程物理研究院總體工程研究所，四川綿陽621900）

HHT方法作為一種新興的非平穩信號處理方法，因其高度的自適應性和高效性，近年來受到學者的廣泛研究。在某炸藥件銷毀的過程中，本文借助振動測試系統開展了炸藥件爆炸振動波的測試，并且采用希爾伯特黃變換對振動信號進行分析處理，從數據分析的結果中探究爆炸波的特點，可為大藥量炸藥爆炸的安全風險評估提供參考。

希爾伯特黃變換；EMD分解；振動信號分析；爆炸試驗

非線性非平穩信號的處理近年來受到廣泛研究。環境試驗中涉及的振動信號如正弦信號、隨機信號大多作平穩信號處理，但是也有很多振動信號是非平穩的，例如炸藥爆炸產生的振動信號，沖擊信號，機械設備運行過程中的摩擦、裂紋等故障信號都表現出非平穩性。炸藥爆炸產生的振動波在地表和空氣中傳播，由此采集到的信號屬于非平穩信號。經典的振動信號數字化處理方法包括傅里葉頻域分析法和小波變換。傅里葉分析的對象是信號的頻域，不能夠同時獲得時域和頻域的特征，也無法反映出信號在時間局部區域上的頻率特征，因此只適用于線性系統。而小波分析作為一種非平穩數據分析方法，其本質上是窗口可調的Fourier譜分析方法，一旦小波基底選定，所有的數據都是基于這一組基進行分析，因此小波分析不是自適應的。

相較于上述兩種方法，1998年由N.E.Huang提出一種新的時頻分析方法，用以分析非線性非平穩信號，美國宇航局將其命名為Hi1bert_Huang Transform，簡稱HHT，具體包括經驗模態分解（Empirica1 Mode Decomposition，EMD）和希爾伯特（Hi1bert）譜分析［1_2］。與小波變換不同的是，HHT的時頻分析方法是自適應，具有很好的時頻聚集性和時頻分辨率，它可以根據信號的局部時變特征進行自適應的時頻分解，非常適合于非平穩信號的分析。本文首先求解爆炸振動波的傳播速度，其次將HHT運用到炸藥件爆炸振動信號的分析中，通過EMD分解得到信號中包含的各IMF分量，經Hi1bert變換繪制出信號的頻譜圖、邊際譜圖等，并分析得出爆破振動信號的時頻特性，能量分布等特性，從而為大藥量炸藥爆炸的安全風險評估提供參考。

1　HHT方法的基本原理

HHT方法是一種兩步驟的信號處理方法。首先通過對信號進行經驗模態分解（Empirica1 Mode Decomposition，簡稱EMD），產生一組信號數據序列，每個序列稱為一個固有模態函數（Intrinsic mode function，IMF），不同的IMF具有不同特征尺度。之后對每個IMF分量分別進行Hi1bert變換，得到各個IMF分量的瞬時頻譜，將所有IMF分量的瞬時頻譜進行綜合即可得到信號的時頻譜—Hi1bert譜［3_4］。

1　.1　EMD的基本原理

信號可通過以下步驟實現IMF分解：

對于原始信號，首先找到此信號所有的局部極值點，然后用3次樣條函數曲線對極值點進行插值，擬合出原始信號的上下包絡線。假設上下包絡線的均值為，從原始信號中減去均值，得到，

理想情況下，h1（t）可以滿足IMF的判定條件。若不滿足條件，則將h1（t）看作原始信號，重復上述步驟，直至最終得到的分量c1（t）滿足IMF的判定條件，此時記為，

判定c1（t）滿足IMF的條件是［5］：SD=0.2～0.3，其中，

得到第一個IMF分量后，從原始信號中減去c1（t），得到殘差r1（t），并將其視為新的“原始”信號繼續重復上述步驟，依次可以得到第2、3直至第n個IMF分量，因此有，

當殘余量rn（t）成為一個單調函數，不能從中提取滿足IMF條件的分量時，循環結束。至此，原始信號x（t）可由n階IMF分量以及殘差rn（t）表示：

1.2Hllbert譜與Hllbert邊際譜

將EMD分解得到的IMF分量分別進行Hi1bert變換，即可得到每個分量ci（t）的瞬時頻譜，所有的瞬時頻譜綜合之后即可得到Hi1bert譜［6］。具體的過程如下：首先對IMF分量ci（t）作Hi1bert變換得到解析信號為

解析信號的幅值和相位為：

在上述基礎上定義解析信號的瞬時頻率為

通過Hi1bert變換得到的振幅以及頻率均為時間的函數，將振幅放入時間-頻率中，就可得到Hi1bert譜：

將H（w，t）對時間進行積分，就可得到Hi1bert邊際譜：

邊際譜表達了每個頻率在全局上的幅度（或能量），是統計意義上的全部累計幅度。

2　炸藥爆炸實驗

一發炸藥件由于內部出現缺陷，現決定將其進行爆破銷毀，在銷毀過程中布置相應的傳感器，對相關數據進行測量，研究爆炸振動信號的特性，以了解信號特點，為意外事故中大藥量炸藥爆炸帶來的安全風險評估提供參考。

2.1振動測試

爆炸振動測試示意圖如圖1所示，傳感器布置在一條測線上，共布置3個點，分別在距藥球85 m、135 m、185 m處。為了測得振動波的平均傳播速度，需要保證3個測點的同步性。

2.2測試結果

測試系統記錄的數據包括3個測點X、Y、Z 3個方向的加速度數值。由于工程實驗中，通常采用爆炸波的垂直振動數據（Z方向）作為分析依據［7］，因此本實驗只關心Z方向的試驗數據。測得的3個通道Z方向加速度曲線如圖2。從圖2可以看出，在每次爆炸試驗過程中各測點處均出現了兩次振動過程，第一次振動過程量級較低，第二次量級較高。

圖1　傳感器測點布置圖

圖2　3個測點Z方向加速度曲線

3　基于HHT的爆炸振動波分析

3.1爆炸波傳播速度

由于3個測點處于同一直線上，并且各測點之間距離固定，因此只要測量各測點之間起振的時間差就可以計算得到爆炸波的傳播速度。爆炸各測點起振時間及由此計算得到的波速見表1。

表1　各測點起振時間

從表1中數據可以看出，第一次振動過程的波速約為3 000～4 500 m/s，第二次振動過程的波速約為310～384 m/s。由文獻資料及測試經驗知，通過地面傳輸的地震波，其傳播速度通常在3 000～5 000 m/s，因此可以判定導致第一次振動過程的為地面傳輸的地震波；導致第二次振動過程的波速較為接近聲波在空氣中的傳播速度（約340 m/s）。由于第二次振動過程的量級明顯高于第一次，因此下面的時頻分析中，主要針對第二次振動過程作分析研究。

3.2波形信號EMD分解

選取測點3的第二次振動過程作為分析對象，截取0.2 s的振動過程（0.79～1.09 s）。將波形數據導入Mat1ab中，對信號波形進行EMD分解，得到圖3和圖4所示的IMF分量。

從圖3～圖4可以看出：

1）原始信號被分解為6個時間尺度從大到小排列的IMF分量（imf1～imf6）和一個殘差r7。從上兩圖中也可以容易得出，隨著IMF分量編號的增大，分量的頻率逐漸變低，波形也相應變長，直到最后得到頻率非常低的信號殘差r7，EMD分解結束。這同時說明在爆炸波的傳播過程中，信號的高頻部分衰減幅度較大，而低頻部分則衰減緩慢。

2）與小波分解需要提前設定某個基函數不同，EMD分解的過程是完全自適應的，從而避免了因基函數選取不同而導致結果不同的問題。

3）從各IMF的振幅可看出，能量主要集中在imf1～imf3 這3個分量上。

3.3爆炸信號譜圖

對各IMF分量進行Hi1bert變換之后，通過Mat1ab繪制出信號的Hi1bert譜圖（圖5）和信號邊際譜圖（圖6）。

圖5的Hi1bert譜圖表達了每個不同的頻率在整個時間歷程上所累積的能量分布。橫坐標是時間（采樣點），縱坐標為頻率值。從圖中可以看出，信號的主要振動能量集中在頻率300 Hz以內的低頻段，并且在100 Hz左右能量最為集中。高頻段處包含的能量占總能量的極少部分，這一點從高頻段處散落的頻點上可很容易得出。

圖6的邊際譜是Hi1bert譜對時間軸的積分，它是對信號中各個頻率成份幅值的整體測度，是統計意義上的全部累加幅度。邊際譜反映了信號的幅值隨頻率變化的情況，可以指示任何頻率在全局上的幅度。與Fourier頻譜幅值不同的是，邊際譜中某一頻率的累計幅值出現時，就表示一定有該頻率的振動波出現。從圖6可以看出，信號能量主要分布在低頻段，主頻帶分為兩段，分布在0～80 Hz和100～200 Hz處，在頻率為2 Hz處頻率幅值達到最大值1.483 4。

圖3　EMD分解后的IMF分量（imf1～imf3）

圖4　EMD分解后的IMF分量和殘余量（imf4～imf6、r7）

圖5　Hi1bert譜圖

圖6　信號邊際譜圖

4　結論

文中首先通過獲取各測點的起震時間，計算得出振動波的傳播速度，從而得出兩次振動過程波形傳播方式均為地表傳播和空氣傳播。其次針對爆炸試驗數據非線性非平穩的特征，本文將HHT方法引入到第二次振動過程信號處理中去，有效提取了時頻曲線的主要特征信息，并詳細地分析了其物理意義。從分析中可以看出，HHT方法是分析非平穩信號數據非常有效的方法，此方法將頻譜隨時間的演變關系明確表現出來，通過繪制Hi1bert譜圖、邊際譜等一系列圖形，可以有效提取特征信息。

從Hi1bert邊際譜可知，爆炸振動信號的能量主要集中于低頻分量，與普通工程結構的自振頻率（1～5 Hz）相近，因此爆炸振動信號對一般工程結構造成的破壞比較大。本文基于HHT方法得出的分析結論可以有效地揭露爆炸振動信號的特性，對大藥量炸藥的安全風險評估有一定的參考意義。

［1］Huang N E，ZhengShen，Long S R.The Empirica1 Mode Decomposition and the Hi1bert Spectrum for Non1inear and Non-stationary Time Series Ana1ysis［J］.Proceedings of the Roya1 Society of London，1998，454（1971）：903_995.

［2］N E HUANG.The empirica1 mode decomposition and the-Hi1bert spectrum for non1inear and non-stastinary time series ana1ysis.Proceedings of the Roya1 Society of London series A-Mathematica1 Physica1 and Engineering Sciences［C］// London，1998.

［3］李琳.HHT時頻分析方法的研究與應用［D］.吉林：吉林大學，2006.

［4］胡勁松，楊世錫，吳昭同.基于EMD和HT的旋轉機械振動信號時頻分析［J］.振動、測試與診斷，2004，24（2）：106_110.

［5］蔣禮.運用分頻段希爾伯特黃變換進行多分量信號的頻散分析［J］.電訊技術，2012，52（4）：472_477.

［6］謝桂海，李浩，楊磊.非平穩數據處理方法與瞬時頻率［J］.軍械工程學院學報，2006，18（6）：70_74.

［7］朱瑞賡，李錚.爆炸地震波的現場測量及其安全距離［J］.爆炸與沖擊，1982（1）：60_67.

The aPPllcatlon of HHT ln blastlng slgnal Processlng

SU Xiu-hong
（Institute of System Engineering，，CAEP，Mianyang 621900，China）

HHT has been extensive1y researched in recent years because of its high adaptabi1ity and efficiency as a new method of non-stationary signa1 processing.This artic1e first1y carried out the test of b1ast vibration wave based on the testing system，then ana1yse the vibration signa1 using the Hi1bert Huang Transform.We can study the characteristics of b1ast wave from the resu1ts of data ana1ysis，the ana1ysis can a1so offer reference for safety risk of a 1arge dose of exp1osive.

hi1bert-huang transform；EMD；vibration signa1 ana1ysis；b1asting test

TN 911.21

A

1674_6236（2016）10_0156_03

2015_05_19稿件編號：201505168

蘇秀紅（1988—），女，山東臨沂人，碩士，工程師。研究方向：動態測試與數字信號處理。