基于粒子濾波的故障預(yù)測

基于粒子濾波的故障預(yù)測

航空航天飛機(jī)和軍用推進(jìn)器是一項(xiàng)復(fù)雜的技術(shù)和資料，通常由各種功能子系統(tǒng)和設(shè)備組成，傳輸關(guān)系復(fù)雜。在其功能/性能不斷增強(qiáng)的同時(shí),也帶來了可靠性、安全性和可用性等一系列問題。綜合運(yùn)載器健康管理技術(shù)(IVHM)就是在這種背景下提出的一種全新的管理飛行器健康狀態(tài)的解決方案。IVHM技術(shù)的核心和難點(diǎn)是故障預(yù)測,這也是它區(qū)別于傳統(tǒng)的健康監(jiān)測技術(shù)的主要特點(diǎn)。故障預(yù)測使得設(shè)備維護(hù)人員可以提前預(yù)知故障的發(fā)生,從而采取一系列維修或預(yù)防的措施,而不必等到故障真正發(fā)生之后再做出反應(yīng),從而能有效避免一些災(zāi)難性的后果。對故障預(yù)測問題的研究,一種基本的解決方案是采用基于對象故障演化模型狀態(tài)估計(jì)的故障預(yù)測技術(shù)。它的基本思想是:對象系統(tǒng)的故障演化過程可以由一個(gè)狀態(tài)空間模型加以表述,即存在著一些包含故障演化信息的狀態(tài)變量(通常稱為系統(tǒng)故障指征),這些狀態(tài)變量與系統(tǒng)的故障演化過程之間有著某種對應(yīng)關(guān)系,通過對系統(tǒng)狀態(tài)變量在未來一段時(shí)間內(nèi)變化情況的預(yù)測,結(jié)合一定的判別準(zhǔn)則,可以預(yù)知對象系統(tǒng)的故障情況和剩余使用壽命。對于上述問題,貝葉斯估計(jì)方法僅僅從理論上給出了故障預(yù)測方案,但對于大多數(shù)實(shí)際系統(tǒng)卻很難求得貝葉斯估計(jì)的解析解。因此,實(shí)際應(yīng)用過程中出現(xiàn)了以卡爾曼濾波、擴(kuò)展卡爾曼濾波和無跡濾波為代表的具體故障預(yù)測算法。相對于前述方法,作為狀態(tài)估計(jì)領(lǐng)域一種新興發(fā)展方向的粒子濾波算法具有不受模型線性、高斯假設(shè)約束的特點(diǎn),理論上能適用于任意非線性非高斯動(dòng)態(tài)系統(tǒng)的狀態(tài)估計(jì)問題。因此,本文在粒子濾波技術(shù)的基礎(chǔ)上提出了一種故障預(yù)測算法。此算法在狀態(tài)估計(jì)階段,采用聯(lián)合估計(jì)和粒子濾波同時(shí)估計(jì)系統(tǒng)狀態(tài)和未知參數(shù)的后驗(yàn)分布;在算法的狀態(tài)預(yù)測階段,采用了兩種不同的計(jì)算方法。這兩種方法的計(jì)算結(jié)果可以在一個(gè)統(tǒng)一的預(yù)測框架之下進(jìn)行有效交互,從而進(jìn)一步提高了預(yù)測的準(zhǔn)確性和可靠性。1對象系統(tǒng)故障模型假設(shè)對象系統(tǒng)的故障演化過程可以采用如下的狀態(tài)空間模型描述:Xk=fk-1(Xk-1,θk,Uk-1)+wk-1(1)Yk=hk(Xk,θk)+vk(2)式中:Xk和Yk分別為系統(tǒng)的狀態(tài)變量和觀測變量(Xk又稱故障指征,直接與對象故障演化過程相關(guān));Uk為系統(tǒng)的輸入變量,通常用來描述系統(tǒng)工作的外部環(huán)境;θk為模型的未知參數(shù)(通常隨時(shí)間緩慢變化),用于描述由于對象系統(tǒng)故障演化而導(dǎo)致的部分模型參數(shù)隨時(shí)間變化的特性;wk和vk分別為系統(tǒng)噪聲和觀測噪聲。對于上述的對象系統(tǒng),故障預(yù)測問題可以轉(zhuǎn)化為一個(gè)已知系統(tǒng)從某一起始時(shí)刻0至當(dāng)前時(shí)刻k的觀測序列{Y1,Y2,…,Yk},同時(shí)在系統(tǒng)存在未知緩變參數(shù)的前提條件下,對未來某一時(shí)刻k+p(p>0)的系統(tǒng)狀態(tài)變量Xk+p求取其在一定的準(zhǔn)則意義下的估計(jì)值的狀態(tài)估計(jì)問題。2實(shí)驗(yàn)2不同實(shí)驗(yàn)條件下高斯混合模型的離散化粒子濾波算法的核心是使用粒子及其權(quán)值組成的離散隨機(jī)測度近似于變量的分布密度。但是,離散形式的分布密度不便于濾波輸出的交互。由于任意連續(xù)隨機(jī)變量的分布密度均可采用一組高斯混合模型近似描述,即對變量Xk的一組隨機(jī)測度{Xik,Wik},其樣本均值?Xk和方差Pk為?Xk=Ν∑i=1XikWikΡk=Ν∑i=1Wik(Xik-?Xk)[Xik-?Xk]Τ}(3)則其分布密度可以近似為Ρ(Xk|Y1∶k)≈Ν∑i=1WikΝ(Xik,νΡk)(4)式中:ν=0.5N-2/d,d為狀態(tài)變量Xk的維數(shù),N為采樣樣本的粒子數(shù)。并且,由于多組高斯混合模型之間的加權(quán)線性組合仍為高斯混合模型,同時(shí)用蒙特卡羅方法可以方便地將高斯混合模型采樣離散化。為此采用高斯混合模型近似變量分布密度,并且將其嵌入基本粒子濾波的算法中,形成了如下的高斯混合模型粒子濾波算法:假設(shè)對象系統(tǒng)的狀態(tài)方程和觀測方程分別為P(Xk|Xk-1)和P(Yk|Xk)。(1)粒子權(quán)值的認(rèn)定當(dāng)T=0時(shí),生成服從先驗(yàn)分布P(X0)的N個(gè)樣本粒子{X10,X20,…,XΝ0},每個(gè)粒子的權(quán)值Wi0設(shè)為1/N。按照式(3)和式(4)計(jì)算狀態(tài)變量分布密度的初始高斯混合模型。(2)混合模型的形式采樣k-1時(shí)刻N(yùn)個(gè)獨(dú)立同分布的樣本粒子服從P(Xk-1|Y1∶k-1)(采用形如式(4)的高斯混合模型形式)。計(jì)算k時(shí)刻的樣本粒子Xik～P(Xk|Xik-1)。(3)測量更新更新k時(shí)刻樣本對應(yīng)的權(quán)值:Wik=P(Yk|Xik)Wik-1。并進(jìn)行權(quán)值歸一化:Wik=WikΝ∑j=1Wjk。(4)高斯混合模型按照式(3)和式(4)計(jì)算狀態(tài)變量在T=k時(shí)刻狀態(tài)變量后驗(yàn)分布密度的高斯混合模型。由于在上述的濾波過程中,每個(gè)時(shí)刻的初始樣本粒子均為高斯混合模型重采樣所獲取的獨(dú)立同分布粒子,有效緩解了濾波過程的粒子權(quán)值退化,因此基本粒子濾波算法中的再采樣環(huán)節(jié)可以省略。3基于粒子濾波故障預(yù)測算法的設(shè)計(jì)針對上述的故障預(yù)測問題,采用粒子濾波算法實(shí)現(xiàn)其故障預(yù)測的過程可按以下3個(gè)步驟:(1)基于聯(lián)合估計(jì)的增廣系統(tǒng)狀態(tài)空間模型為估計(jì)k時(shí)刻的狀態(tài),需要計(jì)算對象系統(tǒng)狀態(tài)變量的濾波密度P(Xk|Y1∶k)。然而,由于對象系統(tǒng)模型中存在未知參數(shù),它的取值將對狀態(tài)變量估計(jì)產(chǎn)生影響。對此問題的一種解決方案是采用聯(lián)合估計(jì),即將對象系統(tǒng)狀態(tài)變量和未知參數(shù)做狀態(tài)增廣,在增廣系統(tǒng)可觀測的前提下,對增廣系統(tǒng)做狀態(tài)估計(jì)。此時(shí)增廣系統(tǒng)的狀態(tài)空間模型為→Xk=[Xkθk]=[fk-1(Xk-1,θk,Uk-1)+wk-1θk-1+ζk-1]Yk=hk(→Xk,vk)=hk(Xk,θk)+vk}(5)式中:→Xk為增廣狀態(tài)變量;原對象系統(tǒng)時(shí)變參數(shù)θk服從一個(gè)隨機(jī)游走模型;模型參數(shù)ζk～N(0,Σζ)服從均值為0,方差為Σζ的高斯隨機(jī)分布。(2)標(biāo)準(zhǔn)權(quán)值表示方法的輸出交互在采用粒子濾波方法獲取的k時(shí)刻增廣系統(tǒng)狀態(tài)變量濾波密度P(→Xk|Y1∶k)的近似測度基礎(chǔ)上,可以估算出未來k+p時(shí)刻的原系統(tǒng)狀態(tài)的先驗(yàn)概率密度P(Xk+p|Y1∶k)。一種方案參照文獻(xiàn),在一定的假設(shè)條件下,只需要采用k時(shí)刻對象系統(tǒng)的狀態(tài)方程(狀態(tài)方程采用式(1)所表述的原對象系統(tǒng)方程,系統(tǒng)未知參數(shù)θk保持k時(shí)刻期望值不變)以濾波密度P(→Xk|Y1∶k)的樣本粒子為初始值,迭代采樣出k+p時(shí)刻原系統(tǒng)狀態(tài)變量的樣本粒子Xik+p,同時(shí)保持其對應(yīng)的權(quán)值不變。則k+p時(shí)刻的原系統(tǒng)變量的先驗(yàn)概率密度為Ρ(Xk+p|Y1∶k)≈Ν∑i=1Wikδ(Xk+p-Xik+p)(6)式中:δ為狄拉克函數(shù)。該方法的缺點(diǎn)是假設(shè)故障演化模型緩變參數(shù)在時(shí)刻k+1至k+p之間保持不變,這種要求有時(shí)難以符合實(shí)際情況,表現(xiàn)為其預(yù)測結(jié)果對系統(tǒng)未來變化趨勢的預(yù)測通常存在一定程度的滯后。另一種方案采用數(shù)據(jù)驅(qū)動(dòng)的預(yù)測算法(例如ARMA模型、灰模型或支持向量機(jī)回歸模型等)預(yù)測時(shí)刻k+1至k+p的量測信息,將k+p時(shí)刻原系統(tǒng)狀態(tài)變量先驗(yàn)分布的預(yù)測問題轉(zhuǎn)化為k+p時(shí)刻增廣系統(tǒng)狀態(tài)變量后驗(yàn)分布的估計(jì)問題,然后采用本文的改進(jìn)粒子濾波算法求解。該方法能有效克服上一種方法的缺點(diǎn),但是其預(yù)測結(jié)果容易受到系統(tǒng)量測噪聲的影響,有可能會偏離系統(tǒng)未來真實(shí)的變化情況。為了提高預(yù)測算法的準(zhǔn)確性和可靠性,最直接的做法是將上述兩種計(jì)算方法的結(jié)果進(jìn)行信息融合。一種有效的融合算法可以參考交互式多模型算法中的輸出交互環(huán)節(jié),對上述兩種方法的結(jié)果做加權(quán)線性組合,即Ρ(Xk+p|Y1∶k)≈u1Ρ(X1,k+p|Y1∶k)+u2Ρ(X2,k+p|Y1∶k)(7)式中:P(X1,k+p|Y1∶k)和P(X2,k+p|Y1∶k)分別為以上兩種方法的預(yù)測結(jié)果;u1和u2分別為兩種方法對應(yīng)的權(quán)值。P(X1,k+p|Y1∶k)和P(X2,k+p|Y1∶k)分別采用式(4)所示的高斯混合模型近似模擬。在步驟(3)中各種指標(biāo)計(jì)算所需的狀態(tài)變量的樣本粒子均由以上兩種方法輸出的加權(quán)線性組合P(Xk+p|Y1∶k)的高斯混合模型重采樣產(chǎn)生。權(quán)值u1和u2的確定是通過比較兩種方法預(yù)測結(jié)果和量測真值的相似程度(計(jì)算似然度函數(shù))逐漸對其進(jìn)行更新來實(shí)現(xiàn)。權(quán)值更新算法如下:假設(shè)在時(shí)刻k,對其后時(shí)刻k+p的預(yù)測中u1和u2的初值分別為u10和u20(u10和u20的取值由先驗(yàn)知識決定,若無先驗(yàn)知識則取u10=u20=1/2)。當(dāng)系統(tǒng)獲取k+1時(shí)刻的量測信息Yk+1時(shí),分別計(jì)算在量測為Yk+1,兩種方法的似然度flik,1和flik,2。設(shè)兩種方法k+1時(shí)刻預(yù)測結(jié)果樣本粒子所產(chǎn)生的量測誤差以及其方差分別為rij,k+1和Sij,k+1,i=1,2為預(yù)測方法標(biāo)號;j=1,2,…,N為采樣粒子標(biāo)號。rij,k+1=Yk+1-hk+1(→Xij,k+1,vj,k+1)Sij,k+1=Ν∑j=1Wij,k+1{(hk+1(→Xij,k+1,vj,k+1)-ˉYik+1)?[hk+1(→Xij,k+1,vj,k+1)-ˉYik+1]Τ}}(8)式中:→Xij,k+1為第i種方法k+1時(shí)刻的預(yù)測結(jié)果;Wij,k+1為第j個(gè)粒子對應(yīng)的權(quán)值;ˉYik+1=Ν∑j=1Wij,k+1hk+1(→Xj,k+1,vj,k+1)為第i種方法預(yù)測結(jié)果的量測均值。對于第i種方法,其似然函數(shù)為flik,i=Ν∑j=1Wij,k+1Ν(rij,k+1;0,Sij,k+1)(9)式中:N(X;μ,Σ)為高斯分布密度函數(shù),Ν(X;μ,Σ)=1√(2π)ndet(Σ)exp{-12[X-μ]Τ(Σ)-1(X-μ)}(10)式中:n為隨機(jī)變量X的維數(shù)。在獲取flik,1和flik,2的前提下,兩種預(yù)測方法的致信度u1和u2分別更新為u1=u10flik,1u10flik,1+u20flik,2u2=u20flik,2u10flik,1+u20flik,2}(11)隨著新的量測信息Yk+2,…,Yk+q(其中1<q<p)的依次獲取,可以采用上述方法逐漸更新k時(shí)刻對于未來k+p時(shí)刻的預(yù)測結(jié)果(所不同之處,權(quán)值的初值u10和u20為前次更新的結(jié)果u1和u2)。在高斯混合模型粒子濾波算法的基礎(chǔ)上,將上述兩種預(yù)測方法和權(quán)值更新算法相結(jié)合就構(gòu)成了本文所提出的預(yù)測框架。預(yù)測框架的具體實(shí)現(xiàn)步驟如下:①在時(shí)刻k,采用上述兩種預(yù)測方法計(jì)算其后k+1至k+p時(shí)刻的系統(tǒng)狀態(tài)信息(兩種預(yù)測方法權(quán)值初值u10和u20由先驗(yàn)知識決定,若無先驗(yàn)知識則取u10=u20=1/2)。②當(dāng)獲取k時(shí)刻后的系統(tǒng)量測信息Yk+q(q≥1)時(shí),判斷對k+q時(shí)刻的預(yù)測結(jié)果與系統(tǒng)量測的差值是否超出一定的閾值范圍(或q大于p):若是,調(diào)用①從當(dāng)前時(shí)刻開始向其后預(yù)測p步;否則,按照上述的權(quán)值更新算法修正兩種預(yù)測方法的權(quán)值,從而將前次的預(yù)測結(jié)果更新(更新在時(shí)刻k,對其后k+1至k+p時(shí)刻的預(yù)測結(jié)果)。(3)基于全概率回歸的剩余分布預(yù)測在獲取的k+p時(shí)刻狀態(tài)變量的概率密度的預(yù)測結(jié)果P(Xk+p|Y1:k)的基礎(chǔ)上,結(jié)合一定的故障判據(jù),可以近似推算出對象系統(tǒng)的剩余壽命的分布密度。參照文獻(xiàn),系統(tǒng)剩余壽命分布計(jì)算的原理如圖1所示。假設(shè)對象系統(tǒng)的故障判據(jù)為Σfailure={Ηlow≤X≤Ηup}(12)式中:Hup和Hlow分別為系統(tǒng)故障的上下界。系統(tǒng)變量從時(shí)刻k+1至k+p的預(yù)測結(jié)果的樣本粒子在狀態(tài)空間中形成了一片“粒子云”,并且該粒子云與故障區(qū)域的交會形成一個(gè)交集區(qū)域。將交集中相同時(shí)刻的樣本粒子權(quán)值累加,其結(jié)果反映其對應(yīng)時(shí)刻系統(tǒng)故障可能性大小,即通過上述方法可以獲得系統(tǒng)剩余壽命分布從時(shí)刻k+1至k+p中間的一組樣本粒子,對其相應(yīng)的權(quán)值進(jìn)行歸一化,然后依照全概率公式即可以獲得系統(tǒng)剩余壽命分布的概率密度為Ρrul(k)=∑i=1ΝΡr{Ηlow≤Xki≤Ηup}W?ki(13)式中:W?ki為歸一化權(quán)值。系統(tǒng)的剩余使用壽命從理論上來說是一個(gè)從當(dāng)前時(shí)刻(假設(shè)未發(fā)生故障)至無窮時(shí)刻之間的連續(xù)分布,實(shí)際上無法獲取其解析結(jié)果。因此該方法的本質(zhì)是通過隨機(jī)變量的采樣分布來近似其真實(shí)分布。分析上述算法,由于預(yù)測運(yùn)算存在時(shí)間意義上的截?cái)?即預(yù)測結(jié)果只能計(jì)算至當(dāng)前時(shí)刻k之后的p(有限長)個(gè)時(shí)刻。因此,在特定的時(shí)間段內(nèi),算法將遭遇“密度泄漏”問題,進(jìn)而影響結(jié)果的準(zhǔn)確性,即算法存在一定的適用性條件。例如,在系統(tǒng)工作的初始階段,系統(tǒng)未發(fā)生故障或者故障征兆不明顯時(shí),系統(tǒng)剩余壽命分布的特點(diǎn)是各個(gè)時(shí)刻的概率密度較為平均,直觀表現(xiàn)為分布密度函數(shù)比較平緩。若在該情況下做截?cái)?計(jì)算結(jié)果將存在較大誤差;相對來說,當(dāng)系統(tǒng)接近發(fā)生故障或者故障征兆較為明顯時(shí),系統(tǒng)剩余壽命分布將主要集中在未來一個(gè)有限的時(shí)間區(qū)域內(nèi),直觀表現(xiàn)為分布密度靠近當(dāng)前時(shí)刻且存在一個(gè)或多個(gè)尖峰。在該情況下做截?cái)鄬ψ罱K計(jì)算結(jié)果的影響較小。4罐系統(tǒng)模型建立的仿真采用文獻(xiàn)中的三罐系統(tǒng)作為故障預(yù)測仿真實(shí)驗(yàn)的對象,如圖2所示。該系統(tǒng)由3個(gè)水槽T1,T2和T3組成,3個(gè)水槽之間依次通過水管連接;水槽T2連接一個(gè)排水閥門,由它控制整個(gè)系統(tǒng)的輸出;系統(tǒng)的輸入分別由連接在水槽T1和T2的水泵控制,其輸入流量分別為Q1和Q2;系統(tǒng)的狀態(tài)變量為3個(gè)水槽的液位h1,h2和h3;三罐系統(tǒng)的數(shù)學(xué)模型為dXdt=A(X)+BU(t)(14)式中:{X=[x1x2x3]Τ=[h1h2h3]ΤU=[Q1Q2]ΤA(X)=1A[-Q13Q32-Q20Q13-Q32],B=1A在上述模型中:Ci為水槽排水管流量系數(shù);hi為相應(yīng)的水槽液位高度;Qij為相互連接的水槽間水流速率;Q1和Q2分別為系統(tǒng)的輸入流量;A為水槽的截面積;Sn為水槽間連接水管的截面積;i=1,2,3并且(i,j)∈{(1,3);(3,2);(2,0)};未知系數(shù)Q13,Q32,Q20可以采用托里切利原理求解:q=CSnsgn(Δh)(2g|Δh|)1/2(15)式中:g為地球的重力加速度;Δh為相鄰的兩個(gè)水槽的液位之差;q為水槽連接水管的水流速率。因此,模型中的未知系數(shù)分別為{Q13=C1Snsgn(h1-h3)(2g|h1-h3|)1/2Q32=C3Snsgn(h3-h2)(2g|h3-h2|)1/2Q20=C2Sn(2gh2)1/2將系統(tǒng)模型離散化,對式(14)建立形如式(1)的非線性離散時(shí)間系統(tǒng)模型。Δt為單位采樣時(shí)間,則X(k+1)=X(k)+ΔtA(X)+ΔtBU(k)+w(k)(16)然后設(shè)計(jì)形如式(2)的觀測方程:Y(k)=kp[Q13(k)Q32(k)Q30(k)]Τ+v(k)(17)式中:kp為觀測方程系數(shù);三罐系統(tǒng)模型有關(guān)參數(shù)的正常取值見表1。3個(gè)水槽的液位初始值分別為h10=0.40,h02=0.30,h03=0.35;系統(tǒng)的輸入Q1和Q2為Qk=4.5×10-5m3/s,k=1,2。系統(tǒng)的故障狀態(tài)?定義為?={|hk-hk0hk0|≥0.2}(k=1,2,3)(18)整個(gè)仿真實(shí)驗(yàn)時(shí)間為200Δt,在100Δt時(shí),系統(tǒng)由于發(fā)生故障導(dǎo)致管2的流量系數(shù)C2按照以下的規(guī)律變化:C2={C20+(k-100)×0.06100<k<1302.4k>130(19)式中:C20為正常狀態(tài)下的C2的值。假