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

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

引 言

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

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

1 滾動軸承衰退性能指標(biāo)的建立

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

圖1 衰退性能指標(biāo)構(gòu)造流程Fig.1 Construction procedure of fading performance index

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

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

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

1.2 滾動軸承相對特征指標(biāo)的建立

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

(1)

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

表1 頻域指標(biāo)

1.3 基于PCA的滾動軸承衰退性能指標(biāo)的建立

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

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

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

(2)

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

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

(3)

2 滾動軸承壽命預(yù)測

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

2.1 相空間重構(gòu)技術(shù)簡介

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

計算以下3個變量

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

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

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

(9)

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

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

2.2 退化過程的對比

通過歷史退化過程的衰退性能指標(biāo)序列建立歷史退化模型，以該模型為基準(zhǔn)與當(dāng)前退化過程進行對比，如圖3所示。

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

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

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

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

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

(10)

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

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

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

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

(11)

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

(12)

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

(13)

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

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

2.3 平均壽命的估計

為了提高預(yù)測的準(zhǔn)確性，將歷史失效時間和當(dāng)前預(yù)測的失效時間進行匹配作為樣本，采用RBF神經(jīng)網(wǎng)絡(luò)的方法擴充樣本建立概率模型，實現(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組作為訓(xùn)練，1組作為當(dāng)前分析對象進行跟蹤和預(yù)測，該跟蹤軸承的時間長度為238min。

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

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

3.2 衰退性能指標(biāo)的建立

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

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

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

圖8 滾動軸承的頻域指標(biāo)Fig.8 Frequency domain index of rolling bearings

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

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

圖10 6組訓(xùn)練滾動軸承的衰退性能指標(biāo)Fig.10 Fading performance indexes of 6 groups of training rolling bearings

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

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

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

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

3.4 平均壽命的計算

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

(16)

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

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

4 結(jié)束語

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