基于EMD 和分形盒維數(shù)的旋轉(zhuǎn)機(jī)械耦合故障診斷方法研究

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

(1.燕山大學(xué) 車輛與能源學(xué)院，秦皇島 066004;2.燕山大學(xué) 河北省測試計(jì)量技術(shù)及儀器重點(diǎn)實(shí)驗(yàn)室，秦皇島 066004)

旋轉(zhuǎn)機(jī)械轉(zhuǎn)子系統(tǒng)中轉(zhuǎn)軸裂紋、動(dòng)靜碰摩和支撐部件松動(dòng)為常見故障，是導(dǎo)致機(jī)械系統(tǒng)失效甚至造成嚴(yán)重事故的主要原因。目前國內(nèi)外學(xué)者［1－3］對裂紋、碰摩、松動(dòng)故障單獨(dú)存在的系統(tǒng)研究較多，提出了許多有效的診斷方法。

但在實(shí)際轉(zhuǎn)子中，常常出現(xiàn)兩種甚至三種故障同時(shí)存在的情況。這種耦合故障轉(zhuǎn)子的動(dòng)力學(xué)行為較單一故障轉(zhuǎn)子更加復(fù)雜，而且相互影響，不容易診斷。文獻(xiàn)［4］采用有限元方法分析了裂紋碰摩耦合故障轉(zhuǎn)子系統(tǒng)的響應(yīng)，在嚴(yán)重碰摩故障信號(hào)中診斷出了裂紋故障。文獻(xiàn)［5］建立了帶有裂紋－支承松動(dòng)耦合故障的雙跨彈性轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)模型，利用數(shù)值仿真對故障非線性響應(yīng)進(jìn)行研究。文獻(xiàn)［6］通過分岔圖、Poincaré映射圖和幅值譜圖，分析了碰摩－松動(dòng)故障機(jī)組軸系隨轉(zhuǎn)頻比和轉(zhuǎn)子質(zhì)量偏心變化的非線性動(dòng)態(tài)響應(yīng)。

EMD是近年來發(fā)展起來的處理非平穩(wěn)、非線性信號(hào)的時(shí)頻分析方法［7］。該方法克服了傳統(tǒng)時(shí)頻分析方法中的不足，具有很強(qiáng)的自適應(yīng)性，并在機(jī)械故障診斷領(lǐng)域得到了廣泛應(yīng)用［8－10］。分形的本意為不規(guī)則的、分?jǐn)?shù)的、支離破碎的物體，可以理解為局部和整體在某個(gè)方面有相似性［11］。分形維數(shù)是度量分形的重要指標(biāo)，有許多種分形維數(shù)可用來刻畫信號(hào)的復(fù)雜度［12］，如Hausdorff維、盒子維、容量維、信息維、關(guān)聯(lián)維數(shù)等，其中盒維數(shù)由于其計(jì)算簡單、高效而應(yīng)用最為廣泛［13－15］。

針對耦合故障信號(hào)復(fù)雜，具有強(qiáng)非線性的特點(diǎn)，本文提出一種基于EMD和分形盒維數(shù)的耦合故障診斷方法。該方法先利用EMD將故障信號(hào)分解，然后對有效IMF進(jìn)行盒維數(shù)計(jì)算，得到耦合故障信號(hào)的故障特征向量，判斷故障狀態(tài)，達(dá)到故障診斷的目的。

1 系統(tǒng)動(dòng)力學(xué)模型和運(yùn)動(dòng)微分方程

如圖1所示，建立含有裂紋－碰摩－松動(dòng)耦合故障的剛性支承轉(zhuǎn)子－軸承系統(tǒng)模型，轉(zhuǎn)子圓盤與軸承之間為無質(zhì)量的彈性軸。模型左端發(fā)生松動(dòng)，軸承座與基礎(chǔ)之間的松動(dòng)最大間隙為ξ1。轉(zhuǎn)子圓盤左側(cè)有一弓形橫向裂紋，其深度為a，并且圓盤受到碰摩故障影響，其碰摩函數(shù)為P。圖中O1為軸瓦幾何中心，O2為轉(zhuǎn)子幾何中心，O3為轉(zhuǎn)子質(zhì)心，kc為定子剛度，k為彈性軸剛度，m1為兩端軸承處的轉(zhuǎn)子集中質(zhì)量，m2為轉(zhuǎn)子圓盤的等效集中質(zhì)量，m3為軸承支座的等效集中質(zhì)量。模型還考慮了左端滑動(dòng)軸承作用在轉(zhuǎn)軸上的非線性油膜力，為別為 Fx、Fy。

圖1 含裂紋碰摩松動(dòng)故障的轉(zhuǎn)子－軸承系統(tǒng)模型Fig.1 Rotor-bearing system consider in crack rub-impact and pedestal looseness

設(shè)轉(zhuǎn)子右端的徑向位移為x1，y1;轉(zhuǎn)子圓盤的徑向位移為x2，y2;松動(dòng)端軸心位移為x3，y3;軸承支座在豎直方向位移為y4，則具有裂紋碰摩松動(dòng)耦合故障的轉(zhuǎn)子－軸承系統(tǒng)運(yùn)動(dòng)微分方程為:

圖2 裂紋截面模型Fig.2 Crack section model

式中:u為轉(zhuǎn)子的偏心量，c1為轉(zhuǎn)子在軸承處的阻尼系數(shù)，c2為轉(zhuǎn)子圓盤的阻尼系數(shù)，cs為支座松動(dòng)阻尼系數(shù)，ks為支承剛度。ω為轉(zhuǎn)子轉(zhuǎn)速，g為重力加速度，ε、δ為僅與裂紋深度a有關(guān)的相對剛度參數(shù)，由式(3)、(4)所得。如圖2所示，F(xiàn)(Ψ)為裂紋開閉函數(shù)，本文采用余弦波模型來表示裂紋開閉過程，粗略地考慮裂紋半開半閉的過渡過程，忽略了裂紋的全閉和全開是一個(gè)持續(xù)過程。

余弦波模型的數(shù)學(xué)表達(dá)式為:

式中:φ0為初相位，β為裂紋方向與偏心之間的夾角，x，y為轉(zhuǎn)子初始位置。

式(1)中的轉(zhuǎn)定子碰摩函數(shù)P為:

式中:Pn為法向碰摩力，Pτ為切向摩擦力，ξ為靜止時(shí)轉(zhuǎn)子與定子之間的間隙，e為轉(zhuǎn)子軸心位移，由圖 3所示，當(dāng)位移e大于等于間隙 ξ時(shí)系統(tǒng)發(fā)生碰摩。

碰撞發(fā)生時(shí)(e≥ξ)，其正向碰摩力與切向摩擦力分解在x－y坐標(biāo)系可以表示為:

圖3 碰摩力模型Fig.3 Rub － impact model

式中:kc為定子徑向剛度，f為不考慮速度影響時(shí)的摩擦系數(shù)，b為速度影響系數(shù)為轉(zhuǎn)子靜件間的相對滑動(dòng)速度為轉(zhuǎn)子的徑向位移，其中x和y分別對應(yīng)式(1)中的x2和y2。

松動(dòng)故障等效成剛度和阻尼的變化，支承間隙系統(tǒng)在位移條件下ks、cs為分段性，其表達(dá)式為:

式(1)中油膜力沿x和y兩個(gè)方向的分量為:

式中:μ為潤滑油粘度，c為軸承徑向間隙，R為軸承半徑，L為軸承長度。

2 經(jīng)驗(yàn)?zāi)Ｊ椒纸?/h2>

經(jīng)驗(yàn)?zāi)Ｊ椒纸釫MD是一種自適應(yīng)分解方法，可以把復(fù)雜的信號(hào)分解為有限個(gè)IMF分量。IMF信號(hào)一般滿足兩個(gè)條件:① 從全局特性上看，極值點(diǎn)數(shù)必須和過零點(diǎn)數(shù)一致或者至多相差一個(gè);② 在某個(gè)局部點(diǎn)，極大值包絡(luò)和極小值包絡(luò)在該點(diǎn)的算術(shù)平均值是零，即兩條包絡(luò)線關(guān)于時(shí)間軸對稱。

可以把任何信號(hào)x(t)按下面步驟分解:

(1)用三次樣條線將所有的局部極大值點(diǎn)連接起來形成上包絡(luò)線。

(2)用三次樣條線將所有的局部極小值點(diǎn)連接起來形成下包絡(luò)線。

(3)上下包絡(luò)線的平均值記為m1，求出:

理想情況，如果h1是一個(gè)IMF，則h1就是x(t)的第一分量。

(4)如果h1不滿足IMF的條件，把h1作為原始據(jù)，重復(fù)(1)、(2)、(3)，得到上下包絡(luò)線的平均值m1再判斷h11=h1－m11是否滿足IMF的條件，如不滿足，重復(fù)循環(huán) k次，得到 h1k=h1(k－1)－m1k，使得 h1k滿足IMF條件。記c1=h1k，則c1為信號(hào)x(t)的第一個(gè)滿足IMF條件的分量。

(5)將c1從x(t)中分離出來，得到:

將r1作為原始數(shù)據(jù)重復(fù)以上過程，得到x(t)的第二個(gè)滿足IMF條件的分量c2，重復(fù)循環(huán)n次，得到n個(gè)滿足IMF條件的分量。這樣就有:

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

3 分形盒維數(shù)及其計(jì)算

3.1 分形盒維數(shù)

設(shè)X是Rn的非空有界子集，記N表示最大直徑為ε且能覆蓋X集合的最少個(gè)數(shù)，則X的盒維數(shù)定義為:

3.2 盒維數(shù)計(jì)算方法

設(shè)離散信號(hào)y(i)?Y，Y是n維歐氏空間Rn上的閉集。用盡可能細(xì)的ε網(wǎng)格劃分Rn，Nε是集合Y的網(wǎng)格計(jì)數(shù)。由于式(13)中的極限無法按定義求出，所以在計(jì)算時(shí)需要采用近似的方法。以ε網(wǎng)格作為基準(zhǔn)，逐步放大到kε網(wǎng)絡(luò)，其中 k∈Z+。這樣，令 Nkε為離散空間上的集合Y的網(wǎng)格計(jì)數(shù)，則由式(14)和式(15)可以計(jì)算得到。

式中:j=1，2，…，N/k;N 為采樣點(diǎn)數(shù)。k=1，2，…，M，M＜N。

網(wǎng)絡(luò)計(jì)數(shù)Nkε為:

式中:Nkε＞1。

在lgkε－lgNkε圖中確定線性較好的一段為無標(biāo)度區(qū)，設(shè)無標(biāo)度區(qū)的起點(diǎn)和終點(diǎn)分別為k1、k2，則:

最后，用最小二乘法確定該直線的斜率:

盒維數(shù)dB為:

4 數(shù)值仿真和故障診斷

由方程(1)可以看出，含有裂紋碰摩松動(dòng)耦合故障的轉(zhuǎn)子系統(tǒng)是一個(gè)有復(fù)雜外激勵(lì)的非線性系統(tǒng)。目前分析這種系統(tǒng)最有效的方法就是數(shù)值仿真，本文采用變步長四階龍格－庫塔法對方程(1)進(jìn)行數(shù)值求解，系統(tǒng)參數(shù)選取如下:m1=4 kg，m2=32.5 kg，m3=50 kg，R=0.025 m，L=0.012 m，c=0.11 mm，a=0.015 m，μ=0.018 Pa·s，c1=1 050 N·s/m，c2=2 100 N·s/m，cs1=350 N·s/m，cs2=500 N·s/m ，f=0.1，b=0.01，k=2.5 ×107N/m，u=0.05 mm，kc=4.5 ×107N/m，ks1=2.5 ×107N/m，ks2=2.5 ×109N/m，ω =789.3 rad/s，ξ=0.5 mm，ξ1=1 mm，φ0=0，β =0。

考慮到三種耦合故障模型結(jié)構(gòu)特點(diǎn)，轉(zhuǎn)子左端的Y方向位移最能體現(xiàn)耦合故障特性。圖4為數(shù)值解得到的轉(zhuǎn)子左端徑向位移y3的時(shí)域圖。由圖4可以看出由于機(jī)械故障的影響，y3的振動(dòng)有很強(qiáng)的非線性。

圖4 y3時(shí)域圖Fig.4 Time domain chart of y3

EMD方法用于處理非線性、非平穩(wěn)信號(hào)有良好的效果。如圖5所示，把由龍格－庫塔法解出的y3振動(dòng)信號(hào)經(jīng)過EMD方法分解，得到含有故障特征的IMF。由于EMD方法本身原因產(chǎn)生虛假模態(tài)，故只給出IMF1～I(xiàn)MF4。

圖5 y3的EMD處理結(jié)果Fig.5 Empirical mode decomposition(EMD)results of y3

由分解得到的IMF可以看出信號(hào)每隔一段時(shí)間就會(huì)有頻率和幅值上的明顯變化，且這種變化具有周期性。根據(jù)已知的單一故障特征可以判斷這是由于松動(dòng)或者碰摩故障引起的。要對故障進(jìn)一步診斷，傳統(tǒng)方法是對每個(gè)IMF進(jìn)行希爾伯特變換，求出邊界譜，通過邊界譜的倍頻關(guān)系進(jìn)行故障診斷。

圖6是轉(zhuǎn)子左端的Y方向位移y3的邊界譜圖。由圖可以看出振動(dòng)主要是由1/2倍頻組成的，其次是1倍頻，2倍頻等高頻成分則不是主要振動(dòng)頻率。為了更好地進(jìn)行故障診斷，需要把裂紋、松動(dòng)、碰摩兩兩耦合情況下振動(dòng)信號(hào)的邊界譜與三種故障同時(shí)存在的情況作比較。

把三種耦合故障模型進(jìn)行簡化并分析，可以得到兩種耦合故障的邊界譜(如圖7～9)，由于篇幅有限不做詳細(xì)說明。圖7為碰摩松動(dòng)故障信號(hào)的邊界譜，由圖可知信號(hào)振動(dòng)主要由1/2倍頻構(gòu)成，和三種耦合故障情況相似，但1倍頻成分不如圖6明顯，說明1/2倍頻振動(dòng)主要是由于碰摩或松動(dòng)引起的。圖8為裂紋松動(dòng)故障信號(hào)的邊界譜。通過與圖6對比可以明顯看出裂紋松動(dòng)故障信號(hào)不僅含有1/2倍頻、1倍頻成分，也含有3/2倍頻、2倍頻等高頻成分，由此可以診斷出裂紋故障是引起高頻振動(dòng)的主要原因。圖9為裂紋碰摩故障信號(hào)的邊界譜，從圖中可以清晰地看出振動(dòng)信號(hào)同時(shí)具有高頻和低頻振動(dòng)。通過上面的分析可知高頻是有裂紋故障引起的，低頻是有碰摩故障引起的。

通過對四幅邊界譜圖的分析，可以得出裂紋松動(dòng)碰摩三種故障同時(shí)存在時(shí)與其兩兩存在時(shí)故障特征的區(qū)別。但由于圖像抽象、頻率繁雜，不能直觀的把這四種情況區(qū)分開，所以本文在EMD處理的基礎(chǔ)上結(jié)合分形盒維數(shù)的方法，求出每種情況下IMF1～I(xiàn)MF4的盒維數(shù)，準(zhǔn)確、清楚地反映出四種耦合故障特征。

圖6 y3邊界譜Fig.6 Boundary spectrum of y3

圖7 碰摩松動(dòng)故障信號(hào)的邊界譜Fig.7 Boundary spectrum ofrub-impact and pedestal looseness signal

圖8 裂紋松動(dòng)故障信號(hào)的邊界譜Fig.8 Boundary spectrum of crack and pedestal looseness signal

圖9 裂紋碰摩故障信號(hào)的邊界譜Fig.9 Boundary spectrum of crack and rub－impact signal

由表1可知，EMD處理后的故障信號(hào)分形盒維數(shù)有明顯的差異，其中IMF1到IMF4的盒維數(shù)相差最大的為裂紋碰摩故障，最小的為裂紋松動(dòng)故障。三種故障耦合的信號(hào)盒維數(shù)值與碰摩松動(dòng)故障和裂紋碰摩故障差異較大，雖與裂紋松動(dòng)故障差異相對較小，但也能夠清晰分辨，所以分形盒維數(shù)可以作為故障的特征向量，對耦合故障進(jìn)行診斷。

表1 四種故障信號(hào)經(jīng)EMD處理后的盒維數(shù)Tab.1 Processed box dimension by EMD of four kinds of fault signal

5 結(jié)論

(1)提出一種基于EMD和分形盒維數(shù)的旋轉(zhuǎn)機(jī)械耦合故障診斷方法。該方法先用EMD對故障信號(hào)進(jìn)行分解，得到IMF，再根據(jù)分形計(jì)算方法求出IMF的盒維數(shù)，通過對盒維數(shù)的分析比較進(jìn)行故障診斷。

(2)建立含有裂紋－碰摩－松動(dòng)故障耦合的轉(zhuǎn)子－軸承系統(tǒng)動(dòng)力學(xué)模型。在考慮單個(gè)故障特性的基礎(chǔ)上得出模型運(yùn)動(dòng)微分方程，并用龍格－庫塔法解出含有耦合故障特征的振動(dòng)信號(hào)。

(3)把新方法應(yīng)用到所建故障模型的診斷中，得到裂紋、碰摩、松動(dòng)三種故障耦合故障特征向量。通過與傳統(tǒng)診斷方法比較，證明該方法能夠有效、準(zhǔn)確地對耦合故障進(jìn)行診斷。

［1］Hu N Q，Wen X S.The application of Duffing oscillator in characteristic signal detection of early fault［J］.Journal of Sound and Vibration，2003，268:917 －931.

［2］Sekhar A S，Mohanty A R，Prabhakar R.Vibrations of cracked rotor system:transverse crack versus slant crack［J］.Journal of Sound and Vibration，2005，279:1203－1217.

［3］馬 輝，孫 偉，任朝暉，等.多盤懸臂轉(zhuǎn)子系統(tǒng)支座松動(dòng)故障研究［J］.航空動(dòng)力學(xué)報(bào)，2009，24(7):1512－1517.MA Hui，SUN Wei， REN Zhao-hui et al. Research on pedestal looseness fault of multi-disc over-hung rotor system［J］.Journal of Aerospace Power，2009，24(7):1512－1517.

［4］姚紅良，李 鶴，李小彭，等.裂紋碰摩耦合故障轉(zhuǎn)子系統(tǒng)診斷分析［J］.振動(dòng)工程學(xué)報(bào)，2006，19(3):307－212.YAO Hong-liang，LI He，LI Xiao-peng et al.Diagnoses of coupling fault of crack and rub-impact in rotor systems［J］.Journal of Vibration Engineering.2006，19(3):307 －212.

［5］羅躍綱，聞邦椿.雙跨轉(zhuǎn)子系統(tǒng)裂紋－松動(dòng)耦合故障的非線性響應(yīng)［J］.航空動(dòng)力學(xué)報(bào)，2007，22(6):996－1001.LUO Yue-gang，WEN Bang-chun.Nonlinear responses of two-span rotor-bearing system with coupling faults of crack and pedestal looseness［J］.Journal of Aerospace Power，2007，22(6):996－1001.

［6］黃志偉，周建中，張孝遠(yuǎn)，等.水輪發(fā)電機(jī)組軸系松動(dòng)－碰摩耦合故障的動(dòng)態(tài)響應(yīng)［J］.西南交通大學(xué)學(xué)報(bào)，2011，46(1):121－126.HUANG Zhi-wei，ZHOU Jian-zhong，ZHANG Xiao-yuan，et al.Dynamic responses of shaft system of hydraulic generator set with coupling faults of bearing looseness and rotor rubimpact［J］.Journal of Southwest Jiaotong University，2011，46(1):121－126.

［7］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceedings of the Royal Society London A，1998，454(1971):903－995.

［8］Yang Y，He Y G，Chen J S，et al.A gear fault diagnosis using Hilbert spectrum based on MODWPT and a comparison with EMD approach［J］.Measurement，2009，42:542－511.

［9］唐宏賓，吳運(yùn)新，滑廣軍，等.基于EMD包絡(luò)譜分析的液壓泵故障診斷方法［J］.振動(dòng)與沖擊，2012，31(9):44－48.TANG Hong-bin，WU Yun-xin，HUA Guang-jun，et al.Fault diagnosis of pump using EMD and envelope spectrum analysis［J］.Journal of Vibration and Shock，2012，31(9):44－48.

［10］邵忍平，曹精明，李永龍.基于EMD小波閾值去噪和時(shí)頻分析的齒輪故障模式識(shí)別與診斷［J］，振動(dòng)與沖擊，2012，31(8):96 －101，106.SHAO Ren-ping，CAO Jing-ming，LI Yong-long.Gear fault pattern identification and diagnosis using time-requency analysis and wavelet threshold de-noising based on EMD［J］.Journal of Vibration and Shock，2012，31(8):96 －101，106.

［11］Mandelbrot B B.The Fractal Geometry of Nature［M］.San Francisco:Freeman，1982.

［12］陳 奇，趙 韓，黃 康.齒輪結(jié)合面切向接觸剛度分形計(jì)算模型研究［J］.農(nóng)業(yè)機(jī)械學(xué)報(bào)，2011，42(2):203－206.CHEN Qi，ZHAO Han，HUANG Kang.Fractal model of tangential contact stiffness between Gears’s Joint Surfaces［J］.Transactions of the Chinese Society for Agricultural Machinery，2011，42(2):203－206.

［13］鐘明壽，龍 源，謝全民，等.基于分形盒維數(shù)和多重分形的爆破地震波信號(hào)分析［J］.振動(dòng)與沖擊，2010，29(1):7－11.ZHONG Ming-shou，LONG Yuan，XIE Quan-min，et al.Signal analysis for blasting seismic wave based on fractal boxdimension and fractal［J］.Journal of Vibration and Shock，2010，29(1):7－11.

［14］郝 研，王太勇，萬 劍，等.分形盒維數(shù)抗噪研究及其在故障診斷中的應(yīng)用［J］.儀器儀表學(xué)報(bào)，2011，32(3):540－545.HAO Yan，WANG Tai-yong，WAN Jian，et al.Research on fractalbox dimension anti-noise performance and its application in faultdiagnosis［J］. Chinese Journalof Scientific Instrument，2011，32(3):540 －545.

［15］趙春暉，馬 爽，楊偉超.基于分形盒維數(shù)的頻譜感知技術(shù)研究［J］.電子與信息學(xué)報(bào)，2011，33(2):475－478.ZHAO Chun-hui，MA Shuang，YANG Wei-chao.Spectrum sensing in cognitive radios based on fractal box dimension［J］.Journal of Electronics＆ Information Technology，2011，33(2):475－478.