程 擂，韓 焱，王 鑒，杜 娟

(中北大學(xué)電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，山西 太原030051)

水中爆炸沖擊波信號(hào)是一種瞬態(tài)非平穩(wěn)信號(hào)，特點(diǎn)是持續(xù)時(shí)間短、頻帶寬。常用的分析方法有3種:時(shí)域分析、頻域分析和時(shí)頻域分析。希爾伯特-黃變換(HHT)［1］是1998 年美國(guó)黃鍔提出的一種處理非平穩(wěn)信號(hào)、且在實(shí)踐中得到充分驗(yàn)證的時(shí)頻域信號(hào)處理手段，一些研究者將其應(yīng)用到了沖擊信號(hào)等非平穩(wěn)信號(hào)的分析中，并且取得了一定的成果。賴(lài)思邈等［2］利用HHT 對(duì)水下目標(biāo)特征進(jìn)行提取，結(jié)果表明HHT 適用于水聲非平穩(wěn)信號(hào)的分析，并且有很強(qiáng)的自適應(yīng)特性和較好的時(shí)頻聚集性。宋敬利等［3］利用HHT 方法分析了某型艦船在非接觸水下爆炸作用下典型部位沖擊響應(yīng)信號(hào)，計(jì)算了Hilbert 幅值譜、Hilbert 能量譜和Hilbert 瞬時(shí)能量，探討了非接觸水下爆炸作用下艦船沖擊響應(yīng)的時(shí)頻特征。

HHT 處理非平穩(wěn)信號(hào)效果很好，但仍存在一些問(wèn)題:第1，在低頻處會(huì)產(chǎn)生一些幅度很小的錯(cuò)誤固有模態(tài)函數(shù)(IMF);第2，在高頻處的第1 級(jí)固有模態(tài)函數(shù)，其頻帶寬度過(guò)大，并不是單一組分;第3，由于干擾問(wèn)題的存在，相鄰的IMF 分量不能被很好地區(qū)分。另外，HHT 還很容易受到噪聲的干擾［4］。

本文中，提出一種改進(jìn)的HHT 算法，用來(lái)對(duì)水中爆炸沖擊波信號(hào)進(jìn)行時(shí)頻分析:首先，采用小波分層的方法對(duì)信號(hào)進(jìn)行預(yù)處理，將輸入信號(hào)分為多個(gè)窄帶子信號(hào)的結(jié)合，使得每一個(gè)IMF 為單一組分;第2，采用基于相關(guān)系數(shù)的IMF 選取原則，對(duì)所有的IMF 分量進(jìn)行篩選，消除錯(cuò)誤的IMF 分量。

1 改進(jìn)的HHT

1.1 HHT

HHT 是一種新的時(shí)頻分析方法，基本思路是:以瞬時(shí)頻率為出發(fā)點(diǎn)，為了分析信號(hào)的瞬時(shí)頻率，希望把信號(hào)分成由有物理意義的瞬時(shí)頻率的各個(gè)信號(hào)分量的組合，包括經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和希爾伯特變換［1］2 個(gè)部分。

1.1.1 經(jīng)驗(yàn)?zāi)B(tài)分解

經(jīng)驗(yàn)?zāi)B(tài)分解包括原始信號(hào)在連續(xù)模式下的信號(hào)分解。從原始信號(hào)中得到的模式對(duì)于信號(hào)s 的分解可以寫(xiě)為式中:ci為第i 分量，稱(chēng)作固有模態(tài)函數(shù)(IMF)，TM是余數(shù)。

IMF 需要滿足2 個(gè)條件:(1)整個(gè)數(shù)據(jù)區(qū)間內(nèi)，極大值和極小值的個(gè)數(shù)需要相同或者最多相差±1;(2)由每個(gè)極小值和極大值產(chǎn)生的包絡(luò)平均必須為零，而且要求分解到最后時(shí)，剩下的余數(shù)不能有極值。

首先，找到信號(hào)s［n］的極值點(diǎn)，用包絡(luò)線分別將上下極值點(diǎn)連接起來(lái)，得到上包絡(luò)線u1和下包絡(luò)線l1，然后算出上下包絡(luò)線的均值，記作m1。此時(shí)m1與原始信號(hào)之差

驗(yàn)證是否滿足IMF 需滿足的2 個(gè)條件。若是，則令c1=h1;若不是，則將h1作為原始數(shù)據(jù)重復(fù)上述過(guò)程

直到找出符合條件的hk，令其為c1，即求出了IMF 的第1 個(gè)分量c1。其次計(jì)算

r1(n)仍然包含原始數(shù)據(jù)的頻率信息，因此將r1(n)作為新信號(hào)，重復(fù)步驟(1)～(4)，循環(huán)k 次，直到解出的IMF 分量ck或者殘余信號(hào)rk非常小或?yàn)閱握{(diào)時(shí)，循環(huán)結(jié)束。

1.1.2 希爾伯特變換

通過(guò)希爾伯特變換，求出幅度-時(shí)間-頻率的三維表示。定義一解析信號(hào)

式中:SH［n］為s［n］的希爾伯特變換

在極坐標(biāo)下用指數(shù)表示為

根據(jù)式(9)即可做出幅度-時(shí)間-頻率的三維表示圖。

1.1.3 希爾伯特邊緣譜

定義了希爾伯特譜后，如在時(shí)間上對(duì)其進(jìn)行積分，即可求得其希爾伯特邊緣譜

與傅立葉頻譜相似，希爾伯特邊緣譜表示了一個(gè)在各個(gè)頻率處幅值或能量的分布圖。

1.2 改進(jìn)的HHT

針對(duì)HHT 的缺陷，首先通過(guò)小波分析將原始信號(hào)分成不同頻段子信號(hào)的預(yù)處理方法。小波的多分辨分析特性能將信號(hào)在不同尺度下進(jìn)行多分辨率的分解，并將交織在一起的各種不同頻率組成的混合信號(hào)分解成不同頻段的子信號(hào)。對(duì)于一個(gè)N 層分解來(lái)說(shuō)，共有N+1 個(gè)途徑分解信號(hào)。在爆炸沖擊波信號(hào)中，高頻部分為主要分析部分，因此采用如圖1 所示的小波分解。其中A 代表高頻分量，D 代表低頻分量。

通過(guò)N 層小波包分解，原始信號(hào)被分解成為N+1 個(gè)不同頻段的子信號(hào)，且每一個(gè)子信號(hào)具有很窄的頻帶，避免了EMD 分解產(chǎn)生非單一組分IMF 的情況。然而該預(yù)處理會(huì)使IMF 分量增加，其中會(huì)包含許多錯(cuò)誤的IMF 分量。

圖1 基于爆炸沖擊特征的小波分解Fig.1 Wavelet decomposition based on the characteristic of explosion shock wave

理想中正確的IMF 分量將會(huì)包含原始信號(hào)中所具有的頻率成分，因此這些IMF 分量與原始信號(hào)應(yīng)具有良好的相關(guān)性。根據(jù)這一特性，采用基于IMF 分量與原始信號(hào)相關(guān)系數(shù)的選取原則，來(lái)對(duì)產(chǎn)生的IMF 分量進(jìn)行篩選，將低頻處幅度很小的錯(cuò)誤IMF 以及通過(guò)小波包分解產(chǎn)生的錯(cuò)誤IMF 去除。設(shè)產(chǎn)生的第i 級(jí)IMF 分量ci與原始信號(hào)的相關(guān)系數(shù)為ri，根據(jù)相關(guān)系數(shù)的定義，則

設(shè)定一閾值ρ，根據(jù)ρ 對(duì)每一個(gè)IMF 分量與原始信號(hào)的ri進(jìn)行對(duì)比。若ri＞ρ，則將該IMF 分量保留，若ri＜ρ，則認(rèn)為該IMF 分量為錯(cuò)誤分量，將其去除。ρ 可以通過(guò)以下方式確定

式中:η 的選擇取決于所需IMF 分量的程度，ρ 越大，IMF 分量越少，反映的高頻信號(hào)越明顯，但因此可能會(huì)丟失某些頻率的信息;ρ 越小，所包含高頻信號(hào)的信息越豐富，但高頻信號(hào)的反映相對(duì)較差。

2 數(shù)值模擬與實(shí)驗(yàn)

2.1 數(shù)值模擬

在上述討論的基礎(chǔ)上，對(duì)模擬出的爆炸沖擊波信號(hào)進(jìn)行數(shù)值分析。爆炸沖擊波的數(shù)學(xué)表達(dá)式為

式中:Ps為常數(shù)，1≤Ps≤3，b=0.5+Ps;t0可以利用Josef Henrych 的經(jīng)驗(yàn)公式來(lái)計(jì)算，即t0=(0.107+0.444R+0.264R2－0.129R3+0.033 5R4)W1/3，其中0.05 ≤R ≤3，W 為炸藥質(zhì)量，kg。

模擬出的爆炸信號(hào)如圖2(a)所示。圖2(c)為加入海洋噪聲后的模擬水中爆炸沖擊波信號(hào)，噪聲如圖2(b)所示。

對(duì)該信號(hào)進(jìn)行HHT，得出該模擬爆炸沖擊波信號(hào)經(jīng)EMD 分解后的IMF 分量，信號(hào)的HHT 時(shí)頻表示圖如圖3 所示。從圖3 中可以看出，分量c1包含了高頻及低頻信息，其頻帶過(guò)寬，不能成為單一組分;分量c6、c7及c8為幅度很小的錯(cuò)誤IMF 低頻分量。在希爾伯特能量譜的時(shí)頻表示圖中，爆炸沖擊波信號(hào)的高頻部分不能很好地被表現(xiàn)出來(lái)，通過(guò)局部放大后方可觀察到具有高能量的高頻信息(如圖3(b)放大圖中的白點(diǎn)部分)。

圖2 模擬爆炸沖擊波信號(hào)Fig.2 Simulated explosion shock wave

圖3 模擬水中爆炸沖擊波信號(hào)的IMF 及HHT 能量分布時(shí)頻圖Fig.3 IMF and HHT time-frequency representation of simulated underwater explosion shock wave

圖4 改進(jìn)后的模擬水中爆炸沖擊波信號(hào)IMF 及HHT 能量分布時(shí)頻圖Fig.4 IMF and HHT time-frequency representation of simulated underwater explosion shock wave using improved HHT

利用HHT 改進(jìn)算法對(duì)信號(hào)進(jìn)行處理，得出經(jīng)篩選后的IMF 分量(如圖4(a)所示)，對(duì)其進(jìn)行希爾伯特變換，得出改進(jìn)的HHT 時(shí)頻表示圖，如圖4(b)所示。

通過(guò)圖4 可以看出，經(jīng)過(guò)小波分層以及相關(guān)系數(shù)篩選后，爆炸沖擊波信號(hào)的IMF 分量中去除了幅度較小的錯(cuò)誤低頻分量，在高頻部分中頻率信號(hào)也較為單一，更好地反映了爆炸沖擊波信號(hào)的EMD 分解。而在能量分布時(shí)頻表示圖中，沖擊信號(hào)中的高頻部分被有效地體現(xiàn)出來(lái)(如圖4(b)中在10 ～20 ms 處的白點(diǎn))，通過(guò)該圖可以看出改進(jìn)的HHT 算法有效地去除了錯(cuò)誤的低頻干擾。對(duì)比圖3 與圖4 的IMF 分量以及HHT 能量時(shí)頻表示圖可以看出，改進(jìn)的HHT 解決了在原始HHT 變換中存在的問(wèn)題，能夠更好地反映出信號(hào)的時(shí)頻特征圖。

2.2 水中爆炸實(shí)驗(yàn)

圖5 所示為水面下70 cm 處，10 g TNT 炸藥的實(shí)測(cè)爆炸沖擊波信號(hào)，沖擊波傳感器位于距離炸點(diǎn)2.45 m的同一水平面上。該信號(hào)為2008 年11 月于地下目標(biāo)毀傷技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室中北大學(xué)爆炸塔處測(cè)得。

采用改進(jìn)的HHT 對(duì)該爆炸沖擊波信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，選取相關(guān)系數(shù)最高的前10 組IMF 分量，并計(jì)算其能量分布時(shí)頻表示圖，結(jié)果如圖6、7 所示。

圖5 實(shí)測(cè)水中爆炸信號(hào)Fig.5 Experimental underwater explosion shock signal

圖6 實(shí)測(cè)水中爆炸信號(hào)IMF 分量Fig.6 IMF components of experimental underwater explosion signal

圖7 實(shí)測(cè)水中爆炸沖擊波信號(hào)的HHT 能量分布譜時(shí)頻表示Fig.7 HHT time-frequency representation of experimental explosion shock wave using improved HHT

從圖6 可以看出，采用改進(jìn)的HHT 可以有效地得出爆炸沖擊波信號(hào)中感興趣的IMF 分量，低頻的IMF 分量通過(guò)相關(guān)系數(shù)篩選法被去除，僅僅保留了在爆炸沖擊波信號(hào)中的高頻信號(hào);另外，由于在HHT之前先進(jìn)行了小波分解，這樣高頻處的第1 級(jí)固有模態(tài)函數(shù)IMF 分量c1的頻帶寬度被降低，解決了原始EMD 分解第1 級(jí)IMF 分量頻帶過(guò)寬的問(wèn)題，使得生成的信號(hào)能量分布時(shí)頻表示更為精確有效。在圖7的能量分布圖中，采用改進(jìn)的HHT，爆炸沖擊波信號(hào)的高頻部分被更好地反映(如圖7 中5 ms 處的時(shí)頻分布)，由于沒(méi)有低頻IMF 分量以及相鄰IMF 分量之間的干擾，信號(hào)的時(shí)頻分布圖為進(jìn)一步的炸點(diǎn)定位、爆炸識(shí)別等奠定了更有效的基礎(chǔ)。

3 結(jié) 論

利用HHT 變換處理非平穩(wěn)信號(hào)的優(yōu)勢(shì)，結(jié)合小波分解及相關(guān)系數(shù)篩選法對(duì)水中爆炸信號(hào)進(jìn)行處理。結(jié)果表明，該方法能夠有效、準(zhǔn)確地分析爆炸沖擊波信號(hào)的時(shí)頻信息，得到爆炸信號(hào)能量分布的時(shí)頻圖，為進(jìn)一步對(duì)爆炸信號(hào)進(jìn)行時(shí)延估計(jì)、效能分析處理等提供了依據(jù)。但如果數(shù)據(jù)量較大時(shí)，在EMD進(jìn)行分解并得出HHT 時(shí)頻分布圖時(shí)需要很多時(shí)間，因此提高算法的速率將是進(jìn)一步研究的關(guān)鍵。

［1］ Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proceedings of the Royal Society of London.Series A:Mathematical，Physical and Engineering Sciences，1998，454(1971):903-995.

［2］ 賴(lài)思邈，張效民，曹丕.HHT 在水下目標(biāo)信號(hào)特征提取中的應(yīng)用［J］.電聲技術(shù)，2009，33(6):36-40.LAI Si-miao，ZHANG Xiao-min，CAO Pi.Application of HHT in feature extraction of underwater target signal［J］.Audio Engineering，2009，33(6):36-40.

［3］ 宋敬利，賈則，張姝紅.希爾伯特-黃變換在艦船沖擊響應(yīng)分析中的應(yīng)用［J］.水雷戰(zhàn)與艦船防護(hù)，2009，17(3):11-15.SONG Jing-li，JIA Ze，ZHANG Shu-hong.Hilbert-Huang transform and its applications on analysis of ships’shock response［J］.Mine Warfare＆Ship Self-Defence，2009，17(3):11-15.

［4］ Penga Z K，Tseb P W，Chua F L.An improved Hilbert-Huang transform and its application in vibration signal analysis［J］.Journal of Sound and Vibration，2005(286):187-205.