爆破振動(dòng)信號(hào)的HHT時(shí)頻能量譜分析＊

關(guān)曉磊，顏景龍

（1.中北大學(xué)電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051；2.北京理工大學(xué)機(jī)電學(xué)院，北京 100081）

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

1 爆破振動(dòng)信號(hào)HHT分析

爆破振動(dòng)信號(hào)HHT分析的步驟如圖1所示。

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

1.1 EMD分解

為了研究瞬態(tài)和非平穩(wěn)現(xiàn)象，頻率必須是時(shí)間的函數(shù)。非平穩(wěn)隨機(jī)信號(hào)通過(guò)EMD過(guò)程后，得到一系列IMF分量，而IMF使瞬時(shí)頻率具有意義。針對(duì)爆破振動(dòng)信號(hào)的信號(hào)特征，在HHT中對(duì)EMD端點(diǎn)效應(yīng)、擬合的欠包絡(luò)現(xiàn)象、篩選終止準(zhǔn)則等問(wèn)題進(jìn)行了改進(jìn)。

1.1.1 首尾極值點(diǎn)對(duì)稱(chēng)延拓

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

步驟如下：計(jì)算端點(diǎn)處首個(gè)極值點(diǎn)x（m）（m為x（n）首個(gè)極值點(diǎn)索引位置，x（n）為X（t）的采樣值）兩側(cè)數(shù)據(jù)對(duì)稱(chēng)度ξ。若ξ在規(guī)定閾值內(nèi)，則以該極值點(diǎn)為對(duì)稱(chēng)軸對(duì)稱(chēng)4～5個(gè)極值點(diǎn)；若首尾極值點(diǎn)兩側(cè)數(shù)據(jù)對(duì)稱(chēng)性較差，則采用鏡像法［4－5］。定義對(duì)稱(chēng)度為

一般取ξ參考門(mén)限為0.1～0.2（實(shí)際ξ越大，說(shuō)明對(duì)稱(chēng)性越差），即：實(shí)際ξ低于門(mén)限，則采用首尾極值點(diǎn)作為對(duì)稱(chēng)軸對(duì)稱(chēng)延拓一定數(shù)量極值點(diǎn)數(shù)據(jù)；若ξ大于門(mén)限，說(shuō)明對(duì)稱(chēng)性差，則采用端點(diǎn)鏡像對(duì)稱(chēng)延拓。

圖2為實(shí)測(cè)數(shù)據(jù)和首極值點(diǎn)對(duì)稱(chēng)向前延拓后的曲線，對(duì)稱(chēng)延拓后曲線仍保持平滑。延拓后進(jìn)行曲線的上下包絡(luò)擬合，得到合格的IMF后截取原始數(shù)據(jù)對(duì)應(yīng)的數(shù)據(jù)。

圖2 原始曲線和首極值點(diǎn)對(duì)稱(chēng)延拓后曲線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

傳統(tǒng)的EMD插值算法采用三次樣條插值，但是容易出現(xiàn)欠包絡(luò)現(xiàn)象，即上下包絡(luò)線局部無(wú)法完整包裹當(dāng)前信號(hào)。于是，采用了分段三次Hemite保形插值算法。它與三次樣條插值算法的不同之處在于計(jì)算插值點(diǎn)處一階導(dǎo)數(shù)的方法不同。保形插值計(jì)算插值點(diǎn)xk一階導(dǎo)數(shù)dk的方法更簡(jiǎn)單、直接，即：當(dāng)δk和δk－1同號(hào)時(shí)，（w1＋w2）／dk＝w1／δk－1＋w2／δk，其中δk是子區(qū)間xk≤x≤xk＋1的折線斜率，w1＝2hk＋hk－1，w2＝hk＋2hk－1（hk＝xk＋1－xk，即區(qū)間步長(zhǎng)）；當(dāng)δk和δk－1異號(hào)時(shí)，dk＝0。

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

1.1.3 篩分終止準(zhǔn)則

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

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

將σ（t）作為判定篩選終止的判據(jù)，設(shè)定門(mén)限值θ1、θ2和α，規(guī)定當(dāng)σ（t）中數(shù)據(jù)小于θ1的比例達(dá)到α、且不存在大于θ2的數(shù)據(jù)時(shí)，終止篩選。一般默認(rèn)θ1＝0.05，θ2＝0.5，α＝0.95。與標(biāo)準(zhǔn)差σs相比，σ（t）更能反映IMF的均值特性，且兩個(gè)條件相互補(bǔ)充，使信號(hào)只能在局部出現(xiàn)較大波動(dòng)，保證整體均值為零［6］。為防止篩分發(fā)散和誤差累積，一般對(duì)篩分次數(shù)做限制。

圖4中，h1（t）頻率范圍較寬，經(jīng)過(guò)不斷篩選得到頻率較為單一的h121（t），即IMF，且h在各個(gè)時(shí)間段上都具有較好的頻率一致性，說(shuō)明局部終止準(zhǔn)則得到的IMF局部特性較好。

1.1.4 模態(tài)分解終止條件

EMD分解得到一個(gè)IMF后，用原信號(hào)減去IMF便是信號(hào)殘差，若殘差為單調(diào)函數(shù)或整體均值小于預(yù)定誤差則完成了EMD。

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

1.2 瞬時(shí)屬性估計(jì)的改進(jìn)算法

采用Hilbert變換法計(jì)算瞬時(shí)頻率，很容易出現(xiàn)無(wú)意義的負(fù)頻率，而基于經(jīng)驗(yàn)的AM－FM分解可以很好地解決這個(gè)問(wèn)題。過(guò)程中，把IMF唯一地分成包絡(luò)部分AM和載波部分FM［7］。其中，F(xiàn)M可直接計(jì)算正交項(xiàng)DQ，且FM是單位振幅，滿(mǎn)足了Bedrosian定理的要求；同時(shí)，標(biāo)準(zhǔn)化的FM可以給出更加局部化、更加實(shí)用的基于能量的誤差。

1.2.1 AM－FM 分解

標(biāo)準(zhǔn)化的過(guò)程基于經(jīng)驗(yàn)的、通過(guò)擬合數(shù)據(jù)的多次迭代實(shí)現(xiàn)。首先，找到IMF數(shù)據(jù)xi（t）絕對(duì)值的所有局部極大值點(diǎn)，然后利用曲線擬合所有極大值點(diǎn)，得到包絡(luò)e1（t）；同樣在過(guò)程中，使用極值點(diǎn)對(duì)稱(chēng)延拓法補(bǔ)充極值點(diǎn)，然后進(jìn)行標(biāo)準(zhǔn)化處理

由于擬合僅穿過(guò)極大值點(diǎn)，在振幅快速變化的地方包絡(luò)樣條可能會(huì)穿到數(shù)據(jù)點(diǎn)之下，所以IMF一般大于單位振幅。為了得到單位振幅的FM，采用迭代法對(duì)y1（t）再進(jìn)行標(biāo)準(zhǔn)化

經(jīng)過(guò)n次迭代后，yn（t）的所有值都等于或者小于1，標(biāo)準(zhǔn)化就結(jié)束了。這時(shí)，得到FM部分Fi（t），同時(shí)得到包絡(luò)AM部分Ai（t），也即調(diào)頻部分，作為瞬時(shí)振幅。這種方法繞開(kāi)了Bedrosian條件的限制，對(duì)于波形局部零均值且對(duì)稱(chēng)的IMF信號(hào)，能夠有效計(jì)算瞬時(shí)振幅

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

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

1.2.2 瞬時(shí)頻率計(jì)算

基于經(jīng)驗(yàn)的AM－FM分解得到了IMF信號(hào)的正交項(xiàng)，跳過(guò)了采用希爾伯特變換計(jì)算瞬時(shí)頻率的過(guò)程。標(biāo)準(zhǔn)化之后得到FM部分為信號(hào)載波部分F（t），得正交項(xiàng)為

它最大的優(yōu)點(diǎn)是沒(méi)有使用希爾伯特變換，不涉及積分過(guò)程，不受臨近點(diǎn)的影響，且瞬時(shí)頻率通過(guò)求導(dǎo)運(yùn)算得到，所以局部性很好。于是，相位角和瞬時(shí)頻率為

式中反正切和反余弦理論上是等價(jià)的，但計(jì)算方法不一樣，反正切的計(jì)算穩(wěn)定性較高。

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

圖7 瞬時(shí)頻率Fig.7 Instantaneous frequencies

2 能量計(jì)算

HHT的定義為

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

瞬時(shí)能量

式（10）反映了信號(hào)能量隨時(shí)間變化的情況，對(duì)于爆破，它代表了雷管爆炸造成的振動(dòng)能量的激發(fā)，可參照譜圖和起爆時(shí)序估計(jì)雷管是否起爆［8］。

以某礦爆破實(shí)測(cè)信號(hào)（豎直方向）為例，爆破振動(dòng)信號(hào)判據(jù)主要用最大振速、主頻和持續(xù)時(shí)間三個(gè)獨(dú)立參數(shù)對(duì)信號(hào)進(jìn)行衡量。采樣頻率5kHz，最大振速6.858cm／s，持續(xù)時(shí)間約400ms，F(xiàn)FT頻域主頻38Hz。從瞬時(shí)能量圖和實(shí)際數(shù)碼電子雷管爆破時(shí)序?qū)Ρ龋ㄒ?jiàn)圖8）可見(jiàn)，每個(gè)尖峰是雷管爆炸產(chǎn)生的能量激發(fā)，能量在120和270ms時(shí)刻附近較高，說(shuō)明這些時(shí)間內(nèi)的爆破振動(dòng)的能量較大，也即爆炸的能量較大，雷管的裝藥量較大。

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

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

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

3 結(jié) 論

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

［1］張義平.爆破震動(dòng)信號(hào)的HHT分析和應(yīng)用研究［D］.武漢：中南大學(xué)，2006.

［2］Huang N E，WU Zhao－h(huán)ua，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］張進(jìn)林.抑制經(jīng)驗(yàn)?zāi)B(tài)分解端點(diǎn)效應(yīng)的常用方法性能比較研究［D］.昆明：云南大學(xué)，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］王國(guó)棟.信號(hào)瞬時(shí)頻率的估算［D］.廣州：中山大學(xué)，2010.

［8］陳銀魯.爆破振動(dòng)信號(hào)的能量分析方法及其應(yīng)用研究［D］.杭州：浙江大學(xué)，2010.