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

王 巖，羅 倩，鄧 輝

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

0 引言

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

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

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

1 數據預處理和特征提取

1.1 數據LCD預處理

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

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

1.2 特征提取

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

1.3 特征篩選方法

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

類內標準差計算如下：

(1)

類間標準差計算如下：

(2)

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

χ=Si/So

(3)

其中：χ表示敏感度指標，特征對不同狀態的區別度與敏感度成正相關，某一特征的χ值越大，該特征對不同狀態的區別度越高。

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

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

2.1 特征值分布

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

2.2 多維高斯分布參數估計

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

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

(4)

(5)

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

(6)

(7)

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

(8)

(9)

軸承故障第j個狀態由均值μj和協方差矩陣Σj構成的多維高斯分布模型為Vj：

(10)

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

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

(11)

2.3 軸承故障診斷

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

3 仿真及實測數據診斷結果分析

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

表1 特征敏感度排名

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

圖1 診斷正確率與特征數關系

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

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

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

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

表2 多維高斯分布模型參數

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

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

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

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

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

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

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

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

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

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

4 結語

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

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

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

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