時(shí)培明，丁雪娟，李 庚，韓東穎

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

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

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

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

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

1.1 EMD分解方法

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

EMD方法的分解步驟為:

(1)確定信號(hào)所有的局部極值點(diǎn)。

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

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

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

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

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

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

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

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

令:

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

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

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

當(dāng)rn(t)成為一個(gè)單調(diào)函數(shù)不能再?gòu)闹刑崛〕鰸M足IMF條件的分量時(shí)，循環(huán)結(jié)束。這樣由式(5)和式(6)得到:

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

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

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

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

和相位函數(shù):

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

這樣，可以得到:

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

再定義Hilbert邊際譜:

式中，T為信號(hào)的總長(zhǎng)度。

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

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

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

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

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

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

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

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

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

具體步驟為:

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

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

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

(2)計(jì)算匹配誤差Ei

式中，最后一項(xiàng)為匹配波形的趨勢(shì)項(xiàng)，反映特征波形在曲線中相對(duì)極值點(diǎn)的位置。

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

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

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

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

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

在上述延拓過(guò)程中兼顧了數(shù)據(jù)的極值點(diǎn)及非極值點(diǎn)的波形數(shù)據(jù)，使得延拓?cái)?shù)據(jù)與原信號(hào)交界處光滑過(guò)渡，避免了邊界處瞬時(shí)頻率的跳躍。

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

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

式中，L為信號(hào)延拓后的長(zhǎng)度，A為信號(hào)兩端延拓中較長(zhǎng)的延拓長(zhǎng)度。

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

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

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

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

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

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

2 仿真分析

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

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

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

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

圖8 Hilbert時(shí)頻譜對(duì)比圖Fig.8 The composition of Hilbert spectrum

圖9 邊際譜對(duì)比圖Fig.9 The composition of marginal spectrum

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

圖10 不對(duì)中故障的時(shí)域波形Fig.10 Time domain waveform of misalignment fault

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

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

圖12 故障信號(hào)的Hilbert譜Fig.12 The Hilbert spectrum of fault signal

圖13 故障信號(hào)的Hilbert邊際譜Fig.13 The marginal spectrum of fault signal

4 結(jié)論

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

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

(3)將改進(jìn)方法應(yīng)用到旋轉(zhuǎn)機(jī)械故障診斷中，通過(guò)對(duì)含有不對(duì)中故障的信號(hào)進(jìn)行分析，證明了該方法能從非線性故障信號(hào)中得到真實(shí)有用的故障信息，實(shí)現(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中端點(diǎn)效應(yīng)［J］.振動(dòng)、測(cè)試與診斷，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ù)測(cè)的數(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變換中的一種端點(diǎn)延拓方法［J］.振動(dòng)與沖擊，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變換端點(diǎn)效應(yīng)處理新方法［J］.機(jī)械工程學(xué)報(bào)，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變換端點(diǎn)效應(yīng)問(wèn)題的處理方法［J］.機(jī)械工程學(xué)報(bào)，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］于慧君，陳章位，王慶豐.一種加窗重疊信號(hào)平滑連接方法及其在振動(dòng)信號(hào)預(yù)處理中的應(yīng)用［J］.振動(dòng)與沖擊，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.