基于迭代制導的運載火箭入軌級彈道設計及仿真

高家智，陳小前

（1.國防科技大學空天科學學院，長沙，410073；2.太原衛星發射中心，太原，030027）

0 引 言

迭代制導制導精度高、任務適應性強、箭上飛行軟件簡單、對地面諸元準備要求相對較低，在運載火箭上應用越來越廣泛[1～3]。文獻[4]針對主發動機推力為常值的情況，通過優化開關機時間來滿足終端位置約束，引入權重因子提高迭代求解的精度；文獻[5]針對多約束入軌條件，提出一種將迭代制導與數值積分相結合的軌跡預測制導方法，增加的推進劑的消耗在可接受范圍之內；文獻[6]在地心慣性系中建立航天器模型，以目標軌道要素為終端約束條件，得到一種簡單有效的控件變軌迭代制導算法；針對飛行過程中發動機推力下降的重大故障，文獻[7]基于迭代制導算法對彈道進行了重構，在完成故障識別和能力預測后，利用迭代制導將有效載荷送入預先設計好的救援軌道上去。

Boeing 公司針對阿波羅計劃使用的土星5 運載火箭迭代制導方法的不足，開發了線性正切制導律（Linear Tangent Guidance，LTG），并在飛馬火箭上升段入軌控制和離軌段返回控制中得到應用[8]。本文方法利用了其閉路制導律和預估公式，但制導計算的思路與文獻提供的相反，引入了制導迭代變量，采用了根據初值求解兩點邊值問題的Newton-Raphson 方法。最后利用數值仿真驗證了本文所設計方法的有效性。

1 最優控制問題及其制導律

在不考慮稀薄大氣的情況下，火箭在真空段內的運動方程可表示如下[9]：

終端約束條件為

式中V，R分別為火箭在地心發射慣性坐標系中的速度和位置；G為重力矢量，G=-ω2R=-3R，其中，fM= 3.986 × 105km3s2為引力常數；ω為地球自轉角速度；為秒流量；a為加速度的模，a=，其中，P為發動機推力，m為質量；Ue為有效排氣速度；u為推力方向的單位向量。

定義哈密頓函數H和增廣函數F如下：

這里的λ，α和β是拉格朗日乘子。

定義如下變量：

式中φ，ψ分別為俯仰角和偏航角。根據最優控制的必要條件有如下歐拉方程和橫截條件：

由歐拉方程可解出：

式中uR=unit(R)為火箭質心地心矢徑的單位向量。

由橫截條件可得出：

相對于地球半徑而言，火箭因推進引起高度的變化量是小量，并且λ與uR的夾角在90°附近，因此在式（8）中可忽略高階向量，得：

該方程的解為

式中t為飛行時間；λK為定常單位向量；為定常速率向量；K為定常時間，一般情況下K≈Tgo2，Tgo為剩余飛行時間。如果是常值加速度，則K=Tgo2。

可以證明，對于圓軌道入軌時，下式成立：

即兩矢量正交。

由于在迭代過程中，不斷刷新K、λK和，為簡便起見，省去下標K，即有：

考察式（13），存在如下特性：

b）如果采用地球平面假設，因為ω=0，則：

在本文的制導研究中，仍考慮地球是球體，且使用上述有關假設：

b）ω=λ˙；

c）忽略式（8）中的 3ω2(λ·uR)uR；

那么，制導律為

簡化為LTG 制導律：

2 迭代制導的實現

2.1 端點狀態預估

運載火箭的迭代制導屬于復雜的兩點邊值問題，關機時刻火箭位置和速度的預估精度關系到迭代制導的精度、穩定性和魯棒性。數值積分是一種精確的解決辦法，但卻是無法在箭上實時實現的辦法，所以推導簡化的解析代數公式就是一項根本任務。

在制導律式（15）的控制作用下，火箭在真空段的運動方程為

假設剩余飛行時間Tgo是已知的，且存在如下封閉形式的積分：

式中F1，F2，F3為待定系數。

積分式（18）有：

式中VD，RD分別為末端速度和位置矢量。

改寫式（19）、（20），得到：

再積分式（17），得到：

式中Vgo，Rgo分別為待增速度和位置矢量。

因此，如果已經知道λ和，就可以根據積分公式LT、ST、QT得出Vgo和Rgo，從而在當前位置R和速度V下，由制導律式（15）進行制導的端點的預估公式為

式中VP，RP分別為預測的末端速度和位置矢量。

2.2 縱平面迭代制導

運載火箭的歐拉角采用3-2-1 型表示的發慣系動力學方程組，可知只有側平面的運動影響縱平面的運動，但縱平面的運動不會影響側平面的運動，特別當偏航角很小時飛行控制可以解耦。因此可以先進行縱平面的迭代控制，再進行側平面的迭代控制。

縱平面內的迭代控制采用標準入軌參數作為邊界條件。這里的邊界條件是采用極坐標邊界條件而非直角指標邊界條件，因為半長軸和偏心率與極坐標直接相關，所以在縱平面內的運動控制邊界條件是滿足標準入軌點對應的地心矢徑、速度大小和當地速度傾角，這樣就保證了橢圓軌道的半長軸和偏心率，加上側平面的標準軌道跟蹤導引，控制了zK和，理論上也是控制軌道根數的5 個變量。升交點赤經Ω主要由發射窗口的起飛時刻、火箭飛行時間和方位角決定，近地點幅角雖然是放開的，但跟蹤Z通道的標準軌道這種約束本身也對第6 個軌道根數進行了隱式控制，這種方法一定程度約束了近地點幅角偏差的放大。

這里采用極坐標邊界條件的另一個目的是提高迭代制導的魯棒穩定性和收斂性，減小迭代制導結尾階段的過早發散效應，三自由度制導仿真計算表明選擇這種邊界條件與視加速度反饋措施一起可以有效提高迭代制導的精度，使跳出迭代的剩余常值飛行時間小于0.5～1.0 s，幾乎迭代至最后，提高了制導精度。

縱平面內的迭代控制采用 LTG 制導律方法，預測縱平面內關機時刻發慣系內的直角坐標參數Xp、、Yp和，然后進一步預測關機時刻的位置矢徑Rp、速度大小Vp和當地速度傾角Θp，即：

迭代公式采用牛頓迭代法，即：

式中α,β為與俯仰角有關的制導變量。

因為采用了極坐標邊界條件，沒有顯式解，校正計算需要采用3 變量的牛頓快速迭代法，涉及3×3 雅可比逆矩陣的實時計算，3×3 雅可比逆矩陣不需要數值計算，可以有解析表達式，并且因為上一步計算得到了的初值，下一步牛頓迭代計算只需一次就很精確，完全滿足制導收斂精度。

由迭代解α,β,Tgo可以計算LTG 制導律的乘子向量λ和得到：

發慣系縱平面推力方向的單位控制向量：

記u=[u(1),u(2)]T，故當前時刻滿足邊界條件的俯仰姿態角應為

2.3 側平面迭代制導

發射慣性系中的姿態角為俯仰角和偏航角為φ和ψ，初始質量m0，秒流量，軸向推力F在發射慣性系中推力分量表示為

動力學方程為

若偏航角為常值規律，火箭推力通常假設是常值，則推力的Z通道分量為常值，那么側向通道在關機點的狀態可以解析預測，即：

利用縱平面預測得到的縱向關機點狀態和關機時間以及上述側向通道關機點狀態公式預測，可以進行軌道傾角i的解析預測計算，方法為：在縱向通道迭代計算完畢，預估縱向關機速度和位置，再根據上述常值偏航角和常值推力計算關機時的發慣系側向速度和側向位置，然后轉換到發射地慣系計算軌道傾角，則i=i(ψ)，根據軌道傾角期望值ic，容易單變量迭代求解偏航角指令。

那么對偏航角進行迭代滿足所要求的軌道傾角屬于單變量迭代，很容易快速求解。

3 仿真結果及分析

仿真起始條件為：發慣系下起滑點及第3 級關機點參數X=593 125.71 m，Y=184 422.32 m，Z=-77 388 m，VX=6676.607 m/s，Vy=648.986 m/s，Vz=-286.533 m/s，標準點火時刻532.3 s，推力28 000 N，初始質量1206.3 kg，秒耗量10 kg/s，標準滑行時間303.7 s[10]。

目標衛星軌道參數為：半長軸a=6860.788 km，偏心率e=0，軌道傾角i=97.3 °。

制導精度要求：半長軸偏差小于1 km，偏心率偏差小于0.0001，軌道傾角偏差小于0.01°。

經過仿真計算，發慣系位移及速度變化情況見圖1所示。

圖1 發慣系位移及速度變化情況Fig.1 Position and Velocity Changes in Launch Inertial System

俯仰角、偏航角以及其它參數變化曲線如圖2 所示，滑行軌道、迭代制導軌道與運行軌道的放大圖如圖3 所示。輸出結果如下：半長軸6860.768 km，偏心率0.000 012，軌道傾角97.3072°，制導時間49.35 s，升交點角距163.42°。可見長半軸誤差20 m,偏心率偏差0.000 012，軌道傾角偏差0.0072°，均滿足制導精度要求。

圖2 主要飛行變量變化情況Fig.2 Major Flight Variable Variations

圖3 滑行軌道、迭代制導軌道及運行軌道放大圖Fig.3 Enlarged View of Dragging Track,Iterated Guidance Track,and Orbit

4 結束語

本文對多級運載火箭入軌段的迭代制導方法進行了較為深入研究，從仿真結果來看本文所提出的入軌級迭代制導方法制導精度高，且計算簡單、易于實現。