周期應(yīng)力下非線性退化建模與剩余壽命估計(jì)

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

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

0 引 言

近年來(lái)，預(yù)測(cè)與健康管理（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ù)測(cè)是PHM的重要組成部分，其中剩余壽命（ remaining useful lifetime，RUL）的預(yù)測(cè)為優(yōu)化設(shè)備檢測(cè)與維護(hù)策略，合理訂購(gòu)備件提供了科學(xué)的依據(jù)，因此也是當(dāng)前可靠性領(lǐng)域研究的重點(diǎn)和熱點(diǎn)。

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

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

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

1 問(wèn)題描述

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

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

在建立退化模型后，要進(jìn)行剩余壽命估計(jì)，首先需要對(duì)設(shè)備的壽命進(jìn)行定義。將設(shè)備失效的時(shí)間定義為隨機(jī)退化過(guò)程{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è)備的退化過(guò)程描述為一類(lèi)擴(kuò)散過(guò)程{X（ t）；t≥0}。

式中t是模型的時(shí)間尺度，當(dāng)退化過(guò)程為線性退化過(guò)程時(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)定性、測(cè)量設(shè)備測(cè)量誤差及穩(wěn)定性、測(cè)試過(guò)程中的外部干擾等隨機(jī)因素決定；μ（ t；θ）是時(shí)間 t的非線性函數(shù)，θ為未知參數(shù)。式（2）將退化過(guò)程就描述為非時(shí)齊的擴(kuò)散過(guò)程。在設(shè)備的退化過(guò)程中，受到的環(huán)境或應(yīng)力影響可能是不斷變化的，μ（t；θ）是隨時(shí)間的非線性函數(shù)，因此可以描述設(shè)備所處環(huán)境影響隨時(shí)間不斷變化的過(guò)程。 并且，當(dāng) μ（ t；θ）=μ 時(shí)，該模型就是一種線性漂移的擴(kuò)散過(guò)程，這也說(shuō)明了式（2）更加具有一般性和靈活性。

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

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

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

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

2.1 退化建模

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

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

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

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

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

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

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

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

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

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

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

2.3 剩余壽命估計(jì)

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

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

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

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

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

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

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

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

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

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

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

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

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

3 仿真研究

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

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

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

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

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

3.1 模型驗(yàn)證

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

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

其中初始 X（ 0）=0，設(shè)定閾值 l=100，樣本軌跡數(shù)量M=10000，離散化步長(zhǎng)Δt=0.01。為了檢驗(yàn)?zāi)Ｐ偷倪m應(yīng)性，對(duì)退化過(guò)程中隨機(jī)性作用比較大的參數(shù)σB分別取3個(gè)不同的值：1）令σB取值較小，也就是考慮擴(kuò)散系數(shù)對(duì)模型的主導(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ì)無(wú)限趨近于真實(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ù)對(duì)隨機(jī)過(guò)程的影響比較小時(shí)，該模型可以近似于線性模型。

3.2 剩余壽命估計(jì)

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

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

表1 Wiener過(guò)程模型與本文退化模型的比較

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

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

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

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

4 結(jié)束語(yǔ)

長(zhǎng)壽命、高成本、高可靠性設(shè)備在貯存或使用過(guò)程中可能受到周期性的應(yīng)力作用，為了描述這種非線性退化過(guò)程，本文提出了一種基于擴(kuò)散過(guò)程的退化模型，并在首達(dá)時(shí)間的意義下通過(guò)近似計(jì)算推導(dǎo)了相關(guān)的剩余壽命分布。通過(guò)數(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)測(cè)維護(hù)與管理提供更加準(zhǔn)確的決策信息，具有潛在的工程應(yīng)用價(jià)值。

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

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

[1]曾聲奎，Petch M，吳際.故障預(yù)測(cè)與健康管理（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.