張衛(wèi)貞，曾建潮，石 慧，董增壽

(1.太原科技大學 工業(yè)與系統(tǒng)工程研究所,山西 太原 030024;2.中北大學 大數(shù)據(jù)學院，山西 太原 030051)

0 引言

隨著機械設備朝著大型化、自動化、精密化方向的不斷發(fā)展，故障發(fā)生的可能性和故障類型的復雜性也隨之增加，設備的突然故障可能導致整個生產(chǎn)流程中斷，產(chǎn)生重大經(jīng)濟損失，甚至危及人身安全。因此，對機械設備的關鍵部件進行實時監(jiān)測及剩余壽命的預測具有非常重要的現(xiàn)實意義?，F(xiàn)有剩余壽命預測方法主要包括：物理模型、基于專家知識的模型和基于數(shù)據(jù)驅(qū)動的預測模型等。對于復雜的機械設備，建立物理模型往往十分困難，所獲得的專家知識也不完備，因此，基于數(shù)據(jù)驅(qū)動的壽命預測方法逐漸受到重視[1-8]。

Huang等[9]利用支持向量機建立退化模型，用于預測剩余壽命;Chryssaphinou等[10]將部件的退化狀態(tài)建模為離散的半馬爾科夫過程;Si等[11]對關于數(shù)據(jù)驅(qū)動方面的回歸模型、比例風險模型、隨機濾波模型和隱馬爾科夫等剩余壽命預測模型進行了綜述。許多表征系統(tǒng)磨損或裂紋發(fā)展的單調(diào)退化過程被建模為Gamma過程[12-14]，但這些數(shù)據(jù)驅(qū)動的實時剩余壽命預測方法，往往都需要進行退化模型的假設與參數(shù)估計。假設的退化模型與實際模型之間往往存在較大差距，且參數(shù)估計的最優(yōu)化有可能收斂到局部最小卻不能保證全局最優(yōu)，因此預測模型并不能保證最終漸近收斂于真實的樣本模型。

核密度估計方法是一種數(shù)據(jù)驅(qū)動的方法，但該方法對數(shù)據(jù)分布的形式不作任何假定，是從數(shù)據(jù)本身出發(fā)研究數(shù)據(jù)分布特征的非參數(shù)估計方法[15-16]。核密度估計方法常用于分類中，Zhang等[17]提出一種對旋轉(zhuǎn)機械的故障類型進行區(qū)分的方法，對表征類型故障的特征，通過加入測試樣本前后概率密度相對熵之間的比較，來判斷屬于哪種類型的故障，其中利用核密度估計的方法求概率密度;李存華等[18]提出一種基于核密度估計的數(shù)據(jù)聚類分析方法，對基于網(wǎng)格數(shù)據(jù)重心分箱后的數(shù)據(jù)點進行核估計,來構造高效的聚類算法。也有研究將核密度估計用于動態(tài)模型，文獻[19-20]利用核密度估計對風速模型進行構建，有效降低了傳統(tǒng)研究需要假設風速服從某種已知分布導致的誤差。目前，關于核密度估計用于壽命預測的研究相對較少，Xu等[21]在基于Bayes的實時壽命預測中，采用實時退化特征信息的核密度估計來估計參數(shù)的先驗分布;王潔[22]提出了基于核密度估計的氣缸疲勞壽命預測，在已知N個氣缸樣本失效循環(huán)次數(shù)的前提下，利用核密度估計求任意第i個氣缸失效時循環(huán)次數(shù)服從的概率密度，再利用循環(huán)次數(shù)的可靠度和平均剩余壽命結果對核密度估計方法和Weibull分布兩種方法進行了比較。該方法是針對有多個同類設備歷史數(shù)據(jù)的情形，然而實時壽命預測中，越來越多的現(xiàn)代設備并沒有大量同類設備的歷史退化數(shù)據(jù)。

本文針對實時監(jiān)測系統(tǒng)，利用當前時刻的狀態(tài)監(jiān)測信息和歷史信息，提出一種基于核密度估計的實時剩余壽命預測方法。首先，針對基于固定窗寬的核估計求樣本概率密度時，容易造成樣本數(shù)據(jù)少的地方擬合不足，而樣本數(shù)據(jù)多的地方擬合過度的問題，改進為根據(jù)樣本數(shù)據(jù)的密度自適應地選擇窗寬值進行核密度估計，即高密度區(qū)域采用較大的核窗寬，而低密度區(qū)域采用較小的核窗寬，以提高核密度估計的準確性;其次，隨著實時監(jiān)測的進行，監(jiān)測到的樣本數(shù)據(jù)不斷增多，樣本的核密度估計也隨之不斷更新，采用傳統(tǒng)的核密度估計模型時，每新增一個樣本數(shù)據(jù)，基于這些樣本的核密度估計都要重新計算，這樣會造成歷史樣本不斷重復計算，計算量也越來越大，為避免實時監(jiān)測系統(tǒng)中樣本核密度估計不斷重復計算的問題，本文提出了核密度估計模型實時更新 遞推算法，進而實現(xiàn)對特征退化分布和實時剩余壽命的不斷實時更新;最后，采用IEEE PHM 2012的軸承全壽命數(shù)據(jù)對本文方法進行了驗證，并與基于Gamma分布[23]的剩余壽命預測方法進行了比較，驗證了本文預測方法的正確性和有效性。

1 核密度估計模型的構建

在實際應用中，隨著現(xiàn)代傳感技術的發(fā)展，越來越多的設備劣化狀態(tài)可以直接或間接獲得，研究者往往可以通過傳感器監(jiān)測到歷史狀態(tài)數(shù)據(jù)，選擇能表征部件連續(xù)退化的特征量，以更好地揭示退化部件的真實狀態(tài),滿足建模的需要。

設tk為當前監(jiān)測時刻，[0,tk]的監(jiān)測數(shù)據(jù)為當前接收到的設備退化數(shù)據(jù)，則相應設備特征可以隨監(jiān)測時間的退化趨勢得到，如圖1所示。設每單位時間采集一次隨機退化特征增量的樣本，ΔX1,ΔX2,…,ΔXk為k個抽樣于[0,tk]的獨立同分布的隨機退化特征增量樣本，將其服從的概率密度函數(shù)記為fk(Δx)，則未知密度函數(shù)fk(Δx)的核密度估計可表示為

(1)

式中：k為已知的隨機退化特征增量的樣本數(shù);hk為決定每個樣本貢獻度的平滑窗寬;K為核函數(shù)。樣本ΔXi對密度估計的貢獻度取決于核函數(shù)和所選擇的窗寬,即在給定樣本之后，核密度估計性能的好壞，取決于核函數(shù)K及窗寬hk的選取是否適當。

對于核密度估計與真實密度之間誤差的測量方法有許多，積分均方誤差(MISE)作為一種最易處理的全局測量方法被廣泛使用。

(2)

Silverman[24]基于積分均方誤差最小的思想，通過對不同核函數(shù)(Epanechnikov、高斯(Gauss)、三角(Triangle)等)的效率進行比較，認為不同核函數(shù)對積分均方誤差的偏差影響非常小。這里選擇實際中應用最為廣泛的高斯核函數(shù)(Gaussian Kernel)，用于核密度估計的模型。

(3)

(4)

將高斯核函數(shù)K(Δx)代入式(4)，初始最優(yōu)窗寬

(5)

式中σk為k個已知的初始隨機退化特征增量樣本的方差，

(6)

(7)

2 核密度估計模型的實時更新

2.1 自適應窗寬的確定及自適應窗寬下的核密度估計

由上述建模過程可以發(fā)現(xiàn)，假設初始窗寬hk為整個區(qū)間的固定窗寬，則隨著動態(tài)監(jiān)測系統(tǒng)時間的推移，當tk+1時刻獲得新的退化特征增量數(shù)據(jù)ΔXk+1時，k+1個樣本數(shù)據(jù)的固定窗寬核密度估計可表示為

(8)

(9)

對于動態(tài)實時監(jiān)測系統(tǒng)，當監(jiān)測到第k+1個樣本時，第k+1個樣本數(shù)據(jù)點處的自適應窗寬可表示為

(10)

設ΔX1,ΔX2,…,ΔXk是來自未知密度函數(shù)f(Δx)的初始樣本，當新增一個樣本數(shù)據(jù)ΔXk+1時，未知密度函數(shù)f(Δx)基于自適應窗寬的核密度估計表示為：

(11)

式中h(ΔXi)(i=1,2,…,k+1)表示不同樣本點處的自適應窗寬。

2.2 自適應核密度估計的實時更新

由于研究對象是實時監(jiān)測系統(tǒng)，每監(jiān)測到新的樣本數(shù)據(jù)，進行核密度估計時都要重新計算，這就造成每新增一個樣本數(shù)據(jù)，求新增樣本后的核密度估計時，其所有歷史樣本的核密度估計也需要重復計算。隨著樣本數(shù)據(jù)的不斷增多，計算量也會變得越來越大。為提高核密度估計模型的實效性，減少不必要的重復計算，通過ΔX1,ΔX2,…,ΔXk這k個初始歷史樣本的核密度估計遞推得到k+1個樣本，即監(jiān)測到ΔX1,ΔX2,…,ΔXk,ΔXk+1時的核密度估計，從而實現(xiàn)核密度估計的實時更新。具體利用式(11)可得：

(12)

式中h(ΔXk+1)為樣本點ΔXk+1處的的自適應窗寬值，由式(10)計算得出。

因此，當任意tk+j時刻新增j(j=1,2,3,…)個樣本時，k+j個樣本數(shù)據(jù)的核密度估計可通過遞推得到：

j=1,2,3,…。

(13)

這樣，實時監(jiān)測過程中每新增一個退化特征增量樣本，核密度估計都可以由其歷史退化特征增量的樣本自適應遞推得到，從而可有效避免實時監(jiān)測系統(tǒng)中的重復計算問題，提高核密度估計過程的效率。

3 特征退化分布的計算

(14)

當tk+1時刻新增一個樣本時，[0,tk+1]特征退化量的概率密度函數(shù)可表示為：

(15)

當tk+j時刻新增j個樣本時，[0,tk+j]時間累積退化量xk+j的概率密度函數(shù)可表示為：

j=1,2,3,…。

(16)

4 剩余壽命預測模型的建立

對于給定的失效閾值xth(設tk+t時刻累積退化量為xth時，系統(tǒng)失效)，如圖2所示。要實現(xiàn)剩余壽命的預測，首先要基于初始時刻到當前tk時刻的特征退化量X1:k(記Xk=X(tk)，X1:k={X1,…,Xk})，預測tk+t時刻的特征退化量Xk+t。設每單位時間監(jiān)測一次，增加一個新的樣本數(shù)據(jù)，則tk+t時刻有k+t個樣本數(shù)據(jù)。設T為設備在tk時刻的剩余壽命，則預測的剩余壽命的概率分布函數(shù)為

FT(t)=p(T≤t)=p(Xk+t≥xth)

(17)

式中g(xk+t)為[0,tk+t]特征退化量的概率密度，

(18)

將式(18)代入式(17)，可將預測的剩余壽命的概率分布FT(t)通過換元積分法化為

(19)

根據(jù)不斷更新的樣本，將來tk+t時刻k+t個隨機特征退化增量樣本的自適應核密度估計可表示為

(20)

因此，[0,tk+t]特征累積退化量的概率密度為

(21)

將式(21)代入式(19)，則tk時刻預測的設備剩余壽命T的概率密度函數(shù)為

(22)

在獲得新的退化特征增量數(shù)據(jù)后，可以根據(jù)剩余壽命預測的概率分布函數(shù)FT(t)，重新計算下一時刻的概率分布函數(shù)FT(t+1)，從而實現(xiàn)對預測的剩余壽命分布的更新。

5 實例計算與分析

本文利用IEEE PHM2012提供的軸承全壽命數(shù)據(jù)對模型進行驗證。該數(shù)據(jù)來源于FEMTO-ST研究中心PRONOSTIA試驗臺對滾動軸承的加速壽命試驗，其中振動信號的采樣頻率為25.6 kHz，每次采樣時長為0.1 s，采樣間隔為10 s，即每次采樣可得到2 560個樣本數(shù)據(jù)。本文以轉(zhuǎn)速為1 800 rpm,載荷為4 000 N工況下的Bearing 1-1的全壽命振動數(shù)據(jù)為例進行分析。

軸承故障診斷中，由于均方根(RMS)可以較好地反映軸承的磨損退化，在實際中得到了廣泛應用。本文以Bearing 1-1的均方根特征隨監(jiān)測時間的退化趨勢為例進行分析(如圖3)。由圖3可以看出，均方根隨時間基本呈現(xiàn)單調(diào)增加的趨勢，能較好地反映其退化趨勢，該軸承在t=2.749×104s時磨損開始加劇，且在t=2.803×104s時失效，均方根的失效閾值為5.607 mm/s2。

5.1 基于固定窗寬與自適應窗寬剩余壽命預測準確性的比較

窗寬選擇的好壞直接影響核密度估計的準確性。在實時剩余壽命預測中，退化特征增量樣本隨時間隨機變化，如果窗寬在整個區(qū)間上取固定值，核密度估計時易造成樣本數(shù)據(jù)少的地方擬合不足，而樣本數(shù)據(jù)多的地方擬合過度；若根據(jù)樣本數(shù)據(jù)的密度自適應地選擇窗寬值，則高密度區(qū)域采用較大的核窗寬，低密度區(qū)域采用較小的核窗寬，可以更符合實際樣本數(shù)據(jù)的需要，提高核密度估計的準確性，進而提高剩余壽命預測的準確性。

圖4和圖5所示為第2個監(jiān)測點處(tk=1.7×104s)和第6個監(jiān)測點處(tk=2.5×104s)基于固定窗寬與基于自適應窗寬剩余壽命(RUL)的概率密度(PDF)比較。

通過比較可以看出，基于自適應窗寬的剩余壽命預測結果相比較于基于固定窗寬的剩余壽命預測結果更接近實際的剩余壽命，且隨著監(jiān)測時間的增加、監(jiān)測數(shù)據(jù)的增多，兩種預測方法與實際剩余壽命之間的誤差進一步減小，基于自適應窗寬的剩余壽命預測結果的誤差相對更小，說明基于自適應的窗寬能夠?qū)Ω怕拭芏冗M行更好地估計，從而能更準確地對剩余壽命進行預測。

5.2 基于核密度估計的實時剩余壽命預測

設以監(jiān)測時間t∈[0,1.5]×104s單位時間隨機退化特征(RMS)的增量作為初始樣本，隨著系統(tǒng)運行時間的增加，接收到的監(jiān)測樣本不斷增多，基于核密度估計的剩余壽命的概率密度實時更新，剩余壽命的概率密度變窄變高，方差越來越小，說明預測的不確定性不斷減小，如圖6所示。此外，預測的剩余壽命值RUL越來越接近實際的剩余壽命值，其中，剩余壽命預測值通過平均剩余壽命(MTTF)[28]得到：

(23)

軸承的磨損過程是一個連續(xù)累積退化的過程，Gamma分布由于具有非負、增長、獨立增量的屬性，被廣泛用于磨損和裂紋擴展等逐漸累積損傷過程的退化建模中。為進一步驗證本文模型預測的有效性，對相同初始樣本下相同監(jiān)測點處的剩余壽命，采用基于Gamma分布的剩余壽命預測方法進行比較。圖7和圖8所示為在第1個監(jiān)測點(tk=1.5×104s)和第7個監(jiān)測點處(tk=2.7×104s)兩種模型得到的預測剩余壽命的概率密度。由圖7和圖8的比較可以看出，本文模型預測的剩余壽命概率密度的方差相對于基于Gamma分布的預測模型有所變小，對退化數(shù)據(jù)擬合的程度更高，且到第7個監(jiān)測點處時，本文方法預測的剩余壽命期望值已經(jīng)很接近實際的剩余壽命值。

為進一步對本文所提方法的預測效果進行評估，對不同監(jiān)測時間、實際剩余壽命、本文模型預測的平均剩余壽命以及基于Gamma分布預測的平均剩余壽命之間的比較，如表1所示。對比表中數(shù)據(jù)可以發(fā)現(xiàn)，隨著監(jiān)測時間的增加，監(jiān)測信息不斷增多，本文方法預測的剩余壽命值與基于Gamma分布的剩余壽命預測值相比，進一步減小，當監(jiān)測到足夠多的信息時，預測的剩余壽命逐漸變得很接近真實壽命，從而驗證了本文基于核密度估計的實時剩余壽命預測模型用于剩余壽命預測的有效性。另外，通過兩種模型預測所得剩余壽命與實際剩余壽命均方根誤差(RMSE1,RMSE2)的比較也可以看出，隨著監(jiān)測時間的增加，RMSE1,RMSE2均呈現(xiàn)逐漸減小的趨勢，且隨著監(jiān)測數(shù)據(jù)的增多，本文模型預測的結果與實際剩余壽命的誤差更小，說明本文模型預測的剩余壽命值更接近實際的剩余壽命值。

表1 兩種模型平均剩余壽命預測值(RUL)對比

6 結束語

在對許多機械設備剩余壽命預測時，很難得到大量破壞性試驗樣本數(shù)據(jù)，使得預測時采用傳統(tǒng)的退化模型假設、參數(shù)估計方法結果往往不夠準確的問題，本文提出一種基于核密度估計的實時剩余壽命預測方法，該方法不需要對數(shù)據(jù)分布的形式做任何假定，從數(shù)據(jù)本身出發(fā)研究數(shù)據(jù)分布特征。對于樣本的核密度估計，針對傳統(tǒng)的固定窗寬核密度估計會因樣本的疏密程度不同導致擬合不足的問題，改進為自適應窗寬核密度估計，以提高擬合優(yōu)度。此外，實時監(jiān)測過程中，隨著監(jiān)測信息不斷增加，針對核密度估計不斷重復計算的問題，建立了核密度估計的實時更新模型。實例分析結果表明：隨著實時監(jiān)測信息的不斷增多，預測剩余壽命的方差越來越小，預測準確度不斷提高，并通過與基于Gamma分布的剩余壽命預測模型的結果比較，進一步驗證了本文模型的有效性。可見，在不對退化數(shù)據(jù)分布做任何假設的前提下，該模型利用實時監(jiān)測信息可以較準確地預測被監(jiān)測設備的實時剩余壽命，從而可以為設備的實時預測維護提供有力的支撐，有效預防設備突然異常故障的發(fā)生，提高系統(tǒng)運行的安全性和可靠性。下一步，將在傳統(tǒng)核密度估計用于實時剩余壽命預測的基礎上，考慮隨機變量的有界性對核密度估計的影響，并將改進的方法用于實時剩余壽命預測中。