基于線性二次滾動(dòng)時(shí)域法的運(yùn)載火箭發(fā)動(dòng)機(jī)推力故障診斷

葉 松，陳 曦，熊寸平

(北京航天自動(dòng)控制研究所，北京 100854)

0 引言

隨著航空航天技術(shù)的快速發(fā)展，運(yùn)載火箭的組成越來(lái)越復(fù)雜，其可靠性和安全性問(wèn)題日益突出[1]。故障診斷技術(shù)的應(yīng)用，能夠及時(shí)檢測(cè)系統(tǒng)故障，提高系統(tǒng)的可靠性，同時(shí)為后續(xù)容錯(cuò)控制奠定基礎(chǔ)。動(dòng)力系統(tǒng)作為運(yùn)載火箭的一個(gè)關(guān)鍵系統(tǒng)，為運(yùn)載火箭提供動(dòng)力，推動(dòng)運(yùn)載火箭前進(jìn)。推力下降故障是動(dòng)力系統(tǒng)故障類(lèi)型之一，會(huì)導(dǎo)致系統(tǒng)性能退化甚至不穩(wěn)定，因此研究推力下降故障診斷具有重要意義。

故障診斷包含故障檢測(cè)、估計(jì)與故障定位。故障檢測(cè)反映故障的發(fā)生，故障估計(jì)表征故障程度，故障定位表征故障發(fā)生的位置。故障診斷常用基于動(dòng)態(tài)模型的方法[2-3]、基于信號(hào)處理[4-5]的方法以及基于知識(shí)[6-7]的方法。基于成熟的運(yùn)載火箭建模理論，采用基于動(dòng)態(tài)模型的方法更具有工程應(yīng)用價(jià)值。

基于動(dòng)態(tài)模型的非線性系統(tǒng)故障診斷方法主要有狀態(tài)估計(jì)法、參數(shù)估計(jì)法。前者通常對(duì)非線性模型的某些特征點(diǎn)進(jìn)行線性化，利用建模誤差作為未知輸入并設(shè)計(jì)未知輸入觀測(cè)器獲取殘差，進(jìn)行故障診斷，后者通常對(duì)非線性模型利用非線性參數(shù)估計(jì)的方法進(jìn)行故障診斷。王青等[8]采用狀態(tài)估計(jì)法，基于線性化后的飛行器模型，設(shè)計(jì)降階故障檢測(cè)濾波器，在減少計(jì)算量的同時(shí)保證了故障檢測(cè)的快速性。符文星等[9]研究了基于強(qiáng)跟蹤濾波器的參數(shù)估計(jì)方法，對(duì)運(yùn)載火箭的推力參數(shù)進(jìn)行了正確估計(jì)。

本文結(jié)合狀態(tài)估計(jì)法與參數(shù)估計(jì)法，提出了一種基于擴(kuò)展卡爾曼濾波器生成殘差，采用線性二次滾動(dòng)時(shí)域算法[10]估計(jì)故障，并根據(jù)過(guò)載方向定位故障的方法，對(duì)發(fā)動(dòng)機(jī)推力下降故障進(jìn)行在線檢測(cè)與定位。最后將該方法進(jìn)行數(shù)學(xué)仿真，驗(yàn)證了算法的有效性。

1 數(shù)學(xué)模型

本文以運(yùn)載火箭為例，火箭動(dòng)力系統(tǒng)包含4臺(tái)搖擺發(fā)動(dòng)機(jī)，其布局如圖1所示。

圖1 火箭發(fā)動(dòng)機(jī)布局圖Fig.1 Distribution of launch vehicle power system

受地球自轉(zhuǎn)的影響，發(fā)射坐標(biāo)系為動(dòng)坐標(biāo)系，這樣存在附加哥氏力和力矩、附加相對(duì)力和力矩的影響，但本方法中這些影響均可以視作濾波器噪聲，因此將發(fā)射坐標(biāo)系視作慣性系，從而忽略附加哥氏力和力矩、附加相對(duì)力和力矩的影響，以突出研究的重點(diǎn)。簡(jiǎn)化后的運(yùn)載火箭動(dòng)力學(xué)及運(yùn)動(dòng)學(xué)方程[11]

(1)

(2)

(3)

式中，γ為滾轉(zhuǎn)角，?為俯仰角，ψ為偏航角。

取狀態(tài)變量如下

(4)

記

可得

(5)

為簡(jiǎn)化計(jì)算，取量測(cè)變量為除發(fā)動(dòng)機(jī)推力外的12個(gè)狀態(tài)變量

yi=xi+vi

(6)

式中，i=1,2,…,12。v為白噪聲序列。

上述模型進(jìn)行線性化和離散化并寫(xiě)為如下非線性形式

x(k+1)=F(k,u(k),x(k))+Γ(k)v(k)

y(k+1)=H(k+1,x(k+1))+e(k+1)

(7)

式中，整數(shù)k≥0為離散時(shí)間變量，x為狀態(tài)變量，u為輸入向量，y為輸出變量。非線性函數(shù)F，H為狀態(tài)的一階連續(xù)偏導(dǎo)數(shù)，Γ為已知矩陣。系統(tǒng)噪聲v(k)、測(cè)量噪聲e(k)為高斯白噪聲，互不相關(guān)。在此基礎(chǔ)上，設(shè)計(jì)擴(kuò)展?fàn)顟B(tài)觀測(cè)器(Extended Kalman Filter，EKF)如下

X(k+1)=X(k+1,k)+K(k+1)[Z(k+

1)-H(k+1)X(k+1,k)]

(8)

濾波增益陣K為

K(k+1)=P(k+1,k)HT(k+1)·

1)+R(k+1)]-1

(9)

預(yù)報(bào)誤差協(xié)方差陣為

P(k+1,k)=Φ(k+1,k)P(k)ΦT(k+

1,k)+Q(k)

(10)

狀態(tài)估計(jì)誤差協(xié)方差陣為

P(k+1)=[I-K(k+1)H(k+1)]P(k+1,k)·
[I-K(k+1)H(k+1)]T+
K(k+1)R(k+1)K(k+1)T

(11)

式中,X(k+1,k)為狀態(tài)的一步預(yù)測(cè)值如下

X(k+1,k)=f(X(t))|X(t)=X(tk)

(12)

Φ(k+1,k)為狀態(tài)轉(zhuǎn)移矩陣如下

(13)

H(k+1)為量測(cè)轉(zhuǎn)移矩陣，根據(jù)狀態(tài)變量和量測(cè)變量的關(guān)系，可以得到H(k+1)為單位陣。

“我上一次到訪查謨-克什米爾大君的斯里那加宮殿時(shí)， 他們端出了三十個(gè)盤(pán)子，如果我說(shuō)任何一個(gè)盤(pán)子上的寶石都能在市場(chǎng)賣(mài)得到一百萬(wàn)元，恐怕是遠(yuǎn)遠(yuǎn)低估了這些寶物的美及它所代表的財(cái)富?！?/p>

將模型實(shí)時(shí)觀測(cè)到的某狀態(tài)量與EKF方法下?tīng)顟B(tài)量的計(jì)算值作比較，可得到系統(tǒng)殘差預(yù)測(cè)向量

(14)

式中,r為1×n維。

(15)

式中,rωz為1×n維。

2 故障診斷算法設(shè)計(jì)

考慮運(yùn)載火箭推力下降故障，記推力下降程度為f，f∈(0,1)，0表示推力下降程度為0%，1表示推力下降程度為100%。當(dāng)發(fā)動(dòng)機(jī)發(fā)生推力故障時(shí)，推力Pf如下

Pf=(1-f)P

(16)

2.1 殘差表示

殘差數(shù)據(jù)Y為某個(gè)時(shí)間間隔內(nèi)殘差向量rωz的采樣值，以如下形式表示

(17)

式中，tb為數(shù)據(jù)采樣開(kāi)始的時(shí)間，Δ為采樣間隔，N為樣本的數(shù)量，f為故障程度。向量Y的大小為(N+1)×1。

2.2 靈敏度矩陣求解

假定動(dòng)力系統(tǒng)運(yùn)行正常，則零故障記為Y(tb,N,Δ;0)；若發(fā)生故障，獲得的殘差數(shù)據(jù)將進(jìn)一步表示為Y(tb,N,Δ;f)。將故障殘差數(shù)據(jù)線性化處理，分為無(wú)故障殘差部分和故障相關(guān)部分，如下

Y(tb,N,Δ;f)≈Y(tb,N,Δ;0)+S(tb,N,Δ)f

(18)

式中，S為關(guān)于Y(tb,N,Δ;f)的最后一個(gè)參數(shù)f的(N+1)×1雅可比矩陣，稱(chēng)為靈敏度矩陣。

(19)

無(wú)故障情況下的殘差為

r0(t)=r(t;f=0)

(20)

由于實(shí)際中殘差無(wú)法直接對(duì)故障求導(dǎo)，可通過(guò)下式求解靈敏度矩陣參數(shù)

(21)

式中，d是有限差分步長(zhǎng)，與采樣時(shí)間保持一致。e為單位故障向量，表征1%的推力下降程度。通過(guò)下式對(duì)特征數(shù)據(jù)進(jìn)行采樣來(lái)計(jì)算矩陣S。

(22)

2.3 估算算法

前文介紹了殘差向量Y以及靈敏度矩陣S的求解，下文對(duì)故障估算算法進(jìn)行介紹。

對(duì)于實(shí)際系統(tǒng)，由于傳感器噪聲和建模誤差的存在，殘差不為零。將所有這些因素匯總在一個(gè)廣義的“噪聲”ew數(shù)據(jù)向量中。假設(shè)殘差數(shù)據(jù)向量Y的統(tǒng)計(jì)模型如下

Y=Sf+ew

(23)

通過(guò)對(duì)該殘差數(shù)據(jù)模型中的變量進(jìn)行統(tǒng)計(jì)假設(shè)，可以得出故障參數(shù)的估計(jì)值。一種設(shè)計(jì)方法是將統(tǒng)計(jì)參數(shù)視為估算算法的調(diào)整參數(shù)，可以對(duì)這些參數(shù)進(jìn)行調(diào)整，以實(shí)現(xiàn)算法的合理實(shí)用性能?；诖耍僭O(shè)噪聲ew是具有協(xié)方差矩陣V的正態(tài)分布向量。同時(shí)假定要估計(jì)的未知故障向量f是獨(dú)立于(噪聲ew)的隨機(jī)變量，并且為具有協(xié)方差矩陣R的正態(tài)分布向量。在這種情況下，可以從噪聲數(shù)據(jù)中獲得故障向量f的最佳統(tǒng)計(jì)估計(jì)值

(24)

在此過(guò)程中如果忽略傳感器噪聲，唯一建模誤差是殘差與故障矢量關(guān)系式中一階近似的線性化誤差。在實(shí)施估計(jì)算法時(shí)，協(xié)方差矩陣V設(shè)為單位矩陣，這對(duì)應(yīng)于所有觀察通道中存在的相同幅度的白噪聲。

將故障的協(xié)方差矩陣R選擇為

R=F2

(25)

F具有預(yù)期故障強(qiáng)度的含義。可以基于故障性質(zhì)的知識(shí)來(lái)設(shè)置故障強(qiáng)度，實(shí)際中通過(guò)多次測(cè)試獲取。

故障向量可以從整個(gè)飛行過(guò)程中收集的數(shù)據(jù)集作為一個(gè)批處理估計(jì)獲得，但是，可能需要在整個(gè)飛行過(guò)程中計(jì)算和更新估算值，以便在線估算和觀察不斷發(fā)展的故障。為滿足此需求，以滾動(dòng)時(shí)域的形式實(shí)現(xiàn)用于故障參數(shù)估計(jì)的GLS算法。在每個(gè)時(shí)間步，滾動(dòng)時(shí)域算法使用N個(gè)最新數(shù)據(jù)點(diǎn)來(lái)計(jì)算故障參數(shù)矢量的估計(jì)值。此估算形式如下

(26)

由于一種故障類(lèi)型對(duì)應(yīng)一種故障強(qiáng)度，因此針對(duì)單一的推力下降故障，故障協(xié)方差陣維數(shù)恒為1×1。因此選用單一殘差向量ωz保證公式(25)的維數(shù)匹配，滿足對(duì)故障敏感的殘差均可選用。

2.4 故障定位

基于殘差信息的故障檢測(cè)僅能反映整體的故障時(shí)刻以及故障程度，當(dāng)4臺(tái)發(fā)動(dòng)機(jī)的某一臺(tái)發(fā)生故障時(shí)，僅基于殘差無(wú)法區(qū)分4臺(tái)發(fā)動(dòng)機(jī)的故障信息。因此提出一種基于故障時(shí)刻的過(guò)載方向信息定位故障的方法。進(jìn)行仿真分析，當(dāng)動(dòng)力系統(tǒng)4臺(tái)發(fā)動(dòng)機(jī)分別在25 s發(fā)生30%推力下降故障時(shí)，ny(箭體系y1方向過(guò)載),nz(箭體系z(mì)1方向過(guò)載)狀態(tài)如圖2、圖3所示。

圖2 4臺(tái)發(fā)動(dòng)機(jī)分別故障下過(guò)載ny響應(yīng)曲線Fig.2 ny response with each individual engine fault

圖3 4臺(tái)發(fā)動(dòng)機(jī)分別故障下過(guò)載nz響應(yīng)曲線Fig.3 nz response with each individual engine fault

由圖2、圖3可以得知，過(guò)載ny在第2，3臺(tái)發(fā)動(dòng)機(jī)故障響應(yīng)相似，第1，4臺(tái)發(fā)動(dòng)機(jī)故障響應(yīng)相似。過(guò)載nz的第1，2臺(tái)發(fā)動(dòng)機(jī)故障響應(yīng)相似，第3，4臺(tái)發(fā)動(dòng)機(jī)故障響應(yīng)相似。因此檢測(cè)出故障發(fā)生后，可根據(jù)ny，nz的突變方向，綜合判定故障發(fā)動(dòng)機(jī)的編號(hào)。利用這種特性，可以應(yīng)用于故障檢測(cè)方法中。

2.5 算法流程

提出的方法分為生成殘差、估計(jì)算法與故障定位3部分。

第1部分在推力無(wú)故障情況下，以P，vx，vy，vz，ωx，ωy，ωz，x，y，z，φ，ψ，γ作為輸入數(shù)據(jù)，基于預(yù)測(cè)模型，利用EKF方法得到除推力P外的預(yù)測(cè)數(shù)據(jù)計(jì)算值，計(jì)算殘差。

由上述狀態(tài)量計(jì)算的殘差不能完全反映推力下降故障，通過(guò)分析殘差數(shù)據(jù)vx，vy，vz，ωx，ωy，ωz,x,y,z,φ,Ψ,γ的曲線發(fā)現(xiàn)，推力下降故障僅對(duì)vx，vy，vz，ωx，ωy，ωz影響非常比較明顯。其中，vy，ωy與ωz在某些時(shí)刻的殘差有明顯突變，且故障后一段時(shí)間能夠?qū)⑾到y(tǒng)控制到穩(wěn)定狀態(tài)，符合實(shí)際情況，因此采用ωz計(jì)算殘差。

將無(wú)故障的狀態(tài)參數(shù)ωz與預(yù)測(cè)值做差，即得到殘差r0。同理，在1%故障下ωz實(shí)際值與預(yù)測(cè)值做差對(duì)應(yīng)生成殘差rf_1%。

在推力故障情況下，通過(guò)推力P，由EKF計(jì)算故障下的狀態(tài)參數(shù)ωz，最終生成rf。

3 仿真分析

3.1 單臺(tái)分析

(1)不同采樣長(zhǎng)度N對(duì)仿真結(jié)果的影響

設(shè)定50 s時(shí)發(fā)生50%的推力下降，考慮到控制器在10 ms左右即可控制住故障，采樣長(zhǎng)度N分別為1，5，10，采樣間隔Δ為0.01 s,協(xié)方差陣R=202,將歷史緩沖區(qū)取NΔ對(duì)應(yīng)值，分別為0.01，0.05，0.1 s，仿真結(jié)果如圖5、圖6所示。

由圖5、圖6可以看出，采樣長(zhǎng)度N越大，故障估計(jì)值f恢復(fù)正常的拍數(shù)越少，越不利于觀測(cè)到故障；采樣長(zhǎng)度過(guò)小，故障估計(jì)值f有較大損失。為了保證故障檢測(cè)的及時(shí)性及準(zhǔn)確性，采樣長(zhǎng)度為5較為合適，實(shí)際應(yīng)用中可根據(jù)此分析方法確定實(shí)際采樣長(zhǎng)度。

圖4 估算算法流程圖Fig.4 Flow chart of estimation algorithm

圖5 不同采樣長(zhǎng)度下故障檢測(cè)結(jié)果Fig.5 Fault detection results under different sampling lengths

圖6 不同采樣長(zhǎng)度下故障檢測(cè)結(jié)果放大圖Fig.6 Enlarged view of fault detection results under different sampling lengths

(2)不同故障協(xié)方差矩陣R對(duì)仿真結(jié)果的影響

設(shè)定50 s時(shí)發(fā)生50%的推力下降，采樣長(zhǎng)度N為5，采樣間隔Δ為0.01 s,協(xié)方差陣R分別為0.012，12，202，5002，5 0002，10 0002,歷史緩沖區(qū)為0.05 s，圖7給出50 s附近故障估計(jì)結(jié)果。

圖7 50 s附近不同協(xié)方差陣下的故障估計(jì)結(jié)果(大范圍)Fig.7 Fault detection results under different covariance matrices around 50 s(wide range)

由圖7可以看出，協(xié)方差陣取得越大，f反映故障發(fā)生強(qiáng)度越顯著，但取得過(guò)大f會(huì)出現(xiàn)二次增大，可能導(dǎo)致故障檢測(cè)的誤報(bào)。根據(jù)上述結(jié)果進(jìn)一步確定合適的協(xié)方差陣大小，協(xié)方差陣R分別為12，52，202，502，1002,仿真結(jié)果如圖8所示。

圖8 50 s附近不同協(xié)方差陣下的故障估計(jì)結(jié)果(小范圍)Fig.8 Fault detection estimation results under different covariance matrices around 50 s(small scale)

考慮到故障檢測(cè)的準(zhǔn)確度，避免二次誤報(bào)，因此選取R=202至502量級(jí)較為合適。

(3)不同推力下降程度仿真

基于上述仿真，給出同一時(shí)刻(20 s)下，推力下降故障分別為0%、30%、50%時(shí)的仿真結(jié)果，取采樣長(zhǎng)度N=5，采樣間隔Δ為0.01 s,協(xié)方差陣R=202,將前0.05 s數(shù)據(jù)作為歷史緩沖區(qū)，故障殘差仿真結(jié)果如圖9、圖10所示，故障向量f的估計(jì)如圖11所示。

(a) 無(wú)故障

(c) 50%故障圖9 20 s時(shí)刻不同故障程度殘差仿真結(jié)果Fig.9 Simulation results of residuals with different fault degrees at 20 s

(a) 無(wú)故障

(b) 30%故障

(c) 50%故障圖10 20 s時(shí)刻不同故障程度故障檢測(cè)仿真結(jié)果Fig.10 Simulation results of fault detection with different fault degrees at 20 s

由上述結(jié)果可得，同一時(shí)刻發(fā)生不同故障程度，f估計(jì)值不同。在20 s時(shí)刻，以5%故障下降程度為間隔，仿真獲取故障時(shí)刻估計(jì)的峰值為fmax，繪制曲線如圖11所示，可以判斷兩者存在近似的線性關(guān)系。

圖11 故障程度與故障估計(jì)值f關(guān)系Fig.11 Relationship between fault degree and fault estiomate f

由圖11可以看出,兩者之間存在線性關(guān)系，由于該曲線僅針對(duì)50 s時(shí)刻故障，遍歷時(shí)間與故障程度，形成一個(gè)二維表格，根據(jù)辨識(shí)結(jié)果在表格內(nèi)進(jìn)行插值，即可得到推力下降程度。

以推力下降程度lossP為因變量、f估計(jì)峰值為自變量進(jìn)行線性擬合，可得到20 s時(shí)刻下估計(jì)值f與實(shí)際推力下降程度lossP關(guān)系式

(27)

利用上式將故障估計(jì)值與預(yù)設(shè)的故障程度對(duì)應(yīng)，表1給出部分驗(yàn)證結(jié)果。

表1 故障程度的估計(jì)結(jié)果

上述結(jié)果表明，該方法能夠較為準(zhǔn)確地獲得故障損失程度。在檢測(cè)出故障后，由于控制器的控制作用，所以突變恢復(fù)正常。殘差曲線中40，50 s附近也產(chǎn)生了與推力下降故障不相關(guān)的突變。由表1可以看出，這些突變?cè)斐傻膄估計(jì)值的變化遠(yuǎn)小于推力下降故障造成的估計(jì)結(jié)果的變化，表明該方法對(duì)推力下降故障具有較高的靈敏度。

3.2 4臺(tái)發(fā)動(dòng)機(jī)故障檢測(cè)

基于過(guò)載在故障時(shí)刻突變方向與發(fā)動(dòng)機(jī)位置相關(guān)的思路進(jìn)行多次故障檢測(cè)驗(yàn)證，結(jié)果如表2所示。

表2 故障發(fā)動(dòng)機(jī)編號(hào)檢測(cè)

上述結(jié)果表明，利用根據(jù)ny,nz的突變方向能夠準(zhǔn)確判定發(fā)動(dòng)機(jī)故障臺(tái)的編號(hào),表明該思路能夠用于4臺(tái)發(fā)動(dòng)機(jī)的故障定位。

4 結(jié)論

本文針對(duì)運(yùn)載火箭動(dòng)力系統(tǒng)的推力下降故障，提出了一種線性二次滾動(dòng)時(shí)域算法用于估計(jì)故障時(shí)間與故障程度，結(jié)合估計(jì)算法檢測(cè)的故障時(shí)間提出一種故障定位的策略。仿真結(jié)果表明，該方法能夠成功檢測(cè)發(fā)動(dòng)機(jī)故障，并對(duì)該類(lèi)型故障具有較高靈敏度。

此外，給出根據(jù)過(guò)載的突變方向定位故障的策略，與估計(jì)算法檢測(cè)的故障時(shí)間相結(jié)合，實(shí)現(xiàn)了發(fā)動(dòng)機(jī)故障的準(zhǔn)確定位。