動力系統故障下運載火箭推力與質量辨識方法

于海森，譚述君，2，劉 浩，毛玉明

(1. 大連理工大學工業裝備結構分析優化與CAE軟件全國重點實驗室，大連 116024；2. 大連理工大學遼寧省空天飛行器前沿技術重點實驗室，大連 116024；3. 上海宇航系統工程研究所，上海 201109)

0 引 言

運載火箭是航天運輸系統的主要組成部分,是目前人類進入空間的主要工具,是發展空間技術、確保空間安全的基石[1]。運載火箭系統結構復雜,發射成本高,確保其安全可靠非常重要。在國內外失敗的發射案例中,相當一部分由動力系統故障[2]引起。動力系統故障一般表現為推力下降或關機[3]。運載火箭某臺發動機的推力下降會產生推力不平衡干擾力矩[4],而推進劑消耗速度發生變化則導致質量、質心等模型參數與預設值偏差變大,同時推力下降也導致發動機擺角效率的下降,即引起姿控系統伺服機構故障,從而引發嚴重的控制問題。事實上,動力系統發生故障有時并不是災難性的,即發生故障后,火箭仍具有繼續飛行并完成或部分完成發射任務的能力[5-6]。因此,及時地診斷出火箭推力大小和火箭當前質量有助于后續的控制重構和飛行任務調整[7]。

故障診斷常用的方法有基于信號處理的診斷方法[8]、基于知識的診斷方法[9]以及基于模型的故障診斷方法[10]。火箭的自主診斷技術主要采用基于模型的分析方法[11]。基于模型的方法已經得到了廣泛的研究,可細分為參數估計方法[12]、狀態估計方法[13]、等價空間法[14]等。在參數估計領域中,最小二乘法是一種得到最廣泛應用的參數估計方法,可用于動態、靜態、線性和非線性系統等[15]。例如文獻[16]采用漸消記憶的遞推最小二乘法對火箭推力參數進行了辨識。卡爾曼濾波是處理隨機估計應用最廣泛的方法[17],為擴展適用范圍,卡爾曼濾波出現了多種改進形式,包括擴展卡爾曼濾波、無跡卡爾曼濾波[18]、UD分解濾波等[19]。例如文獻[20]基于卡爾曼濾波和故障因子設計了一種針對液體火箭發動機穩態故障的檢測與診斷方法。文獻[21]提出了線性-二次后退時域算法進行故障檢測,該算法利用擴展卡爾曼濾波獲取殘差來估計推力損失程度,同時利用過載故障方向定位方法對故障發動機定位。

針對運載火箭推力參數辨識的方向主要有兩個方面:一種是固體火箭利用遙測的發動機燃燒室壓強和預先建立的內彈道程序估計性能參數[22];另一種是利用視加速度等信息對推力參數進行辨識,從而為后續的控制重構和飛行任務調整提供信息[23]。目前利用火箭視加速度等信息對火箭推力進行辨識的研究中,均假設火箭的質量已知,在這種情況下采用成熟的辨識方法顯然可以估計出推力大小。然而火箭實際飛行過程中質量與預設質量是有偏差的,尤其是當火箭發生推力下降的情況下,燃料的消耗速度跟預先設定速度會偏差很大,導致火箭實際質量與預設質量偏差很大,此時質量不能再看作是已知的。而且測量的視加速度信息中運載火箭的推力和質量耦合在一起,直接采用傳統的最小二乘或卡爾曼濾波方法都無法區分出推力參數和質量參數,給辨識工作帶來了很大的困難。同時推進系統故障將導致推力參數的突變,提高辨識方法的跟蹤能力也是一個難點。

精確辨識火箭動力系統故障下的質量和推力,可以更好地進行動力學模型修正、控制重構以及飛行任務調整等。因此,本文針對運載火箭推力與質量辨識相互耦合的問題,利用火箭慣性敏感器件測量的視加速度信息,提出了一種最小二乘法與卡爾曼濾波相結合的聯合校正辨識方法,可以精確有效地辨識火箭的推力參數與質量參數。進一步,提出了一種基于視加速度向量2范數差分的推力參數故障檢測方法,提高了對推力參數突變的跟蹤能力。

1 推力參數與質量參數的辨識問題

本文以兩級運載火箭二級飛行段為研究對象[24],假設發動機比沖已知。定義地心慣性坐標系OX1Y1Z1：原點在地心,X1軸在赤道平面內指向發射時刻本初子午線方向,Z1軸垂直赤道平面指向北極,Y1軸滿足右手定則。在地心慣性坐標系中建立火箭的二級飛行段動力學方程以及質量消耗方程如下:

(1)

(2)

(3)

本文采用地心慣性坐標系下的視加速度信息作為辨識的觀測量,則辨識的觀測方程可以表示為:

y=Fu/m+d

(4)

其中,d表示視加速度的測量噪聲。

為了后續描述上的方便,記觀測矩陣為

h(m,u)=u/m

(5)

則觀測方程也可以表示成:

y=h(m,u)F+d

(6)

當發動機出現動力系統故障后,運載火箭的推進劑質量消耗率會隨之發生變化,不能再采用質量的預設值;而式(4)或式(6)則表明,測量的視加速度信息中推力與質量耦合在一起的,因此需要將質量和推力同時作為未知量進行辨識,這也是本文研究的辨識問題。

2 聯合校正辨識方法

2.1 基本思路

視加速度測量方程式(6)可以看作是對F的觀測,然而,由于運載火箭的推力和質量均是未知的、待估計的量,它們是耦合在一起的,因此,無法直接采用最小二乘方法進行推力參數辨識。而且,當運載火箭推力發生故障后,推力參數會發生突變,因此需要采用跟蹤性強的辨識方法。

圖1 運載火箭推力參數與質量聯合校正辨識算法框圖

上文的聯合校正遞推辨識算法的關鍵在于可以在每個遞推步中實現對質量和推力參數的有效更新:在狀態(含質量參數)估計模塊中,采用的卡爾曼濾波算法包含基于模型的預估和基于測量信息的校正兩個環節,因此即使推力參數存在一定偏差,也可以綜合測量信息實現對狀態(含質量參數)的校正,這正是卡爾曼濾波算法的優點;在推力參數辨識模塊,采用最小二乘算法具有很強的抗干擾能力,可以有效降低質量參數偏差帶來的影響,同時采用漸消記憶的遞推最小二乘算法,通過降低歷史數據的權重來實現對時變推力參數的跟蹤。

進一步,由于當火箭發生動力系統故障時,推力參數發生突變是重要的特征,以前的觀測數據(視加速度)不再反映突變后的推力信息。因此為了提高對推力參數的跟蹤性能,在上述聯合校正遞推辨識算法框架的基礎上,又提出了推力突變的檢測方法和相應的修正措施。

2.2 聯合校正辨識算法的實現

運載火箭推力與質量聯合校正遞推辨識算法將卡爾曼濾波與漸消記憶的遞推最小二乘結合,本節具體介紹該方法的具體實現過程。

選取狀態變量x=[rxryrzvxvyvzm]T,則有火箭的質量m=x7。根據式(1)、式(2)以及式(6)將運載火箭動力學模型轉化為如下狀態空間形式:

(7)

其中:f(x,u,F)和g(x,u,F)表示動力學方程和觀測方程中的相應函數,具體表達形式如下:

(9)

卡爾曼濾波算法中對狀態的預估方程形式如下:

(10)

狀態預估值的協方差矩陣方程為:

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

(11)

式中:P(k+1|k)表示k+1時刻的狀態預估協方差矩陣;Φ(k+1|k)表示狀態轉移矩陣;P(k)表示k時刻的狀態估計協方差矩陣;Γ(k+1|k)表示離散噪聲分布矩陣;Qw(k)表示離散過程噪聲方差。其中,離散噪聲分布矩陣Γ(k+1|k)的計算公式為:

Γ(k+1|k)=Φ(k+1|k)C(k+1)Δt

(12)

式(11)中狀態轉移矩陣Φ(k+1|k)的計算公式為:

Φ(k+1|k)=eA(k+1)Δt

(13)

式中:A表示函數f(x,u,F)的雅可比矩陣,其具體形式為:

(14)

其中,元素的具體計算公式為:

(15)

卡爾曼濾波增益矩陣K的計算公式為:

K(k+1)=P(k+1|k)HT(k+1)[H(k+1)P(k+1|k)HT(k+1)+Rd(k+1)]-1

(16)

式中:Rd(k+1)表示離散觀測噪聲方差;H表示函數g(x,u,F)的雅可比矩陣,其具體形式表示為:

(17)

狀態估計協方差矩陣的計算公式如下:

P(k+1)=[Is-K(k+1)H(k+1)]P(k+1|k)

(18)

式中:Is表示對應狀態維度單位陣。

卡爾曼濾波的狀態估計值為:

(19)

(20)

根據k+1時刻運載火箭質量的辨識結果,通過遞推最小二乘估計k+1時刻的推力參數。由于推力是變化的,且當發生故障時,推力會產生突變,為了提高參數突變情況下的辨識方法跟蹤性能,采用漸消記憶的遞推最小二乘算法。引入漸消因子ρ,ρ∈(0, 1],通過降低歷史數據的權重進而降低歷史數據對當前時刻參數估計的影響。雖然漸消因子越小,辨識的跟蹤性越好,但是由于舊數據的減少,往往會降低對噪聲的抗干擾能力,因此漸消因子一般選取0.95≤ρ<1。

漸消記憶的遞推最小二乘增益矩陣為:

(21)

式中:Ib表示對應觀測維度單位陣;PLs表示信息矩陣,且有:

(22)

從而推力參數的估計值為:

(23)

運載火箭推力與質量的聯合校正遞推辨識算法流程如圖2所示。可以看出,每個循環周期內包含質量參數估計模塊和推力參數估計模塊。兩個模塊依次交替進行,質量參數估計模塊為推力參數估計的最小二乘算法提供觀測矩陣,而推力參數估計模塊為質量參數估計的卡爾曼濾波算法提供模型參數,二者互相校正,最終實現質量和推力參數的聯合辨識。

圖2 聯合校正辨識方法流程圖

2.3 動力系統故障的檢測

在聯合校正辨識算法中,推力參數辨識采用漸消記憶的遞推最小二乘算法,通過不斷相乘漸消因子降低對歷史數據的權重來實現對推力參數變化的跟蹤。因此當發生推力突變后,漸消因子需要一定的時間才能比較徹底的消除突變前測量數據的影響,這就導致也需要一定的時間才能跟蹤上突變后的推力。當然,通過降低漸消因子也可以提高跟蹤能力,但這同時會降低對噪聲的抗干擾能力。因此,如果能夠檢測出動力系統故障導致的推力突變的時刻,那么就可以直接刪除掉突變前的測量數據,從突變時刻重新開始最小二乘辨識,從而大大提高對推力參數突變的跟蹤能力。

本文的觀測信息是地心慣性坐標系下的視加速度,根據式(4),視加速度向量的2范數與推力成正比,當運載火箭推力發生突變時,由于運載火箭質量不會發生突變,因此推力的突變必然導致視加速度的突變。

理論上,發生突變時視加速度的變化率趨近無窮大。然而工程上只能利用采樣數據的差分來近似變化率,所以提出通過計算視加速度向量2范數的差分來判斷推力參數突變時刻。當視加速度向量2范數的差分在某一時刻的絕對值超過設定閾值時,認為在該時刻運載火箭的推力發生了突變。因此,得到判斷推力參數突變的條件如下:

(24)

不等式左邊表示視加速度向量2范數差分的絕對值,不等式右邊中η為給定變化率的閾值。當不等式成立時,認為推力參數在k+1時刻發生突變。在推力沒有突變的時間內,視加速度向量2范數差分的絕對值非常小,約10-3量級;考慮到實際測量的信號中含有噪聲,因此η的設置需要比10-3量級大一些,不妨設置為η=10-1。

實際上,視加速度變化率閾值η的設定可以根據發動機正常工作下推力的變化率和視加速度采樣周期進行預估。同時,在使用式(24)進行判斷時,如果僅采用某一個時刻的計算值進行判別,可靠性過低,容易引起誤判。因此一方面可將測量信號先通過低通濾波器、濾除測量噪聲的影響;另一方面,利用卡爾曼濾波算法本身具有的濾噪特性,也可以將視加速度測量值替換為對應的視加速度預估值進行判斷,也能在一定程度上提高檢測可靠性。

檢測出突變時刻后,在推力參數辨識的最小二乘遞推算法中就可以刪除突變前的測量數據,從而消除突變前測量數據對突變后推力參數辨識的影響。

有故障檢測的運載火箭質量和推力參數聯合校正遞推辨識算法如框圖3所示。與前面無故障檢測的算法框圖1相比,主要增加了故障檢測模塊和對推力參數辨識模塊的修正。具體修正如下。

圖3 有故障檢測的運載火箭推力參數與質量聯合校正辨識算法框圖

u(k+1)))-1

(25)

并利用k+1時刻測量信息直接進行k+1時刻推力參數的最小二乘估計,如下:

u(k+1))y(k+1)

(26)

3 仿真校驗

本節以火箭二級飛行段出現推力下降故障為研究對象,驗證本文所提出運載火箭推力參數和質量參數聯合校正遞推辨識算法的有效性。

火箭二級飛行段仿真時間選取239～249 s,采樣周期為Δt=0.01 s,假設火箭在243.9 s發生推力突變下降故障。動力學方程式(1)和(2)中的引力常數為μ=3.986 0×1014,火箭比沖為Isp=342.786 4 s,海平面處重力加速度為g0=9.806 7 m/s2。過程噪聲方差選取Qw=diag(0.01,0.01,0.01,0.01,0.01,0.01,0.01)。

通過在視加速度仿真值上施加高斯白噪聲來模擬含觀測噪聲的視加速度測量值[17]。為了驗證本文方法對噪聲不確定性的適應性,適當增大噪聲的方差,噪聲的方差取Rd=diag(0.05,0.05,0.05)。

圖4(a)、圖4(b)和圖4(c)給出了推力下降20%故障下的視加速度仿真模擬數據曲線,可以看出,當火箭推力發生突變后,視加速度也隨之突變。

圖4 推力下降20%故障下視加速度模擬值

-1.6×103,3.4×103,-1.7×103,3.8×104]T

(27)

且有

(28)

3.1 聯合校正辨識算法中的參數影響

本節驗證聯合校正辨識算法中漸消因子的大小對辨識結果的影響,以及驗證有故障檢測的效果。在無故障檢測的情況下,漸消因子分別選擇0.99,0.97與0.95;在有故障檢測情況下,漸消因子選擇0.97。則采用本文聯合校正辨識算法的運載火箭推力與推力誤差辨識結果如圖5、圖6所示。

圖5 不同漸消因子與有無故障檢測的推力參數辨識結果

圖6 不同漸消因子與有無故障檢測的推力參數辨識誤差

從圖5、圖6中的不同漸消因子-無故障檢測曲線可以看出,漸消因子選取的越小,在推力發生故障突變后,辨識的跟蹤性越好;在辨識結果跟蹤上理論推力的這些時間內,推力辨識誤差小于1%,漸消因子選取越小,推力辨識誤差越大。因此,漸消因子選取越小雖然可以有效提高辨識的跟蹤性,但是會增加推力辨識誤差,漸消因子選取不宜過大或過小。從圖中漸消因子均取0.97時,有故障檢測與無故障檢測的兩條曲線可以看出,在推力發生故障突變后,有故障檢測的推力辨識結果可以立刻跟蹤上理論的推力;在243.9 s至246.4 s時間內,有故障檢測的推力辨識誤差與無故障檢測的推力辨識誤差幾乎一致,且有故障檢測的結果在推力辨識結果跟蹤突變后推力的這段時間里推力誤差也只有2.5%。

質量與質量誤差的辨識結果如圖7、圖8所示。從圖7、圖8可以看出,漸消因子越小,在故障下推力突變后,質量跟蹤越好;在漸消因子均取0.97,有故障檢測和無故障檢測的結果可以看出,在故障下推力突變后,有故障檢測的質量誤差更低,小于0.002 5%。這是因為在推力突變后,減小漸消因子以及有故障檢測的推力參數估計值較快得接近理論值,進而在聯合校正辨識算法中通過推力估計結果來辨識的質量結果更精確。

圖7 不同漸消因子與有無故障檢測的質量辨識結果

圖8 不同漸消因子與有無故障檢測的質量辨識誤差

3.2 聯合校正辨識方法的有效性

本節驗證聯合校正辨識方法的有效性,包括不同推力初值對辨識結果的影響以及不同推力下降故障下的辨識結果。本節采用漸消因子取0.97、有故障檢測的方法,推力初值在理論值的基礎上分別取0.9倍理論值和1.1倍理論值,則有不同推力初值下的運載火箭推力與質量辨識結果如圖9、圖10所示。

圖9 不同推力初值下的推力辨識結果

圖10 不同推力初值下的質量辨識結果

從圖9中可以看出,當推力初值有偏差時,聯合校正辨識方法可以較快地跟蹤上理論的推力;從圖10中可以看出,當推力初值有偏差時,聯合校正辨識方法對質量的辨識效果依然理想。

當火箭發生推力突變下降故障后,往往在故障發生后的一段時間并不是保持不變,而是變化的,即二次故障。因此,本節驗證本文的有故障檢測的聯合校正辨識方法對不同故障大小以及二次故障下的適應性。通過假設在推力發生突變下降后經過一段正弦波動變化后保持在某個數值來模擬二次故障,推力初值取理論值。

則有在不同故障大小以及二次故障情況下的辨識結果如圖11至圖14所示。

圖11 不同故障下的推力辨識結果

從圖11、圖12中可以看出,當發生不同大小故障或二次故障時,有故障檢測的聯合校正辨識仍然可以快速跟蹤上變化的推力,在驗證不同大小故障或二次故障時,需獲得不同的觀測值,由于噪聲是隨機的高斯白噪聲,因此即便在發生故障時間前,推力辨識誤差仍不相同。在不同推力大小故障下,除了在故障時刻附近推力辨識誤差有時會增加達到2%左右,其它時間內均小于0.6%。

圖12 不同故障下的推力辨識誤差

圖13、圖14表明,不同推力大小的故障由于質量消耗速度不同,因此故障后質量曲線不同;質量辨識誤差仍然非常小。

圖14 不同故障下的質量辨識誤差

4 結 論

本文以運載火箭二級飛行段為研究對象,研究了動力系統故障下運載火箭推力與質量的高精度辨識問題,主要成果如下:

1) 提出了一種最小二乘與卡爾曼濾波相結合的運載火箭推力與質量聯合校正辨識方法,該方法解決了質量未知情況下的推力辨識問題。

2) 提出了利用視加速度向量2范數差分檢測推力參數突變時刻的方法及其改進辨識策略,顯著提高了對推力參數突變的跟蹤能力。

3) 以運載火箭推力下降不同大小的故障以及二次故障為例,通過深入仿真驗證了本文提出的方法可以有效地辨識出動力系統故障下的運載火箭的推力與質量,具有良好的跟蹤性與精度。