基于粒子濾波的故障預測

基于粒子濾波的故障預測

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