孟文俊 ， 張四聰， 淡紫嫣， 蔣 端， 劉 彈， 徐光華,2

(1. 西安交通大學機械工程學院 西安，710049) (2. 西安交通大學機械制造系統(tǒng)工程國家重點實驗室 西安，710054)

引 言

在機器運行過程中，一旦機械設備的某個零部件發(fā)生故障，往往會影響整個機械設備的性能，有可能造成巨大的經(jīng)濟損失。滾動軸承性能和可靠性對整個機械設備的可靠運行起著至關重要的作用[1]。隨著滾動軸承壽命理論研究的不斷發(fā)展，壽命預測方向已經(jīng)具有非常成熟的理論基礎和較豐富的實踐應用。根據(jù)滾動軸承壽命預測模型的不同，預測方法分為以下3種類型[2]：基于概率論的壽命預測方法、基于斷裂力學的壽命預測方法和基于模型的壽命預測方法。其中基于模型的壽命預測方法具有建模簡單、(易于操作的優(yōu)點，被廣泛應用。Gebraeel等[3]將歷史退化過程的性能數(shù)據(jù)采用反向傳播算法(back propagation，簡稱BP)的神經(jīng)網(wǎng)絡進行訓練，然后把當前運行部件的性能數(shù)據(jù)輸入神經(jīng)網(wǎng)絡中進行學習，通過與歷史數(shù)據(jù)的差異估計當前部件的退化軌跡，最后通過Bayes法對部件的剩余壽命進行預測，該方法對滾動軸承壽命的預測也適用。Sun等[4]利用支持向量機建立了軸承的壽命預測模型，進行了滾動軸承壽命的預測。申中杰等[5]以均方根值為特征指標，結合多變量支持向量機預測滾動軸承剩余壽命。Ali等[6]提出了基于簡化自適應諧振神經(jīng)網(wǎng)絡(simplified fuzzy adaptive resonance theory map，簡稱SFAM)神經(jīng)網(wǎng)絡和Weibull分布的壽命預測模型。文獻[7]將改進的指數(shù)模型用于滾動軸承壽命預測，可以有效減少預測誤差。王奉濤等[8]提出基于基于核主分量分析(kernel principal component analysis，簡稱KPCA)和威布爾比例故障率模型(Weibull proportional hazards model，簡稱WPHM)的方法，可以反映滾動軸承退化過程。

雖然基于模型的壽命預測方法占主導地位，但該方法依據(jù)樣本數(shù)據(jù)變化的特點進行建模，對樣本的完善程度和質(zhì)量有較高要求，且會受到各種假設條件的約束。其通過估計結構和參數(shù)確定的模型不會隨滾動軸承運行時間的增加而發(fā)生變化，是一種靜態(tài)模型，而運行性能可靠性分析是在滾動軸承運行過程中進行分析，是動態(tài)問題。本研究通過PCA方法對滾動軸承的多個性能指標進行融合，構建了能有效反映滾動軸承退化過程的衰退性能指標，進行壽命預測。通過相空間重構技術實現(xiàn)當前退化過程和歷史退化過程的對比，得到壽命預測值。同時，結合歷史失效時間進行統(tǒng)計推斷，得到更準確的平均壽命，并隨著觀測樣本的不斷積累實現(xiàn)平均壽命的動態(tài)更新。

1 滾動軸承衰退性能指標的建立

筆者選用PCA技術進行多特征融合，該衰退性能指標建立的流程如圖1所示。

圖1 衰退性能指標構造流程Fig.1 Construction procedure of fading performance index

1.1 滾動軸承時域和頻域特征提取

時域指標包括有量綱指標和無量綱指標，前者對早期故障較為敏感，但不穩(wěn)定；后者不受滾動軸承轉(zhuǎn)速和載荷的影響，對信號幅值和頻率的變化不敏感，會隨故障的發(fā)展敏感性下降。

頻域特征指標如表1所示，s(k)為信號x的頻譜，k=1,2,…,K，K為譜線數(shù)。fk為第k條譜線的頻率值，其中：頻率特征p1反映了頻域振動幅值變化；p2～p4，p6，p10～p13反映了頻譜的分散或集中程度；p5，p7～p9反映了主頻帶位置的變化。

1.2 滾動軸承相對特征指標的建立

為了減少不同軸承之間的差異，在特征融合之前對特征指標進行標準化和滑移處理。任意選取正常期內(nèi)一段趨勢平穩(wěn)的特征值，計算該段特征值的平均數(shù)，得到原始特征值與平均數(shù)之比，即相對特征指標。對相對特征指標進行7點滑移平均處理

(1)

其中：xRRX為原始特征與平均數(shù)的比值；xMA為最終相對特征指標。

表1 頻域指標

1.3 基于PCA的滾動軸承衰退性能指標的建立

PCA實質(zhì)是揭示多維數(shù)據(jù)之間的內(nèi)在關系，通過一種空間變換方式，實現(xiàn)高維數(shù)據(jù)能夠用少數(shù)主成分來描述[9]。其算法流程如圖2所示。

時域和頻域特征指標經(jīng)標準化和滑移處理后的相對特征指標組成的矩陣表示為

X=(X1,X2,…,Xι,…,Xp)T

(2)

對原始矩陣進行零均值處理，通過式(3)的線性變換得到新的坐標系Y1,Y2,…,Yp，將X的協(xié)方差矩陣的特征值λ由大到小排列，通過線性組合得到第1，第2，…，第p個主成分。將第1主成分作為衰退性能指標，用于滾動軸承的壽命預測

圖2 PCA分析算法流程Fig.2 Algorithm flow of PCA analysis

(3)

2 滾動軸承壽命預測

利用滾動軸承退化過程的相似性，對當前需要預測的退化過程與歷史退化過程進行相似性匹配，實現(xiàn)對未知過程的預測。

2.1 相空間重構技術簡介

相空間重構理論指出：一個n維系統(tǒng)有n個狀態(tài)變量，這n個狀態(tài)變量可以構成一個n維空間，這個空間稱為相空間。對非線性時間序列的建模和預測一般都可以在相空間中進行，在相空間重構中嵌入維數(shù)和時間延遲是影響相空間重構的重要參數(shù)。

計算以下3個變量

當E1(m)在某個特定值m0后不再發(fā)生變化或者變化很小，此時的即為時間序列相空間重構的最小嵌入維數(shù)。

2) 時間延遲的選取。時間延遲τ對時間序列進行相空間重構有著重要的影響，若τ過大，則相空間中的矢量變得隨機，重構后的矢量丟失大部分原動力系統(tǒng)信息；若τ太小，重構矢量之間相關性太強，無法很好地體現(xiàn)系統(tǒng)的演化規(guī)律。

如果S為原始數(shù)據(jù)序列{x(t)}(t=1,2,…,n)，Q為S的時間延遲數(shù)據(jù)序列{x(t+τ)}，則由Q和S可得到一個二維重構圖，互信息熵的計算公式為

(9)

其中：si和qi為二維重構圖中的點；Psq(si,qj)為當S=si,Q=qj時，重構圖中聯(lián)合分布概率；Ps(si)，Pq(qj)為邊緣分布概率，采用網(wǎng)格法對其進行計算。

當互信息熵I第1次達到極小值點時，此時的τ即為所求最佳時間延遲。

2.2 退化過程的對比

通過歷史退化過程的衰退性能指標序列建立歷史退化模型，以該模型為基準與當前退化過程進行對比，如圖3所示。

圖3 當前退化過程與歷史退化過程的對比Fig.3 Comparison between the current degradation process and the historical degradation process

1) 歷史退化過程在相空間中展開。對M組歷史退化過程的衰退性能指標序列進行相空間重構，實現(xiàn)滾動軸承退化過程的動力學軌道在相空間中完全展開。

2) 相空間中退化軌跡函數(shù)的學習。利用徑向基(radial basis function，簡稱RBF)神經(jīng)網(wǎng)絡對非線性函數(shù)優(yōu)良的逼近能力獲得非線性軌跡函數(shù)。假設第i組(長度為n)歷史退化過程的衰退性能指標序列經(jīng)過相空間重構后矢量為

xxik(t)={xi(tk),xi(tk+τi),…,

xi(tk+(mi-1)τi)}

(10)

其中：k=1,2,…，L；L=n-(mi-1)τi。

將重構后的矢量{xxik(t)|k=1,2,…,L}和對應的服役時間{tk|k=1,2,…,L}分別作為輸入和輸出進行RBF神經(jīng)網(wǎng)絡的訓練，得到第i組歷史退化過程的相空間軌跡函數(shù)。

3) 預測失效時間。假設當前需要預測的退化過程為第M+1次退化過程，對應的運行時間為TM+1，對衰退性能時間序列{xM+1(t1),xM+1(t2),…,xM+1(tn)}進行相空間重構，將重構矢量輸入到訓練完的RBF神經(jīng)網(wǎng)絡中進行學習，獲得1組估計的運行時間向量[T1,T2,…，TM]，計算該運行時間與實際運行時間偏差

ei=(Ti-TM+1)2(i=1,2,…，M)

(11)

在時刻TM+1處可以得到一個誤差向量e=[e1,e2,…,eM]，從而得到該時刻退化軌跡與歷史退化軌跡的相似程度

(12)

根據(jù)相似程度和歷史服役壽命可預測滾動軸承在當前時刻的失效時間

(13)

其中：tM+1為當前時刻預測的失效時間；ti為滾動軸承歷史服役壽命。

在下一次預測時，通過伸縮窗擴大性能指標數(shù)據(jù)，重復以上的步驟得到一個新誤差向量，將其與上次計算的誤差向量相加重新賦值給新誤差向量。

2.3 平均壽命的估計

為了提高預測的準確性，將歷史失效時間和當前預測的失效時間進行匹配作為樣本，采用RBF神經(jīng)網(wǎng)絡的方法擴充樣本建立概率模型，實現(xiàn)平均壽命的計算，如圖4所示。

圖4 失效時間的動態(tài)概率模型Fig.4 Dynamic probabilistic model of failure time

(14)

用t代替式(22)中的y，得到平均壽命

(15)

在下一次觀測時，將上次計算的平均壽命加入到歷史失效時間庫中，通過歷史樣本的不斷積累使分析樣本不斷擴大，估計不同觀測時刻的平均壽命。

3 試驗驗證

3.1 數(shù)據(jù)來源及說明

筆者采用“IEEE PHM 2012 Prognostic challenge”提供的滾動軸承壽命數(shù)據(jù)進行驗證。試驗裝置如圖5所示，滾動軸承具體參數(shù)見表2。在轉(zhuǎn)速為1.8kr/min、載荷為4kN的情況下采集7組滾動軸承退化失效過程的振動數(shù)據(jù)，采樣周期為10s，采樣頻率為25Hz，采樣時長0.1s。其中6組作為訓練，1組作為當前分析對象進行跟蹤和預測，該跟蹤軸承的時間長度為238min。

圖5 滾動軸承壽命實驗裝置Fig.5 Life test device of rolling bearings

NSK6804RS13/mm20/mm32/N4 000/(r·min-1)1 800/N2 470

3.2 衰退性能指標的建立

通過前文方法建立一個有效的衰退性能指標序列。7組滾動軸承的特征指標共7×27個。圖6中的均值、圖7中的偏斜度指標、圖8中的p3和p4特征指標穩(wěn)定性較差；圖7中的峰值指標和裕度指標雖然在滾動軸承退化過程中呈現(xiàn)了一定趨勢，但是總體上信息較嘈雜。因此剔除這些指標，將剩下指標進行PCA融合得到衰退性能指標，如圖9所示。

圖6 滾動軸承的有量綱時域指標Fig.6 Dimensional time domain index of rolling bearings

圖7 滾動軸承的無量綱時域指標Fig.7 Dimensionless time domain index for rolling bearings

圖8 滾動軸承的頻域指標Fig.8 Frequency domain index of rolling bearings

圖9 滾動軸承的衰退性能指標Fig.9 Fading performance index of rolling bearings

圖9表示通過PCA融合的衰退性能指標，可以看到融合后的衰退性能指標的退化趨勢更加明顯。圖中性能指標在180min存在一個明顯的波動拐點，一般將這個點對應的時間稱為首次預測時間(first predicting time，簡稱FPT)，并以此將軸承運行階段分為正常運行階段和失效階段[7]。在第1階段，軸承運行平穩(wěn)，一旦故障發(fā)生，軸承進入第2階段?？梢钥吹竭M入第2階段后衰退性能指標有一個先上升后下降的過程，這是因為當表面缺陷剛剛形成時，連續(xù)滾動會造成小的剝落或者裂紋，但是隨后會被連續(xù)滾動接觸撫平，當損傷擴展到更廣的區(qū)域時，振動水平再次上升[11]。這種現(xiàn)象被稱為“愈合”現(xiàn)象[12]。圖10中6組訓練軸承的衰退性能指標在進入失效階段后，都會有一定程度的波動，以軸承3最為明顯，說明軸承3在進入失效階段后，存在較為明顯的連續(xù)剝落或者裂紋損傷。

圖10 6組訓練滾動軸承的衰退性能指標Fig.10 Fading performance indexes of 6 groups of training rolling bearings

3.3 相空間重構參數(shù)的確定

利用Cao方法確定6組訓練滾動軸承衰退性能指標序列的嵌入維數(shù)m，如圖11所示。在應用中認為當E1變化很小時，對應的嵌入維數(shù)即為最小嵌入維數(shù)。利用互信息法確定這6組序列的時間延遲τ如圖12所示，當互信息熵I第1次達到極小值點時，相對應的時間延遲τ即是該序列相空間重構的最佳時間延遲τ。通過圖11和圖12看出，每組的嵌入維數(shù)和時間延遲相差不大。選取時間延遲的均值τ=15和6組中最大的嵌入維數(shù)m=10作為相空間重構的參數(shù)。

圖11 6組滾動軸承衰退性能指標序列的嵌入維數(shù)Fig.11 The embedding dimension of 6 sets of rolling bearing fading performance index sequence

圖12 6組滾動軸承衰退性能指標序列的時間延遲Fig.12 The time delay of 6 sets of rolling bearing fading performance index sequence

3.4 平均壽命的計算

為了驗證本研究通過多特征融合的衰退性能指標預測的壽命比單一特征指標(峭度指標)預測的壽命更加接近真實壽命，將兩組預測結果進行對比，如圖13所示。T為預測的平均壽命，T′為實際服役壽命，T′=238min。整體壽命誤差率為

(16)

圖13 滾動軸承平均壽命誤差對比圖Fig.13 Comparison of average life error of rolling bearings

由圖13可知，通過多個特征融合的衰退性能指標較單一特征指標預測的誤差更小，說明其更接近滾動軸承的真實壽命，同時隨著滾動軸承運行時間的增加，平均壽命越來越接近真實壽命。針對“IEEE PHM 2012 Prognostic challenge”，也有學者對此進行了算法研究，Wang[13]采用了基于軸承頻率特征的包絡分析。Sutrisno等[14]采用了3種方法，效果最好的是采用振動頻率特征異常檢測計算生存時間比。以運行到189min為例，相比本研究中的方法，這兩組方法得到的結果誤差分別為20.2%和19.8%，均大于圖13中的誤差。

4 結束語

本研究通過對滾動軸承在其壽命時間軸上間斷地性能指標的累積，實現(xiàn)在線服役滾動軸承壽命的動態(tài)預測。通過PCA技術對滾動軸承振動信號時域和頻域的多個特征進行融合，建立了一個衰退性能指標序列。利用相空間重構技術和RBF神經(jīng)網(wǎng)絡對歷史滾動軸承衰退性能指標序列進行相空間重構，獲得了歷史衰退性能指標序列的相空間軌跡函數(shù)。通過與當前退化過程進行相似對比，得到了一個預測的失效時間，最終利用這個預測的失效時間和歷史失效時間進行匹配組合作為統(tǒng)計樣本，建立了基于失效時間的概率模型，實現(xiàn)了平均壽命的估計。試驗結果表明，通過多個特征融合的衰退性能指標較單一特征指標預測的誤差更小，說明其更接近滾動軸承的真實壽命，驗證了方法的有效性。