逃逸飛行器軌跡魯棒性快速評估及制導

李 臻，楊 彪，荊武興，高長生，常武權

(1.哈爾濱工業大學航天工程系，哈爾濱150001；2.北京宇航系統工程研究所，北京100076)

1 引言

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

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

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

2 逃逸飛行器模型及軌跡設計

2.1 動力學模型

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

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

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

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

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

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

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

2.2 誤差線性化模型

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

Δxs，i為第i個狀態偏差，aij是與該狀態相關的第j個加速度項。

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

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

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

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

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

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

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

式(11)可推導為式(14)：

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

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

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

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

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

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

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

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

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

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

至此，得到了系統矩陣A的完整表達形式。

2.3 逃逸軌跡設計

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

圖1 指令迎角規律Fig.1 Profile of attack angle

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

利用推導得出的誤差線性化模型，應用于彈道評估與制導兩部分，軌跡設計中聯合使用誤差傳播法與LQR控制律將體現出極高的計算效率。

3 軌跡魯棒性評估及制導

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

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

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

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

該系統為線性時變系統，其狀態轉移矩陣Φ是系統矩陣A的函數，滿足式(27)：

對于簡單的時不變系統，狀態轉移矩陣可通過矩陣指數函數快速得到。而對于時變系統，通過共軛法來求解狀態傳遞函數。

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

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

于是有式(30)。

對方程組式(31)：

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

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

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

系統末態誤差為式(35)：

設X*為系統在無偏差狀態下的末狀態，那么存在偏差時，系統末狀態為式(36)：

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

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

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

3.2 LQR制導模型

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

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

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

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

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

相應的狀態空間變為式(40)：

該系統最優控制性能指標函數為式(41)：

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

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

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

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

相應的，Q1＝1，Q2＝900，Q3＝9，N1＝36。

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

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

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

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

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

4 算例及仿真

相關仿真參數由表1給出，二次分離點指標見表2，參數偏差取值見表3。

表3 參數偏差Table 3 Parameter deviations

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

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

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

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

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

此外，為了驗證這種算法正確性，這里在固定的極限初態偏差χ0和參數偏差Δxs下，對比分析2種算法下仿真結果：χ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°最大迎角下，分別使用解析法和誤差傳播法計算系統末態參數如下：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°指令迎角下散布規律Fig.3 Dispersion characteristic under-4°angle of attack

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

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

無偏差情況下，系統的末態參數為：X0tf＝[447.1818 2475.8783-1.2149 22.1663 2.2454-0.1256]T

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

表5 算法結果比較Table 5 Comparison of algorithms

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

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

表6 多次打靶成功率統計Table 6 Success rate of multiple shooting

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

圖6 LQR制導前后對比Fig.6 Comparison before and after LQR guidance

LQR制導前后對比如圖6所示，仿真結果表明，在存在參數偏差情況下，實際的射高、橫向速度、縱向速度均會偏離標稱值。使用LQR制導以調整指令迎角的方式，能夠實現對標稱彈道跟蹤逼近。圖6中射高、橫向速度的跟蹤效果良好，偏差在1 m、1 m/s以內；縱向速度跟蹤效果稍差，偏差在5 m/s以內。

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

圖7 飛行軌跡Fig.7 Flight trajectory

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

5 結論

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

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

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

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