趙榮珍， 李霽蒲， 鄧林峰

(蘭州理工大學(xué)機(jī)電工程學(xué)院 蘭州，730050)

引 言

作為旋轉(zhuǎn)機(jī)械設(shè)備的核心部件，滾動(dòng)軸承的狀態(tài)是否正常直接影響著機(jī)械設(shè)備的使用。統(tǒng)計(jì)發(fā)現(xiàn)，在旋轉(zhuǎn)機(jī)械發(fā)生異常的原因當(dāng)中，有30%是由軸承故障引起的[1]。因此,如何快速有效地對(duì)滾動(dòng)軸承故障給予辨識(shí)，已經(jīng)成為目前研究的一個(gè)重點(diǎn)。

對(duì)振動(dòng)信號(hào)進(jìn)行特征提取和故障識(shí)別是滾動(dòng)軸承故障診斷研究領(lǐng)域關(guān)注的重要問題之一。源于滾動(dòng)軸承的振動(dòng)信號(hào)呈現(xiàn)出的非線性和非平穩(wěn)性，使得傳統(tǒng)的以傅里葉變換為基礎(chǔ)的方法難以取得較好的分析效果[2]。Huang[3]提出的經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，簡稱EMD)是一種自適應(yīng)信號(hào)處理方法，該方法不受不確定性原理的限制，非常適合非線性和非平穩(wěn)信號(hào)的分析。但目前已發(fā)現(xiàn)這種算法仍存在著以下幾個(gè)問題待解決：a.固有模態(tài)函數(shù)(intrinsic mode function, 簡稱IMF)的正交性無法通過理論證明；b.收斂條件不合理、過包絡(luò)和欠包絡(luò)等問題產(chǎn)生的模態(tài)混淆，易導(dǎo)致IMF階數(shù)增加；c.計(jì)算過程需要多次迭代，得到一個(gè)實(shí)際的IMF分量需要較長時(shí)間等[4]。針對(duì)上述不足，Gilles[5]依據(jù)小波變換和窄帶信號(hào)分析理論，提出了一種新的自適應(yīng)信號(hào)處理方法——經(jīng)驗(yàn)小波變換(empirical wavelet transform, 簡稱EWT)，并把它成功應(yīng)用于心電圖信號(hào)(electrocardiogram,簡稱ECG)信號(hào)分離和圖像降噪分析中。李志農(nóng)等[3]不僅提出了基于EWT的機(jī)械故障診斷方法，而且將這種方法與傳統(tǒng)的EMD方法進(jìn)行的對(duì)比表明，EWT明顯優(yōu)于傳統(tǒng)的EMD方法，這一結(jié)論還被后續(xù)的一系列研究所證實(shí)[6-9]。

排列熵是一種檢測隨機(jī)性和動(dòng)力學(xué)突變的方法，它具有計(jì)算簡單、抗噪能力強(qiáng)等特點(diǎn)[10]。Yan等[11]將排列熵引入到旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)特征提取中的結(jié)果表明，該特征能夠有效地表征出滾動(dòng)軸承在不同狀態(tài)下的工況特征。類似于傳統(tǒng)的單尺度非線性參數(shù)，排列熵仍可在單一尺度上描述時(shí)間序列的不規(guī)則性。Aziz等[12]還提出了多尺度排列熵(multi-scale permutation entropy,簡稱MPE)的概念,用于衡量時(shí)間序列在不同尺度下的復(fù)雜性和隨機(jī)性，使魯棒性得到加強(qiáng)。鄭近德等[13]將多尺度排列熵概念運(yùn)用在軸承故障診斷上，也取得了良好的診斷效果。

聚類分析是模式識(shí)別的重要方法之一。常用的有模糊C均值(fuzzy C-mean, 簡稱FCM)、K-means聚類、GK(Gustafson-Kessel，簡稱GK)聚類、譜聚類等。其中的FCM采用歐式距離計(jì)算樣本之間的相似性，一般它僅適用于球形分布的數(shù)據(jù)集。GK算法引進(jìn)了自適應(yīng)距離范數(shù)和協(xié)方差矩陣，可以反映數(shù)據(jù)沿任意方向或子空間的分散程度，并沒有改變聚類算法產(chǎn)生類似于球體的聚類狀態(tài)[14]。GG聚類算法是FCM聚類算法和GK聚類算法被改進(jìn)的結(jié)果。因其通過引入模糊最大似然估計(jì)方法來度量樣本之間的距離，故可以反映不同形狀和不同方向的數(shù)據(jù)類。文獻(xiàn)[15]已將GG聚類算法成功應(yīng)用于滾動(dòng)軸承的故障診斷中，取得了較好的分類效果。

基于上述原因，本研究將對(duì)經(jīng)驗(yàn)小波變換、多尺度排列熵、GG聚類算法相結(jié)合并應(yīng)用于滾動(dòng)軸承的故障辨識(shí)方法進(jìn)行研究，欲通過實(shí)驗(yàn)驗(yàn)證所提出方法的有效性。

1 基礎(chǔ)理論

1.1 經(jīng)驗(yàn)小波變換EWT

EWT的核心思想是通過對(duì)信號(hào)的傅里葉頻譜進(jìn)行自適應(yīng)分割，構(gòu)建出一組適合于待處理信號(hào)的小波濾波器，以提取具有緊支撐傅里葉頻譜的AM-FM成分，然后對(duì)提取的AM-FM成分進(jìn)行Hilbert變換，最終得到有意義的瞬時(shí)頻率和瞬時(shí)幅值，進(jìn)而得到Hilbert譜。

圖1 傅里葉軸的分割Fig.1 Partitioning of the Fourier axis

(1)

(2)

傅里葉分割的關(guān)鍵是確定N。除去0，π兩個(gè)邊界之外，文獻(xiàn)[5]采用了易于理解的闕值法尋找其余N-1個(gè)邊界。尋找的途徑是通過計(jì)算頻域內(nèi)k個(gè)幅值極大值Mi(i=1,2,…,k)，按遞減排序M1≥M2≥…≥Mk并歸一化它們到[0,1]，規(guī)定Mk均需大于闕值Mk+α(M1-Mm)，其中α為取值在[0,1]之間的相對(duì)振幅比。在α被確定之后，N為大于闕值的極大值點(diǎn)的個(gè)數(shù)，因此可定義大于闕值的前N個(gè)最大值所對(duì)應(yīng)的ωn為N-1個(gè)邊界。

重構(gòu)原始信號(hào)f(t)的公式定義為

(3)

經(jīng)驗(yàn)?zāi)B(tài)fk(t)的定義如下

f0(t) =Wf(0,t)*φ1(t)

fk(t)=Wf(k,t)*ψk(t)

經(jīng)過對(duì)若干個(gè)經(jīng)驗(yàn)?zāi)B(tài)函數(shù)進(jìn)行Hilbert變換，即可得到Hilbert譜。

1.2 多尺度排列熵理論

1.2.1 排列熵算法

對(duì)于一維時(shí)間序列X={x(i)|i=1,2,…,n},設(shè)嵌入維數(shù)和延遲時(shí)間分別為m，τ，則按照Takens定理對(duì)X進(jìn)行相空間重構(gòu)，可得到式(4)所示的重構(gòu)矩陣，即

(4)

其中：j=1,2，…，K；K=n-(m-1)τ。

該矩陣共有K行，其中每一行均為一個(gè)重構(gòu)分量。如果用{j1,j2,…,jm}表示重構(gòu)分量中各個(gè)元素所在列的索引，則可將式(4)中的若干重構(gòu)分量按升序重新排列為式(5)，即

x(i+(j1-1)τ)≤x(i+(j2-1)τ)≤

…≤x(i+(jm-1)τ)

(5)

若重構(gòu)分量中有大小相等的情況，需通過比較j1和j2值進(jìn)行排序；當(dāng)j1

(6)

當(dāng)Pj=1/m!時(shí)，Hp(m)將達(dá)到最大值ln(m!)。對(duì)Hp(m)進(jìn)行歸一化處理，即

Hp(m)=Hp(m)/ln(m!)

(7)

顯然，Hp可代表X的隨機(jī)性；Hp越大，表示X的隨機(jī)程度越大；反之，說明X越規(guī)則。

1.2.2 多尺度排列熵算法

該算法首先將原始時(shí)間序列進(jìn)行粗?；幚?在多個(gè)尺度上計(jì)算時(shí)間序列的排列熵，然后計(jì)算各個(gè)不同尺度下的排列熵。具體的計(jì)算步驟如下。

1) 設(shè)長度為N的原始時(shí)間序列x(i)={x(1),x(2),…,x(N)},建立的新粗粒序列為

(8)

其中：s為尺度因子。

當(dāng)s=1時(shí)，粗?；蛄袨樵紩r(shí)間序列，稱為單尺度排列熵；當(dāng)s=n時(shí)，原始時(shí)間序列可被分割成為N/s個(gè)每段長度為s的粗粒序列。

2) 對(duì)得到的N/s個(gè)粗粒序列，求排列熵。

1.3 GG聚類算法

文獻(xiàn)[16]給出的GG聚類算法步驟如下。

1) 設(shè)聚類樣本集合為X={x1,x2,…,xi,…,xn}。其中：元素xi具有d個(gè)特征指標(biāo)。設(shè)利用隸屬度劃分矩陣U=(ui k)K×n作為判據(jù)，可將X聚成K類(2≤K≤n)。其中：ui k表示第k個(gè)樣本隸屬于第i個(gè)類別的程度。

2) 設(shè)終止容限為ε>0，隨機(jī)初始化隸屬矩陣為U。

3) 用式(9)計(jì)算各聚類中心，即

(9)

其中：h>1為加權(quán)指數(shù)；l為整數(shù)且l>1。

4) 計(jì)算模糊最大似然估計(jì)距離測度

(10)

其中：Ai為第i個(gè)聚類的協(xié)方差矩陣，pi為第i個(gè)聚類被選中的先驗(yàn)概率。

計(jì)算公式如式(11，12)，即

5) 更新用于分類的隸屬度矩陣

(13)

直到‖Ul-Ul-1‖<ε。

2 設(shè)計(jì)的故障辨識(shí)方法

本研究提出的將EWT、多尺度排列熵、GG聚類算法相結(jié)合的新故障診斷法具備如下特點(diǎn)。首先，它充分利用了EWT的自適應(yīng)分解和排列熵計(jì)算簡單、抗噪能力強(qiáng)的特點(diǎn)。其次，利用相關(guān)系數(shù)選取EWT分解后的最優(yōu)模態(tài)分量的處理方式，在減少了一定的計(jì)算量之外，將多個(gè)尺度下計(jì)算出的最優(yōu)模態(tài)其排列熵作為了特征向量。針對(duì)用此方式得到的熵值特征向量存在著高維度和數(shù)據(jù)難以可視化的新問題，故隨后利用主成分分析(principal component analysis, 簡稱PCA)進(jìn)行可視化降維，最后將得到的維數(shù)低、敏感度高且分類誤差率小的主要特征向量輸入至GG聚類算法中去實(shí)施聚類分析。

為上述過程設(shè)計(jì)的算法步驟如下：

1) 對(duì)振動(dòng)信號(hào)進(jìn)行經(jīng)驗(yàn)小波變換，得到若干個(gè)AM-FM分量；

2) 對(duì)若干AM-FM分量進(jìn)行相關(guān)性分析，相關(guān)系數(shù)最大的即為最優(yōu)模態(tài)分量，EMD各模態(tài)與原信號(hào)的相關(guān)性約等于各分量的自相關(guān)系數(shù)[17]，因此對(duì)得到的若干固有模態(tài)分量進(jìn)行相關(guān)性分析[18]并選出最優(yōu)模態(tài)作為下一步故障分類和識(shí)別的樣本，以達(dá)到剔除無關(guān)模態(tài)的目的；

3) 計(jì)算最優(yōu)模態(tài)分量的多尺度排列熵值，經(jīng)實(shí)驗(yàn)決定選取前9個(gè)尺度的排列熵作為特征向量；

4) 利用主成分分析對(duì)熵特征向量進(jìn)行降維；

5) 將降維后的低維特征向量作為GG聚類算法的輸入，并采用聚類評(píng)價(jià)指標(biāo)評(píng)價(jià)聚類效果。

與上述算法對(duì)應(yīng)的數(shù)據(jù)處理流程如圖2所示。

圖2 故障聚類方法流程Fig.2 The flow char of the fault classification method

3 實(shí)例分析

本研究采用美國凱斯西儲(chǔ)大學(xué)電氣工程實(shí)驗(yàn)室的滾動(dòng)軸承數(shù)據(jù)[19]對(duì)圖2算法進(jìn)行驗(yàn)證。軸承型號(hào)為SKF公司6205-2RS深溝球軸承。電機(jī)功率約1 494 W，轉(zhuǎn)速1 730 r/min，采用電火花加工技術(shù)在軸承上布置單點(diǎn)故障，凹坑的直徑×深度=0.177 8 mm×0.297 4 mm。設(shè)置的采樣頻率為12 kHz，采集{滾動(dòng)體故障(rolling element fault, 簡稱REF),內(nèi)圈故障(inner race fault, 簡稱IRF),外圈故障(outer race fault, 簡稱ORF),正常(NORM)}={REF, IRF, ORF, NORM}這4種狀態(tài)下的振動(dòng)信號(hào)各50組，每個(gè)信號(hào)長度為2 048。

隨機(jī)選取滾動(dòng)體故障信號(hào)進(jìn)行分析。圖3是經(jīng)過EWT自適應(yīng)分解后得到的10個(gè)AM-FM分量情況，表1是相關(guān)系數(shù)的計(jì)算結(jié)果。可以看出，C7與原始振動(dòng)信號(hào)的相關(guān)系數(shù)最大為0.856 6，因此選擇C7為最優(yōu)AM-FM分量進(jìn)行下一步的故障分類和識(shí)別。

圖3 滾動(dòng)體故障EWT處理結(jié)果Fig.3 Rolling element failure EWT processing results

Tab.1 AM-FM component and correlation coefficient of the original signal

分量序號(hào)相關(guān)系數(shù)分量序號(hào)相關(guān)系數(shù)C10.121 5C60.462 5C20.153 2C70.856 6C30.135 2C80.064 8C40.120 3C90.037 3C50.068 5C100.029 2

在計(jì)算多尺度排列熵時(shí)，需設(shè)定的參數(shù)有：時(shí)間序列長度N、嵌入維數(shù)m、時(shí)延τ和尺度因子s。其中，推薦的m=3～7[13]。m過大，將增大計(jì)算時(shí)間且無法反映序列的細(xì)微變化；m過小，重構(gòu)序列中可能包含的狀態(tài)會(huì)太少，則難以檢測出時(shí)間序列的動(dòng)態(tài)突變；經(jīng)過權(quán)衡，本研究中選取m=4。圖4為不同τ時(shí)的排列熵變化情況?？煽闯?，τ對(duì)信號(hào)熵值的影響較小，故選取τ=1。取s=12，計(jì)算12個(gè)粗粒向量的排列熵，得到的4種狀態(tài)結(jié)果見圖4。

圖4 滾動(dòng)軸承信號(hào)在不同時(shí)延下的排列熵Fig.4 Arrangement entropy of rolling bearing signals at different times

圖5是4種故障隨尺度因子s增大時(shí)最優(yōu)模態(tài)的多尺度排列熵變化情況。顯然，當(dāng)s=1、即單一尺度排列熵時(shí)，雖然正常狀態(tài)與故障狀態(tài)的熵值區(qū)別明顯，但是3種故障狀態(tài)的熵值很接近，此時(shí)很難區(qū)分3種故障狀態(tài)，因而需要對(duì)最優(yōu)模態(tài)分量進(jìn)行多尺度分析。由于前幾個(gè)熵值表征了EWT最優(yōu)模態(tài)的主要信息，通過實(shí)驗(yàn)決定采用前9個(gè)尺度的排列熵值作為特征向量。

圖5 EWT最優(yōu)模態(tài)的多尺度排列熵Fig.5 The optimal modal of multi-scale EWT permutation entropy

對(duì)滾動(dòng)軸承的4種狀態(tài)分別進(jìn)行EWT分解并通過相關(guān)性分析選取最優(yōu)模態(tài)分量，在9個(gè)尺度下計(jì)算最優(yōu)模態(tài)分量對(duì)應(yīng)的排列熵作為特征向量，可得到4組50×9的排列熵。采用極大似然估計(jì)法，用PCA算法將特征向量從9維降至3維。因?yàn)橛?種狀態(tài)，故聚類中心個(gè)數(shù)初步設(shè)定為4個(gè)，設(shè)加權(quán)指數(shù)h=2、迭代終止容差ε=0.000 1。將經(jīng)PCA降維后的特征向量輸入GG聚類算法，得到的結(jié)果如圖6所示。

圖6 EWT多尺度排列熵的GG聚類等高線圖Fig.6 The GG clustering contour for EWT multi-scale dimension permutation entropy

圖6中，V1～V4分別是{REF, IRF, ORF, NORM}4種狀態(tài)的聚類中心，它們的聚類中心坐標(biāo)值見表2。從圖6可看出，滾動(dòng)軸承的4種不同狀態(tài)不僅被明顯地分開、而且各類別均聚集在聚類中心附近，即聚集緊密、沒有出現(xiàn)不同類別相互混疊的現(xiàn)象，類間距較大。各組樣本集合的隸屬度平均值見表3。從表3可看出，第1組樣本對(duì)于V3的隸屬度遠(yuǎn)大于其他3組，說明第1組樣本隸屬于V3類，辨識(shí)效果良好。同理，第2～4組樣本分別可歸屬于V1，V4和V2類。

表2 4種不同類型信號(hào)的聚類中心

Tab.2 The clustering centers of the four different typessignals

樣本組號(hào)聚類中心xy內(nèi)圈故障V10.103 90.473 2外圈故障V20.667 50.648 1滾動(dòng)體故障V30.961 80.540 0正常狀態(tài)V40.754 30.177 3

表3 4種不同類型信號(hào)的隸屬度

為進(jìn)一步驗(yàn)證節(jié)2所提方法的有效性，本研究還采用直接對(duì)原始振動(dòng)信號(hào)進(jìn)行多尺度排列熵計(jì)算和EWT分解后提取單一尺度排列熵作對(duì)比試驗(yàn)，并把結(jié)果分別輸入FCM，GK，GG聚類算法，得到的聚類結(jié)果如圖7所示。通過對(duì)比分析得到的結(jié)論如下。

1) 將圖7(a)、圖7(b)與圖7(g)、圖7(h)進(jìn)行對(duì)比可以得出，直接對(duì)原始信號(hào)進(jìn)行多尺度排列熵提取得到的特征向量經(jīng)FCM，GK聚類算法處理的聚類效果不理想，而經(jīng)EWT多尺度排列熵提取的特征向量聚類效果較好。顯然，這應(yīng)該是原始信號(hào)經(jīng)EWT消噪、保留了更多故障信息的結(jié)果，故聚類效果較好。

圖7 不同模式組合下的2維聚類效果圖Fig.7 Different mode combination 2 dimension clustering results

2) 圖7(d)與圖7(e)聚類效果不佳。這是因?yàn)閱我怀叨扰帕徐夭⒉荒芎芎玫乇碚鞴收蠣顟B(tài)，而且經(jīng)過EWT后并沒有篩選出最優(yōu)模態(tài)分量，沒有去除冗余信息，這對(duì)聚類效果造成了一定的不利影響。

3) EWT多尺度排列熵組合較其他組合模型，其聚類結(jié)果的類內(nèi)緊致性最好。圖6、圖7(g)和圖7(h)的聚類中心基本吻合，圖7(a)、圖7(d)的聚類中心不一樣，說明聚類方法對(duì)聚類中心的影響較小，而聚類中心與輸入的特征向量有關(guān)。

總體而言，如果從聚類等高線的角度去評(píng)價(jià)，則FCM聚類的等高線接近于圓球形；GK聚類的等高線類似于橢圓；而GG聚類的等高線，應(yīng)該說具有等高線形狀不確定的特點(diǎn)。

使用分類系數(shù)(ppartition coefficient, 簡稱PC)和劃分熵(classification entropy, 簡稱CE)這兩個(gè)測度評(píng)價(jià)[16]，對(duì)圖6、圖7不同聚類模型的聚類效果評(píng)價(jià)結(jié)果見表4。

表4 不同模型FCM，GK，GG的聚類指標(biāo)

根據(jù)表4可得到如下推論：a.相對(duì)于FCM和GK聚類，GG聚類具有明顯的優(yōu)勢(shì)，這是由于FCM聚類的圓球形，它僅反映超球形數(shù)據(jù)結(jié)構(gòu)的標(biāo)準(zhǔn)距離規(guī)范；b.GK聚類的橢球形仍近似于球形；c.GG聚類因?yàn)榫垲愋螤钍侨我庑螤畹模示垲愋Ч梢苑从吵鰯?shù)據(jù)集的分散程度;d.EWT多尺度排列熵組合在GG聚類中，PC值最大達(dá)到了1，CE值最小為NAN。由此看出，EWT多尺度排列熵組合較其他模式組合而言具有一定優(yōu)勢(shì)，因此本研究提出的方法具有良好的滾動(dòng)軸承故障分類與辨識(shí)性能。

4 結(jié)束語

經(jīng)驗(yàn)小波變換是近幾年來興起的一種新自適應(yīng)信號(hào)處理方法。它具有計(jì)算簡單、計(jì)算速度快、理論基礎(chǔ)完備的特點(diǎn)。在此基礎(chǔ)上，提出了一種經(jīng)驗(yàn)小波變換、多尺度排列熵、GG聚類算法相結(jié)合的故障診斷方法，并把它應(yīng)用到滾動(dòng)軸承的故障診斷中。該方法首先采用EWT對(duì)原始信號(hào)進(jìn)行分解、得到若干個(gè)固有模態(tài)分量，通過計(jì)算各個(gè)分量與原始信號(hào)的相關(guān)系數(shù)，選取最優(yōu)模態(tài)分量以剔除冗余信息，對(duì)最優(yōu)模態(tài)分量進(jìn)行多尺度的排列熵計(jì)算，由于得到的熵值特征向量具有高維度和數(shù)據(jù)無法可視化問題，故采用PCA進(jìn)行降維之后再輸入到GG聚類算法中。實(shí)驗(yàn)證明，本研究提出的故障診斷方法能夠較好地區(qū)分滾動(dòng)軸承的不同狀態(tài)，故障類間無重疊，是一種有效的自適應(yīng)故障特征提取和故障數(shù)據(jù)聚類與分類方法。