航天器遠程最優攔截方法研究①

符 俊，蔡 洪，張士峰，丁洪波

(國防科技大學航天與材料工程學院，長沙 410073)

0 引言

迄今有很多學者對空間攔截問題進行了研究，一般來說，它們可分為2類:脈沖推力，如Herrick方法、Godal方法、Gauss方法、普適變量法和 Lambert方法等［1］;有限推力，如文獻［2－4］。對于有限推力攔截情況，往往是在一定的性能指標下尋找最優的攔截軌道，這是一個最優控制問題。

隨著計算機水平的高速發展，求解最優控制問題的數值方法得到了廣泛應用，其中魯棒性和通用性更好的直接法往往能更好地求解動態優化問題。它將函數空間中的最優控制問題轉化為歐式空間中的非線性規劃問題，然后再利用各種非線性算法完成問題的求解［5］。

本文基于非線性化C-W方程，利用Legendre偽譜法(Legendre Pseudospectral Method，LPM)對航天器遠程攔截軌道進行優化，同時基于固定時間飛行定理討論了一種脈沖推力假設下的優化方案，并將這2種方法進行了對比分析。

1 數學模型

1.1 動力學方程

C-W方程常用于描述圓參考軌道下的航天器近距離相對運動，這主要是因為C-W方程在建立的過程中作了以下2個假設:(1)目標航天器的軌道為圓軌道;(2)攔截器與目標航天器之間的距離遠小于它們的地心距。但在攔截問題中，假設(2)極大地限制了C-W方程的應用，因為攔截器和目標航天器之間的距離往往很大，故假設2不再成立。在這種情況下，下面給出適于描述遠程攔截問題的航天器相對運動方程(目標航天器不機動):

式中 ［x y z vxvyvz］T表示攔截器在目標航天器LVLH(Local Vertical，Local Horizontal)坐標系中的相對運動狀態矢量;Tx、Ty、Tz為攔截器的推力在LVLH坐標系三軸上的分量;m為攔截器的質量;Isp為攔截器推進劑比沖;a為目標航天器的地心距;g為重力加速度;μ為地球引力常數;rc為攔截器的地心距，為目標航天器的平均軌道角速度為攔截器推力

在優化問題中，控制量為 U=［Tx，Ty，Tz］T，性能指標應為攔截器消耗的燃料最少，即終端時刻質量最大:

式中 tf表示終端時刻。

1.2 無量綱化參數

采用數值方法進行優化時，必須對物理量進行無量綱化，選取的無量綱化參數如下:

式中 Rs為目標航天器的地心距;m0為攔截器的初始質量。

2 有限推力最優攔截

2.1 Legendre偽譜法簡介

Legendre偽譜法由美國海軍研究院Fahroo Fariba和Ross I Micheal提出［6］，它屬于直接法中的一種。該方法將狀態變量和控制變量在一系列LGL(Legendre-Gauss-Lobatto)點上離散，并以離散點為節點構造Lagrange插值多項式來逼近狀態變量和控制變量，從而將動態最優控制問題轉化為靜態參數優化問題(NLP問題)。Legendre偽譜法通過對全局插值多項式求導來近似狀態變量對時間的導數，從而將微分方程約束轉換為一組代數約束，這些約束條件加上問題本身的約束條件，如邊界約束、路徑約束等，共同構成NLP問題的約束條件。Legendre偽譜法的求解步驟可參考文獻［3 －4，7］。

2.2 優化結果及分析

以一具體的軌道攔截任務為例，研究Legendre偽譜法的應用。仿真條件:攔截器總質量5 00 kg，推進劑質量4 00 kg，比沖300 s，最大推力2 000 N。初始時刻 t0，攔截器的軌道根數:a=16 678.137 km，e=0.2，i=28.5°，ω =0，Ω =0，?=1.54°。目標航天器的軌道根數:a=6 878.137 km，e=0，i=30°，Ω =0，u=2.05°。

經過計算，可得到t0時刻攔截器在目標航天器LVLH坐標系中的相對運動狀態:

其中，前3項表示相對位置，后3項表示相對速度，距離單位為km，速度單位為km/s。

仿真計算中采用60個LGL點，優化結果見圖1～圖4。

圖1 控制變量的優化結果Fig.1 Optimal results of control variables

圖1(a)為推力隨時間的變化曲線，可看出在攔截過程中，攔截器受到的推力主要沿著跡向，在另外2個方向上的分量很小。發動機工作時產生的推力一般沿體系x方向，為了使推力在LVLH坐標系三軸上的分量大小滿足要求，可通過控制攔截器的姿態來實現。在本次仿真中，推力作用的時間為506.9 s，最大值為2 000 N，符合最大推力約束條件。之后的過程發動機停止工作，攔截器在地心引力的作用下飛向目標航天器進行攔截。圖1(b)為攔截器的過載變化曲線圖，從圖1(b)可看出，最大過載出現在y方向，但均不超過1，滿足容許過載較小的航天器的變軌要求。

顏曉晨正在試衣服，一個二十五六歲的長發女子走了進來，看了她幾眼，拿了一套顏曉晨試穿的衣服，翻看價格牌。一個營業員在接電話，另一個營業員正低著頭幫顏曉晨整理褲腳，都沒顧上招呼她，顏曉晨笑著說：“全場五折。”

圖2 狀態變量優化結果Fig.2 Optimal results of state variables

圖2為攔截器的狀態變量優化結果圖。其中圖2(a)為相對距離變化曲線，可看到，經歷8 080.5 s后，攔截器與目標航天器的相對距離趨于0，成功實現攔截。分析圖2(a)可發現，攔截器與目標航天器的相對距離經歷了一個先增大后減小的過程，而不是單調減小，這是因為性能指標為能量最省，在前面給定的初始相對運動狀態條件下，如果直接控制攔截器朝目標航天器飛去，勢必會消耗很大的能量，故攔截器與目標航天器之間的距離不是單調減小。圖2(b)為相對速度變化曲線，圖2(c)為攔截器的質量變化圖。攔截器最終質量為 233.5 kg，燃料消耗為 266.5 kg。

圖3 相對運動軌跡Fig.3 Trajectories of relative motion

圖4 哈密頓函數隨時間的變化曲線Fig.4 Time histories of Hamiltonian

圖3為攔截器相對目標航天器的運動軌跡，從圖3看到，它們之間的相對運動主要發生在x-y平面，是由于初始時刻兩者的軌道面空間位置差異不大造成的。

在應用Legendre偽譜法處理最優控制問題時，往往面臨著一個問題:解的最優性能否得到保證?這可通過求解哈密頓函數的值是否滿足一階最優性條件來進行判斷:作為區別于一般直接數值解法的重要特性之一，Legendre偽譜法將最優控制問題轉化為LGL點上的非線性規劃問題后，可將對該問題的求解轉化為對一個增廣性能指標的優化求解，并能在求解該問題的同時得到Lagrange乘子，即LGL點上的協態變量值，從而計算哈密頓函數值以檢驗是否滿足一階最優性必要條件。如果滿足則可認為結果最優，否則就不是最優的［8－9］。

本次仿真中的性能指標為J=－m(tf)，這是一個末值型的性能指標，且終端時刻自由。根據龐特李亞金極大值原理，在優化過程中，哈密頓函數值應恒為0。圖4為哈密頓函數隨時間的變化曲線，從圖4可看到，該值一直處于0附近，最大相差為0.7%，滿足最優性必要條件，因此解的最優性能夠得到保證。

3 基于脈沖推力的最優攔截方案

前面利用Legendre偽譜法得到了最優攔截問題的解。下面討論在脈沖推力假設下航天器遠程最優攔截方案的研究方法，理論基礎是固定時間攔截定理。

3.1 Battin-Vaughn 方法

求解固定時間攔截問題的經典算法是高斯方法和Lambert飛行時間定理。通過對上述2種算法的改進，派生出了很多算法，如海里克方法、歌德方法、普適變量法、p迭代法和Battin-Vaughan方法。其中，普適變量法算法簡單，具有普適性，適用于所有的圓錐曲線，但在某些情況下，收斂速度太慢;p迭代法用牛頓迭代法來修正p的試探值，收斂速度很快，其不足在于r1、r2共線情況下不收斂;而Battin-Vaughan方法通過引入超幾何函數和連分數進行數值迭代求解，迭代速度與精度均比較理想，適合計算機快速求解對于任意情況都收斂很快，具有幾乎一致收斂的性質。據Klumpp在1986年以來對Battin-Vaughan算法就所有合理情況所做的廣泛實驗研究表明，Battin的算法對于載人和不載人自主制導都能提供所需要的可靠性、快速性和列緊性，并且對于行星軌道確定以及整個天體力學也都提供了所需要的精度和應用范圍［10］。

3.2 優化方案

在空間攔截任務中，關鍵技術是在一個時間窗口(也稱任務窗口)［t0，t1］中攔截器能命中目標并使消耗的能量最小。根據固定時間攔截定理，單圈Lambert問題存在唯一解。因此，攔截器不同的機動時刻tman對應著攔截目標航天器所需的不同速度沖量Δv，這稱為一種攔截方案，以(tman，Δv)表征。根據這種思路，在實際計算過程中，首先計算出一次攔截任務中任務窗口內的所有機動攔截方案，然后在這些方案中選擇一種最優的方案(比如所需特征速度最小或攔截時間最短)或者是滿足特定攔截任務要求(在指定時刻攔截目標)的方案，這種靈活性是該方法的一大優點。

尋求最優攔截方案的關鍵是求解機動時刻tman，可采用遺傳算法進行求解，也可采用枚舉法，這里選擇枚舉法。采用枚舉法尋找最優攔截方案，具體要求是在時間窗口內，選擇機動時刻tman作為迭代變量，求解特征速度最小的機動方案，迭代步長可視具體情況靈活控制。為克服枚舉法計算量大的缺點，可通過變步長控制計算量。求解流程見圖5。

圖5 仿真流程Fig.5 Flow of simulation

3.3 仿真分析

仿真條件與2.2節相同，為與偽譜法的優化結果進行對比，選擇攔截器的任務時間窗口為［t0，t1］，其中t1=8 080.5 s。

圖6為機動時刻與所需速度沖量的關系，橫坐標表示機動時刻tman，縱坐標表示所需的速度增量Δv。t0為零時刻，可看出，當機動時刻tman=4 328 s時，所需的速度增量最小，為1 501.73m/s。

圖6 機動時刻與速度沖量的關系Fig.6 Relationship of maneuver time and velocity impulse

攔截器最優變軌過程如表1所示，變軌所需的速度沖量為［1 091.58，－930.36，－445.03］(J2000 坐標系)，單位為 m/s;變軌點的軌道位置為 ?=241.91°，攔截點的軌道位置為?=86.25°。攔截器在原軌道上運行4 328 s后開始變軌。可以看出，變軌前后攔截器的軌道面并沒有改變，只是軌道的形狀和大小發生了變化，因此所需的速度增量不是很大。根據公式Δv=Ispg ln(m+/m－)計算得到消耗的燃料為199.99 kg。

表1 攔截器變軌前后軌道根數Table 1 Orbital elements of interceptor

4 2種方法的比較

為便于描述，Legendre偽譜法用方法Ⅰ表示，基于脈沖推力的優化方法用方法Ⅱ表示。根據前面的分析，發現方法Ⅱ的燃料消耗要小于方法Ⅰ，其中一個原因為方法Ⅱ通過選擇攔截器與目標航天器相對相位關系最優時進行變軌，如圖7所示。方法I是直接在初始相對相位關系下實施攔截機動變軌，因此使得攔截消耗的能量較大。2種方法的綜合對比詳情見表2。

圖7 攔截器變軌圖Fig.7 Orbit change of interceptor

表2 2種方法的比較Table 2 Comparison of two kinds of methods

5 結論

(1)仿真過程中發現，Legendre偽譜法精度高，對變量初值不敏感，收斂半徑大，適于處理航天器遠程攔截問題。

(2)在有限推力攔截中，Legendre偽譜法用來優化推力的大小和方向，從而得到最優攔截軌道;而在脈沖推力最優攔截方法中，主要是通過調整攔截器和目標航天器的相位關系來搜索最優攔截軌道。

(3)與Legendre偽譜法相比，采用基于脈沖的最優機動攔截方案具有更大的靈活性，既能夠滿足性能指標要求，也能滿足具體攔截任務要求。

［1］任萱.人造地球衛星軌道力學［M］.長沙:國防科技大學出版社.1988:141-157.

［2］涂良輝，袁建平，羅建軍.基于偽光譜方法的有限推力軌道轉移優化設計［J］.宇航學報，2009，29(4):1189-1193.

［3］胡正東，丁洪波，曹淵，等.偽譜法在SGKW軌道快速優化中的應用［J］.航天控制，2009，27(4):3-7.

［4］宣穎，張為華，張育林.基于Legendre偽譜法的固體運載火箭軌跡優化研究［J］.固體火箭技術，2008，31(5):425-429.

［5］羅亞中，唐國金，田蕾.基于模擬退火算法的最優控制問題全局優化［J］.南京理工大學學報，2005，29(2):144-148.

［6］Fahroo F，Ross I M.Costate estimation by a legendre pseudospectral method［J］.Journal of Guidance，Control and Dynamics，2002，24(2):270-277.

［7］Huntington Geoffrey Todd.Advancement and analysis of a gauss pseudospectral transcription for optimal control problems［D］.Cambridge，Massachusetts Institute of Technology，2007:115-143.

［8］Qi G，Ross I M，Wei K，et al.Connections between the convector mapping theorem and convergence of pseudospectral methods for optimal control［J］.Compute Optimal Application，DOI.10.1007/s10589-007-910-2-4，2007(1):1-29.

［9］丁洪波，蔡洪，張士峰，等.高超聲速滑翔式再入飛行器最大航程飛行軌跡分析［J］.國防科技大學學報，2009，31(6):67-72.

［10］趙瑞安.空間武器軌道設計［M］.北京:中國宇航出版社，2008:176-178.

［11］Richard H Battin.An introduction to the mathematics and methods of astrodynamics［M］.AIAA Education Series，AIAA，Reston，VA，1999.