馮 磊， 王宏力， 韓曉霞， 周志杰

0 引 言

陀螺儀是慣性導(dǎo)航系統(tǒng)的敏感部件，其性能的優(yōu)劣直接決定整個(gè)系統(tǒng)的導(dǎo)航精度［1］。隨著使用時(shí)間的增長(zhǎng)，陀螺儀性能不可避免地發(fā)生退化。在這種情況下，及時(shí)有效地預(yù)測(cè)陀螺儀的剩余壽命，適時(shí)采取必要的維護(hù)手段，對(duì)于提高整個(gè)慣性導(dǎo)航系統(tǒng)的可靠性與安全性有著至關(guān)重要的作用［2］。

現(xiàn)有的預(yù)測(cè)方法中，一般假設(shè)設(shè)備的性能退化變量直接可觀測(cè)，并且可以完整地反映設(shè)備的退化狀態(tài)。通過(guò)建模設(shè)備的退化過(guò)程，預(yù)測(cè)退化狀態(tài)達(dá)到給定失效閾值的時(shí)間，進(jìn)而預(yù)測(cè)設(shè)備的剩余壽命［3］。但是，工程實(shí)際中，直接測(cè)量的退化變量因受到各種噪聲與誤差的影響，不能精確反映設(shè)備的退化狀態(tài)，這種測(cè)量被稱為不完整測(cè)量［4］。例如，在測(cè)試慣性導(dǎo)航系統(tǒng)的陀螺儀漂移系數(shù)過(guò)程中，由于受到各種因素的影響，測(cè)量結(jié)果不可避免地存在噪聲干擾與誤差。為解決直接監(jiān)測(cè)量帶噪聲的問(wèn)題，可以采用狀態(tài)空間模型，建模設(shè)備實(shí)際的退化過(guò)程以及退化過(guò)程與不完整測(cè)量之間的關(guān)系，通過(guò)濾波技術(shù)預(yù)測(cè)設(shè)備的剩余壽命。此外，根據(jù)濾波算法的預(yù)測(cè)方程與更新方程可以方便地實(shí)現(xiàn)在線預(yù)測(cè)。因此，基于狀態(tài)空間模型的剩余壽命預(yù)測(cè)方法日益受到研究人員的重視。但是，現(xiàn)有研究中大都直接利用退化狀態(tài)的估計(jì)均值作為設(shè)備實(shí)際的退化狀態(tài)值，不能精確反映狀態(tài)估計(jì)的不確定性對(duì)后續(xù)剩余壽命預(yù)測(cè)的影響［5-7］。

基于上述分析，文中考慮直接監(jiān)測(cè)量含有噪聲的剩余壽命預(yù)測(cè)問(wèn)題，利用基于線性漂移的布朗運(yùn)動(dòng)建模陀螺儀退化過(guò)程，然后構(gòu)建狀態(tài)空間模型表征不完整測(cè)量與實(shí)際退化狀態(tài)的關(guān)系。在利用期望最大化算法與卡爾曼濾波聯(lián)合估計(jì)與更新未知參數(shù)與退化狀態(tài)后，充分考慮狀態(tài)估計(jì)的不確定性，將估計(jì)的退化狀態(tài)分布引入剩余壽命預(yù)測(cè)過(guò)程中，得到剩余壽命概率密度函數(shù)的參數(shù)化解析形式。一旦新的監(jiān)測(cè)值可用，即可更新未知參數(shù)與退化狀態(tài)。相應(yīng)地，可以更新剩余壽命的分布，最終實(shí)現(xiàn)在線剩余壽命預(yù)測(cè)。實(shí)驗(yàn)結(jié)果表明了所提方法具有較小的預(yù)測(cè)不確定性。

1 基于狀態(tài)空間技術(shù)的剩余壽命預(yù)測(cè)

1.1 退化過(guò)程建模

假設(shè)設(shè)備實(shí)際的退化過(guò)程｛X（t），t≥0｝可由基于線性漂移的布朗運(yùn)動(dòng)建模，即X（t）可表示為

基于線性漂移的布朗運(yùn)動(dòng)被廣泛用于退化過(guò)程建模，其最大優(yōu)點(diǎn)是可以根據(jù)首達(dá)時(shí)間（First Hitting Time，F(xiàn)HT）的概念定義設(shè)備的壽命，進(jìn)而得到壽命分布的解析形式，即逆高斯分布［8］。因此，時(shí)刻ti設(shè)備的剩余壽命可以定義為退化過(guò)程｛X（t），t≥ti｝達(dá)到預(yù)先設(shè)定的失效閾值的首達(dá)時(shí)間［9］，即剩余壽命Tr可定義為

若直接監(jiān)測(cè)的退化變量不含有噪聲，則可根據(jù)式（2）直接得到剩余壽命分布的概率密度函數(shù)。但是，正如上節(jié)所述，不完整測(cè)量廣泛存在于工程實(shí)際中，直接監(jiān)測(cè)信息不能精確反映設(shè)備的退化狀態(tài)，因此，不能直接利用式（2）計(jì)算設(shè)備剩余壽命。文中利用狀態(tài)空間技術(shù)解決此問(wèn)題。

1.2 狀態(tài)空間模型的構(gòu)建

為建模直接監(jiān)測(cè)與實(shí)際退化之間的關(guān)系，采用如下加性噪聲模型構(gòu)建觀測(cè)方程:

其中，ε（t）為噪聲項(xiàng)，假設(shè)ε（t）～N（0，σ2）且與B（t）相互獨(dú)立。

為了辨識(shí)實(shí)際的退化狀態(tài)，首先將動(dòng)態(tài)系統(tǒng)（1）和（3）離散化，得到離散時(shí)間點(diǎn)tk＝kΔt，k＝1，2，…的狀態(tài)方程和觀測(cè)方程:

根據(jù)狀態(tài)空間模型（4），應(yīng)用卡爾曼濾波［10］估計(jì)當(dāng)前時(shí)刻實(shí)際的退化狀態(tài)。首先定義E（Xi｜Y1∶i）和＝Var（Xi｜Y1∶i）為利用當(dāng)前時(shí)刻的整個(gè)觀測(cè)序列Y1∶i＝Δ｛Y1，Y2，…，Yi｝得到的狀態(tài)Xi的條件期望與方差，即Xi～N（X∧i｜i，Pi｜i）?，F(xiàn)有研究大都直接將狀態(tài)估計(jì)的均值作為其實(shí)際值，直接可以得到時(shí)刻ti，剩余壽命概率密度函數(shù)為:

式（5）沒(méi)有考慮狀態(tài)估計(jì)的不確定性，為了減少剩余壽命預(yù)測(cè)的不確定性，文中將估計(jì)的狀態(tài)分布引入后續(xù)的剩余壽命預(yù)測(cè)中。根據(jù)全概率定律，可以得到當(dāng)前時(shí)刻ti，剩余壽命概率密度函數(shù)為:

2 參數(shù)估計(jì)算法

由于模型中含有未知參數(shù)，并且狀態(tài)是未知的，因此，采用EM算法解決存在未知狀態(tài)時(shí)的參數(shù)極大似然估計(jì)問(wèn)題。為后續(xù)表示方便，首先定義當(dāng)前時(shí)刻的狀態(tài)序列為Xi｝。通過(guò)迭代計(jì)算且最大化包含完整數(shù)據(jù)集（X1∶i，Y1∶i）的對(duì)數(shù)似然函數(shù)的條件期望，EM 算法能夠產(chǎn)生一個(gè)估計(jì)序列收斂到參數(shù)的極大似然估計(jì)值［11］。令θ∧（j）表示 EM 算法第j次迭代所得的估計(jì)值，則完整數(shù)據(jù)對(duì)數(shù)似然函數(shù)的條件期望可表示為:

其中，L（X1∶i，Y1∶i｜θ）＝logp（X1∶i，Y1∶i｜θ）為完整數(shù)據(jù)對(duì)數(shù)似然函數(shù)?？偨Y(jié)起來(lái)，參數(shù)估計(jì)程序包含如下兩個(gè)步驟:

1）E 步。

計(jì)算 Q（θ，θ∧（j））。

2）M 步。

求解:

具體到文中構(gòu)建的狀態(tài)空間模型（4），為了計(jì)算對(duì)數(shù)似然函數(shù)的條件期望，對(duì)于k＝1，2，…，i，首先定義如下變量:

經(jīng)過(guò)一系列代數(shù)運(yùn)算，可以得到

其中

為降低參數(shù)空間的維數(shù)，減少計(jì)算復(fù)雜度，分別計(jì)算式（8）相對(duì)于參數(shù)λ，σB和σ的偏導(dǎo)數(shù)，并且令這3個(gè)偏導(dǎo)數(shù)為零，可得

根據(jù)一階必要條件可知，最優(yōu)的參數(shù)估計(jì)值必滿足式（9）～式（11）。因此，可通過(guò)直接計(jì)算得到未知參數(shù)的估計(jì)值。

由上述推導(dǎo)可知，一旦獲取新的觀測(cè)值，就可以利用卡爾曼濾波和EM算法，估計(jì)和更新退化狀態(tài)及未知參數(shù)。相應(yīng)地，根據(jù)式（6）可以更新剩余壽命的概率分布，實(shí)現(xiàn)實(shí)時(shí)剩余壽命預(yù)測(cè)。

3 陀螺儀剩余壽命預(yù)測(cè)實(shí)驗(yàn)

3.1 問(wèn)題描述

陀螺儀漂移是表征其性能的一項(xiàng)重要指標(biāo)［12］。運(yùn)行狀態(tài)下，陀螺儀的轉(zhuǎn)子高速旋轉(zhuǎn)必然造成轉(zhuǎn)軸的磨損，引起漂移系數(shù)的增長(zhǎng)，當(dāng)漂移系數(shù)增大到預(yù)先設(shè)定的閾值，即認(rèn)為系統(tǒng)發(fā)生失效。本實(shí)驗(yàn)選取敏感軸方向的一次項(xiàng)漂移系數(shù)作為衡量陀螺儀退化狀態(tài)的性能指標(biāo)。該系數(shù)的值一旦達(dá)到失效閾值0.38°／h，即表明陀螺儀已經(jīng)失效。對(duì)陀螺儀連續(xù)工作測(cè)試，采樣間隔為3h，所得到的測(cè)試數(shù)據(jù)共96組，如圖1所示。

圖1 漂移系數(shù)測(cè)試數(shù)據(jù)

3.2 實(shí)驗(yàn)結(jié)果

為預(yù)測(cè)陀螺儀的剩余壽命，首先利用EM算法對(duì)參數(shù)進(jìn)行估計(jì)與更新。各未知參數(shù)的估計(jì)與更新結(jié)果如圖2所示。

由圖2可知，隨著監(jiān)測(cè)數(shù)據(jù)的增加，各未知參數(shù)很快收斂到各自的穩(wěn)定值，說(shuō)明所提參數(shù)估計(jì)方法的有效性。其中，參數(shù)a最終的估計(jì)結(jié)果為0.004 6，σB的估計(jì)結(jié)果為0.007 8，σ的估計(jì)結(jié)果為0.004 5。

圖2 參數(shù)估計(jì)與更新結(jié)果

為了驗(yàn)證文中方法的有效性，將文中模型與不考慮狀態(tài)估計(jì)不確定性的模型做比較研究。為表示方便，將文中所構(gòu)建的模型稱為模型1，不考慮狀態(tài)估計(jì)不確定性的模型為模型2。以第90個(gè)到96個(gè)數(shù)據(jù)為例，剩余壽命預(yù)測(cè)結(jié)果如圖3所示。

圖3 兩種模型的剩余壽命預(yù)測(cè)結(jié)果的比較

由圖3可知，兩種模型的預(yù)測(cè)曲線都比較平滑，但是文中模型考慮了退化狀態(tài)估計(jì)的不確定性，因此，利用文中方法預(yù)測(cè)得到的剩余壽命的概率密度曲線更加緊致，預(yù)測(cè)結(jié)果的不確定性小。由于后續(xù)的維護(hù)策略優(yōu)化需要預(yù)測(cè)的剩余壽命信息，因此，文中方法可以減少維護(hù)策略建模的不確定性，降低失效風(fēng)險(xiǎn)，提高維護(hù)效率。

4 結(jié) 語(yǔ)

陀螺儀剩余壽命預(yù)測(cè)實(shí)踐中，受各種因素的影響，直接監(jiān)測(cè)的退化變量含有噪聲，不能準(zhǔn)確反映陀螺儀的退化狀態(tài)。為解決此問(wèn)題，構(gòu)建了狀態(tài)空間模型建模實(shí)際退化過(guò)程以及直接監(jiān)測(cè)與退化狀態(tài)間的關(guān)系，并且考慮狀態(tài)估計(jì)的不確定性，將其引入剩余壽命預(yù)測(cè)過(guò)程，最終實(shí)現(xiàn)剩余壽命的在線預(yù)測(cè)。仿真實(shí)驗(yàn)表明，文中所提方法能夠較好地解決陀螺儀剩余壽命預(yù)測(cè)的相關(guān)問(wèn)題，并且可以減少預(yù)測(cè)的不確定性，為后續(xù)的維護(hù)策略安排提供準(zhǔn)確而實(shí)時(shí)的預(yù)測(cè)信息。

［1］ 胡昌華，馬清亮，鄭建飛.導(dǎo)彈測(cè)試與發(fā)射控制技術(shù)［M］.北京:國(guó)防工業(yè)出版社，2010.

［2］ Zhou Z J，Hu C H，Xu D L，et al.A model for re-al-time failure prognosis based on hidden Markov model and belief rule base［J］.European Journal of Operational Research，2010，207（1）:269-283.

［3］ Si X S，Wang W，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）:1-14.

［4］ Sun J Z，Zuo H F，Pecht M G.Advances in sequential monte carlo methods for joint state and parameter estimation applied to prognostics ［C］／／Prognostics and System Health Management Conference.Shenzhen，China:［s.n.］，2011:1-7.

［5］ Christer A H，Wang W，Sharp J M.A state space condition monitoring model for furnace erosion prediction and replacement ［J］.European Journal of Operational Research，1997，101（1）:1-14.

［6］ Carr M J，Wang W.Modeling failure modes for residual life prediction using stochastic filtering theory［J］.IEEE Transactions on Reliability，2010，59（2）:90-96.

［7］ Wang W，Carr M J，Xu W，et al.A model for residual life prediction based on Brownian motion with an adaptive drift［J］.Microelectronics Reliability，2011，51（2）:285-293.

［8］ Chhikara R S，F(xiàn)olks J L.The inverse gaussian distribution as a lifetime model［J］.Technometrics，1977，19（4）:461-468.

［9］ Lee M L T，Whitmore G A.Threshold regression for survival analysis:Modeling event times by a stochastic process reaching a boundary［J］.Statistical Science，2006，21（4）:501-513.

［10］ 王志賢.最優(yōu)狀態(tài)估計(jì)與系統(tǒng)辨識(shí)［M］.西安:西北工業(yè)大學(xué)出版社，2004.

［11］ Schon T B，Wills A，Ninness B.System identification of nonlinear state-space models ［J］.Automatica，2011，47（1）:39-49.

［12］ 姜楠，張雙才，蘇莉蔚.基于 MEMS的慣性測(cè)量系統(tǒng)［J］.長(zhǎng)春工業(yè)大學(xué)學(xué)報(bào):自然科學(xué)版，2010，31（2）:181-184.