張 龍，黃文藝，熊國良

(華東交通大學(xué) 機電工程學(xué)院，南昌 330013)

載荷、轉(zhuǎn)速波動等外在因素以及系統(tǒng)內(nèi)部的非線性、故障的發(fā)生和發(fā)展等常導(dǎo)致機械系統(tǒng)出現(xiàn)非線性、非平穩(wěn)特征。以平穩(wěn)信號為對象的傳統(tǒng)信號處理方法難以處理頻域特征隨時間變化的非平穩(wěn)信號。目前已有多種適于非平穩(wěn)過程分析的時頻分析方法以及各種融合診斷方法如小波分析、經(jīng)驗?zāi)B(tài)分解(EMD)、全息譜、神經(jīng)網(wǎng)絡(luò)、支持向量機(SVM)和粗糙集等。

大量文獻利用上述方法進行故障診斷，但大多屬于故障的定性診斷，即只能識別故障的部位和類型，缺少定量指標(biāo)[1]。事實上，機械設(shè)備故障的演變存在一個由輕微到嚴(yán)重的發(fā)展過程，對于許多設(shè)備來說并不是出現(xiàn)了早期輕微故障就必須更換相應(yīng)的零部件。而應(yīng)準(zhǔn)確及時識別運行過程中故障的萌生和演變過程，這對于保障機械設(shè)備安全運行、避免經(jīng)濟損失和災(zāi)難性事故具有重大意義[2]。故障程度評估是指有效確定故障演變歷程和故障程度，已引起了廣泛的關(guān)注[3-6]。竇東陽等[3]利用EMD和Lempel-Ziv指數(shù)對滾動軸承單點損傷工況下的損傷程度進行識別，實驗表明Lempel-Ziv綜合指標(biāo)隨內(nèi)圈故障加重遞減，隨外圈故障加重而遞增。從飛云等[4]基于AR預(yù)測白噪化的Kolmogorov-Smirnov檢測方法，實現(xiàn)了滾動軸承疲勞實驗的全壽命檢測，該方法不僅能夠顯著地表現(xiàn)前期的微弱退化狀態(tài)，而且還能有條件地更早指示設(shè)備的異常狀態(tài)，對于故障預(yù)測的研究具有一定意義。江瑞龍等[5]提出了一種將模型訓(xùn)練和評估診斷相分離的思路，具體闡述了系統(tǒng)結(jié)構(gòu)及軟硬件實現(xiàn)，并通過軸承全壽命加速疲勞實驗數(shù)據(jù)測試，驗證了所提出方法的有效性和系統(tǒng)的可行性。潘玉娜等[6]提出了一種基于小波包和支持向量數(shù)據(jù)描述的軸承性能退化評估方法，以小波包分解的節(jié)點能量構(gòu)成特征向量，僅需要正常狀態(tài)下的數(shù)據(jù)樣本即可用支持向量數(shù)據(jù)描述建立知識庫，在一定程度上實現(xiàn)了對待測樣本退化程度的定量評估，軸承加速疲勞壽命實驗分析結(jié)果驗證了該方法的可行性和有效性。

上述研究工作在故障程度定量評估方面取得了一定的成效，但大多采用的是線性處理方法且需要冗長的數(shù)據(jù)樣本。本文提出一種基于多尺度熵的非線性定量評估方法，并利用滾動軸承故障仿真信號與全壽命周期疲勞實驗信號驗證了其有效性。

1 多尺度熵及多尺度熵均偏值

熵是一種度量時間序列復(fù)雜性的方法。最初，Pincus提出了近似熵，之后Richman等[7]提出了樣本熵，然而近似熵與樣本熵只反應(yīng)單尺度時間序列信息。Costa等[8]發(fā)現(xiàn)將單尺度熵算法用于從健康與疾病領(lǐng)域得到的數(shù)據(jù)時經(jīng)常會得到相互矛盾的結(jié)果，于是提出了多尺度熵(Multi-scale entropy)的概念。即利用粗粒序列來獲得多尺度時間序列代替原始的單尺度時間序列，然后在每個尺度上分析該生理學(xué)信號。對于滾動軸承而言，隨著故障程度的加深，熵值將逐步減小，即意味著故障信號的復(fù)雜度降低。

1.1 樣本熵

樣本熵(Sample entropy)是近似熵的改進算法，反映兩個長度都為m的相似序列(相似容差為r)，在將隨后的一個數(shù)據(jù)點加入變成長度為m+1的序列后變得不相似的概率。樣本熵具體算法如下[7]。

(1)假設(shè)原始數(shù)據(jù)X={x1，x2,…，xN}，定義X(i)=[xi，xi+1,…,xi+m-1]，1≤i≤N-m，其中m為嵌入維數(shù)；

(2)定義X(i)與X(j)之間的距離dij為兩者對應(yīng)元素之差中絕對值的最大值;

(3)對每一個i，計算X(i)與其余向量X(j)之間的距離dij(1≤j≤N-m，j≠i)，統(tǒng)計dij小于r的數(shù)目，記作num(dij

Bim(r)=num(dij

(1)

其中:r為相似容限，dij為小于r的數(shù)目;

(4)求Bim(r)的平均值Bm(r)；

(5)對于維數(shù)m+1，重復(fù)(1)～(4)，得到Bim+1(r)，進一步得到Bm+1(r)；

(6)最后得到樣本熵的表達式：

SampEn(m，r，N)=lnBm(r)-lnBm+1(r)

(2)

SampEn(m，r，N)的值與嵌入維數(shù)m，相似容限r(nóng)和數(shù)據(jù)長度N有關(guān)系。一般情況m取2；r取0.15*SD(SD是原始數(shù)據(jù)的標(biāo)準(zhǔn)差)；N=500～1000。

1.2 多尺度熵

所謂多尺度熵就是在原始信號的多個尺度上計算其樣本熵，不同尺度通過粗?；^程獲得[8]。

(1)原始數(shù)據(jù)X={x1，x2，...，xN}，長度為N，給定嵌入維數(shù)m和相似容限r(nóng)，建立新的粗粒向量y(s)

(3)

其中:1≤j≤N/s，尺度因子s=[1，2，…，smax]為正整數(shù)，本文取smax=20。

(2)對每一個粗粒序列求其樣本熵，則得到s個粗粒序列的樣本熵值，把多尺度熵繪制為尺度因子s的函數(shù)，稱為多尺度熵分析。

1.3 多尺度熵偏均值

最初Costa將多尺度熵應(yīng)用于生物信號分析，之后該方法在眾多領(lǐng)域獲得應(yīng)用。鄭佳波等[9]將多尺度熵分析應(yīng)用于流型轉(zhuǎn)化動力學(xué)特性分析，提出將多尺度熵變化速率特征作為流型識別的新指標(biāo)。何亮等[10]對金屬互連電遷移噪聲進行多尺度熵分析，通過與傳統(tǒng)表征參量對比，證明多尺度熵能夠?qū)﹄娺w移失效過程進行表征。在故障診斷方向，本文作者[11]曾以多尺度熵作為特征向量輸入到模糊神經(jīng)網(wǎng)絡(luò)進行軸承故障類型定性識別，識別效果較好。Wu等[12]利用多尺度熵分析與支持向量機進行軸承故障診斷，結(jié)果表明比樣本熵分析更為優(yōu)秀。大部分關(guān)于多尺度熵的文獻都沒有從多個熵值中提取一個綜合指標(biāo)以反映信號的復(fù)雜度或非線性度。文獻[11]和[12]也只是將多個尺度上的熵值作為特征向量輸入到分類器進行故障的定性識別，并沒有從多個尺度的熵值中提取一個指標(biāo)進行故障程度的定量化描述?；诖?，本文提出一個新的故障程度定量描述指標(biāo)——多尺度熵偏均值(Partial Mean of Multi-scale Entropy，PMME)。

(1)為了分析一組數(shù)據(jù)的總體趨勢，首先是求其均值。均值計算的目的在于觀察數(shù)據(jù)的集中趨勢，但單憑平均值不能完全表征一組數(shù)據(jù)的總體概況，有必要引入偏斜度這一概念。偏斜度越大，均值的效能越有問題，偏斜度越小，均值越可信賴。在對稱條件下，有:

x′=M0=Me

(4)

其中:x′為均值，M0為眾數(shù)，Me為中位數(shù)，在偏態(tài)分布條件下，三者存在著數(shù)量、位置差異。其中Me位于中間，x′與M0分居兩邊，因此，偏態(tài)a′可用x′與M0的絕對差額(距離)來表示，即:

a′=x′-M0

(5)

x′與M0的絕對差值大小相應(yīng)地表示偏斜程度大小。當(dāng)x′>M0，說明偏斜方向為右(正)偏；x′

由于偏態(tài)以絕對數(shù)表示，具有原數(shù)據(jù)量綱，因此不能直接比較不同數(shù)列的偏態(tài)程度。為了使不同數(shù)列的偏態(tài)值具有可比性，可計算偏態(tài)的相對值，即偏斜度Skewness(簡記為Ske)，又稱偏態(tài)系數(shù)，即偏態(tài)的絕對數(shù)與標(biāo)準(zhǔn)差之比。

Ske=(x′-M0)/SD=3(x′-Me)/SD

(6)

其中SD為數(shù)據(jù)的標(biāo)準(zhǔn)差。

偏斜度Ske的取值范圍一般在0與±3之間，0表示對稱分布，+3或-3分別表示極右或極左偏態(tài)。

(2)根據(jù)以上偏斜度的定義以及其與均值的關(guān)系，本文提出的PMME指標(biāo)定義如下。

假設(shè)求得某一故障信號X的多尺度熵MSE(X)={MSE(1)，MSE(2)，...，MSE(20)} (本文取最大尺度因子smax=20)，則有:

PMME=(1+|Ske(MSE)|/3)·mean(MSE)

(7)

其中:Ske(MSE)和mean(MSE)分別為20個尺度上樣本熵的偏斜度和均值。PMME值可比樣本熵更精確地表征樣本信號的復(fù)雜程度，從而可作為軸承故障程度定量化的一個重要指標(biāo)。

2 仿真信號分析

利用MATLAB程序仿真軸承內(nèi)圈不同故障程度信號。仿真信號X由故障信號xi與白噪聲r兩部分構(gòu)成，即X=axi+br，在b一定的情況下，a值越大表示故障越嚴(yán)重。令采樣頻率fs=20 480 Hz，信號長度L=102 400，轉(zhuǎn)頻fr=15 Hz，內(nèi)圈故障特征頻率fi=90 Hz，軸承結(jié)構(gòu)共振頻率fres=7 000 Hz。指數(shù)型衰減包絡(luò)E=∑iAe-ξ(t-i/fi)，其中ξ=0.048 8為衰減系數(shù)；A為脈沖幅值調(diào)制函數(shù)，對于內(nèi)圈故障A=cos(2πfrt)，當(dāng)0.5π≤2πfrt≤1.5π時A=0。進而得到故障信號xi和X。

xi=Esin(2πfrest)

X=axi+0.22r

(8)

圖1 不同程度內(nèi)圈故障仿真信號

式中：a為故障程度系數(shù)，a越大表示故障程度越嚴(yán)重。取a值分別為1、5、10，得到三種不同故障程度的仿真信號X1、X2和X3，分別表示輕度、中度以及嚴(yán)重故障，波形圖分別如圖1(a)、(b)及(c)所示。三個仿真信號的長度均為102 400，將它們各自分解為長1 024的100個短信號，分別計算這些短信號的PMME值，結(jié)果如圖2所示。顯然相同故障程度信號的PMME值集中在一定的范圍之內(nèi)，而不同故障程度的PMME值存在明顯區(qū)別且隨故障程度的增加而降低。仿真信號分析表明軸承內(nèi)圈信號所對應(yīng)的故障程度可以由該信號的PMME值定量表征。

圖2 三種不同故障程度仿真信號的多尺度熵偏均值

3 實驗分析

為了進一步驗證所提出的的多尺度熵偏均值定量表征軸承故障程度的有效性，下面利用該指標(biāo)分析實驗臺實測信號。

3.1 實驗數(shù)據(jù)

圖3(a) 實驗臺整體結(jié)構(gòu)示意圖；(b)實驗臺局部照片

軸承疲勞試驗數(shù)據(jù)來源于美國國家宇航局(NASA)預(yù)診斷數(shù)據(jù)庫，由美國辛辛那提大學(xué)、美國智能維護中心(IMS)Jay Lee教授課題組提供，實驗臺如圖3所示[13-15]。試驗臺主軸裝有四個Rexnord ZA-2115 雙列滾柱軸承，兩端軸承固定，通過中間兩個軸承由杠桿機構(gòu)向主軸施加272 1.6kg的正常大小的徑向載荷。實驗過程中各軸承均強制潤滑，四個軸承座上各裝有一個熱電偶測量軸承外圈溫度，以監(jiān)視潤滑情況。潤滑系統(tǒng)回油管內(nèi)裝有一個磁塞收集潤滑油液中的金屬磨粒，當(dāng)磨粒數(shù)量達到一定水平時說明軸承故障已較為嚴(yán)重，從而切斷電源，結(jié)束實驗。

Rexnord ZA-2115軸承每列有16個圓柱滾子，軸承節(jié)徑為7.15 cm，滾柱直徑0.84 cm，接觸角15.17°。可知軸承外圈故障特征頻率約為235 Hz。主軸由交流電機以2 000 r/min的轉(zhuǎn)速通過皮帶耦合驅(qū)動。每個軸承座在水平和垂直方向分別裝有一個PCB353B33壓電式加速度傳感器。數(shù)據(jù)采集從2004.2.12 10∶32∶39開始，至2004.2.19 06∶22∶39結(jié)束，期間每隔10分鐘采集振動信號一次，實驗過程總共采集到984個數(shù)據(jù)文件。采樣頻率為20 000 Hz，各傳感器每次采集20 480個數(shù)據(jù)點。本文分析的是數(shù)據(jù)文件第二列數(shù)據(jù)的前8 192點。停機后，解體發(fā)現(xiàn)軸承失效形式為嚴(yán)重外圈剝落故障。

3.2 實驗結(jié)果分析

性能退化評估中表征故障程度的指標(biāo)如峭度、烈度等通常不能反映故障類型[3-6]，但可以利用共振解調(diào)等方法確定故障發(fā)生部位。采用中心頻率為2 338 Hz，有效帶寬為900 Hz的Morlet濾波器對第982個文件的第二列數(shù)據(jù)進行共振解調(diào)。由圖4(b)和(c)的分析結(jié)果及軸承故障特征頻率可以確定軸承最終失效形式為外圈故障，與解體結(jié)果一致。

圖4 共振解調(diào)分析結(jié)果

計算984組實驗數(shù)據(jù)的PMME值，并進行5點平滑處理，結(jié)果如圖5所示。從圖中可以看到，PMME值在第575個樣本處出開始有持續(xù)下降趨勢，在第778個數(shù)據(jù)處開始急劇下降，隨著故障程度的劣化逐漸減小。表明PMME值能夠較早地發(fā)現(xiàn)并且正確跟蹤軸承故障的發(fā)展。

峭度(Kurtosis)是軸承故障程度定量檢測的一個重要指標(biāo)，實驗信號的峭度計算結(jié)果如圖6所示，大約第850個樣本處才開始有明顯的下降趨勢。顯然，和PMME值相比，峭度不能發(fā)現(xiàn)早期故障。可能原因有兩方面，一是只計算了原始信號的的峭度值，即單尺度上的結(jié)果；二是峭度不能像樣本熵一樣挖掘信號內(nèi)在非線性和復(fù)雜性特征。

圖5 (a)實驗數(shù)據(jù)PMME值; (b)局部放大圖

圖6 實驗數(shù)據(jù)的峭度

4 結(jié) 論

實現(xiàn)狀態(tài)維修、預(yù)測維修的基礎(chǔ)是準(zhǔn)確評估設(shè)備當(dāng)前的健康狀況。通常采用的峭度、均方根值及振動烈度等時域統(tǒng)計參數(shù)未能充分挖掘設(shè)備狀態(tài)信號中的特征信息，因此故障程度評估效果不佳。本文在多尺度熵的基礎(chǔ)上，提出了一個新的故障程度定量表征指標(biāo)——多尺度熵偏均值(PMME)，該指標(biāo)綜合了狀態(tài)信號在多個尺度上的非線性信息。滾動軸承仿真實驗和疲勞實驗分析表明該指標(biāo)能夠較早地發(fā)現(xiàn)軸承早期故障并且能與故障發(fā)展趨勢保持較好的一致性。

[1]雷亞國，何正嘉.混合智能故障診斷與預(yù)示技術(shù)的應(yīng)用進展[J].振動與沖擊，2011,30(9): 129-135.

LEI Ya-guo，HE Zheng-jia.Advances in applications of hybrid intelligent fault diagnosis and prognosis technique[J].Journal of Vibration and Shock，2011,30(9): 129-135.

[2]王國彪，何正嘉，陳雪峰，等.機械故障診斷基礎(chǔ)研究“何去何從”[J].機械工程學(xué)報，2013(1): 63-72.

WANG Guo-biao，HE Zheng-jia，CHEN Xue-feng，et al.Basic research on machinery fault diagnosis- what is the prescription [J].Journal of Mechanical Engineering，2013(1): 63-72.

[3]竇東陽，趙英凱.基于EMD和Lempel-Ziv指標(biāo)的滾動軸承損傷程度識別研究[J].振動與沖擊，2010,29(3): 5-8.

DOU Dong-yang，ZHAO Ying-kai.Fault severity assessment for rolling element bearings based on EMD and Lempel-Ziv index[J].Journal of Vibration and Shock，2010,29(3): 5-8.

[4]從飛云，陳進，董廣明.基于AR模型的Kolmogorov-Smirnov檢驗性能退化及預(yù)測研究[J].振動與沖擊，2012,31(10): 79-82.

CONG Fei-yun，CHEN Jin，DONG Guang-ming.Performance degradation assessment by Kolmogorov-Smirnov test and prognosis based on AR model[J].Journal of Vibration and Shock，2012,31(10): 79-82.

[5]江瑞龍，陳進，劉韜，等.設(shè)備性能退化評估在巡檢系統(tǒng)中的應(yīng)用[J].振動與沖擊，2012,31(16): 134-137.

JIANG Rui-long，CHEN Jin，LIU Tao，et al.Performance degradation assessment with the application of patrol inspection system[J].Journal of Vibration and Shock，2012,31(16): 134-137.

[6]潘玉娜，陳進.小波包-支持向量數(shù)據(jù)描述在軸承性能退化評估中的應(yīng)用研究[J].振動與沖擊，2009，28(4): 164-167.

PAN Yu-na，CHEN Jin.Wavelet package-support vector data description applied in bearing performance degradation assessment[J].Journal of Vibration and Shock，2009，28(4): 164-167.

[7]Richman J S，Moorman J R.Physiological time-series analysis using approximate entropy and sample entropy[J].American Journal of Physiology-Heart and Circulatory Physiology，2000，278(6): H2039-H2049.

[8]Costa M，Goldberger A L，Peng C K.Multiscale entropy analysis of complex physiologic time series[J].Physical review letters，2002，89(6): 68102.

[9]鄭桂波，金寧德.兩相流流型多尺度熵及動力學(xué)特性分析[J].物理學(xué)報,2009，58(7): 4485-4492.

ZHENG Gui-bo，JIN Ning-de.Multiscale entropy and dynamic characteristics of two-phase flow patterns[J].Acta Physica Sinica，2009，58(7): 4485-4492.

[10]何亮，杜磊，莊奕琪，等.金屬互連電遷移噪聲的多尺度熵復(fù)雜度分析[J].物理學(xué)報，2008，57(10): 6545-6546.

HE Liang，DU Lei，ZHUANG Yi-qi，et al.Multiscale entropy complexity analysis of metallic interconnection electromigration noise[J].Acta Physica Sinica，2008，57(10): 6545-6546.

[11]Zhang L，Xiong G，Liu H，et al.Bearing fault diagnosis using multi-scale entropy and adaptive neuro-fuzzy inference[J].Expert Systems with Applications.2010，37(8): 6077-6085.

[12]Wu S，Wu P，Wu C，et al.Bearing fault diagnosis based on multiscale permutation entropy and support vector machine[J].Entropy，2012，14(8): 1343-1356.

[13]Qiu H，Lee J，Lin J，et al.Robust performance degradation assessment methods for enhanced rolling element bearing prognostics[J].Advanced Engineering Informatics，2003，17(1): 127-140.

[14]Qiu H，Lee J，Lin J，et al.Wavelet filter-based weak signature detection method and its application on rolling element bearing prognostics[J].Journal of Sound and Vibration，2006，289(1): 1066-1090.

[15]Lee J，Qiu H，Yu G，et al.Rexnord technical services: ‘bearing data set’，IMS，University of Cincinnati.NASA Ames Prognostics Data Repository，2007.[http://ti.arc.nasa.g ov/project/prognostic-data-repository].NASA Ames，Moffett Field，CA.