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

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

0 引言

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

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

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

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

1 核密度估計模型的構(gòu)建

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

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

(1)

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

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

(2)

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

(3)

(4)

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

(5)

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

(6)

(7)

2 核密度估計模型的實(shí)時更新

2.1 自適應(yīng)窗寬的確定及自適應(yīng)窗寬下的核密度估計

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

(8)

(9)

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

(10)

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

(11)

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

2.2 自適應(yīng)核密度估計的實(shí)時更新

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

(12)

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

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

j=1,2,3,…。

(13)

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

3 特征退化分布的計算

(14)

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

(15)

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

j=1,2,3,…。

(16)

4 剩余壽命預(yù)測模型的建立

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

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

(17)

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

(18)

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

(19)

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

(20)

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

(21)

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

(22)

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

5 實(shí)例計算與分析

本文利用IEEE PHM2012提供的軸承全壽命數(shù)據(jù)對模型進(jìn)行驗(yàn)證。該數(shù)據(jù)來源于FEMTO-ST研究中心PRONOSTIA試驗(yàn)臺對滾動軸承的加速壽命試驗(yàn)，其中振動信號的采樣頻率為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ù)為例進(jìn)行分析。

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

5.1 基于固定窗寬與自適應(yīng)窗寬剩余壽命預(yù)測準(zhǔn)確性的比較

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

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

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

5.2 基于核密度估計的實(shí)時剩余壽命預(yù)測

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

(23)

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

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

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

6 結(jié)束語

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