UKF參數估計在三體Lambert問題中的應用

張洪禮，羅欽欽，韓潮*

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

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

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

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

1 問題描述

三體Lambert問題的定義是，在三體系統中，求解經過給定位置矢量R1和 R2、轉移時間為Δt12的轉移軌道.令起始點速度矢量記為V1，經過時間 Δt12后，轉移軌道位置矢量為 R′2，定義函數F:

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

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

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

2 基于二體模型的初值猜測

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

假設起始點位置矢量為R1，起始時刻為t1，終止點位置矢量為R2，終止時刻為t2，則總轉移時間為Δt12=t2－t1，初值猜測算法流程見圖2，詳細步驟如下:

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

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

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

圖2 初值猜測流程圖Fig.2 Initial guess flowchart

3 UKF參數估計及其求解方法

參數估計問題，又被稱為系統辨識或機器學習問題，其目的是確定某個非線性映射:

針對目前國內信息化建設普遍存在的問題，我院進行了大量的調研和學習，基于全院上下對現代醫院信息化建設的深度認同，進行了科學詳細的方案規劃和明確的目標設定。

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

將原始參數估計問題寫成狀態空間表達式:

該式代表一個狀態轉移矩陣為單位陣的靜態過程.rk為過程噪聲;期望輸出dk則與對wk的非線性觀測相對應;ek為觀測噪聲.由此，原始參數估計問題就可以用擴展卡爾曼濾波(EKF)，無損卡爾曼濾波(UKF)，粒子濾波(PF)等濾波器求解.從優化問題的角度看，上述參數估計問題相當于求解性能指標如下以w為優化變量的優化問題:

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

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

圖3 UKF參數估計流程圖Fig.3 UKF parameter estimation flow chart

4 UKF參數估計算法求解精確解

由于三體問題不存在解析解，所以只能通過數值積分對轉移軌道進行精確預報，在猜測初值的基礎上進行改進，最終求得精確解.

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

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

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

表1 UKF參數估計算法中的可調參數推薦值Table1 Reference values of adjustable parameters in the UKF parameter estimation method

5 算例與分析

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

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

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

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

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

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

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

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

表4 各算法的收斂域統計Table4 Statistics of convergence domains of various methods m/s

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

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

6 結論

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

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－halo 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］.北京航空航天大學學報，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.