利用微分算子增強(qiáng)EMD算法頻帶分解能力

劉海波，玄志武

(91550部隊(duì)94分隊(duì)，遼寧 大連 116023)

Hilbert-Huang變換［1］為對非平穩(wěn)、非線性信號進(jìn)行分析的有效方法，且具有自適應(yīng)特點(diǎn)，能揭示非平穩(wěn)信號的局部特性，已廣泛應(yīng)用于機(jī)械故障檢測［2-3］、結(jié)構(gòu)識別［4］、外彈道去噪［5］等領(lǐng)域。但由于 EMD 分解由實(shí)踐中總結(jié)而來，理論基礎(chǔ)尚不完善。實(shí)際應(yīng)用中難以避免出現(xiàn)模態(tài)混疊現(xiàn)象，即一個(gè)IMF分量中包含不同時(shí)間尺度組分情況，削弱了EMD算法的頻帶分解能力，導(dǎo)致分解出的本征模態(tài)函數(shù)失去原有物理意義。

Li等［6］用小波對信號進(jìn)行預(yù)處理，濾除高頻間斷事件，一定程度上改善了分解結(jié)果，但小波基的選擇需人工干預(yù)，削弱了EMD最突出的自適應(yīng)性特點(diǎn)。Ryan等［7］在EMD分解前加入掩膜信號，對區(qū)分頻率較接近的頻帶有一定效果，但選取掩膜信號需信號的先驗(yàn)信息，影響EMD算法的自適應(yīng)性。

本文提出改進(jìn)的EMD算法，先將提高原始信號頻帶分解能力的微分信號進(jìn)行EMD分解，再將所得IMFs進(jìn)行積分處理獲得原始信號有效IMFs。仿真計(jì)算與實(shí)際工程應(yīng)用表明，該方法優(yōu)于傳統(tǒng)的EMD算法，可有效分離頻率較接近的本征模態(tài)函數(shù)，增強(qiáng)EMD方法的頻帶分解能力。

1 EMD方法頻帶分解能力

Rilling等［8］對EMD算法的信號頻帶分解能力給出初步解答。Feldman［9］從理論層面進(jìn)行深入研究得出，假設(shè)兩不同頻帶信號為s1，s2，其振幅與頻率分別為A1，A2與f1，f2，以頻率比f2/f1為橫坐標(biāo)，振幅比A2/A1為縱坐標(biāo)，得EMD算法頻帶分解能力見圖1。圖1中，1區(qū)表示A2/A1≤(f2/f1)-2，或f2/f1≤1.5無法用EMD算法進(jìn)行分解;2 區(qū)表示(f2/f1)-2≤A2/A1≤2.4(f2/f1)-1.75，需幾次平移迭代才能將信號分解;3區(qū)表示A2/A1≥2.4(f2/f1)-1.75，一次迭代即能將信號分解。

由上述得，EMD算法的頻帶分解能力與信號振幅比A2/A1及頻率比f2/f1有關(guān)。若能提高信號的振幅比，使其從1區(qū)進(jìn)入2區(qū)或3區(qū)，則可用EMD進(jìn)行頻帶分解。

圖1 EMD算法頻帶分解能力Fig.1 Frequency decomposition ability of EMD

2 利用微分算子增強(qiáng)EMD頻帶分解能力

設(shè)輸入信號x(t)由兩種模態(tài)構(gòu)成，即:

其中:f2＞f1，A2＜A1。信號振幅比y1可表示為:

對式(1)兩邊進(jìn)行一階微分得:

其振幅比為:

由于f2＞f1故y2＞y1，即通過微分增大信號振幅比，可將1區(qū)部分信號“提升”到2區(qū)或3區(qū)(如圖1中箭頭所示)，增強(qiáng)應(yīng)用EMD算法進(jìn)行頻帶分解能力，圖1演化為圖2。

圖2 經(jīng)微分算子處理的EMD頻帶分解能力Fig.2 Frequency decomposition ability of EMD after differentiation

圖2中，① 區(qū)表示A2/A1≤(f2/f1)-3，或f2/f1≤1.5無法用微分算子處理后EMD算法將其分解;② 區(qū)表示(f2/f1)-3≤A2/A1≤2.4(f2/f1)-2.75，需幾次平移迭代，才能用微分算子處理后的EMD將信號分解;③區(qū)表示A2/A1≥2.4(f2/f1)-2.75，一次迭代即能用微分算子處理后EMD算法將信號分解。

圖2與圖1相比，①區(qū)部分顯著變小，說明EMD分解的盲區(qū)減小，即EMD的頻帶分解能力得以增強(qiáng)。

3 頻帶分解能力判定指標(biāo)

為定量分析EMD算法的頻帶分解能力，引入指標(biāo)Error，表示原始信號成分與分解結(jié)果間在時(shí)域上的誤差:

其中:IMFi(t)表示第i階固有模態(tài)函數(shù)，n為EMD分解階數(shù)。

對同一待分解信號，采用不同分解模式，當(dāng)Error趨于0時(shí)，說明本征函數(shù)IMFi(t)與信號si(t)較接近，頻帶分解正確。

4 改進(jìn)的EMD算法

實(shí)際工程應(yīng)用中，可用差分方法代替微分過程，用累積求和方法代替積分過程。再次用EMD方法提取出累積求和后信號的高頻部分作為新的IMF，避免由替代方法所致誤差，而使累積求和后信號包絡(luò)均值不為零的情況發(fā)生。

圖3 算法流程圖Fig.3 Algorithm flowchart

5 仿真分析

正弦合成信號為:

采樣個(gè)數(shù) 210，圖 4(a)為s(t)高頻部分 0.2 sin(4πt)，圖4(b)為s(t)低頻部分 sin(2πt)。s(t)振幅比A2/A1=0.2，頻率比f2/f1=2，因此A2/A1≤(f2/f1)-2，即圖1中信號s(t)處于1區(qū)，無法用EMD算法分離頻帶，已由圖 5結(jié)果得以證實(shí)，此時(shí) Error為24.88。

對s(t)微分得:

工程應(yīng)用中，信號往往含噪聲信息，直接進(jìn)行微分處理會擴(kuò)大噪聲，影響信號分解效果。當(dāng)信噪比WS-NR為8.12時(shí)，信號s(t)經(jīng)微分再經(jīng)EMD分解的結(jié)果見圖8。對圖8各頻帶積分結(jié)果見圖9。圖8第三、四頻帶與圖6相似;圖9第三、四頻帶與圖7相似但存差異，說明信號s(t)中加入噪聲會對改進(jìn)的EMD分解算法產(chǎn)生影響。

為進(jìn)一步考察疊加噪聲的強(qiáng)度對算法影響，在合成信號s(t)中疊加不同強(qiáng)度白噪聲，利用本文改進(jìn)的EMD算法進(jìn)行處理。結(jié)果見表1，其中 WSNR表示s(t)疊加白噪聲后的信噪比。信噪比由13.64降到2.89低4/5，而 Error只增加2.4倍，說明該算法對噪聲適應(yīng)性較好。

表1 疊加不同強(qiáng)度噪聲對算法結(jié)果影響Tab.1 Different noise strength influence on result

圖4 信號高低頻部分Fig.4 High and low frequency of the signal

圖5 直接EMD分解Fig.5 EMD decomposition directly

圖6 EMD微分后Fig.6 EMD decomposition after differentiation

圖7 積分Fig.7 Integration

圖8 加入噪聲的信號經(jīng)微分再經(jīng)EMD分解Fig.8 EMD after noising signal differentiation

圖9 圖8頻帶積分Fig.9 Frequencies Integration on figure

6 工程應(yīng)用

由于受日、月等天體影響，地球自轉(zhuǎn)不規(guī)律導(dǎo)致地球定向參數(shù)存在多種短、長周期變化。以分析地球定向參數(shù)測量數(shù)據(jù)［10］頻帶周期性規(guī)律為例，檢驗(yàn)本文算法的有效性。

圖10 日長變化數(shù)據(jù)前4 096點(diǎn)Fig.10 First 4 096 points of LOD data

地球定向參數(shù)space2000_midnight.eop中的 LOD(length of day:日長變化)數(shù)據(jù)，記錄了從1976年9月28日至2001年1月6日共8 867個(gè)數(shù)據(jù)點(diǎn)。取前4 096點(diǎn)(圖10)進(jìn)行EMD分解，結(jié)果見圖11。IFM1～I(xiàn)FM7每階固有模態(tài)函數(shù)均有不同振蕩周期，代表歲差、章動(dòng)、極移等因素對LOD影響。

圖11 LOD數(shù)據(jù)EMD分解結(jié)果Fig.11 EMD Result of LOD data

圖12 用本文算法重新分解LOD數(shù)據(jù)Fig.12 Decomposition again using new method

由圖11虛線方框看出，IFM5中至少包含兩種尺度分量，并出現(xiàn)模態(tài)混疊現(xiàn)象，使IFM4亦受到影響。用本文改進(jìn)的EMD算法重新對LOD數(shù)據(jù)進(jìn)行分解，結(jié)果見圖12。與圖11虛線框內(nèi)對比看出，該算法將原混疊模態(tài)成功分離，使各階IMF周期規(guī)律更明顯:IFM1為半月周期，IFM2為整月周期，IFM3為季度周期，IFM4為半年周期，IFM5為整年周期等。

7 結(jié)論

由仿真分析與實(shí)際工程應(yīng)用得出，本文所提算法可增強(qiáng)EMD頻帶分解能力。當(dāng)滿足(f2/f1)-3≤A2/A1，且f2/f1＞1.5條件時(shí)，可成功對混疊模態(tài)進(jì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，1998，454(1971):903-995.

［2］沈 路，楊富春，周曉軍，等.基于改進(jìn) EMD與形態(tài)濾波的齒輪故障特征提?。跩］.振動(dòng)與沖擊，2010，29(3):154-157.

SHEN Lu，YANG Fu-chun，ZHOU Xiao-jun，et al.Gear faultfeature extraction based on improved EMD and morphological filter［J］.Journal of Vibration and Shock，2010，29(3):154-157.

［3］李 慧，劉小峰，夏雨峰.基于諧波小波包變換的齒輪箱包絡(luò)解調(diào)分析［J］.振動(dòng)與沖擊，2012，31(12):129-134.

LI Hui，LIU Xiao-feng，XIA Yu-feng.A method of modal parameters identification using ambient vibration measurements for a large-scale bridge［J］. Journal of Vibration and Shock，2012，31(12):129-134.

［4］付 春，姜紹飛，杜 權(quán).基于改進(jìn)EMD的結(jié)構(gòu)模態(tài)參數(shù)識別方法［J］.武漢理工大學(xué)學(xué)報(bào)，2010，32(9):280-285.

FU Chun，JIANG Shao-fei，DU Quan.Structural modal parameters identification method based on improved EMD［J］.Journal of Wuhan University of Technology，2010，32(9):280-285.

［5］郭小紅，徐小輝，趙樹強(qiáng).基于經(jīng)驗(yàn)?zāi)B(tài)分解的外彈道降噪方法及應(yīng)用［J］.宇航學(xué)報(bào)，2008，29(4):1272-1275.

GUO Xiao-hong，XU Xiao-hui，ZHAO Shu-qiang.Noising reduction algorithm based on empirical mode decomposition(EMD)and applicationsin calculating trajectoriesof spacecraft［J］.Journal of Astronautics，2008，29(4):1272-1275.

［6］Li H，Yang L，Huang D.The study of the intermittency test filtering character of Hilbert-Huang transform ［J］.Mathematics and Computers in Simulation，2005，70:22-32.

［7］ Deering R，Kaiser J F.The use of a masking signal to improve empirical mode decomposition［C］.Philadelphia，PA，USA:Proc. of IEEE International Conference on Acoustics，Speech，and Signal Processing，2005.

［8］ Rilling G，F(xiàn)landrin P.One or two frequency the empirical mode decomposition answers［J］.IEEE Transaction on Signal Processing，2008，56(1):85-95.

［9］Feldman M.Analytical basis of the EMD:two harmonics decomposition［J］. Mechanical Systems and Signals Processing，2009，23(7):2059-2071.

［10］ Gross R S.Combinations of earth orientation measurements:SPACE2000，COMB2000，and POLE2000［M］.California:JPL Publication，2001.