曹鵬飛，李維國(guó)，王俊彥，李海陽(yáng)

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

引 言

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

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

1 動(dòng)力學(xué)模型

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

1.1 圓型限制性三體問(wèn)題

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

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

航天器在B-xyz坐標(biāo)系下的無(wú)量綱動(dòng)力學(xué)方程為

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

式中

1.2 高精度模型

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

航天器的動(dòng)力學(xué)方程為

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

1.3 Halo軌道

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

2 高精度模型下Halo軌道設(shè)計(jì)

2.1 初值計(jì)算

CR3BP 下的閉合Halo 軌道可為高精度模型下Halo 軌道設(shè)計(jì)提供良好初值，但由于動(dòng)力學(xué)模型的差異，二者的轉(zhuǎn)換過(guò)程并不直觀。

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

圖1 坐標(biāo)系幾何構(gòu)型圖Fig. 1 The geometries of the coordinate frames

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

其中：EM為偏近點(diǎn)角。

記歸一化后的L為，即

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

下一步，將rB3轉(zhuǎn)換到地心J2000坐標(biāo)系，該轉(zhuǎn)換由ΩM、iM和uM三個(gè)歐拉角決定。計(jì)算公式為

其中：uM為緯度幅角。

由圖1 可知，B-XYZ到地心J2000 坐標(biāo)系的轉(zhuǎn)換矩陣為

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

在完成坐標(biāo)旋轉(zhuǎn)后，將歸一化單位轉(zhuǎn)換為國(guó)際單位制即可完成位置矢量的轉(zhuǎn)換，具體公式如下

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

其中：nM為月球公轉(zhuǎn)平均角速度。

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

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

進(jìn)一步可得，B-XYZ相對(duì)于OE-XIYIZI的牽連速度，歸一化后的表達(dá)式為

綜合式（17）和（19），在完成坐標(biāo)和單位轉(zhuǎn)換后即可完成速度矢量的轉(zhuǎn)換，具體公式如下

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

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

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

在B-xyz坐標(biāo)系下，Halo軌道關(guān)于xz面對(duì)稱，即在xz面處x或z方向的速度分量為零。基于Halo軌道該特性，本節(jié)將采用SQP 構(gòu)造多層迭代格式，給出高精度模型Halo軌道設(shè)計(jì)方法，具體流程為：

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

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

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

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

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

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

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

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

3 算例驗(yàn)證

給出高精度模型下Halo 軌道設(shè)計(jì)實(shí)例，參考文獻(xiàn)[7]給出參數(shù)配置：Halo軌道法向幅值A(chǔ)z為30 000 km，周期為14.25 天，方向?yàn)槟舷颍挥贚L2 點(diǎn)附近；歷元時(shí)刻為2025 年1 月1 日00:00:0.000UTCG；跳秒與STK11 一致，取37 s；選用RKF7（8）變步長(zhǎng)積分器，相對(duì)誤差和絕對(duì)誤差設(shè)為10-13。

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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