王 巖，羅 倩，鄧 輝

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

0 引言

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

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

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

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

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

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

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

1.2 特征提取

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

1.3 特征篩選方法

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

類內(nèi)標(biāo)準(zhǔn)差計(jì)算如下：

(1)

類間標(biāo)準(zhǔn)差計(jì)算如下：

(2)

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

χ=Si/So

(3)

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

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

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

2.1 特征值分布

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

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

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

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

(4)

(5)

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

(6)

(7)

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

(8)

(9)

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

(10)

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

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

(11)

2.3 軸承故障診斷

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

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

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

表1 特征敏感度排名

從圖1中可以看出，當(dāng)特征集的特征個(gè)數(shù)增加到2時(shí)，診斷正確率就達(dá)到100%。圖1中的診斷正確率隨著特征數(shù)的增加呈下降趨勢(shì)，因?yàn)樘卣靼凑毡疚奶岬降奶卣骱Y選排名方法，逐個(gè)順序加入特征集并計(jì)算正確率。按此排名順序，特征對(duì)于區(qū)分不同軸承故障狀態(tài)的敏感度隨著特征數(shù)增加而下降，第20個(gè)特征之后加入的特征敏感度較差，會(huì)對(duì)整體診斷正確率產(chǎn)生影響，因此圖1中診斷正確率隨著特征數(shù)的增加而下降，也驗(yàn)證了特征篩選排名方法的有效性。

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

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

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

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

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

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

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

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

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

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

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

圖5 軸承振動(dòng)信號(hào)時(shí)間序列

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

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

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

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

4 結(jié)語(yǔ)

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

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

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

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