劉勝蘭，高凌寒，杜劍維，劉 晨

(1.中國艦船研究院，北京 100101；2.重慶大學(xué)，重慶 400004)

0 引 言

艦船動力傳動系統(tǒng)是保障艦船可靠性航行和運(yùn)轉(zhuǎn)的重要動力系統(tǒng)，而滾動軸承是艦船動力傳動系統(tǒng)的主要功能部件，當(dāng)滾動軸承發(fā)生異?；蚬收?，可能會造成船舶運(yùn)行工況異常，所以對滾動軸承的壽命進(jìn)行預(yù)測十分必要。滾動軸承工作異常時，可通過對振動信號采集及分析有效地進(jìn)行滾動軸承壽命預(yù)測，而振動特征提取是滾動軸承性能退化狀態(tài)分類與殘余壽命預(yù)測的前提。振動信號中蘊(yùn)含了豐富的設(shè)備狀態(tài)征兆信息，如何從原始信號中全面有效地提取特征，一直是國內(nèi)外學(xué)者研究的熱點課題。常見的特征包括時域和頻域統(tǒng)計特征，時域統(tǒng)計特征有均值，方差，有效值，峭度，波形因子，脈沖因子等。頻域統(tǒng)計特征如平均頻率，中心頻率，故障特征頻率等。

雖然特征越多越可以全面反映設(shè)備的狀態(tài)信息，但過多的特征也會帶來大量的冗余信息，對后續(xù)計算帶來負(fù)擔(dān)，可能會造成維數(shù)災(zāi)難。怎樣從眾多特征中優(yōu)選出維數(shù)最少且敏感度高的最優(yōu)代表性的特征子集，對提高軸承狀態(tài)識別和殘余壽命預(yù)測準(zhǔn)確率及提升算法的效率都有著十分重要的意義，特征選擇流程如圖1所示。常用的方法有特征降維和特征選擇。

圖1 特征選擇流程Fig.1 Feature selection process

特征降維方法的思想主要是將原始的高維（多維）特征空間中的點向低維空間投影，新的特征空間維數(shù)更低，實現(xiàn)了特征的降維，但投影使原始特征的物理意義消失了。降維方法主要包括主成分分析（Principle Component Analysis，PCA），核主成分分析（Kernel-PCA）及線性判別分析（Linear Discriminant Analysis，LDA）等。

特征選擇方法是從d個特征中選擇個特征出來組成特征子集，舍棄其余的m-d個特征，所以特征的物理意義不會發(fā)生變化。但是在傳統(tǒng)的特征選擇方法中，需要預(yù)先設(shè)定選擇的特征維數(shù)m，并通過對比才能確定最終的優(yōu)選特征子集維數(shù)。因此本文提出一種自適應(yīng)順序最優(yōu)特征選擇方法，可以自動確定m。實驗表明所提方法優(yōu)選出的特征子集在軸承殘余壽命預(yù)測中有良好的效果。

1 振動信號特征提取

在滾動軸承運(yùn)行過程中，如果滾動軸承出現(xiàn)故障或異常情況，那么在振動信號中會有相應(yīng)的特征成分對其進(jìn)行反應(yīng)。基于數(shù)據(jù)驅(qū)動的方法對軸承進(jìn)行殘余壽命預(yù)測，第一步也是最重要的一步，就是從振動信號中提取出能夠全面反映滾動軸承在全生命周期中性能退化趨勢的特征?；谳S承振動信號處理理論，常用的統(tǒng)計特征包括時域特征，頻域特征。單一的或者僅僅是某個域的特征難以全面準(zhǔn)確地反映滾動軸承的運(yùn)行狀態(tài)以及捕捉故障信息。為了更好地反映滾動軸承的變化，引入IMF能量特征[1]，希爾伯特邊際譜特征和熵特征[2]。

由于滾動軸承的全生命周期振動信號通常是按固定的時間周期采集一段固定時長的信號，所以數(shù)據(jù)的特征提取也以時間段為單位進(jìn)行，即在每個時間段上提取相應(yīng)的特征向量，然后將每個時間段上的特征組合成特征矩陣其中i=1,2,…,m，m為時間段的個數(shù)，n為所選擇的特征數(shù)。

2 敏感特征選擇

在提取出來的眾多滾動軸承振動信號特征中，有一些特征相互獨立，對軸承的某種特定故障比較敏感，或者對軸承退化的不同時期比較敏感，可以很好地表征軸承的退化狀態(tài)，應(yīng)該重點利用。但是，還有一些特征是冗余或者與軸承退化狀態(tài)不相關(guān)，不敏感。如果使用這些特征，不僅會使殘余壽命預(yù)測模型更復(fù)雜，增加不必要的計算量，甚至?xí)档皖A(yù)測結(jié)果的可靠性。因此，不能將所有的特征直接合成特征集，輸入到殘余壽命預(yù)測模型中。需要通過特征選擇方法，將最具代表性，最有效的的特征篩選出來，剔除冗余，不相關(guān)的特征，進(jìn)而避免維數(shù)災(zāi)難，提升計算效率，減小預(yù)測結(jié)果的誤差。

2.1 特征性能評估

特征評價是很多模式識別任務(wù)中的重要組成成分，目前有許多可靠，高效，實際可行的方法，如核密度估計、皮爾遜相關(guān)系數(shù)，邊界寬度，F(xiàn)isher判別、基于距離，基于互信息理論等，但它們多用于分類任務(wù)，適用于故障識別領(lǐng)域。以上方法的核心思想是：不同類別的樣本在某種特征描述下具有不同的統(tǒng)計特征，根據(jù)可分離性明顯的原則，選擇出敏感的特征。以上方法是建立在樣本類別已知的情況下，所以也被稱為有監(jiān)督算法。然而，滾動軸承殘余壽命預(yù)測屬于回歸預(yù)測任務(wù)，滾動軸承的退化過程本質(zhì)上是一個連續(xù)變化的隨機(jī)過程。一個優(yōu)良的健康指標(biāo)應(yīng)該滿足下列條件：

1）特征隨滾動軸承性能的退化而變化，即與殘余壽命時間順序有一定的相關(guān)性。

2）退化是一個不可逆的單調(diào)過程，所以理想特征應(yīng)該具有一定的單調(diào)性。

3）對隨機(jī)噪聲具有良好的抗干擾能力，即具有一定的魯棒性。

本文采用文獻(xiàn)[2]提出的3種相關(guān)性，單調(diào)性，魯棒性的算法作為獨立評價函數(shù)，并在此基礎(chǔ)上提出自適應(yīng)順序最優(yōu)搜索算法，選擇出最終的優(yōu)選特征。對于特征信號順序，時間順序，表示在時間tk處取得的特征值，其中k=1,2,…,K，K代表時間長度。首先，使用指數(shù)加權(quán)移動平均法（Exponential Weighted Moving Average, EWMA）將特征順序分為兩部分，分別為平穩(wěn)趨勢項和隨機(jī)余量項:

EWMA的計算公式如下：

F的單調(diào)性評估指標(biāo)記作Mon(F)，公式如下[3]：

F的魯棒性評估指標(biāo)記為Rob(F)，公式如下[3]：

以上3個特征評估指標(biāo)的值域為[0,1]，值越大，說明特征的得分越高。在計算3個特征評估指標(biāo)之前，需要對輸入數(shù)據(jù)做歸一化處理，特征順序F被規(guī)范到[0,1]，時間順序被規(guī)范到[-1,0]。-1代表監(jiān)測起始時刻，0代表最終失效時刻。

2.2 自適應(yīng)順序最優(yōu)特征選擇方法

單獨的指標(biāo)只能片面反映備選特征在某方面的敏感性，為了綜合利用3個評估指標(biāo)，可以構(gòu)建一個加權(quán)的線性組合來融合3個評估指標(biāo)，并作為最終的獨立評價判據(jù)，具體表示如下：

其中：J為綜合判據(jù)；為備選特征集合；wi為每個性能評價指標(biāo)的權(quán)重,本文取w1=0.4,w2=0.3,w3=0.3。需要注意的是，計算出3個性能特征指標(biāo)后，需要根據(jù)最大值法歸一化到[0,1]，再進(jìn)行加權(quán)計算，避免因為量級的影響而偏向某個指標(biāo)。由于綜合判據(jù)J是線性的，和每個指標(biāo)正相關(guān)，所以J的值越大，表明該特征對殘余壽命預(yù)測越有效，應(yīng)該保留。

假如將d個特征獨立使用時的綜合評價指標(biāo)全部都計算一遍，然后按照綜合評價指標(biāo)的大小排序：

通常這種排序可以反映出單個特征的敏感程度，但對于組合特征就不一定了。文獻(xiàn)[5]直接選擇使J較大的前m個特征作為最優(yōu)特征組，但是一般來說，前m個最有效的特征組合成的特征子集并非是最優(yōu)的,排在前面的特征相互之間有可能是冗余的。

從d個特征中挑選m個，根據(jù)排列組合知識可知，所有可能的組合數(shù)為：

在m和d不大的情況下，可以使用窮舉法把所有特征組合的J都算出來加以比較，但如果特征維數(shù)較大，則窮舉法的計算量太大，顯然不太合適。本文在順序搜索算法的基礎(chǔ)上，提出一種自適應(yīng)順序最優(yōu)搜索算法來提高搜索效率，并且可以自動確定特征維數(shù)m，尋找最優(yōu)特征子集。算法從空特征集開始，每一步選擇一個最好的特征（即該特征與已入選的特征組合在一起時J值為最大），并將它添加到該集合中，然后用高斯白噪聲取代該特征，然后繼續(xù)搜索下一個特征，直到有效特征數(shù)不再增加，即最后一步新增的特征為高斯白噪聲。

假設(shè)已入選了k個特征構(gòu)成了一個大小為k的特征組Fk，把未入選的d-k個特征fi，j=1,2,…,d-k，和個高斯白噪聲特征，，按與已入選特征子集組合后的值大小排列：

則下一步的特征組選為

然后用一個高斯白噪聲覆蓋特征f1。在式（11）中，通過下列步驟實現(xiàn)（Fk+f）。

計算

得到

其流程如圖2所示。

圖2 自適應(yīng)順序最優(yōu)特征選擇Fig.2 Adaptive sequential optimal feature selection

綜上，本節(jié)提出了一種自適應(yīng)順序最優(yōu)特征選擇方法，相比直接選取得分最高的m個特征，本文的方法不僅考慮了特征之間的相互影響，去除了冗余特征，還能夠自動確定特征子集的維數(shù)，避免了通過試驗，組合不同維數(shù)的特征，最終確定m的取值，節(jié)省了特征選擇的時間。

3 基于相關(guān)向量機(jī)的軸承剩余壽命預(yù)測

相關(guān)向量機(jī)（Relevance Vector Machine, RVM）是Tipping[5]2000年提出的一種稀疏貝葉斯學(xué)習(xí)方法。它與廣泛使用的支持向量機(jī)具有相似的函數(shù)形式，都可以引入核函數(shù)從而解決非線性問題。但與基于結(jié)構(gòu)風(fēng)險最小化原則構(gòu)建學(xué)習(xí)機(jī)的SVM[6]不同，RVM是基于貝葉斯框架構(gòu)建學(xué)習(xí)機(jī)，所以它不僅輸出函數(shù)值，還同時輸出概率值。

整合模型參數(shù)后，目標(biāo)值的對數(shù)邊緣概率函數(shù)為：

其中：

其中：

4 實驗驗證

4.1 實驗平臺

本文將使用IEEE PHM 2012數(shù)據(jù)挑戰(zhàn)賽[3]的滾動軸承加速壽命試驗數(shù)據(jù)，實驗臺如圖3所示。2個型號為DYTRAN 3035B的單軸加速度傳感器正交放置在被測軸承的軸承座上，采集被測軸承的的振動信號。一個安裝在水平軸，一個安裝在垂直軸。采樣頻率為25.6 kHz，每次采集2 560個數(shù)據(jù)點，每次采樣間隔10 s。

圖3 PRONOSTIA滾動軸承試驗臺Fig.3 PRONOSTIA Rolling bearing test bench

4.2 實驗結(jié)果分析

使用bearing1-1與bearing1-3分別作為訓(xùn)練集與測試集，其水平振動時域如圖4所示。將水平方向和垂直方向的振動均用來提取特征，共計提取124維特征。其中包括36維時域特征，18維頻域特征，6維熵特征，16維IMF能量特征和48維希爾伯特邊際譜特征。

圖4 Bearing1-1和bearing1-3水平振動時域圖Fig.4 Bearing1-1AND bearing1-3 Horizontal vibration time domain diagram

通過自適應(yīng)序列最優(yōu)特征選擇方法優(yōu)選的特征子集共18維特征，記作FS1。從FS1中選取前10維記作FS2，選取前17維記作FS3。從特征中不重復(fù)的選擇綜合評價較高的特征與FS1組合為19維特征子集，分別記作FS4～FS7。

使用RVM模型預(yù)測，預(yù)測結(jié)果通過3個指標(biāo)評價，均方根誤差（RMSE），累積相對準(zhǔn)確度（CRA）和耗時。評價指標(biāo)的公式分別為：

為了直觀比較各特征子集在不同指標(biāo)上的表現(xiàn)，并且避免指標(biāo)尺度的影響，將各指標(biāo)歸一化，如圖5所示。各特征子集預(yù)測結(jié)果如圖6所示。

圖5 各特征子集預(yù)測結(jié)果指標(biāo)比較Fig.5 Comparison of prediction indicators of each feature subset

從圖5中可以發(fā)現(xiàn)，當(dāng)維數(shù)少于18維時，RMSE較小，但CRA較低，耗時較少。隨著特征維數(shù)增加到18維，CRA升到最高，耗時也相應(yīng)增加。當(dāng)特征維數(shù)為19維時，RMSE有所增大，CRA有所降低，耗時繼續(xù)增加。在圖6中，直線代表真實值，曲線代表預(yù)測值，區(qū)間代表95%置信區(qū)間，可以看到FS1的預(yù)測結(jié)果在中后期真實值基本都落入置信區(qū)間，預(yù)測值在真實值附近波動，由此可以說明18維確實是一個有效的最優(yōu)特征維數(shù)，該算法可以自動確定最優(yōu)特征子集及其維數(shù)。

5 結(jié) 語

針對滾動軸承殘余壽命預(yù)測特征提取與特征選擇的問題，本文通過在時域，頻域，時頻域提取大量特征，以便全面反映軸承的狀態(tài)信息。然后在序列搜索算法的基礎(chǔ)上，提出了自適應(yīng)序列最優(yōu)特征選擇方法，通過特征融合以及添加高斯白噪聲特征，可以自動確定特征子集的維數(shù)，尋找最優(yōu)特征子集，并通過實驗驗證了方法的有效性。

圖6 各特征子集預(yù)測結(jié)果Fig.6 Each feature subset prediction result