葉 松，陳 曦，熊寸平

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

0 引言

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

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

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

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

1 數(shù)學模型

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

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

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

(1)

(2)

(3)

式中，γ為滾轉角，?為俯仰角，ψ為偏航角。

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

(4)

記

可得

(5)

為簡化計算，取量測變量為除發(fā)動機推力外的12個狀態(tài)變量

yi=xi+vi

(6)

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

上述模型進行線性化和離散化并寫為如下非線性形式

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為離散時間變量，x為狀態(tài)變量，u為輸入向量，y為輸出變量。非線性函數(shù)F，H為狀態(tài)的一階連續(xù)偏導數(shù)，Γ為已知矩陣。系統(tǒng)噪聲v(k)、測量噪聲e(k)為高斯白噪聲，互不相關。在此基礎上，設計擴展狀態(tài)觀測器(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)

預報誤差協(xié)方差陣為

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

1,k)+Q(k)

(10)

狀態(tài)估計誤差協(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)的一步預測值如下

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

(12)

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

(13)

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

“我上一次到訪查謨-克什米爾大君的斯里那加宮殿時， 他們端出了三十個盤子，如果我說任何一個盤子上的寶石都能在市場賣得到一百萬元，恐怕是遠遠低估了這些寶物的美及它所代表的財富?！?/p>

將模型實時觀測到的某狀態(tài)量與EKF方法下狀態(tài)量的計算值作比較，可得到系統(tǒng)殘差預測向量

(14)

式中,r為1×n維。

(15)

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

2 故障診斷算法設計

考慮運載火箭推力下降故障，記推力下降程度為f，f∈(0,1)，0表示推力下降程度為0%，1表示推力下降程度為100%。當發(fā)動機發(fā)生推力故障時，推力Pf如下

Pf=(1-f)P

(16)

2.1 殘差表示

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

(17)

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

2.2 靈敏度矩陣求解

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

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

(18)

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

(19)

無故障情況下的殘差為

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

(20)

由于實際中殘差無法直接對故障求導，可通過下式求解靈敏度矩陣參數(shù)

(21)

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

(22)

2.3 估算算法

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

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

Y=Sf+ew

(23)

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

(24)

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

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

R=F2

(25)

F具有預期故障強度的含義?？梢曰诠收闲再|的知識來設置故障強度，實際中通過多次測試獲取。

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

(26)

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

2.4 故障定位

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

圖2 4臺發(fā)動機分別故障下過載ny響應曲線Fig.2 ny response with each individual engine fault

圖3 4臺發(fā)動機分別故障下過載nz響應曲線Fig.3 nz response with each individual engine fault

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

2.5 算法流程

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

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

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

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

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

3 仿真分析

3.1 單臺分析

(1)不同采樣長度N對仿真結果的影響

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

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

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

圖5 不同采樣長度下故障檢測結果Fig.5 Fault detection results under different sampling lengths

圖6 不同采樣長度下故障檢測結果放大圖Fig.6 Enlarged view of fault detection results under different sampling lengths

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

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

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

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

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

考慮到故障檢測的準確度，避免二次誤報，因此選取R=202至502量級較為合適。

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

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

(a) 無故障

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

(a) 無故障

(b) 30%故障

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

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

圖11 故障程度與故障估計值f關系Fig.11 Relationship between fault degree and fault estiomate f

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

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

(27)

利用上式將故障估計值與預設的故障程度對應，表1給出部分驗證結果。

表1 故障程度的估計結果

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

3.2 4臺發(fā)動機故障檢測

基于過載在故障時刻突變方向與發(fā)動機位置相關的思路進行多次故障檢測驗證，結果如表2所示。

表2 故障發(fā)動機編號檢測

上述結果表明，利用根據(jù)ny,nz的突變方向能夠準確判定發(fā)動機故障臺的編號,表明該思路能夠用于4臺發(fā)動機的故障定位。

4 結論

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

此外，給出根據(jù)過載的突變方向定位故障的策略，與估計算法檢測的故障時間相結合，實現(xiàn)了發(fā)動機故障的準確定位。