一種軸承核密度估計的剩余壽命預測研究

武沖鋒

（晉能控股陽泉固莊煤業有限公司，山西陽泉 045060）

0 引言

隨著機械設備不斷向智能化和精密化方向發展，軸承作為機械設備的關鍵部件，其性能退化直接影響整個設備的性能[1-4]。為此，對軸承進行實時監測和剩余壽命預測有著十分重要的研究價值和意義。現研究階段可將剩余壽命預測的方法大致分為基于物理失效模型、基于知識表示和數據驅動的方法[5]。面對復雜的設備，獲得物理模型是很困難的，知識表示方法更適合定性推理，不太適合定量計算，而且很難獲得完整的知識。然而，數據驅動的預測方法，通過對軸承性能數據的描述，提取有用的信息并預測剩余壽命[6]。Si等[7]將數據驅動方法的伽馬分布、回歸模型、Wiener過程和隨機濾波模型[8]等RUL模型進行總結分析。Pei等[9]提出雙時間尺度的隨機退化設備的非線性退化模型。Zhai等[10]提出基于高階隱半馬爾科夫模型的剩余壽命預測模型。上述數據驅動的方法大多需要假設退化模型和參數估計，而參數估計方法在模型選擇上存在局限性，過分依賴概率密度函數形式的先驗定義，因此不能保證預測模型的準確性和適用性。

本文提出一種軸承核密度估計的非參數剩余壽命預測模型，該方法不對數據分布附加任何假設，而是從數據本身研究數據分布的特點，避免了大多數數據驅動方法需要模型假設和參數估計的問題。該模型在核估計窗寬的選取上引入局部窗寬因子構建自適應窗寬模型。模型通過計算樣本點的密度來自適應選擇樣本點窗寬，提高了核密度估計窗寬選擇的可靠性，從而提高了預測精度。最后通過軸承磨損試驗驗證了所提模型的適應性和準確性。

1 核密度估計模型的構建

核密度估計方法不利用數據分布的先驗知識，不對數據分布進行假設，它是一種用于從數據本身估計未知變量的概率密度函數的方法[11-12]。其由已知的N個樣本點，通過選擇任意核函數及窗寬得到N個核函數，通過線性疊加得到核密度的估計函數。設Δx1,Δx2,…,Δx n為n個獨立且同分布的退化樣本，f(Δx)為其服從十五概率密度函數，則核密度估計的概率密度函數f^(Δx)表示為：

式中：h為窗寬；K(·)為核函數；n為樣本數。

其中，K(·)和h選用的恰當與否決定了核估計的準確度高低。核函數作為影響核密度估計的一個因素，一般情況下任何函數都可作為核函數，常用的有四次核、均勻核、三角核和高斯核。核函數的選擇雖多，但對核密度估計的準確度作用不大。本文選用廣泛應用的高斯核函數。即：

窗寬h作為影響核密度估計平滑性和核函數寬度的主要因素，當h較小時，核密度估計曲線不平滑且曲折，揭示了更多細節；當h較大時，核密度估計曲線平滑，但會覆蓋更多的細節。因此，選擇合適的窗寬對核密度估計是非常重要的。

2 自適應窗寬模型建立

目前自適應窗寬是窗寬選擇的主流方向，其本質是隨著樣本點的增加能夠自適應的選取窗寬，減少了計算量。在現有研究階段常用的自適應窗寬方法大多采用的是通過式（3）積分均方誤差求其最小值得到初始最優窗寬h n。

核密度估計模型中的初始最優窗寬h n，通過式（3）積分均方誤差求其最小值點得到：

將高斯核函數代入式（4）可求出h n為：

式中：σn為n個初始樣本特征退化增量的方差。

該窗寬方法解決了樣本數據實時變化時選擇窗寬和實際中固定窗寬造成的過度（不足）擬合問題。如果樣本分布是近似正態的，這種方法是最佳選擇。然而，當樣本分布不對稱或有多個峰值時，該方法會導致過度平滑，精度需要提高。

然而，在實時監測軸承的核密度估計過程中，為了提高窗寬對真實樣本數據分布的適應性，引入了一個反映基于最優窗寬h n的樣本點分散度的函數λi來自適應選擇窗寬。λi可通過式（6）各個樣本點處的概率密度得到：

將函數λi與式（5）計算出窗寬h n相乘可得自適應窗寬h i：

從而使得不同樣本的窗寬可以根據樣本數據的稀疏程度進行自適應調整，以更好地滿足實踐中核密度估計的需要。

因此，將自適應窗寬h(Δx i)方法引入核密度估計的模型中，從而構建自適應核密度估計的表達式為：

由于在實際應用中樣本數據都是實時更新的，如果每增加一個樣本都對其從頭開始計算，那么計算量會隨數據的增多變得復雜化。所以，為使核估計的計算性能得到提升，通過用已知的n個樣本的核估計推導第n+1個樣本的核密度估計，來實現核密度估計的實時更新。

3 實時剩余壽命

3.1 特征退化分布的計算

假設t n為軸承退化狀態時間，可得到[0,t n]時間內所有退化樣本的特征退化隨時間變化的趨勢可由圖1所示。為使設備的剩余壽命預測準確，由核密度估計可得[0,t n]隨機退化增量的概率密度，推出[0,t n]上隨機退化增量累計特征退化量的，可由n重卷積得到：

圖1 樣本特征退化隨時間變化的曲線

3.2 剩余壽命預測模型

設t n+t時刻，特征退化量達到（圖1）時軸承失效。剩余壽命預測的實現，要對當前t n時刻的剩余壽命預測，通過初始時刻到當前t n時刻的退化量（記，推出t n+t時刻的。設T為設備的剩余壽命，則剩余壽命的概率密度分布函數F T(t)為：

其中，[0,t n+t]特征退化量的概率密度為：

剩余壽命預測的概率密度為：

從而能夠推出t n時刻剩余壽命的概率密度函數為：

隨著實時監測的進行，監測的樣本數量不斷增加，樣本的核密度估計也不斷更新。當使用非實時壽命預測模型時，基于這些樣本的核密度估計必須為每一個新的樣本數據重新計算，這導致歷史樣本不斷被重新計算，計算量越來越大。為避免實時監測系統中樣本核密度估計的重復計算問題，提出了一種遞推核密度估計模型的實時更新算法，并對模型的特征退化分布和剩余壽命進行了連續的實時更新。

4 實例分析

本文使用IEEE PHM2012提供的全壽命軸承數據對該模型進行了驗證。該數據來自于FEMTO-ST研究中心的PRONOSTIA試驗臺的加速軸承壽命測試，振動信號的采樣頻率為25.6 kHz，采樣時間為0.1 s，采樣間隔為10 s，即每次獲得2 560個數據樣本。本文以Bearing 1-5在1 800 r/min和4 000 N載荷下的全壽命振動數據為例。在軸承故障診斷中，均方根（RMS）在實踐中被廣泛使用，因為它能更好地將軸承的退化狀態反映出來。圖2所示給出在時間變化下軸承的均方根特征退化趨勢。從圖中可看出均方根值基本上呈現出隨時間單調上升的趨勢，可以將軸承的退化趨勢很好地展現出來。該軸承在t=2.463×104s時出現故障，均方根的故障閾值為1.396 mm/s2。

圖2 特征值隨監測時間的退化趨勢

利用本文的模型預測軸承的剩余壽命，從圖3可以看出，在軸承運行時間不斷的變化下，接受到的樣本數據逐漸越來越多，在核密度估計的剩余壽命模型不斷的實時更新下，軸承的剩余壽命概率密度不斷的增大且逐漸變窄，其方差也明顯不斷減小，可看出預測值與實際值相差不斷減小，從而模型的準確性不斷提高。

圖3 不同監測時刻剩余壽命的概率密度

為了進一步評估所提方法的預測效果，表1所示為不同監測時間下實際剩余壽命和模型預測的平均剩余壽命的均方根誤差（Root mean square error,RMSE）比較。從表中可看出，軸承在運行過程的實時預測下，預測值與真實值之間的RMSE在逐漸的減小，根據樣本數據的增多，剩余壽命的預測值在不斷的接近真實值，且誤差也越來越小，從而可以驗證本文模型對軸承壽命預測的有效性。

表1 不同時刻的剩余壽命均方根誤差比較

為了能驗證本文模型在軸承應用上的優越性，將本文模型與Wiener過程的剩余壽命預測模型進行比較，Wiener過程多用于對具有非單調退化過程的設備進行建模。在相同條件下對軸承的退化過程進行壽命預測，得到圖4所示不同時刻下本文核密度估計模型與Wiener過程模型的剩余壽命預測值和真實值之間的比較。從圖中可看出，本文模型的預測結果相比Wiener過程模型，更加接近真實壽命值，隨著時間的變化預測的準確性不斷提高，從而驗證了本模型的適用性和精確性。

圖4 不同時刻兩模型剩余壽命值比較

5 結束語

本文提出了一種軸承核密度估計的非參數剩余壽命預測方法，以解決預測設備剩余壽命時在不同樣本分布下核估計的窗寬不能準確選擇的問題。在該方法中，通過引入局部密度因子來選擇核估計的窗寬，構建了一個自適應窗寬模型。該模型根據樣本點的密度自適應地選擇窗寬值進行核密度估計，即較小的窗寬用于樣本點密集區，而較大的窗寬用于樣本點稀疏區。隨著接收到的實時監測數據的實時更新，核密度估計也自適應地遞歸更新，從而避免每增加一個樣本數據就要進行一次計算帶來計算量復雜問題。實例分析表明，隨著樣本數據的增加，所提模型的軸承剩余壽命預測值越來越接近實際值，并于Wiener過程模型比較，驗證了本文模型的準確性和有效性。