一種EMD改進(jìn)方法及其在旋轉(zhuǎn)機(jī)械故障診斷中的應(yīng)用

時培明，丁雪娟，李 庚，韓東穎

(1.燕山大學(xué) 電氣工程學(xué)院，秦皇島 066004;2.燕山大學(xué) 車輛與能源學(xué)院 ，秦皇島 066004)

旋轉(zhuǎn)機(jī)械故障診斷中，常常遇到諸如瞬變、調(diào)幅或調(diào)頻等非平穩(wěn)、非線性信號，提取和分析這些信號中的特征信息是機(jī)械診斷的關(guān)鍵。EMD是Huang首先提出的一種新型的時頻分析方法［1］，能將復(fù)雜信號分解為有限個內(nèi)稟模態(tài)函數(shù)(Intrinsic Mode Function，IMF)，適合處理非線性、非平穩(wěn)信號。自EMD方法問世以來，它就引起了眾多學(xué)者的廣泛關(guān)注，已被應(yīng)用于結(jié)構(gòu)分析、設(shè)備診斷等領(lǐng)域［2-3］。

EMD方法的分析質(zhì)量很大程度上取決于EMD分解的質(zhì)量。而EMD分解由于采用三次樣條插值來獲取信號的瞬時平均，使得這種方法存在特殊的端點效應(yīng)。另外，在Hilbert變換的過程中也存在端點效應(yīng)，嚴(yán)重時會影響整個信號的分析結(jié)果。針對這一問題，研究人員已經(jīng)提出了一些改進(jìn)方法，如鏡像延拓［4］、神經(jīng)網(wǎng)絡(luò)延拓［5］、相似極值延拓［6］、波形特征匹配延拓［7］、支持矢量回歸機(jī)［8］等。但是，這些延拓法存在一個共同的問題，就是在延拓后，信號的端點仍然不確定，這樣隨著分解過程的進(jìn)行，包絡(luò)線兩端仍可能存在發(fā)散現(xiàn)象，并逐漸向內(nèi)“污染”，導(dǎo)致其端點效應(yīng)問題依然存在。

本文提出一種波形特征匹配延拓和余弦窗函數(shù)結(jié)合的EMD端點效應(yīng)處理方法。該方法在保留了波形匹配延拓方法優(yōu)勢的基礎(chǔ)上，通過窗函數(shù)處理，使延拓部分逐漸減小直到歸零，從而使延拓誤差減小［9-10］。通過仿真分析和不對中故障診斷實例，驗證了該方法的有效性。

1 改進(jìn)的EMD時頻分析方法

1.1 EMD分解方法

EMD時頻分析方法適于處理非線性、非平穩(wěn)信號。該方法包括兩個過程:經(jīng)驗?zāi)B(tài)分解(EMD)和Hilbert變換。

EMD方法的分解步驟為:

(1)確定信號所有的局部極值點。

(2)利用三次樣條線將所有的局部極大值點連接起來形成上包絡(luò)線;將所有的局部極小值點連接起來形成下包絡(luò)線。上、下包絡(luò)線應(yīng)該包絡(luò)所有的數(shù)據(jù)點。

(3)上、下包絡(luò)線的平均值記為m1(t)。

上述分解過程如圖1所示，圖(a)為原始信號，圖(b)中·表示極大值點，?為極小值點。圖(c)中a為上包絡(luò)線，b為下包絡(luò)線，c為上、下包絡(luò)線的均值。

圖1 EMD分解示意圖Fig.1 Diagram of the EMD decomposition

(4)信號x(t)和m1(t)的差值為第一個分量，h1(t)為:

如果h1(t)是一個IMF，那么h1(t)就是第一個 IMF分量。

(5)如果h1(t)不滿足IMF條件，把h1(t)作為原始數(shù)據(jù)，重復(fù)上述步驟，得到:

式中，m11(t)是h1(t)上下包絡(luò)的平均值。反復(fù)篩選k次后，使得h1k(t)為IMF分量，即:

令:

從原始信號中獲得的第一個IMF分量c1(t)應(yīng)該包含信號最好的范圍或者最短的周期成分。

(6)從x(t)分離出c1(t)，我們得到:

將r1(t)看作原始數(shù)據(jù)重復(fù)以上步驟，得到x(t)的第2個IMF分量c2(t)。重復(fù)循環(huán)n次，得到信號x(t)的n個IMF分量。這樣就有:

當(dāng)rn(t)成為一個單調(diào)函數(shù)不能再從中提取出滿足IMF條件的分量時，循環(huán)結(jié)束。這樣由式(5)和式(6)得到:

因此可以將信號分解為n個經(jīng)驗?zāi)B(tài)，殘余函數(shù)rn(t)代表信號的平均趨勢。從而對每一個IMF進(jìn)行Hilbert變換，得到Hilbert譜及其邊際譜。

對式(7)中的每一個內(nèi)稟模態(tài)函數(shù)ci(t)作Hilbert變換得到:

構(gòu)造解析信號:

于是得到幅值函數(shù):

和相位函數(shù):

進(jìn)一步可以求出瞬時頻率:

這樣，可以得到:

這里省略了殘余分量rn，RP表示取實部。展開式(13)稱為Hilbert譜，記作:

再定義Hilbert邊際譜:

式中，T為信號的總長度。

1.2 EMD分解中的端點效應(yīng)

EMD分解中由于無法保證端點處的極值點，導(dǎo)致求包絡(luò)平均過程中，會在樣條插值時產(chǎn)生數(shù)據(jù)的擬合誤差。并且隨著分解的進(jìn)行，誤差不停積累，由端點處向內(nèi)逐漸傳播，嚴(yán)重時會使分解的數(shù)據(jù)失去意義。

圖2為仿真信號x(t)的EMD分解結(jié)果，信號x(t)由兩個正弦信號和一個調(diào)幅信號組成，其表達(dá)式為

分解時沒有對原始數(shù)據(jù)進(jìn)行任何處理。由圖2所示，信號x(t)分解出來的3階IMF都有嚴(yán)重的端點效應(yīng)。

圖2 未經(jīng)過處理的EMD分解結(jié)果Fig.2 The decomposition result of simulative signal without being processed

1.3 基于波形特征匹配延拓與余弦窗函數(shù)的EMD方法改進(jìn)

基于波形特征匹配的數(shù)據(jù)延拓方法是通過采用信號內(nèi)部和邊緣處變化趨勢最為相似的子波來對端點處數(shù)據(jù)進(jìn)行延拓，它是一種自適應(yīng)的方法。在具體實現(xiàn)中，通過計算波形匹配來量化兩端波形的變化趨勢。

如圖3所示，以左邊界第一個極值點為極大值為例，Mi、Ni(n=1，2，3，…)分別為曲線的極大值、極小值點，分別對應(yīng)時間tmi、tni，S1為左端點，波形特征匹配延拓法以S1-M1-N1為邊界特征波形，在全部數(shù)據(jù)中找到與S1-M1-N1構(gòu)成的三角形最接近的波形為匹配波形，如Si-Mi-Ni，從Si的前一點(右邊界為后一點)數(shù)據(jù)開始，向前(右邊界向后)延拓波形數(shù)據(jù)，使延拓數(shù)據(jù)符合信號的自然走向。

圖3 波形特征匹配原理圖Fig.3 Diagram of wave characteristic matching

具體步驟為:

(1)根據(jù)S1-M1-N1邊界三角形的時間坐標(biāo)，尋找與S1對應(yīng)的點Si，其對應(yīng)的時間由式(17)求得:

圖4 信號x(t)的波形匹配延拓結(jié)果Fig.4 Results of wave characteristic matching

由于得到的tsi不一定在采樣點上，采用線性插值求其精確值Si。

(2)計算匹配誤差Ei

式中，最后一項為匹配波形的趨勢項，反映特征波形在曲線中相對極值點的位置。

(3)取Ei最小的波形為匹配波形，若存在多個相等的Ei，為獲得足夠的延拓數(shù)據(jù)供選擇，取與起始點距離最遠(yuǎn)的波形為匹配波。

(4)自匹配波形與Si距離最近的數(shù)據(jù)點的前一點開始，將實際波形數(shù)據(jù)延拓到S1前，延拓點數(shù)據(jù)根據(jù)需要選擇，若信號中Si前數(shù)據(jù)點個數(shù)小于需要延拓的點數(shù)，可反復(fù)延拓此段波形。

(5)采用同樣的原理可以對右邊界數(shù)據(jù)進(jìn)行延拓。

(6)采用延拓后的數(shù)據(jù)完成EMD分解。信號包絡(luò)插值計算采用延拓后的全部數(shù)據(jù)，而在EMD分解結(jié)束條件的計算中，仍采用有效數(shù)據(jù)。

(7)對EMD分解得到的IMF作Hilbert變換，此時數(shù)據(jù)長度為延拓后數(shù)據(jù)總長，按原始信號長度及位置取分析結(jié)果，得到有效數(shù)據(jù)的EMD分析結(jié)果。

在上述延拓過程中兼顧了數(shù)據(jù)的極值點及非極值點的波形數(shù)據(jù)，使得延拓數(shù)據(jù)與原信號交界處光滑過渡，避免了邊界處瞬時頻率的跳躍。

圖4是信號x(t)基于波形匹配延拓的結(jié)果。由圖可以看出延拓部分雖然與實際信號比較接近，但還是存在誤差。對延拓數(shù)據(jù)加余弦窗函數(shù)處理，控制端點效應(yīng)向內(nèi)“污染”，得到更準(zhǔn)確的IMF。

定義余弦窗函數(shù):

式中，L為信號延拓后的長度，A為信號兩端延拓中較長的延拓長度。

圖5 余弦窗函數(shù)Fig.5 The shape of cosine window

圖6 信號延拓加窗處理結(jié)果Fig.6 The results processed with the proposed method

余弦窗的形狀如圖5所示。

余弦函數(shù)窗將延拓誤差“控制”在信號兩端，使其無法(或以較慢速度)向數(shù)據(jù)內(nèi)部發(fā)展，保證信號中部數(shù)據(jù)的正確分解。首先，將延拓信號x(t)與余弦窗函數(shù)進(jìn)行內(nèi)積運算，得到信號y(t)=〈x(t)，ω(t)〉。

然后，對處理后的信號y(t)進(jìn)行EMD分解，再將分解得出的IMF的兩端去掉相應(yīng)的延拓部分A。最后對減去延拓部分后的IMF進(jìn)行邊界譜分析。該方法既考慮到了延拓誤差又考慮到了信號的完整性。

圖6為信號x(t)延拓后再加余弦窗處理的結(jié)果。從圖中可以看出延拓的部分逐漸減小直到歸零，從而使延拓部分誤差減小，為得到更準(zhǔn)確的IMF提供了可靠條件。

2 仿真分析

對經(jīng)過處理的信號x(t)進(jìn)行EMD分解得到3階IMF，并與實際信號作比較如圖7。通過與圖2對比發(fā)現(xiàn)經(jīng)過延拓加窗處理得到的IMF明顯更符合實際值，說明該方法對抑制端點效應(yīng)有良好的效果。

圖7 延拓加窗后的分解結(jié)果Fig.7 The decomposition result of simulative signal with the proposed method

圖8、圖9分別是信號x(t)未經(jīng)過處理和用改進(jìn)方法得到的Hilbert譜與邊界譜。由圖8可以看出，未處理的Hilbert譜在信號兩端有比較嚴(yán)重的發(fā)散現(xiàn)象，而通過本文方法處理后的Hilbert譜效果有明顯的改善。圖9中(a)的幅值顯示明顯有很多微弱振蕩，尤其是在高頻部分尤為明顯。通過端點效應(yīng)處理后得到的邊界譜(b)中，幅值振蕩基本消除，得到了更好的處理效果。

3 EMD改進(jìn)方法在不對中故障診斷中的應(yīng)用研究

圖8 Hilbert時頻譜對比圖Fig.8 The composition of Hilbert spectrum

圖9 邊際譜對比圖Fig.9 The composition of marginal spectrum

將本文提出的EMD改進(jìn)方法應(yīng)用于轉(zhuǎn)子不對中故障信號的特征提取及診斷。圖10是采集信號的時域波形，轉(zhuǎn)速為900 r/min，采樣頻率為768 Hz。圖11是故障信號的EMD分解結(jié)果，圖12、圖13分別是改進(jìn)前后信號的Hilbert譜及邊際譜。

圖10 不對中故障的時域波形Fig.10 Time domain waveform of misalignment fault

圖11 故障信號的EMD分解結(jié)果Fig.11 The decomposition result of fault signal

圖11中的內(nèi)稟模態(tài)函數(shù)IMF1，IMF2，IMF3被依次分解出來，分別對應(yīng)著多倍頻、2倍頻和基頻振動模態(tài)，但是因為噪聲的存在，一定程度地影響了EMD的分解精度［11-12］。由Hilbert時頻譜圖可觀察到，基頻分量和2倍頻左右的分量在分析的時間內(nèi)一直穩(wěn)定存在，除此之外還有多倍頻分量，但并不穩(wěn)定;比較圖12，改進(jìn)后的Hilbert譜的端點發(fā)散現(xiàn)象得到了明顯的改善，尤其是信號的左端更為明顯。在Hilbert邊際譜的頻幅譜中觀察到，除了基頻還有其他倍頻信息存在，基頻分量和2倍頻分量占主導(dǎo)地位，且2倍頻分量強(qiáng)度并未超過基頻分量;比較圖13，改進(jìn)后的邊際譜幅值微弱振蕩基本消除，且2倍頻分量得到了突出。根據(jù)以上分析診斷該轉(zhuǎn)子的故障為不對中故障。

圖12 故障信號的Hilbert譜Fig.12 The Hilbert spectrum of fault signal

圖13 故障信號的Hilbert邊際譜Fig.13 The marginal spectrum of fault signal

4 結(jié)論

(1)提出了一種抑制端點效應(yīng)的新方法，首先利用波形特征匹配延拓對信號兩端進(jìn)行延拓，然后根據(jù)延拓情況對信號加余弦窗進(jìn)行處理，把延拓誤差控制在兩邊。最后在EMD分解后，只取信號的有效長度，提高了分解的精度。

(2)通過仿真信號的分析，證明了該方法能有效地抑制端點效應(yīng)，為得到準(zhǔn)確地邊際譜和Hilbert譜提供了保障。

(3)將改進(jìn)方法應(yīng)用到旋轉(zhuǎn)機(jī)械故障診斷中，通過對含有不對中故障的信號進(jìn)行分析，證明了該方法能從非線性故障信號中得到真實有用的故障信息，實現(xiàn)旋轉(zhuǎn)機(jī)械故障的有效診斷。

［1］Huang N E.The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis［J］. Proc. R. Soc. Lond. A，1998，454:903-995.

［2］ Cheng J S，Yu D J，Tang J S，et al.Local rub-impact fault dianosis of the rotor systems based on EMD［J］.Mechanism and Machine Theory，2009，44:784-791.

［3］Yang B，Suh C S.Interpretation of crack-induced rotor nonlinear response using instantaneous frequency［J］.Mechanical System and Signal Processing，2004，18(3):491-513.

［4］韓建平，錢 炯，董小軍.采用鏡像延拓和RBF神經(jīng)網(wǎng)絡(luò)處理EMD中端點效應(yīng)［J］.振動、測試與診斷，2010，30(4):414-417.

HAN Jian-ping，QIAN Jiong，DONG Xiao-jun.Suppression of end-effect in empirical mode decomposition by mirror extension and radial basis function neural network prediction［J］.Journal of Vibration，Measurement＆Diagnosis，2010，30(4):414-417.

［5］胡勁松，楊世錫.EMD方法基于徑向基神經(jīng)網(wǎng)絡(luò)預(yù)測的數(shù)據(jù)延拓與應(yīng)用［J］.機(jī)械強(qiáng)度，2007，19(6):894-899.

HU Jin-song，YANG Shi-xi.Application of EMD method with data extension technique based on RBF neural network to time-frequency analysis［J］.Journal of Mechanical Strength，2007，19(6):894-899.

［6］沈 路，周曉軍，張志剛，等.Hilbert-Huang變換中的一種端點延拓方法［J］.振動與沖擊，2009，28(8):168-174.

SHEN Lu， ZHOU Xiao-jun， ZHANG Zhi-gang， et al.Boundary-extension method in Hilbert-Huang transform［J］.Journal of Vibration and Shock，2009，28(8):168-174.

［7］胡愛軍，安連鎖，唐貴基.HILBERT-HUANG變換端點效應(yīng)處理新方法［J］.機(jī)械工程學(xué)報，2008，44(4):154-158.

HU Ai-jun， AN Lian-suo， TANG Gui-ji. New process method for end effects of Hilbert-Huang transform［J］.Journal of Mechanical Engineering，2008，44(4):154-158.

［8］程軍圣，于德介，楊 宇.基于支持矢量回歸機(jī)的Hilbert-Huang變換端點效應(yīng)問題的處理方法［J］.機(jī)械工程學(xué)報，2006，42(4):23-31.

CHENG Jun-sheng，YU De-jie，YANG Yu.Process method for end effects of Hilbert-Huang transform based on support vector regression machines［J］.JournalofMechanical Engineering，2006，42(4):23-31.

［9］Qi K Y，He Z J，Zi Y Y.Cosine window-based boundary processing method for EMD and its application in rubbing faultdiagnosis［J］. MechanicalSystems and Signal Processing，2007，21:2750-2760.

［10］于慧君，陳章位，王慶豐.一種加窗重疊信號平滑連接方法及其在振動信號預(yù)處理中的應(yīng)用［J］.振動與沖擊，2007，26(8):39-40.

YU Hui-jun，CHEN Zhang-wei，WANG Qing-feng.OLA signal smooth linking method and its application in vibration signal pre-processing［J］.Journal of Vibration and Shock，2007，26(8):39-40.

［11］ Jalan A K，Mohanty A R.Model based fault diagnosis of a rotor-bearing system for misalignment and unbalance under steady-state condition［J］.Journal of Sound and Vibration，2009，327(3-5):604-622.

［12］ Rybczyński J.The possibility of evaluating turbo-set bearing misalignment defects on the basis of bearing trajectory features［J］.Mechanical Systems and Signal Processing，2011，25(2):521-536.