張少釗，胡昌華，張正新，周志杰，張建勛

（第二炮兵工程大學(xué)，陜西 西安 710025）

0 引 言

近年來，預(yù)測與健康管理（prognostics and health management，PHM）[1]作為一種新興技術(shù)，在智能電網(wǎng)、核電站、電力工業(yè)、航空航天、艦船維護(hù)和公共健康管理等領(lǐng)域得到廣泛應(yīng)用，有效地提高了設(shè)備的安全性和可靠性，降低了設(shè)備的維護(hù)保障費(fèi)用和失效事件發(fā)生的風(fēng)險(xiǎn)。預(yù)測是PHM的重要組成部分，其中剩余壽命（ remaining useful lifetime，RUL）的預(yù)測為優(yōu)化設(shè)備檢測與維護(hù)策略，合理訂購備件提供了科學(xué)的依據(jù)，因此也是當(dāng)前可靠性領(lǐng)域研究的重點(diǎn)和熱點(diǎn)。

剩余壽命預(yù)測的方法主要包括基于物理模型的方法、數(shù)據(jù)驅(qū)動(dòng)的方法及兩種方法的融合[2]。當(dāng)系統(tǒng)的物理失效模型容易建立時(shí)，基于物理失效模型的方法最為準(zhǔn)確；但更多的情況下，系統(tǒng)較為復(fù)雜，物理失效模型難以建立，其應(yīng)用也受制于此?；跀?shù)據(jù)驅(qū)動(dòng)的剩余壽命估計(jì)方法由于其適用性更廣，得到了廣泛的研究和關(guān)注。近年來，由于技術(shù)的進(jìn)步和產(chǎn)品開發(fā)周期的縮短，設(shè)備的可靠性明顯提高，失效數(shù)據(jù)也越來越少，使得基于退化數(shù)據(jù)的RUL估計(jì)方法得到了更多的青睞。其中，基于隨機(jī)過程的方法由于更能體現(xiàn)退化過程的隨機(jī)性和動(dòng)態(tài)特征，逐漸成為退化建模的主流[3]，這類方法主要有Wiener過程[4]、Markov 鏈[5]、Gamma 過程[6]、逆高斯過程[7]等。 Wiener過程是一種特殊的擴(kuò)散過程，能夠?qū)Ψ菃握{(diào)退化過程進(jìn)行建模[8]，其在首達(dá)時(shí)間意義下具有解析形式的壽命分布，得到了廣泛的應(yīng)用。選擇不同形式的漂移系數(shù)，過程可以對線性和非線性退化過程進(jìn)行建模。Lu和Peng等[9-10]對線性情況下的維納退化過程做了研究。Si等在文獻(xiàn)[11]中提出了一種基于擴(kuò)散過程的非線性建模方法，將線性退化下的結(jié)果包含為一種特例，并在首達(dá)時(shí)間意義下，導(dǎo)出了剩余壽命的分布函數(shù)。

在實(shí)際的退化過程中，設(shè)備普遍受到周期性變化應(yīng)力的影響。如設(shè)備存貯的環(huán)境溫度應(yīng)力在一年中會(huì)出現(xiàn)周期性的變化，某些設(shè)備的工作載荷也具有周期性。周期應(yīng)力是汽輪機(jī)軸累積損傷的重要因素。然而，目前對這類情況下的退化建模和RUL估計(jì)的研究較少。

基于此，本文針對設(shè)備的退化過程中應(yīng)力的作用符合周期性變化的情況，通過在擴(kuò)散過程的漂移系數(shù)中引入周期函數(shù)，考慮這種情況下設(shè)備的退化建模和RUL估計(jì)問題。首先在考慮周期應(yīng)力作用下，基于擴(kuò)散過程建立了非線性退化模型，在首達(dá)時(shí)間的意義下，推導(dǎo)出了剩余壽命的分布，并給出了參數(shù)的極大似然估計(jì)。通過數(shù)值仿真和統(tǒng)計(jì)分析驗(yàn)證了剩余壽命分布的正確性和本文所提剩余壽命估計(jì)方法的有效性。

1 問題描述

圖1是某型電容器耐壓強(qiáng)度隨使用次數(shù)下降的退化曲線，可以看出，電容器的工作電壓作為其退化過程中的主要應(yīng)力，具有明顯的周期性特征。

利用退化過程建模是將退化過程描述為一種隨機(jī)過程{X（ t）；t≥0}，其中 X（ t）是 t時(shí)刻的隨機(jī)退化信號(hào)，由確定部分和隨機(jī)部分所組成。確定部分用來描述系統(tǒng)的固有特征和共性物理特征，而隨機(jī)部分用來描述系統(tǒng)退化過程中存在的個(gè)體差異以及固有的動(dòng)態(tài)性和差異性。

在建立退化模型后，要進(jìn)行剩余壽命估計(jì)，首先需要對設(shè)備的壽命進(jìn)行定義。將設(shè)備失效的時(shí)間定義為隨機(jī)退化過程{X（t）；t≥0}首次達(dá)到失效閾值l的時(shí)間（首達(dá)時(shí)間），因此，設(shè)備壽命T定義為

圖1 某型電容器樣品試驗(yàn)數(shù)據(jù)示意圖

假設(shè) X（ 0）=0，失效閾值 l＞0。 因此，進(jìn)行剩余壽命估計(jì)的關(guān)鍵在于確定首達(dá)時(shí)間，將設(shè)備的退化過程描述為一類擴(kuò)散過程{X（ t）；t≥0}。

式中t是模型的時(shí)間尺度，當(dāng)退化過程為線性退化過程時(shí)，t是采樣時(shí)間點(diǎn)；B（ t）為標(biāo)準(zhǔn)布朗運(yùn)動(dòng)，即 B（ t）～N（ 0，t），N（ 0，t）是均值為 0，方差為 t的正態(tài)分布；μ（ t；θ）為漂移系數(shù)，與設(shè)備所受應(yīng)力有關(guān)；σB為擴(kuò)散系數(shù)，由設(shè)備本身存在的不一致性與不穩(wěn)定性、測量設(shè)備測量誤差及穩(wěn)定性、測試過程中的外部干擾等隨機(jī)因素決定；μ（ t；θ）是時(shí)間 t的非線性函數(shù)，θ為未知參數(shù)。式（2）將退化過程就描述為非時(shí)齊的擴(kuò)散過程。在設(shè)備的退化過程中，受到的環(huán)境或應(yīng)力影響可能是不斷變化的，μ（t；θ）是隨時(shí)間的非線性函數(shù)，因此可以描述設(shè)備所處環(huán)境影響隨時(shí)間不斷變化的過程。 并且，當(dāng) μ（ t；θ）=μ 時(shí)，該模型就是一種線性漂移的擴(kuò)散過程，這也說明了式（2）更加具有一般性和靈活性。

考慮用式（2）描述周期性應(yīng)力作用下的退化過程。因?yàn)閄（t）是隨機(jī)過程，顯然壽命T也是一個(gè)隨機(jī)變量，設(shè)壽命T的壽命分布函數(shù)為FT（t），概率密度函數(shù)（probability density function，PDF）為fT（t）。為了實(shí)現(xiàn)壽命估計(jì)，關(guān)鍵需要解決以下兩個(gè)問題：

1）退化模型的參數(shù)估計(jì)；

2）基于退化模型的剩余壽命估計(jì)。

2 周期應(yīng)力作用下的非線性退化模型與壽命估計(jì)

2.1 退化建模

如何對周期性的應(yīng)力進(jìn)行描述是退化建模過程中的關(guān)鍵問題。在對具有線性特征的退化數(shù)據(jù)進(jìn)行剩余壽命估計(jì)時(shí)，通常采用線性漂移的擴(kuò)散過程-Wiener過程的退化模型進(jìn)行建模。但是，周期應(yīng)力作用下的退化是一種非線性的過程；并且，在工程實(shí)際中非線性的退化過程更具有普遍性，因此考慮將基于Wiener過程的退化模型推廣到非線性的情況。

Wiener過程是帶漂移的布朗運(yùn)動(dòng)。在實(shí)際的退化過程中，真實(shí)的漂移系數(shù)和擴(kuò)散系數(shù)往往不是常數(shù)，不能通過直接觀測得到；因此，通過擴(kuò)散過程來描述退化過程中的潛在波動(dòng)在理論上是可行的。

由于μ（t；θ）表示退化過程中的漂移系數(shù)，是一種應(yīng)力作用的結(jié)果。因此，當(dāng)應(yīng)力的作用出現(xiàn)周期性特征時(shí)，可以將一類滿足式（2）的擴(kuò)散過程X（t）描述為

式中，參數(shù)A用來描述退化過程的線性趨勢，（B/ω）sin（ωt）用來描述一種周期性的波動(dòng)，反映了設(shè)備受到的環(huán)境等周期應(yīng)力作用下的退化軌跡。

2.2 參數(shù)估計(jì)

本文采用一種基于歷史退化數(shù)據(jù)的參數(shù)估計(jì)方法，通過參數(shù)極大似然估計(jì)對模型的參數(shù)進(jìn)行估計(jì)。

假設(shè)對n個(gè)樣品進(jìn)行性能退化測試，對第i個(gè)樣品的測量時(shí)間為ti1，…，tim，并且在初始時(shí)刻ti0時(shí)的性能退化量Xi0=0，樣品退化量的測量值為Xi1，…，Xim。記ΔXij=Xij-Xi（j-1）為樣品i在時(shí)刻ti（j-1），tij之間的性能退化量,測量時(shí)間間隔Δtij=tij-ti（j-1）。

由Wiener過程的性質(zhì)可知：

假設(shè)各不同退化設(shè)備之間的退化量測量值具有獨(dú)立性，由式（ 4）可以得到未知參數(shù) Θ=[A，B，σB2，ω]′的似然函數(shù)

未知參數(shù)的對數(shù)似然函數(shù)可以表示為

因此，Θ的極大似然估計(jì)可通過對式（7）取最大化得到。

2.3 剩余壽命估計(jì)

確定退化過程的描述形式后，通過推導(dǎo)首達(dá)時(shí)間概率密度函數(shù)，就可以得到剩余壽命的估計(jì)值；所以，首達(dá)時(shí)間概率密度函數(shù)在剩余壽命估計(jì)的過程中，起到連接退化過程和壽命分布的橋梁作用，推導(dǎo)首達(dá)時(shí)間概率密度函數(shù)是剩余壽命估計(jì)的核心問題。

但是在式（ 2）的情況下，當(dāng) μ（ t；θ）不是常數(shù)時(shí)，要通過推導(dǎo)得到首達(dá)時(shí)間分布的解析式比較困難。如果采用數(shù)值仿真的方法，計(jì)算量很大，所得結(jié)果也不能滿足PHM的決策優(yōu)化要求。為了能夠得到首達(dá)時(shí)間概率密度函數(shù)的解析解，Si等在文獻(xiàn)[11]中通過以下合理假設(shè)，得到了首達(dá)時(shí)間的漸進(jìn)解。

假設(shè)1：在式（1）下如果設(shè)備在t時(shí)刻沒有失效，那么認(rèn)為設(shè)備在t時(shí)刻之前沒有發(fā)生過失效。

假設(shè) 2：潛在的退化過程{X（ t）；t≥0}在 t時(shí)刻達(dá)到失效閾值，即 X（ t）=l，那么{X（ t）；t≥0}在 t之前穿越失效閾值l的概率可以忽略。

假設(shè)1是不考慮設(shè)備維修的影響，因?yàn)槿绻O(shè)備在運(yùn)行過程中進(jìn)行了維修，設(shè)備的退化數(shù)據(jù)就出現(xiàn)了更新，使退化數(shù)據(jù)從另一個(gè)初始狀態(tài)開始退化，本文考慮的過程是設(shè)備在一個(gè)退化周期內(nèi)發(fā)生的退化。假設(shè)2是為了能夠推導(dǎo)出首達(dá)時(shí)間概率密度函數(shù)的解析解。在設(shè)備的運(yùn)行過程中，一般情況下，當(dāng)退化達(dá)到失效閾值時(shí)，認(rèn)為設(shè)備繼續(xù)運(yùn)行下去存在風(fēng)險(xiǎn)，因此會(huì)對設(shè)備進(jìn)行維修檢查或者停止使用。然而，在工程實(shí)際中，也可能存在設(shè)備的退化量在某一時(shí)刻達(dá)到了失效閾值，又由于某種原因在短時(shí)間內(nèi)回到正常范圍的情況；但是，這種事件只是一種小概率的事件，假設(shè)2的目的是在計(jì)算過程中將這種小概率事件忽略，這在現(xiàn)實(shí)中是合理可行的。

基于以上兩種假設(shè)，可以推導(dǎo)出剩余壽命概率密度函數(shù)[12]為

具體的，對于式（3）描述的退化過程，由式（2）可得 μ（ t；θ）=A+Bcos（ ωt）。因此在滿足假設(shè) 1，2 的情況下，可由式（ 3）推導(dǎo)出{X（ t），t≥0}穿越首達(dá)失效閾值l的首達(dá)時(shí)間概率密度函數(shù)為

已知設(shè)備的退化過程，假設(shè)設(shè)備運(yùn)行到τ時(shí)刻仍未失效，且在τ時(shí)刻的退化量為xτ（xτ＜l），則設(shè)備的剩余壽命Tl可以定義為

基于此，可以推導(dǎo)出剩余壽命密度函數(shù)的解析式

2.4 剩余壽命估計(jì)算法步驟

基于以上分析，現(xiàn)將周期應(yīng)力作用下的非線性退化過程剩余壽命估計(jì)步驟總結(jié)如下：

1）根據(jù)退化數(shù)據(jù)建立式（3）的退化模型；

2）利用退化數(shù)據(jù)根據(jù)式（7）進(jìn)行模型參數(shù)的極大似然估計(jì)；

3）確定參數(shù)值后利用式（11）進(jìn)行剩余壽命預(yù)測，得出結(jié)論。

3 仿真研究

通過數(shù)值仿真主要解決兩個(gè)問題：

1）首達(dá)時(shí)間概率密度函數(shù)的推導(dǎo)是進(jìn)行壽命預(yù)測的一個(gè)難點(diǎn)，為了得到式（9），在推導(dǎo)的過程中使用了近似的算法，因此有必要通過數(shù)值仿真驗(yàn)證式（9）的正確性和近似準(zhǔn)確度。

2）對于數(shù)據(jù)中失效時(shí)間的確定問題，考慮在實(shí)際的情況中，退化數(shù)據(jù)超過閾值的真實(shí)時(shí)間往往無法獲得；因?yàn)樵趯?shí)際的數(shù)據(jù)的采樣中，退化超過閾值的時(shí)刻不可能剛好在采樣點(diǎn)上，而是在兩個(gè)采樣點(diǎn)之間的某個(gè)時(shí)刻?；诖?，本文將數(shù)據(jù)首次超過失效閾值時(shí)所對應(yīng)的測試時(shí)間近似為失效時(shí)間。

采用Euler離散化[13]的方法來近似退化過程{X（ t），t＞0}：

其中 Y～N（ 0，1），Δt為離散化步長，θ是未知參數(shù)，為了比較首達(dá)時(shí)間數(shù)值解和本文的解析解，在仿真過程中，假設(shè)模型的參數(shù)已知。

3.1 模型驗(yàn)證

通過數(shù)值仿真的方法來近似退化過程{X（ t），t≥0}，通過統(tǒng)計(jì)仿真數(shù)據(jù)中的首達(dá)時(shí)間數(shù)據(jù)，繪制首達(dá)時(shí)間分布直方圖，然后通過近似計(jì)算得到的解析式（9）繪制的首達(dá)時(shí)間分布圖進(jìn)行類比。

為了驗(yàn)證模型，首先設(shè)定退化過程為

其中初始 X（ 0）=0，設(shè)定閾值 l=100，樣本軌跡數(shù)量M=10000，離散化步長Δt=0.01。為了檢驗(yàn)?zāi)Ｐ偷倪m應(yīng)性，對退化過程中隨機(jī)性作用比較大的參數(shù)σB分別取3個(gè)不同的值：1）令σB取值較小，也就是考慮擴(kuò)散系數(shù)對模型的主導(dǎo)性較小的情況，令σB=2；2）擴(kuò)散系數(shù)適中的情況，令σB=6；3）擴(kuò)散系數(shù)較大的情況，令σB=10。分別繪制這3種情況的首達(dá)時(shí)間分布直方圖和解析結(jié)果曲線，如圖2所示。

可以看出，數(shù)值仿真的首達(dá)時(shí)間統(tǒng)計(jì)結(jié)果和解析漸進(jìn)解的結(jié)果基本符合。當(dāng)仿真的樣本數(shù)量不斷增大，達(dá)到一定程度時(shí)，仿真數(shù)據(jù)的首達(dá)時(shí)間統(tǒng)計(jì)結(jié)果直方圖就會(huì)無限趨近于真實(shí)的首達(dá)時(shí)間PDF，基于此，可以驗(yàn)證本文的方法在一定的前提下可以進(jìn)行首達(dá)時(shí)間的估計(jì)。值得一提的是，當(dāng)σB=2時(shí)，也就是圖2（a）的情況下，可以看出，首達(dá)時(shí)間的PDF近似于逆高斯分布，這也驗(yàn)證了當(dāng)擴(kuò)散系數(shù)對隨機(jī)過程的影響比較小時(shí)，該模型可以近似于線性模型。

3.2 剩余壽命估計(jì)

鑒于Wiener過程在剩余壽命估計(jì)中使用較為廣泛，將Wiener過程（記為模型M1）同本文提出的模型（記為模型M2）進(jìn)行比較，這里將失效閾值設(shè)為1000。分別采用模型M1和M2對6組仿真數(shù)據(jù)的退化過程進(jìn)行建模分析。

根據(jù)失效時(shí)間的確定方法，平均失效時(shí)間（mean time to failure，MTTF）為 226.3h。具體數(shù)據(jù)如圖 3 所示，可以看出，退化數(shù)據(jù)具有一定的非線性特征。利用上述數(shù)據(jù)，進(jìn)行參數(shù)估計(jì)，通過比較退化模型的平均失效時(shí)間和總體均方誤差（total mean square error，TMSE）的值來比較模型的擬合效果。具體的結(jié)果如表1所示。

表1 Wiener過程模型與本文退化模型的比較

圖2 數(shù)值仿真首達(dá)時(shí)間統(tǒng)計(jì)和解析漸進(jìn)解的比較結(jié)果

可以看出，本文提出的方法所得到的MTTF值更為接近真實(shí)結(jié)果，并且本文方法得到的TMSE值也好于線性模型。

圖4為不同退化檢測點(diǎn)上估計(jì)的剩余壽命概率密度函數(shù)曲線、實(shí)際的剩余壽命曲線和估計(jì)的平均剩余壽命曲線?？梢钥闯?，本文采用的方法同Wiener過程建模方法得到的估計(jì)結(jié)果顯然有很大的差異，本文方法得到的剩余壽命估計(jì)值更接近實(shí)際的剩余壽命值。這說明對于受到周期應(yīng)力作用的非線性退化過程，本文方法所得到的結(jié)果更準(zhǔn)確。

圖3 仿真退化數(shù)據(jù)

4 結(jié)束語

長壽命、高成本、高可靠性設(shè)備在貯存或使用過程中可能受到周期性的應(yīng)力作用，為了描述這種非線性退化過程，本文提出了一種基于擴(kuò)散過程的退化模型，并在首達(dá)時(shí)間的意義下通過近似計(jì)算推導(dǎo)了相關(guān)的剩余壽命分布。通過數(shù)值仿真和統(tǒng)計(jì)分析驗(yàn)證壽命分布的正確性。可以看出，當(dāng)退化數(shù)據(jù)由于周期應(yīng)力作用而出現(xiàn)非線性波動(dòng)時(shí)，本文提出的模型得到的實(shí)驗(yàn)結(jié)果優(yōu)于現(xiàn)有方法，能夠?yàn)樵O(shè)備的監(jiān)測維護(hù)與管理提供更加準(zhǔn)確的決策信息，具有潛在的工程應(yīng)用價(jià)值。

圖4 模型M1、M2剩余壽命估計(jì)結(jié)果比較

本文方法對于周期應(yīng)力作用的描述更為準(zhǔn)確，相比之下得到的估計(jì)結(jié)果也更加準(zhǔn)確，但同時(shí)由于壽命分布函數(shù)相對復(fù)雜，在參數(shù)估計(jì)時(shí)難以得到解析的結(jié)果，使得參數(shù)估計(jì)相比Wiener過程更加復(fù)雜，這是下一步亟需解決的問題。

[1]曾聲奎，Petch M，吳際.故障預(yù)測與健康管理（PHM）技術(shù)的現(xiàn)狀與發(fā)展[J].航空學(xué)報(bào)，2005，26（ 5）：626-632.

[2]Jardine A K S， Lin D， Banjevic D.A review on machinery diagnostics and prognostics implementing conditionbased maintenance[J].Mechanical Systems and Signal Processing，2006，20（ 7）：1483-1510.

[3]Si X S， Wang W B， Hu C H， et al.Remaining useful life estimation：A review on the statistical data driven approaches[J].European Journal of Operational Research，2011（ 213）：1-14.

[4]Doksum K A，Hoyland A.Models for variable-stress accelerated life testing experimentsbased on wiener processes and the inverse gaussian distribution[J].Theory of Probability and its Applications，1993，37（ 1）：137-139.

[5]Kharoufeh J P.Explicit results for wear processes in A markovian environment[J].Operations Research Letters，2003，31（ 3）：237-244.

[6]Abdel-Hameed M.A gamma wear process[J].IEEE Transactions on Reliability，1975，24（ 2）：152-153.

[7]Wang X，Xu D.An inverse gaussian process model for degradation data[J].Technometrics，2010，52（ 2）：188-197.

[8]Whitmore G A F.Schenkelberg，modelling accelerated degradation data using wienerdiffusion with a time scale transformation[J].Lifetime Data Analysis，1997（ 3）：27-45.

[9]Lu C J，Meeker W Q.Using degradation measures to estimate a time-to-failure distribution[J].Technom-etrics，1993，35（ 2）：161-174.

[10]Peng C Y，Tseng S T.Mis-specification analysis of linear degradation models[J].IEEE Transactions on Relia bility，2009，58（ 3）：444-455.

[11]Si X S， Wang W B， Hu C H， et al.Remaining useful life estimation based on a nonlinear diffusion degradation process[J].IEEE Transactionson Reliability，2012，61（ 1）：50-67.

[12]Kloeden P，Platen E.Numerical solution of stochastic differential equations[M].New York：Springer，1995：71.

[13]Beskos A，Papaspiliopoulos O，Roberts G O，et al.Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes[J].Journal of the Royal Statistical Association：Series B，2006，68（ 3）：333-382.