謝金松，薛林，高慶豐

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3 導(dǎo)彈簡化歐拉梁模型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ā)動機開機的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)驗公式：

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

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

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

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

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

圖4 彈上加表量測到的振動信號Fig.4 Vibration signal measured by accelerometer

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

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

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

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

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

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

由圖6 可以看出，極值鏡像延拓有效抑制了各階IMF 左側(cè)的端點效應(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 可以看出，對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作為精度評價指標(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)混疊和端點效應(yīng)，EEMD 雖然提高了辨識精度，但仍無法完全抑制模態(tài)混疊，改進(jìn)后的EEMD 篩分精度最高。

3.3 振動參數(shù)辨識

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

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