張洪禮，羅欽欽，韓潮*

(1.北京航空航天大學(xué) 宇航學(xué)院，北京100191;2.北京空天技術(shù)研究所，北京100074)

在航天器軌道動(dòng)力學(xué)中，最著名的兩點(diǎn)邊值問題是二體Lambert問題，即在一中心引力場(chǎng)中，求解經(jīng)過兩位置矢量R1和R2、轉(zhuǎn)移時(shí)間為Δt12的開普勒軌道.到目前為止，已有大量學(xué)者對(duì)二體Lambert問題進(jìn)行了深入廣泛的研究，文獻(xiàn)［1－3］等都給出了經(jīng)典的求解方法.二體Lambert問題的特點(diǎn)是中心引力場(chǎng)為主要影響因素，其余影響可作為攝動(dòng)考慮，主要應(yīng)用于近地軌道設(shè)計(jì).但是，在深空軌道設(shè)計(jì)中，需要考慮第三體引力影響，不能單純將其作為攝動(dòng)處理，否則會(huì)導(dǎo)致很大誤差，這就引入了三體Lambert問題.與二體Lambert問題相比，三體Lambert問題的邊值條件是相同的，但由于第三體引力不可作為攝動(dòng)考慮，增加了動(dòng)力學(xué)模型的非線性程度，使得問題求解更加困難.由于三體問題不存在解析解，三體Lambert問題必須依靠數(shù)值方法求解，其本質(zhì)上是一個(gè)兩點(diǎn)邊值問題，一般的求解過程分為兩步:首先根據(jù)邊值約束猜測(cè)初值，然后通過數(shù)值預(yù)報(bào)對(duì)軌道進(jìn)行精確預(yù)報(bào)，采用微分修正等數(shù)值方法對(duì)初值進(jìn)行修正，最終得到收斂的精確解.

本文選取存在一次相對(duì)次天體引力輔助變軌的三體Lambert問題，作為一般三體Lambert問題的一類特例，此類問題可作為多次引力輔助轉(zhuǎn)移軌道設(shè)計(jì)［4－5］的基礎(chǔ)模塊，可用于地月間［6－7］或地火間［8］軌道設(shè)計(jì)中，還可用于平動(dòng)點(diǎn)軌道設(shè)計(jì)［9－10］中，有較大的研究?jī)r(jià)值.針對(duì)三體 Lambert問題的求解，文獻(xiàn)［11］提出了一種簡(jiǎn)單迭代算法，但文中沒有給出搜索變量初值的估計(jì)方法.文獻(xiàn)［12］利用偽狀態(tài)理論，將相對(duì)于次天體的近拱點(diǎn)以及相應(yīng)的位置速度作為設(shè)計(jì)變量，構(gòu)建一個(gè)7維的微分修正，迭代搜索飛掠次天體的轉(zhuǎn)移軌道，但文中同樣沒有給出猜測(cè)初值的方法.文獻(xiàn)［13］ 將轉(zhuǎn)移軌道按其形態(tài)上的特點(diǎn)進(jìn)行歸納分類，建立初值樣本庫，然后根據(jù)給定兩點(diǎn)的相對(duì)位置關(guān)系，從樣本庫中選擇合適的初值，進(jìn)而進(jìn)行精確軌道搜索，但這種方法需要設(shè)計(jì)者具有十分豐富的深空軌道設(shè)計(jì)經(jīng)驗(yàn)，對(duì)特定的約束選擇合適的設(shè)計(jì)初值，人工干預(yù)較大.文獻(xiàn)［14］在微分修正算法的基礎(chǔ)上，推導(dǎo)并利用二階微分修正算法，通過偽狀態(tài)理論給出的初值，可獲得收斂的精確解，但二階微分修正算法的收斂性能仍然有限，尤其對(duì)于轉(zhuǎn)移時(shí)間較長(zhǎng)、飛掠高度較低的情況，轉(zhuǎn)移軌道對(duì)于近拱點(diǎn)狀態(tài)是十分敏感的，可能會(huì)導(dǎo)致搜索過程不斷振蕩，甚至發(fā)散.

針對(duì)以上研究現(xiàn)狀，本文提出一種基于無損卡爾曼濾波(UKF)參數(shù)估計(jì)算法的三體Lambert問題求解方法.首先，只需基于二體Lambert問題在慣性空間中求解出初始設(shè)計(jì)結(jié)果，然后將原問題的求解轉(zhuǎn)化為參數(shù)估計(jì)問題，利用UKF參數(shù)估計(jì)算法求解真實(shí)攝動(dòng)環(huán)境下的精確轉(zhuǎn)移軌道，就可以得到收斂的精確解.UKF濾波算法是基于概率估計(jì)理論的，不依賴于梯度信息，在初值精度不高的情況下仍然可以獲得收斂的最終解.同時(shí)，該方法避免了微分修正算法等傳統(tǒng)數(shù)值方法推導(dǎo)相關(guān)梯度矩陣的復(fù)雜性，從而大大降低了問題求解的難度.除此之外，該方法又具有較好的收斂性和準(zhǔn)確性，可用來有效地求解高非線性程度、高敏感性的三體Lambert問題.

1 問題描述

三體Lambert問題的定義是，在三體系統(tǒng)中，求解經(jīng)過給定位置矢量R1和 R2、轉(zhuǎn)移時(shí)間為Δt12的轉(zhuǎn)移軌道.令起始點(diǎn)速度矢量記為V1，經(jīng)過時(shí)間 Δt12后，轉(zhuǎn)移軌道位置矢量為 R′2，定義函數(shù)F:

其中F是V1的函數(shù)，則三體Lambert問題的數(shù)學(xué)描述為求解下述方程的根:

上述方程沒有解析解，只能通過數(shù)值方法求解，但其非線性程度很高，一般情況存在多個(gè)解，很難確定其解得個(gè)數(shù)并將其全部求出.本文以地月系為例，地月系中兩天體間距相對(duì)較小，質(zhì)量比相對(duì)較大，故非線性程度和敏感性相對(duì)較高，考慮R1和R2都位于月球影響球之外，并且要求轉(zhuǎn)移軌道經(jīng)過月球影響球(否則可采用二體Lambert近似)，如圖1所示.

圖1 三體Lambert問題示意圖Fig.1 The three－body Lambert problem

2 基于二體模型的初值猜測(cè)

在三體Lambert問題的初始設(shè)計(jì)階段，不需要在整個(gè)轉(zhuǎn)移軌道上都考慮月球的影響，而僅僅在月球飛掠時(shí)，考慮由于月球引力引起的入射速度和出射速度之間的方向改變，其余都是地心圓錐曲線軌道.

假設(shè)起始點(diǎn)位置矢量為R1，起始時(shí)刻為t1，終止點(diǎn)位置矢量為R2，終止時(shí)刻為t2，則總轉(zhuǎn)移時(shí)間為Δt12=t2－t1，初值猜測(cè)算法流程見圖2，詳細(xì)步驟如下:

1)將近月點(diǎn)時(shí)刻記為tp，令tp=(t1+t2)/2，計(jì)算tp時(shí)刻月球的位置和速度向量Rm，Vm.

2)求解從R1到Rm、轉(zhuǎn)移時(shí)間為tp－t1和從Rm到R2、轉(zhuǎn)移時(shí)間為t2－tp的兩段地心圓錐曲線軌道，由此可以得到Rm處的地心速度矢量V1，V2，相應(yīng)的月心速度矢量為 v1=V1－Vm，v2=V2－Vm.

4)計(jì)算R1處速度矢量 V1，即求得猜測(cè)初值，作為下一步精確解求解的基礎(chǔ).

圖2 初值猜測(cè)流程圖Fig.2 Initial guess flowchart

3 UKF參數(shù)估計(jì)及其求解方法

參數(shù)估計(jì)問題，又被稱為系統(tǒng)辨識(shí)或機(jī)器學(xué)習(xí)問題，其目的是確定某個(gè)非線性映射:

針對(duì)目前國(guó)內(nèi)信息化建設(shè)普遍存在的問題，我院進(jìn)行了大量的調(diào)研和學(xué)習(xí)，基于全院上下對(duì)現(xiàn)代醫(yī)院信息化建設(shè)的深度認(rèn)同，進(jìn)行了科學(xué)詳細(xì)的方案規(guī)劃和明確的目標(biāo)設(shè)定。

該映射的輸入量為xk，輸出量為yk，w為非線性映射的參數(shù).一般來說，映射的輸入量xk和期望輸出量dk是不變的，輸出誤差定義為ek=dk－G(xk，w).求解參數(shù)估計(jì)問題就是要估計(jì)w的均值，使映射G(xk，w)的輸出誤差最小.

將原始參數(shù)估計(jì)問題寫成狀態(tài)空間表達(dá)式:

該式代表一個(gè)狀態(tài)轉(zhuǎn)移矩陣為單位陣的靜態(tài)過程.rk為過程噪聲;期望輸出dk則與對(duì)wk的非線性觀測(cè)相對(duì)應(yīng);ek為觀測(cè)噪聲.由此，原始參數(shù)估計(jì)問題就可以用擴(kuò)展卡爾曼濾波(EKF)，無損卡爾曼濾波(UKF)，粒子濾波(PF)等濾波器求解.從優(yōu)化問題的角度看，上述參數(shù)估計(jì)問題相當(dāng)于求解性能指標(biāo)如下以w為優(yōu)化變量的優(yōu)化問題:

其中，Re為觀測(cè)噪聲的方差陣，若Re為常值對(duì)角陣，則可以提到求和符號(hào)之外而不影響問題的求解，因此可以任意設(shè)定.過程噪聲r(shí)k的方差陣則影響濾波器的收斂速度和跟蹤性能.一般來說，Rrk越大，當(dāng)前狀態(tài)的濾波值中早期數(shù)據(jù)所占的比例就衰減得越快，越突出新息的作用.用濾波器求解參數(shù)估計(jì)問題的相關(guān)理論，文獻(xiàn)［15］作了詳細(xì)闡述，本文只作簡(jiǎn)要介紹.這里直接給出基于UKF濾波器的參數(shù)估計(jì)算法流程，如圖3所示.其中:

N是w的維數(shù);η是尺度參數(shù);常量ε決定了無損變換(UT)的σ點(diǎn)相對(duì)于w當(dāng)前均值的分布范圍，一般設(shè)為小量，取值范圍為［10－4，1］;常量 κ一般取為0或者3－N;β是與w的先驗(yàn)分布相關(guān)的常量，對(duì)于高斯分布，β=2是最優(yōu)的.ρRLS是遺忘因子，用于防止因模型誤差較大造成的濾波發(fā)散，其取值范圍為(0，1］.α是權(quán)重因子，取值范圍為［0，1］.

圖3 UKF參數(shù)估計(jì)流程圖Fig.3 UKF parameter estimation flow chart

4 UKF參數(shù)估計(jì)算法求解精確解

由于三體問題不存在解析解，所以只能通過數(shù)值積分對(duì)轉(zhuǎn)移軌道進(jìn)行精確預(yù)報(bào)，在猜測(cè)初值的基礎(chǔ)上進(jìn)行改進(jìn)，最終求得精確解.

將式(2)代表的三體Lambert問題改寫為參數(shù)估計(jì)問題，選擇待估計(jì)參數(shù)wk為起始點(diǎn)地心速度矢量V1，輸出dk為R2，輸入xk包括起始點(diǎn)位置矢量R1、起始時(shí)刻t1、終止時(shí)刻t2，則該問題可以表示為

其中rk和ek分別為系統(tǒng)噪聲和觀測(cè)噪聲.顯然，式(4)與式(7)形式上一一對(duì)應(yīng)，因而三體Lam－bert問題的精確解求解已轉(zhuǎn)化為參數(shù)估計(jì)問題，可以通過第3節(jié)的UKF參數(shù)估計(jì)算法求解.另外，待估計(jì)參數(shù)wk的各分量需要進(jìn)行單位化，以利于精確解的搜索過程.

值得注意的是，在求解三體Lambert問題的精確解的過程中，UKF濾波收斂的過程受算法中若干可調(diào)參數(shù)的影響很大，如系統(tǒng)噪聲矩陣Rr、尺度參數(shù)常量ε、遺忘因子ρRLS和權(quán)重因子α等，如選取不合適，則會(huì)導(dǎo)致濾波收斂過程的振蕩幅度很大，不能較快收斂，甚至發(fā)散.這些量的選取和更新方法屬于UKF濾波器算法的改進(jìn)范疇，不是本文的討論重點(diǎn)，可以作為下一步工作.本文通過大量數(shù)值仿真，總結(jié)算法中各個(gè)可調(diào)參數(shù)的推薦值，以供參考，如表1所示.

表1 UKF參數(shù)估計(jì)算法中的可調(diào)參數(shù)推薦值Table1 Reference values of adjustable parameters in the UKF parameter estimation method

5 算例與分析

本節(jié)給出一個(gè)包含引力輔助變軌的三體Lambert問題的求解算例.在J2000坐標(biāo)系中，起始點(diǎn)R1坐標(biāo)為［5048258，893447，－33213306］，終止點(diǎn)R2坐標(biāo)為［9472144，－7816649，31557762］，單位均為 m，轉(zhuǎn)移時(shí)間為 2014－01－01—2014－01－07，整個(gè)轉(zhuǎn)移時(shí)間為6 d.此算例起始點(diǎn)和終止點(diǎn)的位置在地球附近，類似于地月自由返回軌道的邊值條件，終端狀態(tài)相對(duì)于起始狀態(tài)具有很高的敏感度.

首先，利用基于二體模型的初值猜測(cè)方法，可得到該三體Lambert問題的初值為V1=［2924.54，－2100.25，－3 000.33］，單位均為m/s.為利于精確解的搜索，單位化后的V1作為待估計(jì)參數(shù)wk.

接下來，在初值的基礎(chǔ)上，進(jìn)一步對(duì)初值進(jìn)行修正，以獲得三體Lambert問題的收斂的精確解.在這里對(duì)微分修正算法和UKF參數(shù)估計(jì)算法進(jìn)行對(duì)比，精確解搜索的收斂標(biāo)準(zhǔn)為終點(diǎn)位置誤差在1 m以內(nèi).從以上的初值出發(fā)，精確解的搜索過程見表2和表3.其中，表2為微分修正算法的搜索過程，從表中可得出，其搜索過程是不斷震蕩甚至是發(fā)散的，表3為UKF參數(shù)估計(jì)算法的搜索過程，從表中可得出，經(jīng)過12次迭代后，其搜索過程最終收斂.

表2 微分修正算法精確解搜索過程Table2 Iteration of searching the final solution using the differential-correction method m

表3 UKF參數(shù)估計(jì)算法精確解搜索過程Table3 Iteration of searching the final solution using the UKF parameter estimation method m

由此，UKF參數(shù)估計(jì)算法在解決三體Lambert問題中的有效性得以驗(yàn)證，并且UKF參數(shù)估計(jì)算法比微分修正算法具有更大的收斂域.該三體Lambert問題最終的飛行軌跡見圖4.

圖4 三體Lambert問題的飛行軌跡Fig.4 Trajectory of the three－body Lambert problem

為了詳細(xì)研究UKF參數(shù)估計(jì)算法的收斂域，并與微分修正算法、二階微分修正算法［14］進(jìn)行對(duì)比，可以采取下面方法.對(duì)于上述算例最終解的某一個(gè)分量添加擾動(dòng)，而另外兩個(gè)分量保持不變.從這個(gè)擾動(dòng)點(diǎn)出發(fā)，分別使用微分修正算法、二階微分修正算法和UKF參數(shù)估計(jì)算法來搜索轉(zhuǎn)移軌道的精確解，不斷增加擾動(dòng)量，一直到精確解搜索過程發(fā)散，由此得到這3種算法對(duì)于特性的擾動(dòng)分量的收斂域.雖然這不是該問題收斂域的完整描述，但是也部分揭示了各種算法收斂域的基本特性，也可以體現(xiàn)各種算法的優(yōu)劣.3種精確解搜索算法的收斂域統(tǒng)計(jì)信息見表4.

表4 各算法的收斂域統(tǒng)計(jì)Table4 Statistics of convergence domains of various methods m/s

上述結(jié)果可以得出結(jié)論:①設(shè)計(jì)變量V1的單個(gè)分量收斂域與約束條件之間沒有特定的規(guī)律，而僅僅有很大的變化區(qū)間，體現(xiàn)了該問題收斂域具有十分復(fù)雜的幾何結(jié)構(gòu)，這也是由于三體Lambert問題的高度非線性特性導(dǎo)致的;②在保證相當(dāng)精度的情況下，平均水平上看，UKF參數(shù)估計(jì)算法的收斂范圍是微分修正算法、二階微分修正算法收斂域的3～5倍.另外，對(duì)于3種算法均收斂的算例，在Intel Core 2.53 GHz，3 GB RAM的計(jì)算條件下，微分修正算法、二階微分修正算法、UKF參數(shù)估計(jì)算法的平均計(jì)算時(shí)間分別為1.96，2.70，14.95 s.

為了進(jìn)一步研究UKF參數(shù)估計(jì)算法搜索精確解的整體性能，采用更多的隨機(jī)數(shù)值算例來驗(yàn)證.令起始時(shí)間在 2014－01－01—2014－01－30(1 個(gè)月球周期)之間隨機(jī)變化，轉(zhuǎn)移時(shí)間在5～7 d之間隨機(jī)變化，起始點(diǎn)和終止點(diǎn)位置類似于地月自由返回軌道的邊值條件，計(jì)算100個(gè)算例.在基于二體模型猜測(cè)初值的基礎(chǔ)上，選擇可調(diào)尺度參數(shù)ε 為5×10－4或8×10－4，UKF 參數(shù)估計(jì)算法收斂概率為98%，收斂次數(shù)在7～18次.然后，只需稍微更改尺度參數(shù)常量ε(如1×10－4)，可使得余下的2%算例收斂，且具有相當(dāng)?shù)氖諗看螖?shù).經(jīng)過進(jìn)一步的算例驗(yàn)證，若采用更精確的初值猜測(cè)方法，如偽狀態(tài)方法，也可獲得相當(dāng)?shù)氖諗啃阅?由此可見，采用UKF參數(shù)估計(jì)算法求解三體Lambert問題的精確解具有良好的收斂性能.

6 結(jié)論

本文提出了一種基于UKF參數(shù)估計(jì)的從初步設(shè)計(jì)到精確設(shè)計(jì)的三體Lambert問題求解方法.通過數(shù)值算例驗(yàn)證，該方法收斂次數(shù)較少，具有較好的魯棒性，而且降低了對(duì)初值精確度的要求，即使利用二體模型給出的初值，也可以收斂得到精度較高的精確解，同時(shí)避免了傳統(tǒng)數(shù)值方法對(duì)相關(guān)梯度矩陣的推導(dǎo)，因此顯著降低了三體Lambert問題求解的難度，可以有效地解決高非線性、高敏感度的三體Lambert問題.另外，由于該方法適用于各種非線性映射的參數(shù)估計(jì)，可以在三體Lambert問題的基礎(chǔ)上，進(jìn)一步研究星際間引力輔助飛行等問題，具有廣泛的應(yīng)用前景.

References)

［1］ Bate R，Mueller D，White J.Fundamentals of astrodynamics［M］.New York:Dover Publications，1971:177－275.

［2］ Battin R H，Vaughan R M.An elegant Lambert algorithm［J］.Journal of Guidance，Control and Dynamics，1984，7(6):662－670.

［3］ Gooding R H.A procedure for the solution of Lambert’s orbital boundary－value problem［J］.Celestial Mechanics ＆ Dynamical Astronomy，1990，48(2):145－165.

［4］ D’Amarion L，Byrnes D，Sackett L.Optimization of multiple flyby trajectories［C］//AAS/AIAA Astrodynamics Specialists Conference.Provincetown:AIAA Paper 1979:79－162.

［5］ Armellin R，Di Lizia P，Topputo F，et al.Gravity assist space pruning based on differential algebra［J］.Celestial Mechanics and Dynamical Astronomy，2010，106(1):1－24.

［6］ Jesicak M，Ocampo C.Automated generation of symmetric lunar free－return trajectories［J］.Journal of Guidance，Control and Dynamics，2011，34(1):98－106.

［7］ Luo Q，Yin J，Han C.Design of earth－moon free－return trajectories［J］.Journal of Guidance，Control，and Dynamics，2012，36(1):263－271.

［8］ Okutsu M，Longuski J.Mars free returns via gravity assist from Venus［J］.Journal of Spacecraft and Rockets，2002，39(1):31－36.

［9］ Prado A F B A.Traveling between the Lagrangian points and the Earth［J］.Acta Astronautica，1996，39(7):483－486.

［10］ Lian Y J，Jiang X Y，Tang G J.Halo－to－h(huán)alo cost optimal transfer based on CMA－ES［C］//Proceedings of the 32nd Chinese Control Conference，CCC 2013.Piscataway，NJ:IEEE，2013:2468－2473.

［11］ Zazzera F B，Topputo F，Massari M.Assessment of mission design including utilization of libration points and weak stability boundaries，18147/04/NL/mv［R］.Frascati，Italy:ESA，2003.

［12］ Byrnes D V.Application of the pseudostate theory to the threebody Lambert problem［J］.Journal of the Astronautical Sciences，1989，37:221－232.

［13］ Sukhanov A，Prado A F B A.Lambert problem solution in the Hill model of motion［J］.Celestial Mechanics ＆ Dynamical Astronomy，2004，90(3):331－354.

［14］ 羅欽欽，韓潮.包含引力輔助變軌的三體Lambert問題求解算法［J］.北京航空航天大學(xué)學(xué)報(bào)，2013，39(5):679－687.Luo Q Q，Han C.Solution algorithm of the three－body Lambert problem with gravity assist maneuver［J］.Journal of Beijing University of Aeronautics and Astronautics，2013，39(5):679－687(in Chinese).

［15］ Haykin S.Kalman filtering and neural networks［M］.New York:John Wiley ＆ Sons Inc，2002.