綜合多階段退化和多不確定軸承剩余壽命預測

范樂賢，劉淑杰，張洪潮，2

(1.大連理工大學機械工程學院，遼寧 大連 116024；2.德克薩斯理工大學工業工程系，美國 德克薩斯 盧伯克 79409)

0 引言

滾動軸承是機械設備的關鍵部件之一，廣泛應用于各種設備中，它的退化或失效將直接影響設備的性能和可靠性，這使得盡早檢測和預測任何類型的潛在異常和故障并實施實時容錯操作以最小化性能降級和避免危險情況成為必要[1]。RUL預測作為預后與健康管理的關鍵環節，近年來有大量的學者研究。目前主流的RUL預測方法主要有三類：基于物理原理的方法、數據驅動的方法和混合方法[2]。其中，有限的失效機理研究、復雜的物理模型和實時的參數估計方法限制了基于物理原理方法的進一步應用。數據驅動方法可分為兩類：人工智能方法和基于退化趨勢的方法[2]。前者從大量同類設備的狀態監測(condition monitoring，CM)數據中提取深層特征，使用智能算法學習深層特征和使用壽命之間聯系，進一步實現RUL預測[3-4]。后者從監測信號中提取健康指標(Health Indicator，HI)，根據服役環境、退化過程、專家知識和可用數據對退化趨勢進行建模，實現對退化趨勢的預測[5-6]，預測結果可以表示為隨機、概率或確定性變量，近年來基于退化趨勢的數據驅動方法在滾動軸承RUL預測中得到了廣泛的應用。

隨著故障程度的發展，軸承的HI一般呈現出變化的退化趨勢，單一的退化模型很難描述隨時間變化的退化過程。因此，有必要在RUL預測之前根據趨勢變化將退化過程劃分為不同階段并建立對應的退化模型[7]。Lim等使用三個動態模型描述軸承的三種退化過程，并通過貝葉斯算法計算不同模型的概率劃分退化階段[8]。Cui等用三個多項式函數表示軸承的正常、緩慢退化和加速退化，通過設置相對誤差閾值和計算狀態概率判斷退化階段[9]。Wang等將退化過程分為正常和退化階段，分別用不同的函數形式描述退化階段，并用基于3σ區間的檢驗方法確定退化起始[10]。Ahmad等使用線性和二次動態回歸模型描述健康指標的兩種趨勢，并用警報界限技術確定退化起始點[11]。此外，HI的退化趨勢不僅具備多階段的特點，往往還受制作工藝、服役時間和工況環境等因素的影響，導致每個階段的退化路徑具備不確定性變化，因此對退化路徑不確定度進行有效管理同樣是廣泛關注的問題。Hu等建立了退化速度隨運行時間和外部因素變化的線性維納模型[12]。Lei等建立一種同時考慮退化時變性、個體差異性、非線性可變性和測量可變性四種不確定性源的隨機過程模型[13]。Wen等在冪函數模型中引入導致退化路徑不確定性的三種因素[14]。Gao等提出一種能夠擬合多級線性退化同時考慮個體差異的隨機模型[15]。

綜合以上基于退化趨勢的軸承RUL預測方法，可以發現考慮多階段退化和趨勢不確定性的預測方法是相互獨立的。然而，根據軸承的實際退化情況，HI趨勢是同時具備多階段和多種不確定性特征的綜合過程，目前還沒有將這兩種特征結合的預測方法。因此，本文提出一種能夠擬合多階段退化趨勢，同時考慮趨勢多不確定性的RUL預測方法。

1 退化建模和RUL預測

1.1 退化建模

軸承時域振動信號由隨機噪聲和周期性波動組成。在振動數據的時域分析中，可以提取具有良好演化趨勢的時域特征作為HI，而均方根(Root Mean Square，RMS)是具有良好上升趨勢特征的代表，因此本文采用RMS作為HI表示軸承的性能退化。根據文獻[7]的分析，軸承在全壽命周期會依次經歷健康、緩慢退化和加速退化三個階段，如圖1所示。

圖1 滾動軸承的退化過程

因此基于維納過程采用常數、線性和非線性形式描述三個階段HI的變化。即：

x(t)=x0+σB(t)

(1)

x(t)=x0+ηt+σB(t)

(2)

(3)

式中，η是漂移系數，表示軸承的退化率；σ是擴散系數，表示同類軸承的共同特性；x0是每個階段的初始狀態；Λ(t；θ)是非線性退化函數，本文采用冪函數形式Λ(t；θ)=tθ。

y(t)=x(t)+ε

(4)

實際中，CM數據通常為離散形式Y1：k={y1，y2，…，yk}。同樣地，X1：k={x1，x2，…，xk}表示對應CM數據的退化狀態。因此，基于上述多階段退化過程及不確定性的定義，建立具有三種不確定性的隨機退化模型。

常數模型：

(5)

式中，zk=[xk]，A=[1]，Q～N(0，σ2)，H=[1]，R～N(0，δ2)。

線性和非線性模型：

(6)

1.2 RUL預測

在工程應用中，將軸承HI首次超過失效閾值的使用時間視為有效壽命是很自然的，則軸承RUL可以看作是當前時間和首達時之間的時間間隔，其中，首達時是隨機過程首次達到某值的時刻。根據有效壽命定義和首達時的概念，tk時刻的RULLk表示如下：

Lk=inf{x(l+tk)≥ω|x(tx)<ω}

(7)

則上文提出的退化模型對應的RUL的條件概率密度函數(Probability density function，PDF)為：

(8)

Γk=Π(lk+tk;θ)-Π(tk;θ)，Π(t;θ)為HI退化形式。上式詳細證明參見文獻[18-19]，此處省略。然而，利用上述分布在線預測RUL時，不僅需要失效閾值ω，還需要估計狀態zk|k的期望、方差以及SSM中的相關參數。

2 自適應參數估計和階段劃分

2.1 參數估計

然后構造參數向量Θ的聯合似然函數關于CM數據集Y和狀態集Z的條件期望：

[(yi-Hzi|k)(yi-Hzi|k)T]}

(9)

(10)

假設(10)式中

(11)

2.2 階段劃分

SPC是一種利用數理統計控制過程和減少波動的統計技術，可以實現對過程進行及時有效的調整和決策。由于退化過程的不穩定性，CM數據偶爾會出現偏離正常范圍的“奇異值”，這種現象會對階段劃分產生一定的干擾，從而導致階段的錯誤估計。為了準確識別退化階段，避免“奇異值”引發階段錯誤估計，提出基于SPC的階段劃分方法，通過監測殘差變化識別退化階段。

由于退化時變和測量噪聲等不確定性因素的存在，即使在同一階段內，退化狀態和HI也并不保持恒定，而是在一定范圍內波動。這里將殘差定義為退化狀態的期望E(x1：k)和HIY1：k之間的差值，即k=yk-E(xk)。通過對建立的退化模型分析可知，殘差k應服從則使用EWMA控制圖監測k的動態變化，其上下控制限(Up/Low control limit，UCL/LCL)為：

(12)

式中，K為控制限寬度；λ為EWMA權重；μ為EWMA統計量的均值。

若Ei在控制限內，則認為退化階段不發生變化是合理的。否則，有理由判斷退化階段發生改變。然而，僅用一個超過控制限的殘差進行階段劃分仍然是不準確的。因此，根據統計過程控制的判穩準則，計算不同控制限下的最小采樣次數，如表1所示。

表1 不同控制限下滿足穩定判據的最小采樣次數

3 實例研究

圖2 軸承生命周期的振動數據

為了驗證該方法的有效性，本節給出一個滾動軸承的實例，CM數據來自XJTU-SY軸承數據集[20]的Bearing3_1。圖2顯示出了所采集的全壽命振動信號，本文從500min開始對軸承進行狀態監測，并使用相對法確定軸承失效閾值[20]。

在整個監測過程中，從CM數據提取的RMS變化如圖3中的左小圖曲線所示。將提取的RMS帶入退化模型估計退化狀態，結果如圖3所示。圖3左側子圖為t=1000～1500min的局部放大圖，圖3右側子圖為2300min后的局部放大圖。從數據的總體趨勢來看，前2300min軸承處于健康階段，估計狀態在一定范圍內波動，沒有明顯的退化趨勢。2300min后，軸承進入退化階段，估計狀態在2340～2450min間呈現線性退化趨勢，在2450min后呈現非線性退化趨勢。當軸承運行到2490min時，退化狀態首次達到失效閾值，其壽命被認為終止。在500～2490min監測期間共采集1990組CM數據，使用提出的參數估計方法根據這些數據估計不同階段下的參數向量θ，如圖4所示。

圖3 基于RMS的狀態估計結果

使用提出的參數估計方法在每個CM點獲取參數向量θ。為保證模型參數預測精度，待處理數據累積到5組再進行估計，然后每次采樣新數據后更新估計結果。不同階段的估計結果如圖4所示，圖4(a)、圖4(b)、圖4(c)分別為健康階段、緩慢退化階段和加速退化階段的結果。從圖4(a)、圖4(b)可以看出，隨著CM數據的積累，模型參數可以快速收斂并保持穩定。在每個階段結束時(t=2342min和t=2462min)，模型參數都發生一定程度的變化，這與相應時刻發生階段變化的實際情況一致。從圖4(c)可以看出，估計參數在一定范圍內波動，這也符合加速退化階段內軸承退化嚴重、運行狀態不穩定的實際情況。這些結果表明了本文提出的參數估計方法具有自適應估計能力。

圖4 每個階段的參數估計過程

圖5 階段劃分結果

計算估計狀態和RMS的殘差，得到每個監測點EWMA統計量Ei，結果如圖5實線所示，其中(a)、(b)分別表示健康階段和退化階段，虛線為控制限K=3求得的上下控制限。

圖5中的統計量Ei大部分在控制限內，少數(小圓圈)發生“越界”，由穩定判據可知上述“越界點”是該階段內的奇異值，表明建立的模型符合退化階段的實際情況。但在t=2342min和t=2462min后，短時間內出現大量的“越界點”，相鄰“越界點”間的距離不足穩定判據中規定的最小采樣次數，此時有充分的理由認為退化階段發生了變化。將改變點后的監測數據帶入相應的隨機退化模型中，統計量Ei又恢復到控制限內，說明所建模型又符合變化后的退化情況。上述現象表明，所提階段劃分方法可以有效劃分不同階段的監測數據，并及時確定階段變化點。

圖6 退化階段估計RUL的PDF

基于上述狀態、參數及階段估計結果，根據(8)式計算退化階段內各監測點的RUL分布，結果如圖6所示。在每個退化階段初期，由于缺少足夠的監測數據，估計參數不確定性較大。因此，PDF曲線多為“扁平”形狀，可能的RUL范圍較廣，但PDF分布仍覆蓋了真實的RUL。隨著監測數據的積累，特別是在估計參數收斂后，PDF曲線越來越“尖銳”，可能的RUL范圍變窄，這表明預測結果的不確定性降低，可以為后續決策提供有效的信息。

以PDF最大值代表各時刻的預測RUL，預測結果如圖7所示?？梢钥闯?，退化初期CM數據較少，導致估計參數不確定性較大，預測結果誤差較大。隨著CM數據的不斷積累，估計參數逐漸收斂并接近真實值，RUL預測結果也逐漸趨于穩定，預測精度上升。但此時根據CM數據只能判斷出軸承具有線性退化趨勢，不能預測未來的非線性退化趨勢。因此只能預測緩慢退化對應的RUL，如圖7(a)虛線所示，導致真實RUL和預測RUL有一定的差別。倘若只考慮軸承緩慢退化，即只考慮軸承線性退化部分，如文獻[9]。圖8中顯示了只考慮軸承線性退化時，文獻[9]和本文預測的均方根誤差比較。

圖7 退化階段的預測RUL

可以看出，如果只考慮軸承線性退化部分，本文預測的均方根誤差與文獻[9]相比更小，表明該方法的預測精度比文獻[9]有所提高。此時，文獻[9]認為軸承失效，而根據“相對法”設定的失效閾值更接近真實失效閾值。因此，本文進一步預測下一個時期的RUL，結果如圖7(b)所示?？梢钥闯霎斴S承進入最后的加速退化階段，由于退化階段不再發生變化，根據CM數據預測的退化趨勢即是最終趨勢，因此預測精度大大提高，預測RUL皆在30%的誤差限度內。隨著CM數據的積累，模型參數逐漸收斂，RUL預測也進一步靠近真實值。這些結果表明，所提方法在預定失效閾值下，能夠預測全壽命周期內的RUL。

圖8 只考慮軸承緩慢退化時預測RUL的均方根誤差

4 結論

為了描述軸承健康指標的退化趨勢同時具有多階段變化和多不確定性的特點，本文提出一種新的隨機退化模型，該模型同時考慮個體差異性、退化時變性、測量可變性和多階段退化。為了準確劃分退化階段，提出一種基于EWMA控制圖的階段劃分方法，用于自適應切換退化模型。同時，在缺少同類軸承先驗信息的情況下利用基于EM算法的參數估計方法在線更新模型參數。最后，通過對仿真數據和XJTU-SY數據集的實證研究，驗證該方法的有效性。結果表明，該方法能夠準確劃分滾動軸承的不同階段，并根據監測數據提供不同退化階段的RUL預測。