曹鵬飛，李維國，王俊彥，李海陽

（1.北京航天飛行控制中心，北京100094；2.中國衛(wèi)星發(fā)射測控系統(tǒng)部，北京100120；3.國防科學技術大學空天科學學院，長沙410073）

引 言

早在20 世紀Apollo 時代，美國的Farquhar 提出了利用地月L2 點（以下簡稱LL2）附近的Halo 軌道實現(xiàn)月球背面與地球持續(xù)通信[1]。近年來，隨著平動點任務的不斷提出與實施，攝動模型下的Halo 軌道設計問題一直是研究熱點。Farquhar 等[2]較早研究了如何控制航天器使其運行在Halo軌道的近似軌道上。Gomez等[3]基于Floquet理論和不變流形理論，給出了一種攝動模型下的軌道設計方法。Howell 等[4-5]基于標稱軌跡提出了Halo 軌道靶點法控制策略，該方法具有適用范圍廣、簡單可靠等優(yōu)點，但由于其控制精度和代價函數(shù)中的權重矩陣與保持間隔密切相關，有一定局限性。Qi等[6]在Howell研究基礎上，進一步提出了Halo 軌道正切靶點法軌道設計方法。Mehrdad等[7]考慮了太陽引力攝動與月球軌道偏心率攝動，采用線性二次調節(jié)器（Linear Quadratic Regulator，LQR）和多重打靶策略（Multiple Shooting，MS）設計了地月平動點Halo 軌道，但過程較為復雜，且燃料成本偏高。徐明等[8]基于Halo 軌道偏差動力學方程，提出了Halo 軌道線性周期保持控制策略。彭海軍等[9]基于圓型限制性三體問題（Circular Restricted Three Body Problem，CR3BP），采用SDRE（State Dependent Riccati Equation，SDRE）方法設計了用于平動點軌道維持控制的非線性次優(yōu)跟蹤控制器。車征等[10]在CR3BP 的基礎上引入了太陽引力攝動，研究了LL2 點Halo 軌道設計問題。劉磊等[11]基于CR3BP下的平動點軌道運動特性和微分修正策略提出了連續(xù)環(huán)繞控制方法，該方法中平動點軌道初值取自CR3BP一階近似結果。

本文在文獻[11]研究思路基礎上，提高了初值精度，進一步提出了高精度多層序列二次規(guī)劃（Se‐quence Quadratic Program，SQP）修正的求解方法。首先，推導了CR3BP 質心會合坐標系與地心J2000坐標系之間的轉換關系，并在此基礎上，將CR3BP下的閉合Halo 軌道轉換到地心J2000 坐標系得到了高精度模型下Halo 軌道迭代初值。其次，采用SQP構造多層迭代格式，在高精度模型下對初值進行逐層修正。通過仿真測試，文章所述方法的可行性與有效性得到了驗證。

1 動力學模型

本文采用兩種動力學模型，即CR3BP 模型和高精度模型，模擬地月空間飛行器的運動規(guī)律。

1.1 圓型限制性三體問題

CR3BP，即假設大天體和小天體繞其公共質心做勻速圓周運動，研究第三體在兩主天體共同引力下的運動問題。在CR3BP 中，常用的坐標系為質心會合坐標系B-xyz，對于地月系統(tǒng)，其原點位于地月公共質心B，x軸由地球指向月球，z軸指向地月系統(tǒng)角動量方向，y軸與其他兩軸構成右手坐標系。

為簡化動力學方程，通常引入歸一化單位，地月系統(tǒng)長度單位為地月質心平均距離DU=384 400 km，質量單位MU 為地球質量M1和月球質量M2之和，時間單位可有上述兩參數(shù)導出。引入質量比μ，即

航天器在B-xyz坐標系下的無量綱動力學方程為

其中：U為與位置相關的偽勢能函數(shù)，即

式中

1.2 高精度模型

高精度模型考慮了地球中心引力、日月引力攝動、地月非球形引力攝動、太陽光壓攝動以及金星、木星等三體攝動等。在地心J2000坐標系下，忽略地球潮汐和相對論效應等微小攝動量的影響。

航天器的動力學方程為

其中：r為地心位置矢量；μE為地球引力常數(shù)；AN為N體引力攝動加速度，星體間空間幾何關系通過DE405/LE405星歷求解；ANSE為地球非球形攝動，取WGS84引力場模型8×8階計算；ANSM為月球非球形攝動，取LP165P引力場模型8×8階計算；AR為太陽光壓攝動；AP為推進加速度。

1.3 Halo軌道

Halo 軌道是共線平動點附近的三維周期軌道，在平動點任務中被廣泛應用。CR3BP會合坐標系下，Halo 軌道關于xz面對稱，與xz面交于兩點，通常取其中距離x軸較遠的點與x軸的距離作為表征Halo軌道大小的參數(shù)，稱之為法向幅值Az。CR3BP 下，Halo軌道的計算通常先采用Richardson三階近似解析解[12]獲取初值，然后再利用微分修正法對初值進行修正，詳細推到過程可參考文獻[13]。

2 高精度模型下Halo軌道設計

2.1 初值計算

CR3BP 下的閉合Halo 軌道可為高精度模型下Halo 軌道設計提供良好初值，但由于動力學模型的差異，二者的轉換過程并不直觀。

本節(jié)以地月系統(tǒng)為例，詳細推導了CR3BP 質心會合坐標系B-xyz與地心J2000 坐標系OE-XJYJZJ之間的轉換關系。如圖1所示，給出了兩坐標系的幾何構型圖。假設航天器在OE-XJYJZJ中的位置矢量為R=[xJ,yJ,zJ]T， 單 位 為 km， 速 度 矢 量 為V=，單位為km/s；在B-xyz中的位置矢量為r=[x,y,z]T，速度矢量為v=，單位采用歸一化單位。引入自定義坐標系——高精度質心會合坐標系B-XYZ，其原點位于地月質心，坐標軸指向與坐標系B-xyz相同，區(qū)別在于該坐標系中的地月和平動點的位置并非固定，而是實時變化。

圖1 坐標系幾何構型圖Fig. 1 The geometries of the coordinate frames

首先，推導位置矢量轉換關系。假設歷元時刻為t0，月球軌道半長軸為aM，偏心率為eM，升交點赤經(jīng)為ΩM，白赤交角為iM，近地點幅角為ωM，真近點角為fM。由二體問題可知，地月距離為

其中：EM為偏近點角。

記歸一化后的L為，即

由坐標系定義可知，在B-XYZ中，地球位置矢量為rB1=， 航 天 器 位 置 矢 量 為rB2=。因此，航天器相對于地球的位置矢量為

下一步，將rB3轉換到地心J2000坐標系，該轉換由ΩM、iM和uM三個歐拉角決定。計算公式為

其中：uM為緯度幅角。

由圖1 可知，B-XYZ到地心J2000 坐標系的轉換矩陣為

矩陣M1[λ]與M3[λ]滿足的條件為

在完成坐標旋轉后，將歸一化單位轉換為國際單位制即可完成位置矢量的轉換，具體公式如下

速度轉換思路與位置轉換類似，區(qū)別在于前者需要考慮牽連加速度項。在高精度模型下，地月距離實時變化，其徑向變化速率為

其中：nM為月球公轉平均角速度。

由式（9）可知，在B-XYZ中航天器相對地心的速度矢量為

引入地心白道慣性系OE-XIYIZI，其原點位于地心，坐標軸指向與某一瞬時的B-XYZ相同。B-XYZ相對于地心的角速度矢量在OE-XIYIZI中可以表示為ωΜ=[0,0,ω]Τ，其中ω為瞬時旋轉速率。記歸一化后的ω為ωnorm，表達式為

進一步可得，B-XYZ相對于OE-XIYIZI的牽連速度，歸一化后的表達式為

綜合式（17）和（19），在完成坐標和單位轉換后即可完成速度矢量的轉換，具體公式如下

2.2 多層序列二次規(guī)劃修正

由于復雜攝動力的影響，高精度模型下的Halo軌道將不再閉合。為了便于分析，本節(jié)將參照CR3BP 下的Halo 軌道特性，給出高精度模型下Halo軌道圈數(shù)的定義。如圖2 所示，在坐標系B-XYZ中，Halo 軌道由xz面出發(fā)，第2 次到達xz面時軌道為1/2圈，第3 次到達xz面為1 圈，第4 次到達xz面時為2圈，依次類推。

圖2 高精度模型下Halo軌道的圈數(shù)定義Fig. 2 The definition of the number of Halo orbits in the high precision model

在B-xyz坐標系下，Halo軌道關于xz面對稱，即在xz面處x或z方向的速度分量為零?；贖alo軌道該特性，本節(jié)將采用SQP 構造多層迭代格式，給出高精度模型Halo軌道設計方法，具體流程為：

第1步：參考文獻[13]，計算CR3BP質心會合坐標系下Halo軌道在xz處的積分狀態(tài)量x0；

第2 步：采用2.1 小節(jié)方法，將x0轉換到地心J2000坐標系，轉換結果記為X0，將其作為高精度模型下Halo軌道修正初值；

第3步：在高精度模型下采用SQP對X0的3個速度分量進行多層修正，具體流程如圖3所示。

圖3 初值修正流程圖Fig. 3 The flow chart of initial value correction

第1層：優(yōu)化變量為施加的速度增量ΔV1，即

約束條件為軌道第1 次到達xz面時Vx= 0，優(yōu)化結果記為ΔV2。高精度模型下，軌道在xz面處x或z方向的速度分量接近于零但不嚴格為零，因此ΔV1并非最優(yōu)解，需要進一步修正。

第2 層：優(yōu)化變量為ΔV2，約束條件為軌道第2次到達xz面時Vx= 0，優(yōu)化結果記為ΔV3。依次類推，當修正層數(shù)n=4時，修正脈沖的量級已非常小。當n> 4 時，由于修正脈沖的量級接近攝動力量級，程序將終止。因此對于地月系統(tǒng)，修正層數(shù)n不宜大于4次，該結論與文獻[11]結果相吻合。

高精度模型下為了得到飛行多圈的光滑Halo 軌道，需要定期對軌道進行維持控制，單次維持所消耗的速度增量的求解流程與圖3 類似，這里不再贅述。且由前面的分析可知，速度增量的施加位置為xz面處，維持間隔應不宜大于2圈。

3 算例驗證

給出高精度模型下Halo 軌道設計實例，參考文獻[7]給出參數(shù)配置：Halo軌道法向幅值Az為30 000 km，周期為14.25 天，方向為南向，位于LL2 點附近；歷元時刻為2025 年1 月1 日00:00:0.000UTCG；跳秒與STK11 一致，取37 s；選用RKF7（8）變步長積分器，相對誤差和絕對誤差設為10-13。

首先，通過采用文獻[13]提供的方法，計算CR3BP質心會合坐標系下Halo軌道，并利用式（13）和式（20）將其初始點參數(shù)轉換到地心J2000 坐標系。其次，采用2.2 節(jié)所述方法對初始點速度分量進行修正，為加快計算收斂速度，這里只修正Y方向的速度分量，修正層數(shù)n設為4。地心J2000坐標系下修正前后的初始點參數(shù)見表1。

表1 Halo軌道初始點參數(shù)Table 1 The initial point parameters of Halo orbit

直接將修正前的初始點參數(shù)代入高精度模型積分一個周期，并將軌道數(shù)據(jù)實時轉換到高精度LL2點會合坐標系，得到的軌道如圖4所示。將修正后的初始點參數(shù)代入高精度模型，積分1圈，即終止條件為軌道第2次到達xz面時Vx= 0，并將軌道數(shù)據(jù)實時轉換到高精度LL2 點會合坐標系，得到的軌道如圖5所示。

圖4 修正前初始點參數(shù)遞推得到的軌跡Fig. 4 The trajectory derived from the unrevised initial point parameters

圖5 修正后初始點參數(shù)遞推得到的軌跡Fig. 5 The trajectory derived from the revised initial point parameters

對比圖4和圖5可知，CR3BP下的閉合Halo軌道在高精度模型下飛行約半個周期后將急劇發(fā)散，而修正后軌道在高精度模型下周期性保持較好，但仍未閉合。因此，對于長期飛行于Halo 軌道的航天器，需要定期進行軌道維持控制。

仍以上述算例為基礎，給出高精度模型下Halo軌道維持實例。采用圖3的計算流程，SQP修正層數(shù)n設為4，維持間隔設為1圈，對圖5中的軌道進行為期1年的維持，得到的飛行軌跡在高精度LL2點會合坐標系下的投影如圖6所示，維持總次數(shù)和消耗的總速度增量見表2。

表2 維持間隔與維持總速度增量的關系Table 2 The relationship between the maintenance interval and the maintenance impulsive consumption

圖6 高精度地LL2點會合坐標系下飛行一年的Halo軌道軌跡Fig. 6 The flying trajectories of the Halo orbit for one year in the high precision LL2 point synodic coordinate system

此外，表2還給出了不同維持間隔方案時的維持總次數(shù)和總的速度增量消耗信息。由表2可知：①當維持間隔為0.5圈、1圈和1.5圈時，所消耗的總速度增量差異并不大；②當維持間隔由1.5 圈增加為2 圈時，所消耗的速度增量急劇增加，這是由于維持間隔過大造成累積偏差較大，從而導致速度增量消耗增多；③對于LL2 點Halo 軌道，綜合考慮速度增量消耗與維持頻率，最佳維持間隔應為1.5 圈。對比文獻[7]的計算結果，發(fā)現(xiàn)多層SQP 修正方法在高精度Halo 軌道維持方面具有維持間隔長和燃料成本低等優(yōu)點。

4 結 論

本文研究了高精度模型下共線平動點附近Halo軌道的設計問題。推導了CR3BP 質心會合坐標系與高精度模型地心J2000坐標系之間的轉換關系，并在此基礎上，將CR3BP 下的閉合Halo 軌道轉換到地心J2000 坐標系得到了高精度模型下Halo 軌道迭代初值。通過采用SQP 對初值進行多層修正，得到了在高精度模型下的Halo 軌道。進一步研究發(fā)現(xiàn)，該方法還可用于攝動模型下Halo軌道的維持設計。最后，仿真算例驗證了方法的可行性與正確性，研究結果可為未來平動點任務的標稱軌道設計提供參考。