改進(jìn)HHT 算法在導(dǎo)彈工作模態(tài)辨識(shí)中的應(yīng)用*

謝金松，薛林，高慶豐

（1.北京電子工程總體研究所，北京 100854；2.中國航天科工集團(tuán)有限公司 第二研究院，北京 100854）

0 引言

導(dǎo)彈飛行時(shí)彈體受到周身分布的復(fù)雜激勵(lì)，會(huì)產(chǎn)生彈性振動(dòng)。橫向振動(dòng)被彈上加表和陀螺儀量測(cè)，使回路引入額外輸入，嚴(yán)重時(shí)會(huì)造成導(dǎo)彈失穩(wěn)，如美國的Lark 和印度ASLV-D2 火箭，因此需對(duì)導(dǎo)彈彈性振動(dòng)進(jìn)行分析。工程中常通過地面模態(tài)試驗(yàn)確定導(dǎo)彈模態(tài)，依據(jù)模態(tài)參數(shù)設(shè)計(jì)凹陷濾波器以抑制振動(dòng)影響［1-2］。導(dǎo)彈實(shí)際飛行工況與地面振動(dòng)試驗(yàn)環(huán)境存在差異：一方面，導(dǎo)彈飛行中噪聲源（環(huán)境雜波、熱噪聲等）復(fù)雜且激勵(lì)（氣動(dòng)力與擾流）無法量測(cè)；另一方面，導(dǎo)彈彈性頻率會(huì)隨燃料燃燒造成的質(zhì)量衰減而改變，是典型的時(shí)變過程。較之地面振動(dòng)試驗(yàn)，辨識(shí)導(dǎo)彈工作模態(tài)更有意義。

利用傳統(tǒng)算法辨識(shí)導(dǎo)彈工作模態(tài)存在缺陷：時(shí)域或頻域辨識(shí)方法都難以對(duì)時(shí)變模態(tài)進(jìn)行分析；時(shí)域辨識(shí)對(duì)噪聲敏感，辨識(shí)精度低；頻域辨識(shí)存在傅里葉變換截?cái)嗾`差，且大都要求激勵(lì)已知［3］。Hilbert-Huang 變換（HHT）是一種非平穩(wěn)信號(hào)處理技術(shù)［4］，與傳統(tǒng)的數(shù)據(jù)分析方法相比，HHT 具有分析精度高、自適應(yīng)性強(qiáng)、物理意義明確、可處理非平穩(wěn)信號(hào)等優(yōu)勢(shì)。HHT 存在端點(diǎn)效應(yīng)和模態(tài)混疊2 個(gè)固有缺陷，國內(nèi)外學(xué)界為解決此二類問題提出了大量算法。如針對(duì)端點(diǎn)效應(yīng)的鏡像延拓方法［5］、ARMA 時(shí)間序列線性預(yù)測(cè)方法［6］、多項(xiàng)式擬合延拓方法［7］、邊界波形匹配延拓法［8］等，均可有效抑制端點(diǎn)漂移現(xiàn)象。Huang 等［9］提 出 的 集 合 經(jīng) 驗(yàn) 模 態(tài) 分 解（en?semble empirical mode decomposition，EEMD）是抑制模態(tài)混疊的經(jīng)典算法，后人對(duì)模態(tài)混疊的抑制大都為該算法的改進(jìn)；Yeh 等［10］提出的互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解方法（complementary ensemble empirical mode decomposition，CEEMD）有效減少了EEMD 算法中輔助噪聲的影響，練繼建等［11］將該方法應(yīng)用到水壩模態(tài)識(shí)別中；陳永高等［12］參考肖瑛等［13］的思路，認(rèn)為各階IMF 的非確定正交性是模態(tài)混疊的來源，通過在EEMD 算法中嵌入“解相關(guān)算法”和“譜系聚類”方法使各階IMF 強(qiáng)制正交，對(duì)模態(tài)混疊加以抑制，并對(duì)斜拉橋振動(dòng)模態(tài)進(jìn)行辨識(shí)。

本文提出一種多次自適應(yīng)加噪的改進(jìn)HHT 模態(tài)辨識(shí)方法，以兩端自由的變質(zhì)量歐拉梁作為飛行導(dǎo)彈模型，仿真出飛行狀態(tài)下的導(dǎo)彈前三階振動(dòng)模態(tài)，使用改進(jìn)后的HHT 算法辨識(shí)其振動(dòng)參數(shù)。與傳統(tǒng)算法相比，改進(jìn)HHT 方法有效提高了模態(tài)篩分精度，更適用于具有時(shí)變性的導(dǎo)彈工作模態(tài)辨識(shí)。

1 HHT 模態(tài)辨識(shí)方法

HHT 算法包含經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition，EMD）和希爾伯特變換2 個(gè)步驟。經(jīng)驗(yàn)?zāi)B(tài)分解是將不同時(shí)間尺度的信號(hào)分量提取出來，即篩分出振動(dòng)信號(hào)中的各階固有模態(tài)函數(shù)（in?trinsic mode function，IMF）。EMD 結(jié) 束 后，對(duì) 各 階IMF 進(jìn)行希爾伯特變換，擬合變換后的解析信號(hào)即可解算出各階頻率與阻尼比。

利用EMD 篩分振動(dòng)信號(hào)各階模態(tài)的流程如下：

（1）對(duì)原始信號(hào)x(t)的極大值點(diǎn)和極小值點(diǎn)分別進(jìn)行3 次樣條插值，求出信號(hào)上下包絡(luò)線xmax(t)和xmin(t)。

（2）得到上下包絡(luò)線xmax(t)和xmin(t)的均值函數(shù)m(t) = (xmax(t) +xmin(t))/2。

（3）用原始信號(hào)x(t)減去m(t)，得到待判別信號(hào)h1(t)。

（4）若h1(t)滿足IMF 條件，即認(rèn)為h1(t)為原始信號(hào)的一階IMF，用原始信號(hào)減去h1(t)得到的新信號(hào)c1(t)作為原始信號(hào)，重復(fù)以上步驟，直至把所有的IMF 全部找出。

（5）若h1(t)不滿足IMF所需條件，用h1(t)作為原始信號(hào)重復(fù)以上步驟，直至把所有的IMF 全部找出。

Huang 給出的IMF 判定條件為：信號(hào)零極點(diǎn)數(shù)目相等或差1；信號(hào)局部包絡(luò)線均值為0。此條件使各階IMF 具有均勻的波動(dòng)性與對(duì)稱性，滿足各階振動(dòng)模態(tài)的要求。

篩分完成后，對(duì)各階模態(tài)分別進(jìn)行希爾伯特變換得到其在復(fù)域上的解析信號(hào)：

式中：P.V.表示柯西主值積分。

各階模態(tài)信號(hào)的瞬時(shí)幅相函數(shù)為

已知單自由度振動(dòng)表達(dá)式為

式中：A為初始振幅；ωdi(t)為第i階模態(tài)有阻尼自振角速度；ωi(t)為該階模態(tài)無阻尼自振角速度；ξi為該階模態(tài)阻尼比。

聯(lián)立式（2）與式（3），可得各階振動(dòng)模態(tài)幅相函數(shù)的表達(dá)式：

利用式（2）求單階信號(hào)的幅相函數(shù)后，對(duì)幅值函數(shù)取對(duì)數(shù)，得出的曲線斜率即為ξi ωi(t)，相位曲線各個(gè)時(shí)刻的導(dǎo)數(shù)為該時(shí)刻的有阻尼瞬時(shí)頻率。二式聯(lián)立，即可解算出該階模態(tài)的阻尼比ξi與角速度ωi(t)。在導(dǎo)彈飛行過程中，由于質(zhì)量減小，各階模態(tài)振動(dòng)頻率不斷改變，希爾伯特變換提出的瞬時(shí)頻率的概念可以準(zhǔn)確反應(yīng)模態(tài)頻率隨時(shí)間變化的趨勢(shì)。

2 改進(jìn)的HHT 模態(tài)辨識(shí)方法

EMD 存在端點(diǎn)效應(yīng)和模態(tài)混疊2 個(gè)缺陷。端點(diǎn)效應(yīng)是IMF 在信號(hào)初末位置的漂移現(xiàn)象，原因?yàn)樾盘?hào)的兩端點(diǎn)在篩分過程中被視作極點(diǎn)，然而其未必為真實(shí)極值，造成極值擬合的包絡(luò)線在端點(diǎn)處出現(xiàn)漂移。模態(tài)混疊現(xiàn)象是EMD 過程難以避免的固有缺陷：由于EMD 算法在數(shù)學(xué)上無法保證各階IMF 具有嚴(yán)格正交性，因此篩分出的某階IMF 中可能包含不同特征尺度的信號(hào)。

2.1 極值鏡像延拓法抑制端點(diǎn)效應(yīng)

學(xué)界關(guān)于端點(diǎn)效應(yīng)抑制研究已經(jīng)頗有成果：鏡像延拓、極值延拓、數(shù)據(jù)預(yù)測(cè)、波形匹配延拓均能有效抑制端點(diǎn)效應(yīng)。本文通過對(duì)信號(hào)兩端數(shù)據(jù)的極值點(diǎn)進(jìn)行鏡像延拓消除端點(diǎn)效應(yīng)的影響。在信號(hào)左端，摒棄第一采樣點(diǎn)到第一極大值或第一極小值之間的數(shù)據(jù)，利用第一極大值或極小值作為鏡面，對(duì)其后5 個(gè)極值進(jìn)行反射，在信號(hào)末端做相同處理，構(gòu)成新的極值信號(hào)段。通過極值鏡像延拓，可以避免采樣初末點(diǎn)參與EMD 過程，消除端點(diǎn)效應(yīng)。

2.2 集合經(jīng)驗(yàn)?zāi)B(tài)分解抑制模態(tài)混疊

噪聲與高頻信號(hào)的干擾、各階模態(tài)頻率接近都可能造成模態(tài)混疊。根據(jù)經(jīng)驗(yàn)，當(dāng)某模態(tài)瞬時(shí)頻率小于另一模態(tài)瞬時(shí)頻率的2 倍，或者不同階模態(tài)的瞬時(shí)角速度與幅度之積不能始終保持Ai ωi≥Aj ωj時(shí)便會(huì)出現(xiàn)模態(tài)混疊現(xiàn)象。

EEMD 算法利用白噪聲的零均值性，在原始信號(hào)中施加輔助白噪聲后再進(jìn)行EMD，通過多次加噪后加權(quán)平均來中和高頻噪聲和施加噪聲的影響，其具體算法流程如圖1 所示。

圖1 EEMD 算法流程Fig.1 Flowchart of EEMD algorithm

圖1 中，ni(t)為輔助白噪聲，k為輔助白噪聲標(biāo)準(zhǔn)差與原始信號(hào)標(biāo)準(zhǔn)差的比值，σ(x)為原始信號(hào)標(biāo)準(zhǔn)差。經(jīng)過實(shí)驗(yàn)驗(yàn)證，k取0.2，施加300 次噪聲時(shí)篩分效果最好，此時(shí)輔助噪聲對(duì)篩分結(jié)果的影響小于1%。

2.3 改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解方法

通過分析EEMD 分解結(jié)果發(fā)現(xiàn)，系統(tǒng)對(duì)高頻IMF 的篩分精度最高，對(duì)低頻IMF 篩分精度較低；在振動(dòng)能量較大時(shí)段篩分精度較高，在振動(dòng)能量較小時(shí)段篩分精度較低。該現(xiàn)象原因?yàn)檩o助噪聲僅符合高頻IMF 功率比特性，且輔助白噪聲功率時(shí)不變，與逐漸衰減的原始信號(hào)功率變化趨勢(shì)不匹配，故而只對(duì)高頻高能信號(hào)篩分效果好?；诖朔治?，作者對(duì)EEMD 算法提出2 種改進(jìn)方案：

（1）每次EMD 過程中多次添加輔助噪聲：篩分出某階IMF 后，在殘余數(shù)據(jù)中繼續(xù)加噪，且該噪聲應(yīng)符合待分解信號(hào)的功率比要求。此舉可以提高高頻高能IMF 之外的各階IMF 的篩分精度。

（2）為避免隨著原始信號(hào)衰減，輔助噪聲與其功率比逐漸增加，輔助噪聲與原始信號(hào)應(yīng)具備相似的功率衰減率。通過實(shí)驗(yàn)驗(yàn)證，輔助噪聲變?yōu)槿缡剑?）形式時(shí)篩分效果最好：

式中：|x(1)(t)| 為 歸 一化的待分解信 號(hào)；k，σ(x)，ni(t)與EEMD 算法中取值相同。改進(jìn)后的算法如圖2 所示。

圖2 改進(jìn)的EEMD 算法流程Fig.2 Flowchart of improved EEMD algorithm

3 模態(tài)辨識(shí)分析

3.1 導(dǎo)彈振動(dòng)仿真

Euler-Bernoulli 梁兩端自由，且忽略其內(nèi)剪力作用，常用以分析細(xì)長體彈性振動(dòng)問題。本文采用變質(zhì)量Euler-Bernoulli 梁模擬飛行導(dǎo)彈彈性振動(dòng)。彈體模型如圖3 所示。彈體參數(shù)由表1 給出。

表1 導(dǎo)彈彈體參數(shù)Table 1 Missile body parameters

圖3 導(dǎo)彈簡(jiǎn)化歐拉梁模型Fig.3 The missile is simplified to Euler-Bernoulli beam

設(shè)導(dǎo)彈前三階阻尼比ξ1= 0.008 5，ξ2= 0.022 1，ξ3= 0.033 2。在發(fā)動(dòng)機(jī)開機(jī)的1 s 內(nèi)，導(dǎo)彈壓心位置受到了持續(xù)0.3 s 大小1 800 N 的側(cè)向力，由于燃料消耗，導(dǎo)彈體密度以60 kg/m3/s 的速度持續(xù)衰減。根據(jù)Euler-Bernoulli 梁的固有頻率經(jīng)驗(yàn)公式：

可解得該秒內(nèi)導(dǎo)彈一階頻率由41.24 Hz 增加至41.96 Hz；二階頻率由114.56 Hz 增加至116.58 Hz；三階頻率由224.54 Hz 增加至228.49 Hz。

已 知Euler-Bernoulli 梁 受 迫振動(dòng)方 程（7）和Euler-Bernoulli 梁各位置振型方程（8）：

聯(lián)立式（7），（8）可解得各階模態(tài)廣義坐標(biāo)：

根據(jù)振型疊加公式可得出彈上加表由于導(dǎo)彈彈性振動(dòng)測(cè)得的相應(yīng)加速度：

設(shè)彈上加表距離彈頭距離1 264 mm，量測(cè)信號(hào)的信噪比為10，該秒內(nèi)加表量測(cè)到的振動(dòng)信號(hào)如圖4 所示。

圖4 彈上加表量測(cè)到的振動(dòng)信號(hào)Fig.4 Vibration signal measured by accelerometer

3.2 各階模態(tài)篩分

利用圖4 數(shù)據(jù)模擬導(dǎo)彈飛行過程中加表量測(cè)到的振動(dòng)信號(hào)，截取0.3~0.7 s 內(nèi)的數(shù)據(jù)，利用EMD 方法篩分出的各階IMF 如圖5 所示。

由圖5 可以看出，EMD 篩分后各階IMF 左端均出現(xiàn)了端點(diǎn)漂移，各階IMF在0.35 s前后以及IMF3中0.45 s 前后出現(xiàn)了模態(tài)混疊，其原因?yàn)楦麟AIMF 幅值與頻率的乘積Ai ωi在這個(gè)時(shí)間點(diǎn)出現(xiàn)了大小交替。

圖5 EMD 算法分解結(jié)果Fig.5 Decomposition results of EMD

利用極值鏡像延拓方法和EEMD 算法對(duì)端點(diǎn)漂移和模態(tài)混疊進(jìn)行抑制，分解結(jié)果如圖6 所示。

圖6 EEMD 算法分解結(jié)果Fig.6 Decomposition results of EEMD

由圖6 可以看出，極值鏡像延拓有效抑制了各階IMF 左側(cè)的端點(diǎn)效應(yīng)，EEMD 有效抑制了各階IMF 中的模態(tài)混疊，但I(xiàn)MF2 中0.45 s 之后仍出現(xiàn)了模態(tài)混疊。

利用改進(jìn)后的EEMD 算法進(jìn)行篩分，分解結(jié)果如圖7 所示。

由圖7 可以看出，對(duì)EEMD 進(jìn)行改進(jìn)后，IMF2中0.45 s 之后出現(xiàn)的模態(tài)混疊也被抑制。

圖7 改進(jìn)EEMD 算法分解結(jié)果Fig.7 Decomposition results of improved EEMD

為了定量分析各算法篩分精度，引入?yún)?shù)v作為精度評(píng)價(jià)指標(biāo)：

式中：下標(biāo)2L表示向量的2 范數(shù)。

分解精度越高，參數(shù)v取值越小。以上3 種算法得出的各階IMFv值如表2 所示。

表2 各算法分解精度vTable 2 Decomposition accuracy of each algorithm

結(jié)合本節(jié)圖表仿真結(jié)果可知，EMD 幾乎無法抑制模態(tài)混疊和端點(diǎn)效應(yīng)，EEMD 雖然提高了辨識(shí)精度，但仍無法完全抑制模態(tài)混疊，改進(jìn)后的EEMD 篩分精度最高。

3.3 振動(dòng)參數(shù)辨識(shí)

傳統(tǒng)工程中常使用峰值對(duì)數(shù)衰減法進(jìn)行模態(tài)辨識(shí)：提取各階衰減信號(hào)的峰值點(diǎn)的對(duì)數(shù)：

對(duì)lnci(t) 進(jìn)行一階擬合，求得斜率即為ξi ωi(t)，再根據(jù)峰值點(diǎn)的時(shí)間間隔確定振動(dòng)角速率ωi(t)，從而可以解算出阻尼比ξi。

以IMF1 為例，對(duì)其峰值對(duì)數(shù)擬合結(jié)果如圖8所示。

圖8 峰值對(duì)數(shù)衰減法辨識(shí)過程Fig.8 Identification process of peak logarithmic attenuation

峰值對(duì)數(shù)擬合得到的直線斜率ξ1ω1為2.045 8。根據(jù)峰值時(shí)間差的均值則可得出采樣時(shí)間內(nèi)平均周期0.022 5 s，可算出利用該算法辨識(shí)出的頻率為44.4 Hz，阻尼比為0.007 4。

依據(jù)式（4），利用希爾伯特變換得出一階IMF 幅相函數(shù)后，對(duì)幅值對(duì)數(shù)進(jìn)行擬合，如圖9 所示。

圖9 一階IMF 幅值對(duì)數(shù)變化趨勢(shì)Fig.9 Variation trend of the logarithm of the first oder IMF amplitude

對(duì)相位進(jìn)行求導(dǎo)，可以得出瞬時(shí)頻率變化趨勢(shì)如圖10 所示。

圖10 一階IMF 瞬時(shí)頻率變化趨勢(shì)Fig.10 Variation trend of the instantaneous frequency of first oder IMF

對(duì)該時(shí)段內(nèi)所有信號(hào)幅值對(duì)數(shù)進(jìn)行擬合，得到的直線斜率ξ1ω1為2.212 2。由圖10 中的擬合直線可以看出，該時(shí)段內(nèi)一階IMF 瞬時(shí)頻率有增加的趨勢(shì)，瞬時(shí)頻率由260.68 rad/s 增加到了264.51 rad/s，正是這一時(shí)段內(nèi)由于導(dǎo)彈質(zhì)量衰減造成的頻率增加。將瞬時(shí)頻率均值與ξ1ω1取值進(jìn)行聯(lián)立，可得利用該算法辨識(shí)出的阻尼比為0.008。

分別利用以上2 種方法對(duì)三階模態(tài)參數(shù)進(jìn)行辨識(shí)，結(jié)果如表3 所示。

希爾伯特辨識(shí)用到采樣周期內(nèi)所有數(shù)據(jù)點(diǎn)，比僅利用的峰值數(shù)據(jù)的對(duì)數(shù)衰減法對(duì)噪聲影響的容錯(cuò)率更高，也避免掉了采樣頻率不足造成的誤差。由于引入瞬時(shí)頻率的概念，希爾伯特辨識(shí)方法可以反應(yīng)出即導(dǎo)彈的頻率的時(shí)變性，這是峰值對(duì)數(shù)衰減法無法做到的。

分析表3 所示結(jié)果可知，無論哪種辨識(shí)方法，IMF2 與IMF3 的辨識(shí)精度不如IMF1，其原因?yàn)楦哳l信號(hào)頻率與采樣頻率更近，峰值處的采樣值與真實(shí)極值的偏差更大，且高頻信號(hào)衰減迅速，衰減至低能后噪聲影響增大，故而辨識(shí)精度較低。考慮到高頻信號(hào)衰減迅速且頻率高于制導(dǎo)回路帶寬，工程中常采用低通濾波方法將其濾去或忽略其影響。

表3 各算法辨識(shí)結(jié)果Table 3 Identification results of each algorithm

4 結(jié)束語

本文對(duì)于HHT 方法進(jìn)行了改進(jìn)，通過極值鏡像延拓方法抑制了模態(tài)分解時(shí)的端點(diǎn)效應(yīng)，通過多次施加自適應(yīng)輔助白噪聲的改進(jìn)EEMD 方法抑制模態(tài)混疊，并對(duì)導(dǎo)彈振動(dòng)的仿真信號(hào)進(jìn)行處理，成功辨識(shí)出各階振動(dòng)頻率和阻尼比。

實(shí)際飛行之中，考慮到導(dǎo)彈的體密度變化較大，發(fā)動(dòng)機(jī)初始點(diǎn)火階段振動(dòng)頻率約為空彈振頻的0.75 倍左右，利用該方法對(duì)導(dǎo)彈在不同空域內(nèi)彈道飛行的加表量測(cè)信號(hào)進(jìn)行處理，可以得出導(dǎo)彈振動(dòng)模態(tài)的全彈道變化趨勢(shì)，為導(dǎo)彈穩(wěn)定性分析和凹陷濾波器的設(shè)計(jì)提供更加精確的參考。