基于MMFE和可拓k-medoids聚類的軸承性能退化評估

趙聰聰，劉玉梅，趙穎慧，白 楊，施繼紅

(1.吉林農業大學 工程技術學院,長春 130118；2.吉林大學 交通學院,長春 130022；3.一汽研發總院 智能網聯開發院,長春 130011；4.一汽大眾汽車有限公司 技術開發部,長春 130011)

軸承是機械系統重要的零件之一，其運行狀態和使用壽命直接影響機械系統的運行可靠性[1]。軸承性能退化評估與故障診斷也一直是機械故障診斷領域的研究熱點與難點。軸承實際運行狀態受多種因素的影響，若能建立有效的性能退化評估模型，構建狀態參數與軸承性能退化之間的映射關系，則可實現對軸承性能退化程度的實時評估。這對提高軸承的運行安全性和機械系統的工作可靠性均具有重要意義[2]。

軸承性能退化評估涉及兩方面內容，即軸承狀態參數的提取和性能退化評估模型的構建。目前，在軸承性能退化評估和故障診斷領域，大多以軸承垂向振動信號為研究對象，忽略了其他方向的振動信息以及不同方向振動信息之間的關聯性，導致所提取的特征無法全面反應軸承的運行狀態。隨著非線性科學理論的發展，許多非線性動力學方法如樣本熵、近似熵、模糊熵、多尺度熵等被廣泛應用于機械故障診斷領域，并取得了較好效果[3]。但上述熵為單變量分析方法，不能描述不同方向振動信息之間的關聯性。Ahmed等[4]將多元樣本熵的概念拓展到多尺度，提出了多元多尺度樣本熵(Multivariate Multi-Scale Entropy，MMSE)，以從復雜性、互預測性和長時相關性的角度來描述多通道時間序列的動態相關關系[5]。MMSE自提出以來，就在機械故障診斷領域得到了廣泛應用。

為實現對軸承性能退化程度的分析，需構建定量評估模型[6]。退化評估模型包括概率模型和距離模型，前者如高斯混合模型、隱馬爾可夫模型等，后者如支持向量機、支持向量描述及聚類算法等。k-medoids算法作為聚類分析方法的一種，改進了k-means算法初始中心點可能取到非樣本點的缺陷，提高了k-means算法的聚類效果。Rai等[7]將EMD分解與k-medoids聚類算法相結合，利用無故障樣本和失效樣本建立聚類模型，實現了對軸承性能退化的有效評估。利用k-medoids算法對樣本進行聚類可以得到聚類中心，但需進一步根據樣本與聚類中心的距離來定量分析軸承的性能退化程度。為此，將可拓學關聯函數的概念引入到軸承性能退化模型的構建，使用可拓關聯度實現對軸承性能退化程度的定量評估。

基于上述分析，本文提出了一種結合多元多尺度熵、k-medoids算法和可拓學理論的軸承性能退化評估方法。以多尺度模糊熵為基礎，利用軸承多向振動信號構建多元變量，提取多元多尺度模糊熵(multivariate multi-scale fuzzy entropy，MMFE)信息。對軸承正常狀態樣本進行k-medoids聚類得到聚類中心，根據樣本與聚類中心之間的歐式距離確定可拓集合的經典域和節域，同時分析可拓關聯函數的適用范圍。最后根據待測樣本到聚類中心的距離與軸承正常狀態距離區間的關聯度實現對其性能退化程度的定性定量評估，并利用軸承全壽命疲勞試驗進行驗證。

1 MMFE算法

1.1 多元模糊熵

為計算多元模糊熵，首先需根據Takens嵌入理論生成多元嵌入向量。

Xm(i)=[x1,i,x1,i+λ1,…,x1,i+(m1-1)λ1,x2,i,

x2,i+λ2,…,x2,i+(m2-1)λ2,…,xp,i,

xp,i+λp,…,xp,i+(mp-1)λp]=

[zi,zi+1,…,zi+m-1](i=1,2,…,N-n)

(1)

經過多元嵌入重構后，得到N-n個復合延時向量Xm(i)。

max{|zi+l-1-zj+l-1|,l=1,2,…,m}

(2)

定義模糊隸屬度函數。

(3)

式中，n1和r分別為模糊隸屬度函數的邊界梯度和寬度。

(4)

(5)

對所有i求平均值，即得到嵌入維度m時的條件概率Bm(r)。

(6)

保持其他變量不變，將多元復合延遲向量由m維嵌入到m+1維。對于一個嵌入向量M=[m1,m2,…,mp]，系統可以轉換成嵌入向量為M=[m1,m2,…,mk+1,…,mp]的任意空間。由于包含p個時間序列，通過分別拓展mk+1(k=1,2,…,p)，可獲得p×(N-n)個重構向量Xm+1(i)。

(7)

(8)

多元模糊熵定義如下。

(9)

1.2 多元多尺度模糊熵MMFE

MMFE考慮了多通道時間序列中各序列的自相似性和序列之間的互預測性[8]，計算過程如下。

(10)

式中，τ=1,2,…為尺度因子。

2 可拓k-mediods軸承性能評估模型

2.1 k-mediods算法

為實現上述過程，在使用k-medoids算法進行聚類時，首先在數據集A中隨機選擇k個點作為初始中心，即初始代表點，將剩余數據點依據就近原則分配到各類中；然后在各類中按順序遍歷選取非代表點代替初始代表點，重新計算聚類效果；通過反復迭代改進聚類質量，直至獲得穩定的聚類中心[10]。

2.2 可拓距及關聯函數

可拓學是我國學者蔡文教授創立的一門新興學科，主要研究事物拓展的可能性，為解決矛盾問題和不相容問題提供了新的方法和途徑[11]。可拓距作為可拓數學的基本概念，是計算可拓關聯函數的基礎。對于區間X0=[a,b]內的任一點x，經典數學認為x與X0的距離為0，而可拓學利用距描述x與X0的相對位置關系，將經典數學“類內即同”的思想拓展到“類內尚可分為不同層次”。設實軸上任意點x，其與區間X0的可拓距表示如下。

(11)

由式(11)知，x?X0時，ρ(x,X0)>0；x=a或b時，ρ(x,X0)=0；x∈X0時，ρ(x,X0)<0。

可拓距的概念說明了點與區間的位置關系，而關聯函數k(x)能夠描述點與區間的隸屬關系，且關聯函數以可拓距為基礎。設x0表示關聯函數取最大值時的點，即最優值點；x為實域上的任意元素，區間X0=〈a,b〉，X=〈c,d〉，X0?X，且X0、X無公共端點，則元素x關于區間X0的關聯函數與最優值點x0的位置有關。當x0在X0中點處取得時，初等關聯函數k(x)表示如下。

(12)

當x0不在X0中點時，關聯函數k(x)表示如下。

(13)

式(12)～(13)中，X0=〈a,b〉和X=〈c,d〉分別代表可拓集的經典域和節域；k(x)的正負和大小表明了x屬于或不屬于X0的程度；ρ(x,x0,X0)稱為側距。若最優值點x0∈[a,(a+b)/2]，則ρ(x,x0,X0)稱為左側距，記為ρl(x,x0,X0)；若x0∈[(a+b)/2,b]，則ρ(x,x0,X0)稱為右側距，記為ρr(x,x0,X0)。

(14)

(15)

關聯度k(x)描述了事物具有某種屬性的程度，能夠對事物屬性進行定性定量評價。由式(12)～(15)知，k(x)∈(-∞,+∞)，k(x)>0表示事物具有某種特征或滿足某種屬性；k(x)<0表示事物不具有該特征或不滿足該屬性；k(x)=0表示事物即具有該性質，又不具有該性質，屬于性質變化的中間階段。

2.3 軸承性能退化評估模型

由于可拓學關聯函數能定量描述點與區間的關聯程度，故將其與k-mediods算法相結合構建軸承性能退化評估模型。具體思想如下：利用k-mediods算法對軸承正常狀態樣本進行聚類，獲得聚類中心Cnormal；通過軸承狀態樣本與Cnormal之間的距離確定經典域X0和節域X；對軸承任一狀態樣本而言，計算該樣本到Cnormal的歐氏距離與經典域X0之間的關聯度，根據關聯度完成軸承性能退化程度的定性定量評估。具體流程如圖1所示。

圖1 軸承性能退化評估流程Fig.1 Assessment process of bearing performance degradation

在利用上述方法構建軸承性能退化評估模型時，應明確兩點問題：①可拓集合經典域X0和節域X的確定；②關聯函數的形式選擇。可拓集合的確定過程如圖2所示：首先利用k-mediods算法對軸承正常狀態樣本進行聚類，得到聚類中心Cnormal。設n1和n2是正常樣本中與Cnormal之間歐式距離的最小值點和最大值點，將歐式距離記為d(Cnormal,n1)和d(Cnormal,n2)，以此作為經典域邊界值的近似參考，即令a=d(Cnormal,n1)，b=d(Cnormal,n2)。為確定節域的邊界值，計算全部樣本到Cnormal的歐式距離，并以最近點n3和最遠點n4與Cnormal之間的距離d(Cnormal,n3)和d(Cnormal,n4)作為節域的近似參考，即令c=d(Cnormal,n3)，d=d(Cnormal,n4)。本文中，樣本與Cnormal之間的距離(記為dis)越趨于a，則樣本越靠近聚類中心，表明其所代表的軸承運行狀態屬于正常的程度越大；反之，dis越趨于b，樣本越遠離Cnormal，其所代表的軸承運行狀態屬于正常的程度越弱。因此，最優點在區間X0的左端點處取得，故采用式(13)和式(14)進行關聯度的計算。

圖2 可拓集合參考點的選取Fig.2 Selection of reference points to extension sets

3 MMFE對軸承運行狀態的表征能力分析

3.1 試驗裝置

為說明MMFE在表征軸承運行狀態方面的優越性，現利用美國辛辛那提大學智能維護中心提供的滾動軸承全壽命疲勞試驗數據進行分析[12]。試驗臺結構如圖3所示，主軸裝有四個Rexnord ZA-2115雙列滾柱軸承。軸承節圓直徑7.15 cm，滾動體直徑0.84 cm，接觸角15.17°，利用杠桿機構通過軸承座2和3向主軸施加2 721.6 kg的徑向載荷。

圖3 滾動軸承全壽命疲勞試驗臺Fig.3 Test bench for rolling bearing run-to-failure

交流電機通過皮帶耦合驅動主軸，電機轉速2 000 r/min。每個軸承座的水平和垂直方向裝有PCB353B33壓電式加速度傳感器，采集8通道時間序列。軸承疲勞試驗全過程共采集2 156個數據文件，采樣頻率20 kHz。前43個文件的采樣時間間隔為5 min，剩余文件采樣間隔為10 min。試驗結束時，軸承3發生內圈故障，軸承4發生滾動體故障。本文以軸承3為分析對象，利用垂向和橫向的振動加速度信號構建二元多尺度模糊熵。

由于所獲得的樣本文件較多，為便于分析，每半小時提取一次軸承振動數據作為樣本，共獲得690個樣本。軸承在試驗初期處于磨合狀態，振動幅值較大，會影響分析結果的準確性。因此，本文剔除前40組樣本，利用剩余的650組樣本進行分析。

3.2 MMFE參數選擇

由MMFE定義知，嵌入維數向量M、多維時延向量λ、時間序列長度N、模糊隸屬度函數的梯度n1和寬度r都會影響MMFE的計算結果。一般嵌入維數mk=2或3、λk=1(k=1,2,…,p)時，對MMFE的影響較小。樣本長度N與嵌入維數m滿足N=10m～30m。由于采樣頻率為20 kHz，為保證頻率分辨率和信息的完備性，樣本長度N應大于1 000點，故嵌入維數mk=3。依據文獻[13]，設定n1=2，r=0.15SD(SD為多通道數據的標準差)。現利用試驗分析的方法確定樣本長度N。首先利用db3小波分解方法對軸承振動信號進行降噪，而后在軸承全壽命疲勞試驗的初期、中期和末期各取10組數據(初期選擇第21～30組，中期為331～340組，末期為651～660組)，以代表軸承性能退化的不同階段。分別計算樣本長度N=1 000、2 000、3 000、5 000和7 000時的MMFE均值，其中尺度因子τ=20，結果如圖4所示。

由圖4可知，隨著樣本長度的增加，各尺度下MMFE波動減小，趨于平穩。這說明樣本長度越大，MMFE對軸承運行狀態的表征能力越穩定。但應注意，樣本長度增加，計算量隨之增加。當樣本長度N=1 000時，滾動軸承任一運行狀態下振動信號的二元MMFE波動較大，此時MMFE對軸承運行狀態的表征能力較差。軸承二元振動信號各尺度下的MMFE在N=5 000和7 000時相差不大，幾近相等。當N=3 000時，各尺度下MMFE波動較小，且與N=5 000或7 000時的MMFE計算結果相差較小，特別是在疲勞試驗末期。綜合上述分析，本文取樣本長度N=3 000，得到滾動軸承在不同運行狀態下的MMFE變化情況，如圖5所示。

(a) 滾動軸承初期振動信號的MMFE

圖5 N=3 000時MMFE變化Fig.5 The change curve of MMFE with N=3 000

由圖5知，在軸承性能退化的不同階段，各尺度下的MMFE不同，說明MMFE可作為表征軸承運行狀態的特征參數。隨著疲勞試驗的進行，軸承從正常狀態經一系列性能退化過程到最終失效，振動信號的隨機性減弱，確定性增強，信號構成由復雜逐漸變得簡單，故MMFE逐漸減小，與圖5結果相符。

3.3 MMFE與MFE性能對比

為說明軸承垂向振動信號與橫向振動信號之間的關聯性，表明MMFE比多尺度模糊熵(MFE)能更有效地表征軸承運行狀態，現分別對前文滾動軸承三種狀態下的30組垂向和橫向樣本進行分析，計算每種狀態下10組樣本的MFE均值。參數設置同MMFE算法，數據點數N=3 000，結果如圖6所示。

(a) 滾動軸承橫向振動信號MFE

由圖6可以看出，單向振動信號的MFE對軸承運行狀態的區分能力弱于MMFE，尤其在小尺度和高尺度因子的情況下。此外，對比圖5和圖6可知，MFE難以區分軸承的初期和中期狀態，而MMFE在小尺度因子下能有效區分這兩種狀態，說明MMFE比MFE能更有效地識別軸承早期性能退化。為進一步說明軸承橫向與垂向振動信息之間的相關性，現以滾動軸承運行中期的10組樣本為例，分析兩向振動信號MMFE和單向振動信號MFE的均值，結果如圖7所示。

圖7 MMFE與MFE特征對比Fig.7 Comparison of MMFE and MFE

根據相關性理論知，兩組數據之間的相關性越高，復雜度越低，熵值越小；反之亦然。由于MMFE均值在整個尺度因子上小于MFE，表明MMFE考慮了軸承兩向振動信息之間的相關性。相比于軸承單向振動信號的MFE特征，各尺度下MMFE變化更為平穩，且在較高尺度下MMFE趨于平穩。因此，利用軸承兩向振動信號所構建的MMFE特征能更準確的描述軸承運行狀態。

4 軸承全壽命疲勞試驗分析

4.1 最優尺度MMFE選取

在利用20個尺度下的MMFE構建高維特征向量進行聚類時，會導致計算量過大，聚類精度降低。此外，由圖5知，不同尺度下的MMFE所包含的信息量不同：1≤τ≤9時，MMFE能較好地區分滾動軸承的不同運行狀態；10≤τ≤20時，MMFE的區分能力減弱。因此，需選取最優尺度下的MMFE來構建軸承運行狀態特征向量。主成分分析(principal component analysis，PCA)具有去除冗余信息、降低特征向量維數、提高識別率、保留主要特征的特點[14]。因此，本文利用PCA方法對20維MMFE特征向量進行降維。前10個尺度的MMFE累計貢獻率見表1。

表1 MMFE特征主成分分析Tab.1 Principal component analysis of MMFE

由表1知，前7個尺度下的MMFE累計貢獻率大于90%，涵蓋了絕大部分的有用信息。由圖4和圖5知，軸承不同運行狀態的MMFE在第8個尺度取得最大值。故為保證有用信息的完備性，利用前8個尺度下的MMFE構建滾動軸承的運行狀態特征向量。

4.2 可拓k-medoids聚類分析

在利用k-medoids算法獲得表征軸承正常運行狀態的聚類中心Cnormal時，所采用的樣本全部來自于軸承正常運行狀態，故樣本類別k=1。軸承正常狀態樣本數量Num的多少會影響Cnormal的準確性。樣本量過少，無法覆蓋軸承正常運行狀態的可能范圍，Cnormal與實際值存在較大偏差。由于無法確定軸承性能退化的初始時刻，若樣本量過大，部分樣本可能來源于軸承早期性退化或故障狀態，導致樣本分布變得分散，同樣會使Cnormal偏離實際值。本文采用不同的樣本數量進行分析。令樣本量Num=20～440，每次分析時使樣本量增加20，計算不同樣本數量下各樣本點與Cnormal之間歐式距離的平均值，記為Dis，結果如圖8所示。

由圖8知，當樣本量Num=20時，平均距離Dis取得最小值。但由于樣本量過小，所獲得的Cnormal為局部中心，不能代表軸承全部正常狀態樣本的聚類中心。隨著Num的增加，樣本分布先分散后集中，Dis先增大后減小。當Num=180時，Dis取得極小值。隨著樣本量的進一步增加，樣本的聚類性減弱，平均距離再次增大，表明樣本分布變得分散。故利用k-dedoids聚類方法確定Cnormal時，樣本數Num=180，由此得到Cnormal=[0.371 6,0.415 4,0.478 2,0.496 2,0.541 0,0.572 7,0.559 6,0.615 2]。

圖8 樣本與Cnormal之間的平均距離Fig.8 The average distance between the samples and Cnormal

為表明Cnormal與軸承全部運行狀態下聚類中心Call的位置區別，現利用k-dedoids方法對代表軸承全部狀態的650組樣本進行聚類。此時不需要劃分軸承的性能退化階段，樣本代表了軸承運行中的任一狀態，故樣本類別k=1。聚類結果如圖9所示，聚類中心Call=[0.350 0,0.400 5,0.473 8,0.485 4,0.542 7,0.566 3,0.546 8,0.603 5]。

(a) 尺度因子τ=1 ～3

在獲得聚類中心Cnormal后，結合圖2利用樣本點與Cnormal之間的歐式距離來確定經典域X0和節域X。經典域X0代表了軸承正常狀態樣本與Cnormal之間歐式距離的可能范圍，故計算軸承180組正常狀態樣本與Cnormal之間的歐式距離，并以其最小值和最大值作為X0的下限和上限。節域代表了軸承任一狀態樣本與Cnormal之間歐氏距離的可能范圍。同理，根據軸承全部的650組樣本與Cnormal之間的歐式距離來確定節域X。通過計算，180組正常狀態樣本與Cnormal之間歐式距離的最小值和最大值分別為0.030 5和0.090 6，650組樣本與Cnormal之間歐式距離的最小值和最大值分別為0.029 1和0.162 4。因此，a=0.030 5，b=0.090 6，c=0.029 1，d=0.162 4，即X0=〈0.030 5,0.090 6〉，X=〈0.029 1,0.162 4〉。

在確定經典域X0和節域X后，進一步根據式(13)和式(14)計算待測樣本到Cnormal的歐氏距離與軸承正常狀態距離區間的關聯度，由此獲得軸承性能退化曲線。由于軸承兩向振動信號的獲取不可避免的受環境影響，其二元MMFE變化不平穩，導致關聯度的計算存在波動，如圖10中實線所示。為減小關聯度波動對軸承性能退化分析的影響，現利用四次多項式對關聯度曲線進行擬合，獲得關聯度的變化趨勢。擬合關系見式(16)，擬合結果如圖10中虛線所示。

圖10 軸承性能退化曲線Fig.10 The curve of bearing performance degradation

f(k)=-0.064 91×k4-0.194 9×k3-

0.113 5×k2-0.002 5×k+0.551 5

(16)

式中，k為關聯度的計算值。

結合關聯度定義和圖10關聯度擬合曲線知，從試驗開始到第440個樣本時，關聯度大于0.5且其擬合曲線變化平穩，表明待測樣本與Cnormal之間的歐式距離較小，待測樣本較靠近Cnormal，其所代表的軸承運行狀態良好。第440個樣本到第540個樣本時，關聯度擬合曲線出現下降趨勢，表明待測樣本與Cnormal之間的歐式距離增大，樣本逐漸遠離Cnormal，軸承開始出現早期性能退化。由于關聯度大于0，軸承仍處于正常運行狀態。試驗進行到第540個樣本時，關聯度為0，說明待測樣本與Cnormal之間的歐式距離達到經典域的上限，即達到正常樣本與Cnormal之間歐式距離的最大值。第540個樣本后，擬合曲線的關聯度值小于0，說明待測樣本與Cnormal之間的歐式距離已超出其正常運行范圍，軸承性能出現惡化。

4.3 時域特征對比分析

為說明本文所提方法對軸承性能退化程度的表征能力，提取軸承3兩向振動信號的典型有量綱和無量綱時域特征，結果如圖11所示。

(a) 峭度指標

由圖11知，隨著疲勞試驗的進行，軸承故障程度逐漸加重，時域指標隨之發生變化。峭度在第540個樣本處發生明顯波動，說明軸承發生故障，與圖10所示結果一致，但峭度無法有效檢測軸承的早期性能退化。平均值在軸承疲勞試驗中變化較為平緩，僅在試驗結束時才出現明顯波動，說明平均值不能有效表征軸承的性能退化過程。波形因子的變化趨勢與軸承性能退化趨勢相一致，但也僅在試驗末期才出現較大波動，無法檢測軸承的早期性能退化。

綜上分析，以上時域指標不能有效檢測軸承的早期性能退化，且無法定量表述軸承的性能退化程度。相比而言，本文所提方法可以描述軸承的性能退化趨勢，并能夠定量表示其性能退化程度。由圖11還可看出，上述時域指標在垂向和橫向存在明顯的相關關系，說明軸承兩向振動信號之間存在較大的相關性，進行兩向振動信號的綜合分析能獲得更準確的特征參數。

4.4 軸承早期性能退化分析

由前文分析知，軸承在第440個樣本處開始出現早期性能退化，在第540個樣本處發生內圈故障。通過計算得到軸承內圈故障的特征頻率理論值為296.9 Hz，利用頻譜分析法對以上兩個樣本的垂向振動信號進行功率譜分析，結果如圖12所示。

(a) 第440個樣本的降噪信號

由圖12可知，第440個樣本的特征頻率為271 Hz，表明軸承開始出現早期性能退化。第540個樣本的功率譜中出現了與軸承內圈故障頻率相近的291 Hz及其倍頻成分，說明軸承發生內圈故障，與本文分析結果相符。

5 結 論

(1) 將MMFE引入到軸承運行狀態的特征提取。利用軸承垂向和橫向振動信號構建了二元變量，并提取二元多尺度模糊熵。與軸承單向振動信號的MFE特征相比，二元MMFE能更有效的表征軸承運行狀態。

(2)將可拓學關聯函數與k-medoids聚類算法相結合構建了軸承性能退化評估模型。在對可拓集合和關聯函數形式進行分析的基礎上，利用待測樣本到Cnormal的歐式距離與軸承正常狀態距離區間之間的關聯度實現了軸承性能退化的定量評估。與常用的時域指標相比，本文所提方法能有效檢測軸承的早期性能退化，并能實現對軸承性能退化程度的定量表述。