爆破振動信號的HHT時頻能量譜分析＊

關曉磊，顏景龍

（1.中北大學電子測試技術國家重點實驗室，山西 太原 030051；2.北京理工大學機電學院，北京 100081）

隨著信號分析技術的發展，對爆破的了解已不止于時域和頻域的獨立分析，短時傅立葉變換和小波變換在一定程度上可對信號進行時頻分析，但都基于傅立葉變換，要求被分析信號必須平穩。爆破振動信號是典型的非平穩隨機信號，嚴格意義上基于傅立葉變換的分析方法不具有適用性［1］。希爾伯特－黃變換（Hilbert－Huang transform，HHT）是近幾年興起的非平穩隨機信號分析新方法［2－3］，它對信號進行平穩化處理，即將信號經過經驗模態分解（empirical mode decomposition，EMD），產生一系列具有不同特征尺度的固有模態函數（intrinsic mode function，IMF），然后計算IMF信號的瞬時屬性。用時頻域聯合能量法分析振動信號，能夠較好反映振動強度、頻譜特征，特別是持續時間對建筑物的影響，比現行的獨立域參數更科學。

1 爆破振動信號HHT分析

爆破振動信號HHT分析的步驟如圖1所示。

圖1 HHT分析步驟Fig.1 The HHT analysis steps

1.1 EMD分解

為了研究瞬態和非平穩現象，頻率必須是時間的函數。非平穩隨機信號通過EMD過程后，得到一系列IMF分量，而IMF使瞬時頻率具有意義。針對爆破振動信號的信號特征，在HHT中對EMD端點效應、擬合的欠包絡現象、篩選終止準則等問題進行了改進。

1.1.1 首尾極值點對稱延拓

對于原始信號X（t），先要找出所有極值點進行上下包絡的擬合。由于數據首尾段數據缺少極值點，擬合時就會造成端點飛翼現象。爆破振動信號屬于低頻信號，極值點之間的時間跨度大，端部的端點效應容易影響信號內部［4］。通常在數據兩端加入特征波抑制端點效應。微差爆破（毫秒延期爆破）能夠很好地實現降震，所以現有工程爆破多為微差爆破，這使得爆破振動的首尾極值點兩側數據多具有對稱特點。利用這個特性進行對稱延拓，消除爆破振動信號HHT變換的端點效應。

步驟如下：計算端點處首個極值點x（m）（m為x（n）首個極值點索引位置，x（n）為X（t）的采樣值）兩側數據對稱度ξ。若ξ在規定閾值內，則以該極值點為對稱軸對稱4～5個極值點；若首尾極值點兩側數據對稱性較差，則采用鏡像法［4－5］。定義對稱度為

一般取ξ參考門限為0.1～0.2（實際ξ越大，說明對稱性越差），即：實際ξ低于門限，則采用首尾極值點作為對稱軸對稱延拓一定數量極值點數據；若ξ大于門限，說明對稱性差，則采用端點鏡像對稱延拓。

圖2為實測數據和首極值點對稱向前延拓后的曲線，對稱延拓后曲線仍保持平滑。延拓后進行曲線的上下包絡擬合，得到合格的IMF后截取原始數據對應的數據。

圖2 原始曲線和首極值點對稱延拓后曲線Fig.2 The original curve and the curve after extension

1.1.2 插值算法

圖3 三次樣條插值與分段三次Hemite保形插值Fig.3 The cubic spline interpolation and the Hermite interpolation of picecwise cubic

傳統的EMD插值算法采用三次樣條插值，但是容易出現欠包絡現象，即上下包絡線局部無法完整包裹當前信號。于是，采用了分段三次Hemite保形插值算法。它與三次樣條插值算法的不同之處在于計算插值點處一階導數的方法不同。保形插值計算插值點xk一階導數dk的方法更簡單、直接，即：當δk和δk－1同號時，（w1＋w2）／dk＝w1／δk－1＋w2／δk，其中δk是子區間xk≤x≤xk＋1的折線斜率，w1＝2hk＋hk－1，w2＝hk＋2hk－1（hk＝xk＋1－xk，即區間步長）；當δk和δk－1異號時，dk＝0。

由圖3可以看出，分段三次Hemite保形插值較好地保留了原始曲線的形狀，而三次樣條插值的效果欠佳，在300和320ms時刻的包絡發生形狀變化，出現過包絡。

1.1.3 篩分終止準則

擬合得到上下包絡線Xmax（t）和Xmin（t）后，能夠得到一條均值線，然后得到h1（t）

需要對h1（t）進行判斷，判斷h1（t）是否為合格的IMF分量，若不是則將h1（t）作為原信號X（t）重復上述步驟，直到得到合格的h1k（t），即IMF0。一般地，采用基于整體數據標準差σs必然忽略了數據的局部特性，所以出現了基于局部的終止條件。設局部均值線為q（t）＝［Xmax（t）－Xmin（t）］／2，幅值函數為a（t）＝［Xmax（t）＋Xmin（t）］／2，定義函數

將σ（t）作為判定篩選終止的判據，設定門限值θ1、θ2和α，規定當σ（t）中數據小于θ1的比例達到α、且不存在大于θ2的數據時，終止篩選。一般默認θ1＝0.05，θ2＝0.5，α＝0.95。與標準差σs相比，σ（t）更能反映IMF的均值特性，且兩個條件相互補充，使信號只能在局部出現較大波動，保證整體均值為零［6］。為防止篩分發散和誤差累積，一般對篩分次數做限制。

圖4中，h1（t）頻率范圍較寬，經過不斷篩選得到頻率較為單一的h121（t），即IMF，且h在各個時間段上都具有較好的頻率一致性，說明局部終止準則得到的IMF局部特性較好。

1.1.4 模態分解終止條件

EMD分解得到一個IMF后，用原信號減去IMF便是信號殘差，若殘差為單調函數或整體均值小于預定誤差則完成了EMD。

圖5為信號分解得到的IMF及其頻譜。由圖5可以看出，從IFM0～IFM7頻率逐漸降低，說明了EMD具有自適應性，且主要分布在50Hz以下，到最后幾乎成為一條直線；IMF0頻帶較寬，含有較多的高頻量，說明含有高頻噪聲；IMF1～IMF3的幅度比其他分量高，說明振動有用信號主要分布在這些分量中。

1.2 瞬時屬性估計的改進算法

采用Hilbert變換法計算瞬時頻率，很容易出現無意義的負頻率，而基于經驗的AM－FM分解可以很好地解決這個問題。過程中，把IMF唯一地分成包絡部分AM和載波部分FM［7］。其中，FM可直接計算正交項DQ，且FM是單位振幅，滿足了Bedrosian定理的要求；同時，標準化的FM可以給出更加局部化、更加實用的基于能量的誤差。

1.2.1 AM－FM 分解

標準化的過程基于經驗的、通過擬合數據的多次迭代實現。首先，找到IMF數據xi（t）絕對值的所有局部極大值點，然后利用曲線擬合所有極大值點，得到包絡e1（t）；同樣在過程中，使用極值點對稱延拓法補充極值點，然后進行標準化處理

由于擬合僅穿過極大值點，在振幅快速變化的地方包絡樣條可能會穿到數據點之下，所以IMF一般大于單位振幅。為了得到單位振幅的FM，采用迭代法對y1（t）再進行標準化

經過n次迭代后，yn（t）的所有值都等于或者小于1，標準化就結束了。這時，得到FM部分Fi（t），同時得到包絡AM部分Ai（t），也即調頻部分，作為瞬時振幅。這種方法繞開了Bedrosian條件的限制，對于波形局部零均值且對稱的IMF信號，能夠有效計算瞬時振幅

圖6 迭代算法的過程Fig.6 The process of the iterative algorithm

由圖6可以看出，FM0曲線存在幅度大于1的數據段，一次迭代后FM1接近單位振幅，最后FM2完全為單位振幅，即IMF0的FM部分F0（t），然后通過計算得到AM部分。通過對FM部分進行計算得到正交項，進而計算相位角，最后得到瞬時頻率。

1.2.2 瞬時頻率計算

基于經驗的AM－FM分解得到了IMF信號的正交項，跳過了采用希爾伯特變換計算瞬時頻率的過程。標準化之后得到FM部分為信號載波部分F（t），得正交項為

它最大的優點是沒有使用希爾伯特變換，不涉及積分過程，不受臨近點的影響，且瞬時頻率通過求導運算得到，所以局部性很好。于是，相位角和瞬時頻率為

式中反正切和反余弦理論上是等價的，但計算方法不一樣，反正切的計算穩定性較高。

由圖7可以看出，所有IMF分量的瞬時頻率主要集中在60Hz以下，60Hz以上的為IMF0的部分時刻，說明高頻噪聲主要集中在IMF0中。

圖7 瞬時頻率Fig.7 Instantaneous frequencies

2 能量計算

HHT的定義為

式中：瞬時頻率為各IMF經過AM－FM分解后的FM部分經過直接正交計算所得的ωi（t），H（ω，t）為相應IMF的AM部分Ai（t）幅值。三維時頻能量譜為H2（ω，t）隨時間和瞬時頻率的分布。

瞬時能量

式（10）反映了信號能量隨時間變化的情況，對于爆破，它代表了雷管爆炸造成的振動能量的激發，可參照譜圖和起爆時序估計雷管是否起爆［8］。

以某礦爆破實測信號（豎直方向）為例，爆破振動信號判據主要用最大振速、主頻和持續時間三個獨立參數對信號進行衡量。采樣頻率5kHz，最大振速6.858cm／s，持續時間約400ms，FFT頻域主頻38Hz。從瞬時能量圖和實際數碼電子雷管爆破時序對比（見圖8）可見，每個尖峰是雷管爆炸產生的能量激發，能量在120和270ms時刻附近較高，說明這些時間內的爆破振動的能量較大，也即爆炸的能量較大，雷管的裝藥量較大。

三維時頻能量譜如圖9所示。HHT時頻能量譜能夠直觀地從時域（持續時間）和頻域（主振頻帶）角度反映爆破振動信號的能量分布。相對于單純的時域分析，HHT能夠表明任意頻段和幅值范圍的振動的持續時間；相對于FFT分析，HHT能夠表明任意時間段和范圍的振動能量的頻率分布。如圖9所示，能量主要集中在10～60Hz，時間主要集中在120和270ms左右，能量集中范圍中心頻率40Hz左右。按照GB6722－2004《爆破安全規程》，建筑物的允許頻率在10Hz以內，因此爆破不會對建筑物構成影響。

圖8 原始數據曲線和瞬時能量譜Fig.8 The original data curves and the instantaneous energy spectrum

圖9 FFT幅度譜和三維時頻能量譜Fig.9 The FFT amplitude spectrum and the three－dimensional spectrum of time－frequency power

3 結 論

采用HHT變換可以對爆破振動信號進行更加科學的分析。對于爆破振動信號，按照信號特征采用了首尾極值點對稱法對信號進行延拓，較好地抑制了端點效應；采用分段三次Hemite保形插值對極值點進行擬合，比三次樣條插值擬合更能保留信號的原始外形；局部的終止準則能夠較好地保留信號的局部特征，但會出現篩分次數較多、篩分發散的問題；采用AM－FM分解計算IMF瞬時屬性，比Hilbert變換更準確，并且不會出現負頻率現象。最后，HHT時頻能量譜由時頻角度分析振動信號能量分布，可更好地參照爆破安全規程就爆破振動對建筑物是否有影響作出判斷。

［1］張義平.爆破震動信號的HHT分析和應用研究［D］.武漢：中南大學，2006.

［2］Huang N E，WU Zhao－hua，Long S R，et al.On instantaneous frequency［J］.Advances in Adaptive Data Analysis，2009，1（2）：177－229.

［3］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition method and the Hilbert spectrum for nonstationary time series analysis［J］.Proceedings of the Royal Society of London：A，1998，454：903－995.

［4］張進林.抑制經驗模態分解端點效應的常用方法性能比較研究［D］.昆明：云南大學，2010.

［5］WU Fang－ji，QU Liang－sheng.An improved method for restraining the end effect in empirical mode decomposition and its applications to the fault diagnosis of large rotating machinery［J］.Journal of Sound and Vibration，2008，314（3）：586－602.

［6］Niazy R K，Beckmann C F，Brady J M，et al.Performance evaluation of ensemble empirical mode decomposition［J］.Advances in Adaptive Data Analysis，2009，1（2）：231－242.

［7］王國棟.信號瞬時頻率的估算［D］.廣州：中山大學，2010.

［8］陳銀魯.爆破振動信號的能量分析方法及其應用研究［D］.杭州：浙江大學，2010.