楊創(chuàng)艷，馬 軍，王曉東，羅 亭，李卓睿

（1.昆明理工大學(xué)信息工程與自動化學(xué)院，云南昆明 650500；2.昆明理工大學(xué)云南省人工智能重點實驗室，云南昆明 650500）

1 引言

滾動軸承是旋轉(zhuǎn)機械的關(guān)鍵零部件之一，其運行質(zhì)量決定設(shè)備運行的整體穩(wěn)定性和高效安全性，這就使利用退化指標(biāo)監(jiān)測軸承的運行狀態(tài)顯得至關(guān)重要［1～3］.但因軸承的振動特性受多種因素影響，傳統(tǒng)退化指標(biāo)無法準(zhǔn)確評估軸承退化狀態(tài)［4，5］.因此，如何構(gòu)建能有效反映滾動軸承運行狀態(tài)的退化指標(biāo)，已成為當(dāng)前設(shè)備狀態(tài)監(jiān)測領(lǐng)域的研究熱點之一.

近年來涌現(xiàn)出大量的退化指標(biāo)構(gòu)建方法，如：鄭小霞等人［6］提出一種基于隱馬爾科夫模型對數(shù)似然概率的退化指標(biāo)；Xu 等人［7］利用一個滑動平均磨損退化指標(biāo)來描述機械設(shè)備的退化過程；SHAO 等人［8］提出一種基于小波包能量矩熵的軸承性能退化指標(biāo).但以上單一退化指標(biāo)所包含的故障信息量少且抗干擾能力差，嚴重影響滾動軸承性能退化評估的準(zhǔn)確性［9］.所以，構(gòu)建一個能全面反映軸承退化信息的高維特征集D是軸承退化狀態(tài)評估的關(guān)鍵.

高維特征集D可體現(xiàn)不同特征間的差異互補性，但特征間可能存在重復(fù)性和冗余性［10～12］.因此，需要建立性能退化指標(biāo)評價準(zhǔn)則，有效剔除冗余特征.針對這一問題，鄭近德等人［10］利用Fisher得分評價準(zhǔn)則進行特征選擇，剔除冗余特征；戴豪民等人［11］采用加權(quán)最大相關(guān)最小冗余的特征評價準(zhǔn)則，選擇有效特征；Zhu 等人［12］引入多類特征評價準(zhǔn)則，選擇信息量最大的特征.然而，合理的退化指標(biāo)應(yīng)具有良好的單調(diào)性（Monotonicity，Mon）、相關(guān)性（Correlation，Corr）和魯棒性（Robust，Rob）［13，14］，應(yīng)綜合考慮Mon、Corr、Rob 的評價準(zhǔn)則，選擇能敏感表征軸承退化趨勢的特征，構(gòu)建敏感指標(biāo)集F.此外，如何綜合利用敏感退化指標(biāo)集F的特征優(yōu)勢是軸承退化狀態(tài)評估需要考慮的另一個核心問題.Tian 等人［15］將主成分分析（Principal Component Analysis，PCA）和K 近鄰距離融合成一個退化指標(biāo)；WU等人［16］通過度量PCA 的第一主元分量與其余分量之間的馬氏距離形成了退化指標(biāo)，但是經(jīng)過PCA 分解分量不能保證相互獨立；Sun 等人［17］使用KICA 非線性核函數(shù)將原始特征空間映射到核特征空間，分離得到相互獨立的分量.因此，結(jié)合KICA和馬氏距離融合特征，可以獲得敏感退化指標(biāo)KICAMD.

目前，利用敏感退化指標(biāo)進行軸承退化性能評估還存在一些尚未解決的問題，如敏感退化指標(biāo)存在虛假波動，且難以自適應(yīng)確定SFT.軸承退化狀態(tài)主要分以下階段：正常狀態(tài)，退化指標(biāo)保持穩(wěn)定；初始故障，退化指標(biāo)開始線性增長；嚴重退化，退化指標(biāo)開始斜率線性增長［18］.但是隨著時間的推移，軸承故障點會被磨平，導(dǎo)致退化指標(biāo)存在虛假波動［19］.因此，修復(fù)特征的虛假波動對于軸承性能退化的準(zhǔn)確評估至關(guān)重要［20～23］.近年來許多研究者在此方面進行了研究，如：Ahmad 等人［19］提出了一種線性回歸技術(shù)處理退化指標(biāo)的虛假波動；高彩霞等人［20］使用最小二乘法確定線性回歸方程來修復(fù)虛假波動；Li 等人［23］利用指數(shù)模型修復(fù)退化指標(biāo)的虛假波動.但上述方法無法實現(xiàn)軸承退化指標(biāo)的在線監(jiān)測，因此，本文提出了基于灰色回歸模型（Grey-regression Model，GM）修復(fù)的虛假波動.此外，確定軸承SFT，可提高性能退化評估的準(zhǔn)確性.確定SFT 的方法有：利用工程標(biāo)準(zhǔn)ISO-10816 的μ+3σ方法［21］，以及通過大量機器統(tǒng)計［22］.但是，這些方法計算結(jié)果誤差較大.因此，準(zhǔn)確確定SFT 及定量評估軸承的退化階段也直接關(guān)系到軸承運行狀態(tài)評估的有效性.

綜上所述，本文主要工作在于提出了一種特征融合與灰色回歸的滾動軸承性能退化評估方法.本文通過構(gòu)建敏感退化指標(biāo)，克服了傳統(tǒng)方法依靠經(jīng)驗評估軸承退化狀態(tài)的局限性；同時使用GM 修復(fù)退化指標(biāo)的虛假波動，基于時間序列轉(zhuǎn)折突變點自適應(yīng)確定SFT，提高了對軸承退化性能的評估能力；最后，通過美國國家航空航天局（National Aeronautics and Space Administration，NASA）和西安交通大學(xué)-長興昇陽科技有限公司（Xi’an Jiaotong University-Changxing Sumyoung Technology，XJTU-SY）的滾動軸承全壽命周期數(shù)據(jù)庫的數(shù)據(jù)驗證了所提方法的準(zhǔn)確性和有效性.

2 構(gòu)建敏感退化指標(biāo)

2.1 特征提取

2.1.1 時域統(tǒng)計特征

選擇文獻［16，17］設(shè)計的有量綱特征和無量綱特征來表征軸承振動信號xi（N為采樣點）的變化趨勢.有量綱特征包括均值（MV）、方差值（VARV）、均方根值（RMSV）、偏度值（SKV）；無量綱特征包括峭度因子（KF）、脈沖因子（PF）、形狀因子（SF）、裕度因子（MF）.特征計算方程如表1所示.

表1 時域統(tǒng)計特征

2.1.2 能量特征

軸承發(fā)生故障時，故障頻率范圍內(nèi)的信號發(fā)生變化，從而引起能量的變化.因此，引入能量E表征振動信號的沖擊性故障.計算公式為

2.1.3 熵特征

為反映故障深層次的信息，通過有效度量熵值大小來提取軸承故障狀態(tài)的本質(zhì)信息.因此，選用幾種常見熵，計算公式為

Shannon熵：

Renyi熵：

Tsallis熵：

其中，P(xi)代表軸承隨機退化事件xi的概率.

2.2 特征選擇

為了選擇高維特征集D中的有效退化特征，需要建立相關(guān)退化特征的評價準(zhǔn)則.一個優(yōu)良的退化指標(biāo)應(yīng)該滿足下列條件［16］：

（1）單調(diào)性（Mon），退化是一個不可逆的單調(diào)過程；

（2）相關(guān)性（Corr），特征與時間順序具有良好的相關(guān)性；

（3）魯棒性（Rob），對隨機噪聲具有良好的抗干擾性.

因此，基于Mon、Corr、Rob的綜合評價準(zhǔn)則，來選擇合理的退化特征，具體步驟如下.

Step1對高維特征集D中的特征使用平滑法分解為平滑趨勢項和隨機項，即

其中，f(tk)是tk時刻的特征值；fT(tk)是趨勢項；fR(tk)是隨機項.

Step2計算特征的Mon、Corr、Rob，即

其中，δ(·)是單位脈沖函數(shù).

Step3使用熵值法確定加權(quán)線性組合系數(shù)ωi，得到綜合評價指標(biāo)Z，即

Z和信號正相關(guān)，值越大說明該特征越應(yīng)該保留.保留大于閾值A(chǔ)（取Z的均值）［15］的特征，構(gòu)建敏感指標(biāo)集F；否則予以剔除.

2.3 特征融合

為了綜合F中所有特征的優(yōu)點，利用KICA 把輸入空間中的數(shù)據(jù)映射到特征空間，然后在特征空間中進行分析和處理，具體步驟如下.

Step1輸入數(shù)據(jù)集x=(x1,x2,…,xn)和核函數(shù)kij=(xi,xj),其中i,j=1,2,…,n（n為數(shù)據(jù)樣本）.

Step2構(gòu)造核矩陣K=k-Ink-kIn+InkIn，其中In為單位矩陣.

Step3數(shù)據(jù)集x白化處理，Χ=-1HTK，其中Λ為核矩陣K最大特征值，H為Λ對應(yīng)的特征向量.

Step4采用ICA 算法得到正交分離矩陣Α和獨立元T：T=ΑTΧ.

Step5計算第一獨立元與其余獨立元之間的馬氏距離，并選取其中最大值作為最佳投影向量，最終得到敏感指標(biāo)KICAMD，即

其中，μ和S-1分別是獨立成分的均值和協(xié)方差.

3 回歸修復(fù)虛假波動及狀態(tài)定量評估

在KICAMD 的基礎(chǔ)上，還需解決以下問題：①判定KICAMD是否存在虛假波動并修復(fù)；②基于退化指標(biāo)確定SFT和定量評估軸承的退化狀態(tài).

3.1 灰色回歸模型

GM 是由一個單變量的一階微分方程所構(gòu)成，其建立過程如下：對原始數(shù)據(jù)X={x(t),t=0,1,2,···,n}平滑處理，得到X={x0(t),t=0,1,2,···,n}，然后對x一階累加得X={x1(t),t=0,1,2,···,n}.

回歸方程為

方程的解為

離散化形式為

還原數(shù)據(jù)，得到GM，即

3.2 修復(fù)虛假波動

在軸承性能退化評估中，若出現(xiàn)如圖1所示的虛假波動，將嚴重影響退化指標(biāo)對軸承退化狀態(tài)的刻畫，降低軸承性能退化評估的準(zhǔn)確性［19～21］.

圖1 退化指標(biāo)虛假波動

因此，若存在虛假波動，需要修復(fù)才能作為后續(xù)評估的退化指標(biāo)，具體步驟如下.

Step1基于敏感退化指標(biāo)KICAMD建立滑動窗（n=50，滑動步長LS=1）［20］.

Step2在滑動窗上建立GM 得到預(yù)測值，同時，以滑動窗內(nèi)數(shù)據(jù)的3δ原則設(shè)置閾值B.

Step3判斷預(yù)測值是否大于閾值B.若大于則判定為虛假波動，并用預(yù)測值替代原始值；反之則用式（13）所示的相對指標(biāo)KICAMDr替代.

Step4更新滑動窗，進而得到退化指標(biāo)HI，即

3.3 確定SFT及退化狀態(tài)定量評估

引入T 分布隨機鄰接嵌入（T-Distributed Stochastic Neighbor Embedding，T-SNE）的高維數(shù)據(jù)可視化方法呈現(xiàn)軸承全壽命周期數(shù)據(jù)的全局退化趨勢.該方法實現(xiàn)了高維特征集D到（X，Y）兩個方向的特征的約簡［24，25］，如圖2所示.

圖2 T-SNE分布散點圖

根據(jù)圖2 可將滾動軸承運行狀態(tài)分為6 個階段，但不能直接確定SFT 和各狀態(tài)變化的起止點［25］.為此，引入時間序列突變點檢測算法計算狀態(tài)突變點［26］.對給定時間序列x=(x1,x2,…,xn)，按式（15）計算突變點k值［27，28］，即

4 特征融合與灰色回歸的滾動軸承性能退化評估過程

所提方法實現(xiàn)流程如圖3所示，具體包括如下三大過程.

圖3 方法實現(xiàn)流程

（1）構(gòu)建敏感指標(biāo)集

step1構(gòu)建高維特征集D=［MV，VARV，RMSV，SKV，KF，PF，SF，MF，E，SE，RE，TE］.

step2計算高維特征集D中每個特征的評價指標(biāo)：Mon、Corr、Rob.

step3計算加權(quán)綜合評價準(zhǔn)則Z及篩選閾值A(chǔ).

step4若Z大于A則作為有效特征，構(gòu)建敏感指標(biāo)集F；否則，剔除該特征.

（2）修復(fù)虛假波動

step1使用KICA 對敏感指標(biāo)集F進行映射處理獲得最佳投影向量，建立敏感退化指標(biāo)KICAMD.

step2基于KICAMD 建立滑動窗（n=50，LS=1），在滑動窗上建立GM 模型，得到預(yù)測值；同時，以滑動窗內(nèi)數(shù)據(jù)的3δ原則設(shè)置閾值B.

step3若預(yù)測值大于閾值B則判定為虛假波動，用預(yù)測值取代原始值；否則，取相對指標(biāo)ICAMDr.

step4更新滑動窗，得到修復(fù)虛假波動后的退化指標(biāo)HI.

（3）確定SFT和退化狀態(tài)定量評估

step1計算敏感退化曲線的轉(zhuǎn)折突變點.

step2最初的轉(zhuǎn)折突變點即為SFT.若退化指標(biāo)曲線處于SFT 前，數(shù)據(jù)處于正常狀態(tài)；否則，處于退化階段.

step3以SFT 為起點，基于退化時間序列的轉(zhuǎn)折突變點完成軸承的運行狀態(tài)的定量評估.

5 實驗結(jié)果分析

5.1 NASA數(shù)據(jù)實驗分析

實驗數(shù)據(jù)采用IMS 中心軸承試驗平臺［26～28］，如圖4所示.軸承轉(zhuǎn)速為2000 r/min，采樣頻率為20 kHz，采樣點數(shù)為20480，采樣時間間隔為10 min，共采集到984組樣本.圖5是其全壽命周期數(shù)據(jù)時域波形.

圖4 全壽命周期數(shù)據(jù)實驗平臺

圖5 全壽命周期數(shù)據(jù)時域波形圖

從圖5僅能看到趨勢的變化，無法清楚劃分軸承實時運行狀態(tài).為此，根據(jù)所提方法構(gòu)建高維特征集D，并繪制其歸一化特征圖，如圖6所示.

圖6 歸一化高維特征集

從圖6可看出，不同特征從不同維度表征軸承的運行狀態(tài).進一步計算各特征的Mon、Corr和Rob值，如圖7所示.計算綜合評價指標(biāo)Z，并按升序排序，選擇Z≥A（A取Z的均值為0.1092）的8 個有效特征構(gòu)建敏感指標(biāo)集F=［VARV，RMSV，TE，SE，RE，E，MF，SF］，如圖8所示.

圖7 退化特征Mon、Corr和Rob

圖8 有效退化特征選擇

為了綜合敏感指標(biāo)集F中的8 個特征的優(yōu)勢，首先，分別使用PCA 和KICA 對高維特征集D（樣本數(shù)為984 的8 維特征集）進行映射處理，分別得到得分矩陣S和轉(zhuǎn)換矩陣T（樣本數(shù)為984的8維矩陣）；然后，分別計算得分矩陣S和轉(zhuǎn)換矩陣T的第一分量與其余分量之間的馬氏距離，從而將984×8 的矩陣投影為高維空間984×984 的重構(gòu)矩陣；最后，選擇每個樣本空間的最大值（樣本數(shù)為984的1維向量）作為能全面反映軸承退化趨勢的退化指標(biāo)PCAMD和KICAMD.并利用K-medoids聚類方法將綜合敏感指標(biāo)集F融合為退化指標(biāo)K-medoids［29］，如圖9所示.

圖9 敏感退化指標(biāo)變化趨勢

從圖9不能直觀看出3個融合指標(biāo)之間的差別.因此，分別計算3個融合指標(biāo)與原始信號的信噪比和相關(guān)系數(shù)，如圖10 所示.從 圖10 可看出，KICAMD 相 比PCAMD 和K-medoids 信噪比分別提高了75.27% 和62.83%，相關(guān)系數(shù)分別增加了57.75%和43.58%，表明了KICAMD退化指標(biāo)的有效性.

圖10 融合指標(biāo)性能對比

但將KICAMD 退化階段進行局部放大（圖11）發(fā)現(xiàn)，其也存在明顯虛假波動.為此，分別引入線性回歸［19］、指數(shù)回歸［23］和GM 方法對其進行修復(fù)，得到修復(fù)后的退化指標(biāo)HI，如圖11中紅色曲線所示.

圖11 退化特征虛假波動修復(fù)

進一步計算HI 和KICAMD 的均方根誤差（RMSE）和相關(guān)系數(shù)，如表2 所示.對比表2 可知，GM 方法的RMSE最小，相關(guān)系數(shù)最大，論證了GM方法的有效性.

表2 HI和KICAMD的RMSE和相關(guān)系數(shù)

為定量評估軸承的退化狀態(tài)，將完全失效時間（第984 組）作為為終止故障時間（End Failure Time，EFT），進而計算得到HI 時間序列的5 個轉(zhuǎn)折突變點（532.5，699.5，800.5，900.5，965.5），如圖12所示.

從圖12 可看出，5 個轉(zhuǎn)折突變點將HI 分為呈梯度上升的6 個狀態(tài).由于轉(zhuǎn)折突變點在第532.5 組（非整數(shù)），為確定SFT 的取值，獲取第532 組和第533 組的Teager 能量算子解調(diào)頻譜圖，如圖13 所示.從圖13 可知，僅在第533 組的頻譜中發(fā)現(xiàn)與軸承外圈故障頻率234.4Hz 相近的頻率，且存在明顯的2～4 倍頻，由此確定SFT=533.

圖12 時間序列突變點監(jiān)測

圖13 Teager能量算子解調(diào)譜

確定SFT 后，根據(jù)HI 的轉(zhuǎn)折突變點繪制滾動軸承性能退化階段（后續(xù)狀態(tài)點均選擇了轉(zhuǎn)折突變點前一組數(shù)據(jù)），如圖14所示.從圖14可清楚觀察到軸承退化的階段過程：正常狀態(tài)［1，533）；早期故障［533，699）；輕微衰退［699，800）；衰退增長［800，900）；嚴重衰退［900，965）；失效狀態(tài)［965，984］.

圖14 軸承性能退化階段

為驗證所提方法的有效性，表3 對比了HI 及修復(fù)虛假波動后單一退化指標(biāo)的初始退化時間.從表3 可以看出，HI 的初始退化時間更早出現(xiàn)，結(jié)合圖13 也充分說明了其合理性.

表3 性能退化指標(biāo)SFT對比

5.2 XJTU-SY軸承數(shù)據(jù)實驗驗證

實驗所用的軸承加速退化測試平臺XJTU-SY 的平臺圖如圖15所示［30，31］.其采樣頻率為25.6 kHz，采樣間隔1 min，每次采樣時長為1.28 s，轉(zhuǎn)速為2400 r/min，外圈理論故障頻率為123.36 Hz.

圖15 加速退化測試平臺

圖16 是其外圈全壽命周期數(shù)據(jù)時域波形.從圖16不能直接觀察到2410 組之前是否發(fā)生故障.為此，構(gòu)建高維特征集D，并基于Mon、Corr、Rob 評價準(zhǔn)則得到綜合評價指標(biāo)Z，如圖17 所示.Z按升序排序后選擇了Z≥A（A=0.1082）的6 個有效特征，構(gòu)建敏感特征集F=［SE，VARV，RMSV，RE，E，MF］，如圖18所示.

圖16 全壽命周期數(shù)據(jù)時域波形圖

圖17 退化特征的Mon、Corr和Rob

圖18 有效退化特征選擇

采用與5.1 節(jié)相同的實現(xiàn)過程，獲得了經(jīng)GM 修復(fù)后的退化指標(biāo)HI，如圖19 所示.進一步計算得到退化指標(biāo)HI 時間序列的4 個轉(zhuǎn)折突變點（2344，2392，2464，2511），如圖20 所示.同時，為驗證SFT 的準(zhǔn)確性，呈現(xiàn)了如圖21 所示的第2344 組的Teager 能量算子解調(diào)頻譜.從圖21 的頻譜可知軸承發(fā)生外圈故障.因此，選擇SFT=2344是合理的.

圖19 GM虛假波動修復(fù)

圖20 時間序列突變點

圖21 第2344組Teager能量算子解調(diào)譜

在圖20 中，4 個轉(zhuǎn)折突變點將HI 分為呈梯度上升的5 個狀態(tài)，對應(yīng)軸承的5 個退化過程（圖22）：正常階段［1，2344）；早期故障［2344，2392）；輕微衰退［2392，2464）；衰退增長［2464，2511）；失效狀態(tài)［2511，2538］.

圖22 性能退化階段

6 結(jié)論

針對傳統(tǒng)退化指標(biāo)無法準(zhǔn)確反映滾動軸承全壽命周期內(nèi)退化狀態(tài)的問題，引入KICA 和馬氏距離融合多退化指標(biāo)，提出了一種特征融合與灰色回歸的滾動軸承性能退化評估方法，并通過實驗驗證了所提方法的準(zhǔn)確性和有效性，具體總結(jié)如下：

（1）基于Mon、Corr、Rob 的綜合評價準(zhǔn)則選擇有效退化特征，構(gòu)建了有效且全面地表征軸承性能退化狀態(tài)的敏感指標(biāo)集，為性能退化的準(zhǔn)確評估奠定基礎(chǔ)；

（2）利用GM 方法修復(fù)退化指標(biāo)的虛假波動，提高了滾動軸承狀態(tài)定量評估的可靠性；

（3）自適應(yīng)確定退化指標(biāo)HI的時間序列的突變點，實現(xiàn)了滾動軸承性能退化定量評估，克服了傳統(tǒng)初始故障時刻的確認方法依靠經(jīng)驗和實驗觀察的主觀性和局限性.

所提方法可以準(zhǔn)確且有效地評估軸承的退化過程，但進一步討論退化趨勢失效閾值自適應(yīng)判定方法，將更利于實現(xiàn)滾動軸承在線評估和制定更有針對性的維護策略.