基于Gibbs抽樣的軸承故障診斷方法

王 巖，羅 倩，鄧 輝

(北京信息科技大學 信息與通信工程學院，北京 100101)(*通信作者電子郵箱wyleo7@qq.com)

0 引言

在旋轉(zhuǎn)機械中，滾動軸承是易于損壞的零部件之一。據(jù)有關(guān)資料統(tǒng)計，旋轉(zhuǎn)機械的故障中振動故障占70%，而30%的振動故障是由滾動軸承故障引起的[1]，因此滾動軸承故障診斷的理論和應用研究一直是一個重點。

隨著大數(shù)據(jù)和人工智能時代的到來，越來越多的機器學習算法被用于軸承故障診斷中的模式識別，如Ali等[2]采用神經(jīng)網(wǎng)絡，焦衛(wèi)東等[3]使用支持向量機(Support Vector Machine, SVM)進行軸承故障診斷，但它們或多或少地存在一些缺陷，比如支持向量機中針對每一種軸承故障都設置一個判別函數(shù)進行診斷，造成輸入空間中存在無法分類的區(qū)域，即診斷片面性的問題[4]。為解決此問題，本文提出了一種新的基于Gibbs抽樣的軸承故障診斷算法，該方法將概率統(tǒng)計學應用到軸承故障診斷中。Gibbs抽樣是馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo， MCMC)理論的一個分支，用來獲取一系列近似等于指定多維概率分布(比如2個或者多個隨機變量的聯(lián)合概率分布)觀察樣本的算法。使用Gibbs抽樣建立每一種軸承狀態(tài)的分布模型，針對一個待診斷的樣本，分別計算出屬于每一種軸承狀態(tài)的概率，選取最大的判定為該類別，從而達到故障診斷的目的，解決了軸承故障診斷中存在的診斷片面性問題。

基于Gibbs抽樣的滾動軸承故障診斷思路如下：首先對軸承故障振動信號進行局部特征尺度分解(Local Characteristic scale Decomposition， LCD)得到若干個內(nèi)稟尺度分量(Intrinsic Scale Component， ISC)，再對軸承故障振動信號和ISC分別提取時域特征，對所有時域特征使用特征的類間標準差和類內(nèi)標準差的比值篩選出敏感特征并組成特征集。針對不同種類的軸承故障，使用基于Gibbs抽樣算法訓練特征集產(chǎn)生不同的多維高斯分布模型。對于待診斷的故障數(shù)據(jù)，分別計算在每個模型中的概率密度，通過后驗分析得到概率，選取概率最大的結(jié)果將待診斷數(shù)據(jù)診斷為該故障類別。

1 數(shù)據(jù)預處理和特征提取

1.1 數(shù)據(jù)LCD預處理

準確提取軸承振動信號特征是軸承故障診斷中重要的一步，如果直接從這些非平穩(wěn)或非線性振動信號中提取特征勢必影響故障診斷的效果，因此特征提取前必須進行信號預處理。本文采用局部特征尺度分解進行信號預處理，它能夠同時在時域和頻域提供非平穩(wěn)信號的局部化信息，且避免了經(jīng)驗模態(tài)分解(Empirical Mode Decomposition, EMD)的頻率混淆問題[5]和局部均值分解(Local Mean Decomposition， LMD)的信號突變問題[6]，所以它對非平穩(wěn)信號處理效果比較好。

LCD是一種基于極值點的局部尺度參數(shù)自適應的非平穩(wěn)信號分解方法。該方法以任意相鄰的兩個極值點為跨度，通過分段的形式對信號進行線性變換來構(gòu)造基線函數(shù)，從而將原始信號自適應地分解為若干個ISC分量。其具體步驟為:將原始信號通過兩個循環(huán)分解為具有物理意義的內(nèi)稟尺度分量ISC和一個剩余項。第一個循環(huán)是尋找基線信號，并經(jīng)循環(huán)迭代求出ISC分量的過程，ISC分量應滿足ISC分量條件；第二個循環(huán)是將原始信號減掉已獲得的ISC分量后的剩余信號作為原始信號再次經(jīng)循環(huán)迭代得到ISC分量。重復上述步驟直至所得剩余項為一個單調(diào)函數(shù)為止[7]。對于任一時間序列，通過LCD分解可以得到若干個ISC和一個剩余項。LCD分解本質(zhì)是對非平穩(wěn)信號進行平穩(wěn)化處理，將信號中不同尺度的波動或趨勢逐級分解。因為信號的能量主要集中在前幾個ISC分量中，因此本文只取前4個ISC分量。

1.2 特征提取

在LCD預處理之后，對原始振動信號和4個ISC分量分別進行時域特征提取。故障特征可以分為無量綱特征和有量綱特征兩大類[8-9]。無量綱指標因其不受工況影響的特點，對剝落、壓痕等故障有很好的診斷效果；有量綱指標對磨損類故障診斷效果比較好。因此，為了能夠更加全面、準確地診斷故障，本文對原始信號和4個ISC信號分別提取了6個無量綱特征(偏斜度、峭度、峰值、波形指標、脈沖、裕度)和5個有量綱特征(有效值、均值、標準差、方根幅值、最大值)，共得到55個特征。

1.3 特征篩選方法

如果使用上述55個特征作為一個特征集，計算量過大且診斷效果不一定最佳，所以需要從中挑選出敏感特征用于模型訓練。本文使用了基于故障類間、類內(nèi)標準差比值的特征篩選方法。

類內(nèi)標準差計算如下：

(1)

類間標準差計算如下：

(2)

對于某一特征，其類間標準差與類內(nèi)標準差的比值，稱為該特征的敏感度指標。計算公式如下：

χ=Si/So

(3)

其中：χ表示敏感度指標，特征對不同狀態(tài)的區(qū)別度與敏感度成正相關(guān)，某一特征的χ值越大，該特征對不同狀態(tài)的區(qū)別度越高。

2 基于Gibbs抽樣的軸承故障診斷

本文只涉及到軸承運行的4種狀態(tài)：內(nèi)圈故障狀態(tài)、外圈故障狀態(tài)、滾動體故障狀態(tài)和正常狀態(tài)?；贕ibbs抽樣的軸承故障診斷基本思路是：從上述55個特征中選取排名靠前的若干敏感特征組成特征集，針對軸承的每一個狀態(tài)，對這些特征建立多維概率分布，使用Gibbs算法對多維概率分布進行參數(shù)估計，最終得到4個多維概率分布模型。對于一個待診斷軸承數(shù)據(jù)，分別計算在4種模型中的概率密度，通過后驗分析得到概率，選取最大概率判定待診斷數(shù)據(jù)為該軸承狀態(tài)類別。

2.1 特征值分布

建立特征的多維概率分布，首先需要分析每個特征在不同狀態(tài)下的一維分布，然后構(gòu)建出合適的多維概率分布。生成軸承不同狀態(tài)下各特征的分布直方圖，通過直觀觀察發(fā)現(xiàn)絕大多數(shù)直方圖服從高斯分布。進一步，使用Quantile-Quantile正態(tài)檢驗方法[10]和Michael擬合優(yōu)度檢驗方法[11]檢驗各特征分布，發(fā)現(xiàn)除信號峭度特征在任何軸承狀態(tài)下不滿足高斯分布外，其他特征都服從高斯分布，這些特征都落在90%接受區(qū)間內(nèi)，可以認為這些特征都服從高斯分布，因此，針對軸承的每一個狀態(tài)，可以建立特征集的多維高斯分布。

2.2 多維高斯分布參數(shù)估計

使用若干個排名靠前的特征組成特征集構(gòu)建多維高斯分布，對于軸承故障的第j個狀態(tài)的多維高斯分布模型Vj的參數(shù)有均值μj、協(xié)方差矩陣Σj，對這兩個參數(shù)進行參數(shù)估計。高斯分布均值μj的先驗概率服從高斯分布，協(xié)方差矩陣Σj的先驗概率服從逆Wishart分布，通過對參數(shù)先驗分布進行Gibbs抽樣，可以得到這兩個參數(shù)具體的值[12]。故首先對每個多維高斯分布分別建立參數(shù)的先驗分布。

1)從逆Wishart先驗分布中隨機抽樣協(xié)方差矩陣Σj，下標j表示軸承第j個狀態(tài)：

(4)

(5)

2)從高斯先驗分布中抽樣協(xié)方差矩陣μj，如式(6)：

(6)

(7)

第j個軸承狀態(tài)的多維高斯分布參數(shù)估計的Gibbs算法通過將上述兩個步驟運算N次，可以分別得到每個參數(shù)的N個抽樣結(jié)果，最終μj和Σj參數(shù)值由式(8)～(9)得到，其中舍棄了前N/2的抽樣結(jié)果[13]。

(8)

(9)

軸承故障第j個狀態(tài)由均值μj和協(xié)方差矩陣Σj構(gòu)成的多維高斯分布模型為Vj：

(10)

其中：X是待診斷數(shù)據(jù)根據(jù)訓練特征集提取的時域特征；D是特征集中包含特征的個數(shù)，即多維高斯分布的維數(shù)；uj表示X在第j個軸承狀態(tài)模型中計算得到概率密度。

得到概率密度后，需要進行后驗概率分析。在軸承故障診斷之前，假設該軸承被診斷為4種狀態(tài)的先驗概率Zj是相等的，即Zj=0.25。由貝葉斯定理[14]可得，X被診斷為第j個狀態(tài)的后驗概率pj為：

(11)

2.3 軸承故障診斷

使用軸承不同狀態(tài)的振動信號樣本數(shù)據(jù)經(jīng)由LCD并提取時域特征后，選擇D個敏感特征組成特征集，訓練產(chǎn)生4個多維高斯分布模型。對于待診斷的數(shù)據(jù)，先根據(jù)訓練特征集提取的時域特征分別計算式(10)得到概率密度u1、u2、u3和u4，再使用式(11)分別計算概率p1、p2、p3和p4，選取最大概率診斷為該待診斷數(shù)據(jù)的類別。

3 仿真及實測數(shù)據(jù)診斷結(jié)果分析

使用了凱斯西儲大學的滾動軸承振動數(shù)據(jù)進行仿真(數(shù)據(jù)來源：http://csegroups.case.edu/Bearingdatacenter/pages/download-datafile)。振動信號由驅(qū)動端加速度傳感器采集，軸承型號為SKF6205，電機轉(zhuǎn)速為1 772 r/s，使用電火花加工技術(shù)在驅(qū)動端軸承上布置了單點故障，內(nèi)圈、滾動體及外圈故障直徑都為0.007 inch(1 inch=2.54 cm)，外圈故障位于6點鐘位置，采樣頻率為12 kHz，實際故障頻率：內(nèi)圈為159.93 Hz，滾動體為139.20 Hz，外圈為105.88 Hz。訓練模型時，將軸承振動信號若干個采樣點作為一組，從凱斯西儲大學原始數(shù)據(jù)中提取一共240組組成訓練樣本，其中軸承4種狀態(tài)各占60組。測試模型時，從原始數(shù)據(jù)中另外提取一共240組組成測試數(shù)據(jù)，其中軸承4種狀態(tài)也各占60組。將訓練數(shù)據(jù)的時域特征按照敏感度從高到低排列，其中排名前10的特征如表1所示，逐個將特征加入特征集進行模型訓練，可以得到55組不同特征集的模型。測試數(shù)據(jù)按照訓練特征集提取特征，分別計算其在55組不同模型中的診斷正確率，可以得到特征集取不同數(shù)量的特征時該算法的診斷正確率。

表1 特征敏感度排名

從圖1中可以看出，當特征集的特征個數(shù)增加到2時，診斷正確率就達到100%。圖1中的診斷正確率隨著特征數(shù)的增加呈下降趨勢，因為特征按照本文提到的特征篩選排名方法，逐個順序加入特征集并計算正確率。按此排名順序，特征對于區(qū)分不同軸承故障狀態(tài)的敏感度隨著特征數(shù)增加而下降，第20個特征之后加入的特征敏感度較差，會對整體診斷正確率產(chǎn)生影響，因此圖1中診斷正確率隨著特征數(shù)的增加而下降，也驗證了特征篩選排名方法的有效性。

圖1 診斷正確率與特征數(shù)關(guān)系

將240組測試集數(shù)據(jù)分別代入使用4個特征作為特征集的訓練模型中，按照診斷結(jié)果的標簽對4維數(shù)據(jù)使用了多維標度法(Multi-Dimensional Scaling， MDS)[15]映射在二維空間中，如圖2所示。

圖2 高維測試集數(shù)據(jù)的二維映射(240組測試集數(shù)據(jù))

從圖2中可以看出，4種狀態(tài)數(shù)據(jù)特征分布較為鮮明，說明本文提出的算法可以很好地判別4種不同的軸承狀態(tài)。

對前4個特征組成的特征集進一步分析，使用上述240組訓練集產(chǎn)生軸承4種不同狀態(tài)的模型：內(nèi)圈故障V1、外圈故障V2、滾動體故障V3和正常狀態(tài)V4。每個模型的參數(shù)μj和Σj計算結(jié)果列于表2，在圖3中使用多維標度法將不同模型的參數(shù)μj映射在二維空間中。

表2 多維高斯分布模型參數(shù)

圖3 高維高斯分布均值的二維映射

為了驗證本文提出的算法優(yōu)于基于支持向量機的軸承故障診斷算法，使用凱斯西儲大學軸承數(shù)據(jù)中采樣頻率為48 kHz的數(shù)據(jù)，振動信號由風扇端加速度傳感器采集，軸承型號為SKF6205，電機轉(zhuǎn)速為1 797 r/s，使用電火花加工技術(shù)在驅(qū)動端軸承上布置了單點故障，內(nèi)圈、滾動體及外圈故障直徑都為0.021英寸，外圈故障位于3點鐘位置，實際故障頻率：內(nèi)圈為162.20 Hz，滾動體為141.09 Hz，外圈為107.30 Hz。使用本文中的特征提取、篩選及排名的方法，在軸承故障診斷中分別使用了基于Gibbs抽樣的診斷方法和基于支持向量機的診斷方法。原始數(shù)據(jù)中缺少采樣頻率為48 kHz的軸承正常狀態(tài)數(shù)據(jù)，因此只建立了有關(guān)內(nèi)圈故障、外圈故障及滾動體故障的模型，從凱斯西儲大學原始數(shù)據(jù)中提取一共180組數(shù)據(jù)組成訓練樣本，其中軸承3種狀態(tài)各占60組。測試模型時，從原始數(shù)據(jù)中另外提取一共180組組成測試數(shù)據(jù)，軸承3種狀態(tài)也各占60組。特征按照本文中的方法進行排名并依次順序加入特征集，計算每個特征集有關(guān)模型的診斷正確率。

如圖4所示：在特征數(shù)比較少時，兩種診斷方法都能達到最高100%的診斷正確率；隨著特征數(shù)的增加，特征對軸承不同狀態(tài)的敏感度降低，二者的診斷正確率都下降，但基于Gibbs抽樣診斷算法下降的幅度明顯比支持向量機?。惶貏e在特征數(shù)達到43時，基于Gibbs抽樣方法診斷正確率為82.8%，而基于支持向量機方法診斷正確率為71.7%。與基于支持向量機診斷算法相比，基于Gibbs抽樣診斷算法正確率提升了11.1個百分點；當特征數(shù)大于43后，兩種算法的診斷正確率差值進一步增大。說明Gibbs抽樣診斷算法對于敏感度低的特征診斷效果好于支持向量機，特別對于高維復雜數(shù)據(jù)的離群點及支持向量機輸入空間中存在無法分類的區(qū)域[4]，Gibbs抽樣診斷算法通過概率計算依然可以實現(xiàn)對軸承故障的診斷。

圖4 兩種算法的診斷正確率比較

凱斯西儲大學的數(shù)據(jù)是實驗室條件下得到的數(shù)據(jù)，與真實的數(shù)據(jù)還有一定差別，使用該方法對中國鐵路某局實際的滾動軸承數(shù)據(jù)進行了診斷。數(shù)據(jù)來自中國鐵路某局甲2014年某月。該數(shù)據(jù)只有軸承故障和正常兩個狀態(tài)，并沒有對故障進行細分，數(shù)據(jù)是經(jīng)過共振解調(diào)[15]預處理過的。圖5是其振動信號200個抽樣點的時間序列波形。

圖5 軸承振動信號時間序列

使用60組故障數(shù)據(jù)和60組正常數(shù)據(jù)作為訓練樣本，另外60組故障數(shù)據(jù)和60組正常數(shù)據(jù)作為測試數(shù)據(jù)。圖6是該局使用本文提出的診斷方法在不同數(shù)量特征集下的診斷正確率，由圖可知取14個特征組成特征集時診斷正確率最大，為95.8%。

圖6 某局甲的故障診斷正確率

使用14個特征組成特征集對測試集數(shù)據(jù)進行診斷，診斷結(jié)果的標簽使用了多維標度法(Multi-Dimensional Scaling， MDS)[15]顯示二維空間中，如圖7所示。從實際數(shù)據(jù)仿真結(jié)果可以看出，本文方法依然有較高的故障診斷正確率。

圖7 高維測試集數(shù)據(jù)的二維映射(14個特征組成特征集)

4 結(jié)語

本文提出了一種基于Gibbs抽樣的軸承故障診斷算法，通過理論分析及數(shù)據(jù)仿真結(jié)果分析表明:

1)本文算法將概率統(tǒng)計學運用到軸承故障診斷中，采用概率模型代替支持向量機中的判別函數(shù)進行故障診斷，參考每種模型概率密度的計算結(jié)果，通過后驗分析得到概率，以判別概率最大為準則，對于支持向量機輸入空間中無法分類的區(qū)域[4]也能有效判別，避免了判別偶然性和片面性的出現(xiàn)，使得識別結(jié)果更加客觀。與基于SVM的軸承診斷方法相比，在特征數(shù)為43時診斷正確率提升了11.1個百分點。

2)提出了將Gibbs抽樣算法應用于滾動軸承故障診斷，可以實現(xiàn)從信號處理到分類的智能故障診斷。

3)本文分別使用凱斯西儲大學實驗室條件下的數(shù)據(jù)和中國鐵路機車軸承的實際采集數(shù)據(jù)進行了算法驗證，結(jié)果表明本文方法能夠有效地診斷軸承故障狀態(tài)。