趙榮珍, 李霽蒲, 鄧林峰
(蘭州理工大學(xué)機電工程學(xué)院 蘭州,730050)
作為旋轉(zhuǎn)機械設(shè)備的核心部件,滾動軸承的狀態(tài)是否正常直接影響著機械設(shè)備的使用。統(tǒng)計發(fā)現(xiàn),在旋轉(zhuǎn)機械發(fā)生異常的原因當中,有30%是由軸承故障引起的[1]。因此,如何快速有效地對滾動軸承故障給予辨識,已經(jīng)成為目前研究的一個重點。
對振動信號進行特征提取和故障識別是滾動軸承故障診斷研究領(lǐng)域關(guān)注的重要問題之一。源于滾動軸承的振動信號呈現(xiàn)出的非線性和非平穩(wěn)性,使得傳統(tǒng)的以傅里葉變換為基礎(chǔ)的方法難以取得較好的分析效果[2]。Huang[3]提出的經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,簡稱EMD)是一種自適應(yīng)信號處理方法,該方法不受不確定性原理的限制,非常適合非線性和非平穩(wěn)信號的分析。但目前已發(fā)現(xiàn)這種算法仍存在著以下幾個問題待解決:a.固有模態(tài)函數(shù)(intrinsic mode function, 簡稱IMF)的正交性無法通過理論證明;b.收斂條件不合理、過包絡(luò)和欠包絡(luò)等問題產(chǎn)生的模態(tài)混淆,易導(dǎo)致IMF階數(shù)增加;c.計算過程需要多次迭代,得到一個實際的IMF分量需要較長時間等[4]。針對上述不足,Gilles[5]依據(jù)小波變換和窄帶信號分析理論,提出了一種新的自適應(yīng)信號處理方法——經(jīng)驗小波變換(empirical wavelet transform, 簡稱EWT),并把它成功應(yīng)用于心電圖信號(electrocardiogram,簡稱ECG)信號分離和圖像降噪分析中。李志農(nóng)等[3]不僅提出了基于EWT的機械故障診斷方法,而且將這種方法與傳統(tǒng)的EMD方法進行的對比表明,EWT明顯優(yōu)于傳統(tǒng)的EMD方法,這一結(jié)論還被后續(xù)的一系列研究所證實[6-9]。
排列熵是一種檢測隨機性和動力學(xué)突變的方法,它具有計算簡單、抗噪能力強等特點[10]。Yan等[11]將排列熵引入到旋轉(zhuǎn)機械振動信號特征提取中的結(jié)果表明,該特征能夠有效地表征出滾動軸承在不同狀態(tài)下的工況特征。類似于傳統(tǒng)的單尺度非線性參數(shù),排列熵仍可在單一尺度上描述時間序列的不規(guī)則性。Aziz等[12]還提出了多尺度排列熵(multi-scale permutation entropy,簡稱MPE)的概念,用于衡量時間序列在不同尺度下的復(fù)雜性和隨機性,使魯棒性得到加強。鄭近德等[13]將多尺度排列熵概念運用在軸承故障診斷上,也取得了良好的診斷效果。
聚類分析是模式識別的重要方法之一。常用的有模糊C均值(fuzzy C-mean, 簡稱FCM)、K-means聚類、GK(Gustafson-Kessel,簡稱GK)聚類、譜聚類等。其中的FCM采用歐式距離計算樣本之間的相似性,一般它僅適用于球形分布的數(shù)據(jù)集。GK算法引進了自適應(yīng)距離范數(shù)和協(xié)方差矩陣,可以反映數(shù)據(jù)沿任意方向或子空間的分散程度,并沒有改變聚類算法產(chǎn)生類似于球體的聚類狀態(tài)[14]。GG聚類算法是FCM聚類算法和GK聚類算法被改進的結(jié)果。因其通過引入模糊最大似然估計方法來度量樣本之間的距離,故可以反映不同形狀和不同方向的數(shù)據(jù)類。文獻[15]已將GG聚類算法成功應(yīng)用于滾動軸承的故障診斷中,取得了較好的分類效果。
基于上述原因,本研究將對經(jīng)驗小波變換、多尺度排列熵、GG聚類算法相結(jié)合并應(yīng)用于滾動軸承的故障辨識方法進行研究,欲通過實驗驗證所提出方法的有效性。
EWT的核心思想是通過對信號的傅里葉頻譜進行自適應(yīng)分割,構(gòu)建出一組適合于待處理信號的小波濾波器,以提取具有緊支撐傅里葉頻譜的AM-FM成分,然后對提取的AM-FM成分進行Hilbert變換,最終得到有意義的瞬時頻率和瞬時幅值,進而得到Hilbert譜。
圖1 傅里葉軸的分割Fig.1 Partitioning of the Fourier axis
(1)
(2)
傅里葉分割的關(guān)鍵是確定N。除去0,π兩個邊界之外,文獻[5]采用了易于理解的闕值法尋找其余N-1個邊界。尋找的途徑是通過計算頻域內(nèi)k個幅值極大值Mi(i=1,2,…,k),按遞減排序M1≥M2≥…≥Mk并歸一化它們到[0,1],規(guī)定Mk均需大于闕值Mk+α(M1-Mm),其中α為取值在[0,1]之間的相對振幅比。在α被確定之后,N為大于闕值的極大值點的個數(shù),因此可定義大于闕值的前N個最大值所對應(yīng)的ωn為N-1個邊界。
重構(gòu)原始信號f(t)的公式定義為
(3)
經(jīng)驗?zāi)B(tài)fk(t)的定義如下
f0(t) =Wf(0,t)*φ1(t)
fk(t)=Wf(k,t)*ψk(t)
經(jīng)過對若干個經(jīng)驗?zāi)B(tài)函數(shù)進行Hilbert變換,即可得到Hilbert譜。
1.2.1 排列熵算法
對于一維時間序列X={x(i)|i=1,2,…,n},設(shè)嵌入維數(shù)和延遲時間分別為m,τ,則按照Takens定理對X進行相空間重構(gòu),可得到式(4)所示的重構(gòu)矩陣,即
(4)
其中:j=1,2,…,K;K=n-(m-1)τ。
該矩陣共有K行,其中每一行均為一個重構(gòu)分量。如果用{j1,j2,…,jm}表示重構(gòu)分量中各個元素所在列的索引,則可將式(4)中的若干重構(gòu)分量按升序重新排列為式(5),即
x(i+(j1-1)τ)≤x(i+(j2-1)τ)≤
…≤x(i+(jm-1)τ)
(5)
若重構(gòu)分量中有大小相等的情況,需通過比較j1和j2值進行排序;當j1 (6) 當Pj=1/m!時,Hp(m)將達到最大值ln(m!)。對Hp(m)進行歸一化處理,即 Hp(m)=Hp(m)/ln(m!) (7) 顯然,Hp可代表X的隨機性;Hp越大,表示X的隨機程度越大;反之,說明X越規(guī)則。 1.2.2 多尺度排列熵算法 該算法首先將原始時間序列進行粗?;幚?在多個尺度上計算時間序列的排列熵,然后計算各個不同尺度下的排列熵。具體的計算步驟如下。 1) 設(shè)長度為N的原始時間序列x(i)={x(1),x(2),…,x(N)},建立的新粗粒序列為 (8) 其中:s為尺度因子。 當s=1時,粗?;蛄袨樵紩r間序列,稱為單尺度排列熵;當s=n時,原始時間序列可被分割成為N/s個每段長度為s的粗粒序列。 2) 對得到的N/s個粗粒序列,求排列熵。 文獻[16]給出的GG聚類算法步驟如下。 1) 設(shè)聚類樣本集合為X={x1,x2,…,xi,…,xn}。其中:元素xi具有d個特征指標。設(shè)利用隸屬度劃分矩陣U=(ui k)K×n作為判據(jù),可將X聚成K類(2≤K≤n)。其中:ui k表示第k個樣本隸屬于第i個類別的程度。 2) 設(shè)終止容限為ε>0,隨機初始化隸屬矩陣為U。 3) 用式(9)計算各聚類中心,即 (9) 其中:h>1為加權(quán)指數(shù);l為整數(shù)且l>1。 4) 計算模糊最大似然估計距離測度 (10) 其中:Ai為第i個聚類的協(xié)方差矩陣,pi為第i個聚類被選中的先驗概率。 計算公式如式(11,12),即 5) 更新用于分類的隸屬度矩陣 (13) 直到‖Ul-Ul-1‖<ε。 本研究提出的將EWT、多尺度排列熵、GG聚類算法相結(jié)合的新故障診斷法具備如下特點。首先,它充分利用了EWT的自適應(yīng)分解和排列熵計算簡單、抗噪能力強的特點。其次,利用相關(guān)系數(shù)選取EWT分解后的最優(yōu)模態(tài)分量的處理方式,在減少了一定的計算量之外,將多個尺度下計算出的最優(yōu)模態(tài)其排列熵作為了特征向量。針對用此方式得到的熵值特征向量存在著高維度和數(shù)據(jù)難以可視化的新問題,故隨后利用主成分分析(principal component analysis, 簡稱PCA)進行可視化降維,最后將得到的維數(shù)低、敏感度高且分類誤差率小的主要特征向量輸入至GG聚類算法中去實施聚類分析。 為上述過程設(shè)計的算法步驟如下: 1) 對振動信號進行經(jīng)驗小波變換,得到若干個AM-FM分量; 2) 對若干AM-FM分量進行相關(guān)性分析,相關(guān)系數(shù)最大的即為最優(yōu)模態(tài)分量,EMD各模態(tài)與原信號的相關(guān)性約等于各分量的自相關(guān)系數(shù)[17],因此對得到的若干固有模態(tài)分量進行相關(guān)性分析[18]并選出最優(yōu)模態(tài)作為下一步故障分類和識別的樣本,以達到剔除無關(guān)模態(tài)的目的; 3) 計算最優(yōu)模態(tài)分量的多尺度排列熵值,經(jīng)實驗決定選取前9個尺度的排列熵作為特征向量; 4) 利用主成分分析對熵特征向量進行降維; 5) 將降維后的低維特征向量作為GG聚類算法的輸入,并采用聚類評價指標評價聚類效果。 與上述算法對應(yīng)的數(shù)據(jù)處理流程如圖2所示。 圖2 故障聚類方法流程Fig.2 The flow char of the fault classification method 本研究采用美國凱斯西儲大學(xué)電氣工程實驗室的滾動軸承數(shù)據(jù)[19]對圖2算法進行驗證。軸承型號為SKF公司6205-2RS深溝球軸承。電機功率約1 494 W,轉(zhuǎn)速1 730 r/min,采用電火花加工技術(shù)在軸承上布置單點故障,凹坑的直徑×深度=0.177 8 mm×0.297 4 mm。設(shè)置的采樣頻率為12 kHz,采集{滾動體故障(rolling element fault, 簡稱REF),內(nèi)圈故障(inner race fault, 簡稱IRF),外圈故障(outer race fault, 簡稱ORF),正常(NORM)}={REF, IRF, ORF, NORM}這4種狀態(tài)下的振動信號各50組,每個信號長度為2 048。 隨機選取滾動體故障信號進行分析。圖3是經(jīng)過EWT自適應(yīng)分解后得到的10個AM-FM分量情況,表1是相關(guān)系數(shù)的計算結(jié)果。可以看出,C7與原始振動信號的相關(guān)系數(shù)最大為0.856 6,因此選擇C7為最優(yōu)AM-FM分量進行下一步的故障分類和識別。 圖3 滾動體故障EWT處理結(jié)果Fig.3 Rolling element failure EWT processing results Tab.1 AM-FM component and correlation coefficient of the original signal 分量序號相關(guān)系數(shù)分量序號相關(guān)系數(shù)C10.121 5C60.462 5C20.153 2C70.856 6C30.135 2C80.064 8C40.120 3C90.037 3C50.068 5C100.029 2 在計算多尺度排列熵時,需設(shè)定的參數(shù)有:時間序列長度N、嵌入維數(shù)m、時延τ和尺度因子s。其中,推薦的m=3~7[13]。m過大,將增大計算時間且無法反映序列的細微變化;m過小,重構(gòu)序列中可能包含的狀態(tài)會太少,則難以檢測出時間序列的動態(tài)突變;經(jīng)過權(quán)衡,本研究中選取m=4。圖4為不同τ時的排列熵變化情況??煽闯?,τ對信號熵值的影響較小,故選取τ=1。取s=12,計算12個粗粒向量的排列熵,得到的4種狀態(tài)結(jié)果見圖4。 圖4 滾動軸承信號在不同時延下的排列熵Fig.4 Arrangement entropy of rolling bearing signals at different times 圖5是4種故障隨尺度因子s增大時最優(yōu)模態(tài)的多尺度排列熵變化情況。顯然,當s=1、即單一尺度排列熵時,雖然正常狀態(tài)與故障狀態(tài)的熵值區(qū)別明顯,但是3種故障狀態(tài)的熵值很接近,此時很難區(qū)分3種故障狀態(tài),因而需要對最優(yōu)模態(tài)分量進行多尺度分析。由于前幾個熵值表征了EWT最優(yōu)模態(tài)的主要信息,通過實驗決定采用前9個尺度的排列熵值作為特征向量。 圖5 EWT最優(yōu)模態(tài)的多尺度排列熵Fig.5 The optimal modal of multi-scale EWT permutation entropy 對滾動軸承的4種狀態(tài)分別進行EWT分解并通過相關(guān)性分析選取最優(yōu)模態(tài)分量,在9個尺度下計算最優(yōu)模態(tài)分量對應(yīng)的排列熵作為特征向量,可得到4組50×9的排列熵。采用極大似然估計法,用PCA算法將特征向量從9維降至3維。因為有4種狀態(tài),故聚類中心個數(shù)初步設(shè)定為4個,設(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)的聚類中心,它們的聚類中心坐標值見表2。從圖6可看出,滾動軸承的4種不同狀態(tài)不僅被明顯地分開、而且各類別均聚集在聚類中心附近,即聚集緊密、沒有出現(xiàn)不同類別相互混疊的現(xiàn)象,類間距較大。各組樣本集合的隸屬度平均值見表3。從表3可看出,第1組樣本對于V3的隸屬度遠大于其他3組,說明第1組樣本隸屬于V3類,辨識效果良好。同理,第2~4組樣本分別可歸屬于V1,V4和V2類。 表2 4種不同類型信號的聚類中心 Tab.2 The clustering centers of the four different typessignals 樣本組號聚類中心xy內(nèi)圈故障V10.103 90.473 2外圈故障V20.667 50.648 1滾動體故障V30.961 80.540 0正常狀態(tài)V40.754 30.177 3 表3 4種不同類型信號的隸屬度 為進一步驗證節(jié)2所提方法的有效性,本研究還采用直接對原始振動信號進行多尺度排列熵計算和EWT分解后提取單一尺度排列熵作對比試驗,并把結(jié)果分別輸入FCM,GK,GG聚類算法,得到的聚類結(jié)果如圖7所示。通過對比分析得到的結(jié)論如下。 1) 將圖7(a)、圖7(b)與圖7(g)、圖7(h)進行對比可以得出,直接對原始信號進行多尺度排列熵提取得到的特征向量經(jīng)FCM,GK聚類算法處理的聚類效果不理想,而經(jīng)EWT多尺度排列熵提取的特征向量聚類效果較好。顯然,這應(yīng)該是原始信號經(jīng)EWT消噪、保留了更多故障信息的結(jié)果,故聚類效果較好。 圖7 不同模式組合下的2維聚類效果圖Fig.7 Different mode combination 2 dimension clustering results 2) 圖7(d)與圖7(e)聚類效果不佳。這是因為單一尺度排列熵并不能很好地表征故障狀態(tài),而且經(jīng)過EWT后并沒有篩選出最優(yōu)模態(tài)分量,沒有去除冗余信息,這對聚類效果造成了一定的不利影響。 3) EWT多尺度排列熵組合較其他組合模型,其聚類結(jié)果的類內(nèi)緊致性最好。圖6、圖7(g)和圖7(h)的聚類中心基本吻合,圖7(a)、圖7(d)的聚類中心不一樣,說明聚類方法對聚類中心的影響較小,而聚類中心與輸入的特征向量有關(guān)。 總體而言,如果從聚類等高線的角度去評價,則FCM聚類的等高線接近于圓球形;GK聚類的等高線類似于橢圓;而GG聚類的等高線,應(yīng)該說具有等高線形狀不確定的特點。 使用分類系數(shù)(ppartition coefficient, 簡稱PC)和劃分熵(classification entropy, 簡稱CE)這兩個測度評價[16],對圖6、圖7不同聚類模型的聚類效果評價結(jié)果見表4。 表4 不同模型FCM,GK,GG的聚類指標 根據(jù)表4可得到如下推論:a.相對于FCM和GK聚類,GG聚類具有明顯的優(yōu)勢,這是由于FCM聚類的圓球形,它僅反映超球形數(shù)據(jù)結(jié)構(gòu)的標準距離規(guī)范;b.GK聚類的橢球形仍近似于球形;c.GG聚類因為聚類形狀是任意形狀的,故聚類效果可以反映出數(shù)據(jù)集的分散程度;d.EWT多尺度排列熵組合在GG聚類中,PC值最大達到了1,CE值最小為NAN。由此看出,EWT多尺度排列熵組合較其他模式組合而言具有一定優(yōu)勢,因此本研究提出的方法具有良好的滾動軸承故障分類與辨識性能。 經(jīng)驗小波變換是近幾年來興起的一種新自適應(yīng)信號處理方法。它具有計算簡單、計算速度快、理論基礎(chǔ)完備的特點。在此基礎(chǔ)上,提出了一種經(jīng)驗小波變換、多尺度排列熵、GG聚類算法相結(jié)合的故障診斷方法,并把它應(yīng)用到滾動軸承的故障診斷中。該方法首先采用EWT對原始信號進行分解、得到若干個固有模態(tài)分量,通過計算各個分量與原始信號的相關(guān)系數(shù),選取最優(yōu)模態(tài)分量以剔除冗余信息,對最優(yōu)模態(tài)分量進行多尺度的排列熵計算,由于得到的熵值特征向量具有高維度和數(shù)據(jù)無法可視化問題,故采用PCA進行降維之后再輸入到GG聚類算法中。實驗證明,本研究提出的故障診斷方法能夠較好地區(qū)分滾動軸承的不同狀態(tài),故障類間無重疊,是一種有效的自適應(yīng)故障特征提取和故障數(shù)據(jù)聚類與分類方法。1.3 GG聚類算法
2 設(shè)計的故障辨識方法
3 實例分析
4 結(jié)束語