逃逸飛行器軌跡魯棒性快速評估及制導(dǎo)

李 臻，楊 彪，荊武興，高長生，常武權(quán)

(1.哈爾濱工業(yè)大學(xué)航天工程系，哈爾濱150001；2.北京宇航系統(tǒng)工程研究所，北京100076)

1 引言

逃逸飛行器是載人航天器上用以保障航天員生命安全的一個重要組成部分。國外最早在開展阿波羅計劃時便已經(jīng)進行相關(guān)研究[1]；國內(nèi)在設(shè)計CZ-2F載人運載火箭時也開展過相關(guān)試驗工作[2-4]。逃逸飛行器設(shè)計過程涉及逃逸模式選定、逃逸能力評估等各個方面。朱仁璋等[5]進行過對神舟系列飛船的分離動力學(xué)分析，許鋒等[6]研究了逃逸主推力及結(jié)構(gòu)彈性變形對安全逃逸距離的影響，李家文等[7]分析了各種工況下爆炸沖擊波對逃逸飛行器的損壞情況。而關(guān)于逃逸飛行器軌跡設(shè)計與制導(dǎo)方法，目前鮮見文章詳細論述。

考慮到逃逸飛行器逃逸過程中，各種偏差對軌跡魯棒性和制導(dǎo)會帶來較大影響。針對該類飛行器，本文利用誤差傳播法與線性二次調(diào)節(jié)器(Linear Quadratic Regulator，LQR)控制方法，提出一種工程實用性強的軌跡設(shè)計與制導(dǎo)方法。其中，誤差傳播法用于軌跡魯棒性快速評估，難點在于確定時變系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣。陳國強[8]最早在研究引力異常對慣性制導(dǎo)影響時，給出過只考慮重力場作用時，狀態(tài)轉(zhuǎn)移矩陣的近似解析解；鄭偉[9]在研究地球物理攝動對導(dǎo)彈命中精度影響時，采用伴隨矩陣近似求解該矩陣解析解。而LQR理論最早用于再入飛行器的彈道跟蹤中，Dukeman[10]，Zhou等[11]設(shè)計的狀態(tài)調(diào)節(jié)器取得了良好效果；張大元等[12]針對防空導(dǎo)彈，也設(shè)計過基于LQR的彈道跟蹤制導(dǎo)律。

本文針對逃逸過程中的偏差影響，先利用誤差傳播法快速評估軌跡魯棒性得到標稱軌跡，進而采用LQR方法進行軌跡跟蹤制導(dǎo)。本文在設(shè)計誤差傳播法時，系統(tǒng)模型復(fù)雜度增加，區(qū)別于文獻[8]、[9]，不再能給出一個狀態(tài)轉(zhuǎn)移矩陣的解析解形式，因此采用數(shù)值解求解辦法，以較短時間獲得較高計算精度；此外，在不同背景下，LQR應(yīng)用方式略有區(qū)別，通過合理減少狀態(tài)參數(shù)個數(shù)，成功實現(xiàn)多變量跟蹤。

2 逃逸飛行器模型及軌跡設(shè)計

2.1 動力學(xué)模型

本文研究對象為具有軸對稱布局的逃逸飛行器。其直到救生塔分離前的這一飛行過程，可分為主動段和被動段，被動段結(jié)束時需要滿足分離約束。其工作時序如下：在接收到箭體分離指令后，逃逸飛行器從運載火箭主體脫離，逃逸主發(fā)動機工作，經(jīng)過短時間延遲t0后，開始調(diào)整飛行器姿態(tài)。主發(fā)動機只能工作小段時間，t2時間后進入無動力被動飛行段。飛行器的制導(dǎo)只在主動段進行，經(jīng)過一段時間的被動飛行后，救生塔準備分離。

地面發(fā)射系下建立逃逸飛行器的精確動力學(xué)模型見式(1)。

式中，為地面發(fā)射系下速度、加速度；aT為推力加速度，主動段t2時間內(nèi)沿彈體方向為一固定值，被動段內(nèi)大小為0；aR為氣動加速度，為簡化分析本文不考慮側(cè)向運動；g為重力加速度，考慮到J2擾動項；ak、aω分別為科氏加速度與離心加速度，由地球自轉(zhuǎn)引起。各項在地面發(fā)射系下具體表達式如式(2)。

其中引力相關(guān)項表達式見式(3)[9]。

以上方程中，P為體系下軸向推力矢量；C為氣流系下氣動系數(shù)矩陣，包含阻力系數(shù)Cx和升力系數(shù)Cy，可當作迎角α和馬赫數(shù)Ma的函數(shù)，用fxMa，α()與gyMa，α()擬合；q為動壓，SM為逃逸飛行器特征面積；gr與gω分別為重力加速度在徑向與自轉(zhuǎn)方向上的分量；R為地心至飛行器矢量，μ為地球引力常量，J為地球扁率修正項；ae為地球赤道半徑；φ為飛行器處地心緯度。

Γ1與Γ2分別為彈體系與氣流系到發(fā)射系的坐標轉(zhuǎn)換陣。為簡化研究，可假定飛行器沿射面內(nèi)飛行不考慮滾轉(zhuǎn)、偏航，則有式(4)：

φ為彈道傾角，近似有φ＝α＋θ；σ與θ分別為速度偏角與速度傾角，見式(5)。

2.2 誤差線性化模型

在小擾動假設(shè)條件下，每時刻飛行器狀態(tài)誤差為小量。針對逃逸飛行器動力學(xué)模型式(1)，以時間為自變量，可在標準彈道上特征點附近線性化展開，見式(6)。

Δxs，i為第i個狀態(tài)偏差，aij是與該狀態(tài)相關(guān)的第j個加速度項。

取系統(tǒng)狀態(tài)量χ ＝[ΔrΔv]T，上式的狀態(tài)空間表達式見式(7)：

用A來表示上式中狀態(tài)量前的系統(tǒng)矩陣，見式(8)。

根據(jù)不同任務(wù)需求，狀態(tài)偏差Δxs，i可以取不同值。例如，在進行制導(dǎo)系統(tǒng)設(shè)計時，需要得到迎角指令修正量，取為式(9)：

此時，式(7)中最后一項可看作控制輸入，在進行誤差傳播分析時，考慮質(zhì)量偏差、2種氣動系數(shù)偏差、大氣密度偏差影響，則取式(10)：

最后一項又可看作狀態(tài)攝動項。在不同的任務(wù)中，系統(tǒng)矩陣A形式相同，下面將依次分析A中各項具體表達形式。

1)引力加速度偏差。引力加速度與位矢r相關(guān)，用R0表示地心至發(fā)射系原點矢量，則有式(11)。

由于J2擾動導(dǎo)致的誤差是微小的，因此可忽略其影響，此時有式(13)。

式(11)可推導(dǎo)為式(14)：

2)科氏加速度偏差。記ωe×為ωe的反對稱陣，見式(15)：

于是科氏加速度又可以表示為式(16)。

對速度偏導(dǎo)可得式(17)。

3)離心加速度偏差。利用反對稱陣，離心加速度可寫為式(18)：

由于R＝r＋R0，則有式(19)。

4)推力加速度偏差。推力加速度矢量形式在式(2)中可見，其中彈道傾角φ與速度v相關(guān)，見式(20)。

5)氣動加速度偏差。氣動加速度與位矢r、速度v均相關(guān)，不考慮偏航運動下得式(21)：

在飛行中大氣密度模型見式(22)[12]：

則氣動加速度對位矢的偏導(dǎo)可以表示為式(23)。

注意到σ、θ、q、C均是速度v的函數(shù)，則氣動加速度對速度的偏導(dǎo)可以表示為式(24)。

至此，得到了系統(tǒng)矩陣A的完整表達形式。

2.3 逃逸軌跡設(shè)計

在工程任務(wù)中，軌跡設(shè)計之初依據(jù)經(jīng)驗一般給出如圖1所示的指令迎角，迎角極值在區(qū)間-1°到-10°內(nèi)待定。自逃逸開始到t0段，迎角保持為0；t0至t1段，設(shè)計迎角以固定斜率4°/s線性減少直到極值α；t2為飛行器逃逸主推發(fā)動機關(guān)機點時刻，設(shè)計迎角在t1至t2段保持不變；t3為主推發(fā)動機關(guān)機后5 s，在此處迎角回歸至0并一直持續(xù)到被動段飛行結(jié)束時刻t4。

圖1 指令迎角規(guī)律Fig.1 Profile of attack angle

為確定合適最大指令迎角，需要先在不考慮偏差時，遍歷搜索出末態(tài)滿足二次分離指標的迎角極值子區(qū)間。然后進行彈道魯棒性評估，在考慮偏差干擾下，通過打靶篩選出仍能滿足二次分離點指標且裕度較大的迎角極值。選定合適值后，在主動段設(shè)計制導(dǎo)律，保證在偏差作用下，能跟蹤上主動段部分標稱軌跡，這樣在被動段自由飛行后末態(tài)才有可能滿足末態(tài)二次分離要求。

利用推導(dǎo)得出的誤差線性化模型，應(yīng)用于彈道評估與制導(dǎo)兩部分，軌跡設(shè)計中聯(lián)合使用誤差傳播法與LQR控制律將體現(xiàn)出極高的計算效率。

3 軌跡魯棒性評估及制導(dǎo)

3.1 基于誤差傳播模型的魯棒性分析

誤差傳播模型用來快速計算有偏差情況下末端狀態(tài)。基本原理是根據(jù)以位置偏差、速度偏差為狀態(tài)量的系統(tǒng)狀態(tài)空間表達式，求解其狀態(tài)轉(zhuǎn)移矩陣的解析解，進而得出狀態(tài)偏差在任意時刻的通解。考慮在質(zhì)量偏差、氣動系數(shù)偏差、大氣密度偏差下，系統(tǒng)產(chǎn)生的偏差。按照式(7)，取：Δxs＝

可得攝動狀態(tài)方程式(25)。

系統(tǒng)矩陣A由式(8)得出，而攝動項V在此處為式(26)：

該系統(tǒng)為線性時變系統(tǒng)，其狀態(tài)轉(zhuǎn)移矩陣Φ是系統(tǒng)矩陣A的函數(shù)，滿足式(27)：

對于簡單的時不變系統(tǒng)，狀態(tài)轉(zhuǎn)移矩陣可通過矩陣指數(shù)函數(shù)快速得到。而對于時變系統(tǒng)，通過共軛法來求解狀態(tài)傳遞函數(shù)。

式(27)所代表的攝動系統(tǒng)中，共軛方程可寫為式(28)。

G (t，τ)為引入的共軛矩陣，可得式(29)。

于是有式(30)。

對方程組式(31)：

從tf到t0一次積分，即可得到該時間區(qū)間內(nèi)G t，tf( ) 值，由式(32)：

在初始偏差條件χ(t0)＝χ0下，根據(jù)微分方程理論，可得偏差狀態(tài)系統(tǒng)通解為式(33)：

若攝動項V中狀態(tài)偏差源Δxs不是時間的函數(shù)，記為式(34)：

系統(tǒng)末態(tài)誤差為式(35)：

設(shè)X*為系統(tǒng)在無偏差狀態(tài)下的末狀態(tài)，那么存在偏差時，系統(tǒng)末狀態(tài)為式(36)：

通過已經(jīng)計算好的系數(shù)矩陣M1與M2，在誤差源不同取值下可利用式(36)快速計算新的末狀態(tài)參數(shù)。

彈道魯棒性分析過程中，一般直接采用解析法來完成打靶，需要將偏差項直接代入原系統(tǒng)動力學(xué)模型中積分，多次計算速度受限于積分過程。使用誤差傳播法，在確定一種模型的系數(shù)矩陣后，通過矩陣乘法能直接得出末態(tài)參數(shù)，這種方式打靶效率顯然高于傳統(tǒng)的解析法。具體計算流程如圖2所示。

圖2 誤差傳播法計算流程Fig.2 Flowchart of the error propagation method

3.2 LQR制導(dǎo)模型

實際軌跡在各種誤差影響下不可避免地會偏離標稱軌跡，采用LQR方式進行制導(dǎo)，通過以調(diào)整指令迎角的方式，使得實際軌跡能夠跟蹤上標準軌跡。采用式(9)表達形式，并引入控制量u＝Δxs＝Δα。

這里迎角偏差作為控制量，于是偏差系統(tǒng)狀態(tài)方程寫為式(37)：

其中矩陣B滿足式(38)：

對于上述系統(tǒng)，求解控制量u，使得χ＝O即可實現(xiàn)對標準軌跡的精確跟蹤。

系統(tǒng)狀態(tài)χ描述了飛行軌跡位置偏差、速度偏差共計6個量。但是考慮到，飛行器無側(cè)滑運動在射面內(nèi)飛行，因此狀態(tài)χ中z方向運動量，不關(guān)心也不可控。此外，x方向位置代表了射程，在之后計算實際軌跡上各點的反饋矩陣時，會利用其在標稱軌跡上插值，使得每一時刻下Δx＝0。從系統(tǒng)狀態(tài)空間中去除x、z、vz，新的狀態(tài)量為式(39)：

相應(yīng)的狀態(tài)空間變?yōu)槭?40)：

該系統(tǒng)最優(yōu)控制性能指標函數(shù)為式(41)：

其中Q與N分別為狀態(tài)向量與控制向量的加權(quán)矩陣。在工程實踐中，Q、N陣常取對角矩陣，這里χ′為三維量，u為一維量，所以可以令：

性能指標改寫為式(42)：

根據(jù)式(43)所示Bryson法則：

這里取：Δymax＝30 m，Δvxmax＝1 m/s；Δvymax＝10 m/s，Δαmax＝5°。

相應(yīng)的，Q1＝1，Q2＝900，Q3＝9，N1＝36。

為了讓性能指標F最小，最優(yōu)控制量應(yīng)符合式(44)：

其中，時變矩陣P是式(45)所示Riccati方程的解：

在得到u*即Δα*后，進一步算得當前需要的指令迎角，見式(46)：

式中αref為標稱軌跡下指令迎角。

在標稱軌跡上，每一個制導(dǎo)周期內(nèi)取一個特征點，在各特征點附近認為系統(tǒng)矩陣A、B保持不變。計算得到各特征點下的反饋陣K后，預(yù)先裝訂成反饋增益參數(shù)。對于逃逸飛行器，在飛行全過程中，每一個制導(dǎo)周期來臨時，根據(jù)當前射程x，找到標稱軌跡上相鄰的兩個特征點，利用這兩點的反饋增益系數(shù)，插值計算得到當前射程下對應(yīng)的反饋增益系數(shù)作為該制導(dǎo)周期內(nèi)通用的反饋增益系數(shù)。

4 算例及仿真

相關(guān)仿真參數(shù)由表1給出，二次分離點指標見表2，參數(shù)偏差取值見表3。

表3 參數(shù)偏差Table 3 Parameter deviations

氣動系數(shù)擬合函數(shù)為式(47)：

第一步，在不考慮偏差下遍歷搜索-1°到-10°區(qū)間，找到滿足末態(tài)指標的子區(qū)間。該問題中，最大指令迎角在區(qū)間 [-6°，-4°]時，能滿足二次分離指標。 -4°、-5°、-6°情況下的末態(tài)參數(shù)在表 4 中可見，4種主要參數(shù)滿足了表2要求的約束條件。

表4 3種指令迎角下仿真結(jié)果Table 4 Simulation results under 3 angles of attack

第二步，分別計算以上3種情況在被動段的誤差傳遞系數(shù)。被動段初值使用各迎角在不考慮偏差下主動段飛行結(jié)束時刻值。以指令迎角為-4°時為例，被動段初始狀態(tài)為：[r0v0]T＝[132.06 1605.93-0.25 39.00 182.97-0.047]T，相應(yīng)誤差傳播系數(shù)為：

第三步，利用誤差傳遞系數(shù)，由式(36)快速計算末態(tài)參數(shù)。考慮到主動段在進行制導(dǎo)后，被動段的初始參數(shù)相比標準狀態(tài)不可避免地會有差異。χ0取被動段標準狀態(tài)下初值10%范圍內(nèi)隨機大小的偏差。Δxs取參數(shù)偏差區(qū)間內(nèi)隨機大小的偏差。進行N＝10 000次打靶，繪制各指令迎角下的射程射高散布、橫向速度縱向速度散布，見圖3～圖5。

此外，為了驗證這種算法正確性，這里在固定的極限初態(tài)偏差χ0和參數(shù)偏差Δxs下，對比分析2種算法下仿真結(jié)果：χ0＝[Δr0Δv0]T＝[-13.21 160.59-0.03-3.90 18.30 0.00]T；Δxs＝[ΔmΔfxΔgyΔρ]T＝[-1000 0.2 4×10-50.05]T。

在-4°最大迎角下，分別使用解析法和誤差傳播法計算系統(tǒng)末態(tài)參數(shù)如下：X1tf＝[370.0270 2618.6789-1.1965 16.9327-1.8686-0.1214]T；[365.5944 2633.1293-1.1561 16.3831-0.9445-0.1183]T。

圖3 -4°指令迎角下散布規(guī)律Fig.3 Dispersion characteristic under-4°angle of attack

圖4 -5°指令迎角下散布規(guī)律Fig.4 Dispersion characteristic under-5°angle of attack

圖5 -6°指令迎角下散布規(guī)律Fig.5 Dispersion characteristic under-6°angle of attack

無偏差情況下，系統(tǒng)的末態(tài)參數(shù)為：X0tf＝[447.1818 2475.8783-1.2149 22.1663 2.2454-0.1256]T

誤差傳播法相對解析法的計算誤差與計算時間已記錄在表5中。

表5 算法結(jié)果比較Table 5 Comparison of algorithms

解析法下計算結(jié)果可以視為標準值，采用誤差傳播得到的結(jié)果基本與之接近。從表5中數(shù)據(jù)可以看出，x、z方向上各項相對誤差均在4%以內(nèi)；y方向速度相對誤差較大，但絕對誤差在可接受范圍內(nèi)。這種情況可能來自于模型的缺陷，誤差傳播法模型基于誤差線性化模型得到，在基準量小時，引入一定誤差源后，若狀態(tài)偏差大，則線性化效果不理想；相反，基準量越大，線性化效果越好。因此，對比無偏差下末態(tài)參數(shù)可以發(fā)現(xiàn)，該系統(tǒng)中y方向位置相對誤差最小，只有0.5%；z方向位置、速度基準值較小，但在該模型中z方向運動可以忽略，誤差源并不會使其狀態(tài)偏差過大，故相對誤差也較小，約為3%；y方向速度基準值小，受誤差源影響，故相對誤差較大。在用于多次打靶時，采用誤差傳播法計算速度明顯快于解析法。以上結(jié)果，驗證了誤差傳播法的準確性與快速性。

第四步，在存在偏差時仍要滿足二次分離指標，打靶后各散點應(yīng)該分布在散布規(guī)律圖中虛線分割后的右上側(cè)區(qū)域內(nèi)。進行多輪仿真試驗，綜合各迎角下散布情況，選擇滿足末態(tài)分離要求概率更高的一種迎角作為標稱軌跡的指令最大迎角。以3輪N＝10 000次的快速打靶為例，在隨機大小參數(shù)偏差影響下，3種指令迎角工況中，能夠滿足末態(tài)約束要求的概率如表6所示。可以看出，在偏差源影響下，選擇-4°作為標稱軌跡的指令最大迎角時存在55%的成功率完成末態(tài)分離指標，相比其他角度下軌跡，抗干擾能力更強，魯棒性更高。因此，通過這種方式就可以篩選出更理想的標稱軌跡。

表6 多次打靶成功率統(tǒng)計Table 6 Success rate of multiple shooting

第五步，在存在偏差情況下，使用LQR制導(dǎo)律在主動段跟蹤標稱軌跡，被動段自由飛行。這里取如下偏差進行仿真： [ΔmΔfΔgΔρ]Txy＝[800 0.08 4×10-50.01]T

圖6 LQR制導(dǎo)前后對比Fig.6 Comparison before and after LQR guidance

LQR制導(dǎo)前后對比如圖6所示，仿真結(jié)果表明，在存在參數(shù)偏差情況下，實際的射高、橫向速度、縱向速度均會偏離標稱值。使用LQR制導(dǎo)以調(diào)整指令迎角的方式，能夠?qū)崿F(xiàn)對標稱彈道跟蹤逼近。圖6中射高、橫向速度的跟蹤效果良好，偏差在1 m、1 m/s以內(nèi)；縱向速度跟蹤效果稍差，偏差在5 m/s以內(nèi)。

全段軌跡如圖7所示，末態(tài)參數(shù)見表7。在20 s的飛行時間內(nèi)，主動段采用LQR制導(dǎo)律后基本實現(xiàn)了對射高、橫向速度、縱向速度的多變量跟蹤，被動段結(jié)束時射程、射高等末態(tài)參數(shù)也滿足了二次分離點指標要求。

圖7 飛行軌跡Fig.7 Flight trajectory

表7 全段飛行仿真結(jié)果Table 7 Simulation result of the full flight

5 結(jié)論

本文針對逃逸飛行器展開了在分離點約束條件下的軌跡設(shè)計與制導(dǎo)方法研究。區(qū)別于一般彈箭，逃逸飛行器只能在主動段短時間內(nèi)實施制導(dǎo)指令，這就對標稱軌跡在偏差作用下的魯棒性提出了一定要求。主要結(jié)論如下：

1)基于誤差線性化模型，推導(dǎo)建立了誤差傳播模型和LQR制導(dǎo)模型。相較于傳統(tǒng)的解析法，采用誤差傳播法，在快速性上有著巨大的優(yōu)勢，通過具體算例的仿真驗證了這一算法的準確性；

2)基于LQR設(shè)計了逃逸火箭的制導(dǎo)模型，實現(xiàn)了多變量跟蹤，通過仿真驗證了方法的可行性。

3)依據(jù)以上2種模型，詳細介紹了一種可行逃逸飛行器軌跡的設(shè)計方法，在考慮偏差情況下，仍盡可能滿足末態(tài)二次分離點指標，具一定有魯棒性。