上海微小衛星工程中心,上海 201203
探測器在地月轉移過程中,由于各種誤差的存在,實際軌道偏離標稱軌道。為確保探測器順利進入目標軌道,探測器在飛行期間需要進行多次修正。中途修正解決的問題有兩方面:如何確定修正時機使修正速度總量最少;終端誤差滿足任務要求及如何確定單次修正速度增量。
針對第一個方面的研究相對成熟[1-5],這里只做簡單的介紹。國外研究中有間距比策略[1],即盡快進行第一次中途修正,以后各次修正時機的選擇是按照修正效果與前一次修正效果的比值為常數這一條件確定的。國內研究主要圍繞嫦娥探測器展開,文獻[2-3]給出了典型的兩次中途修正策略,即第一次修正在飛行的第一天實施,第二次修正在到達月球的前一天實施。
本文主要針對第二個問題展開研究。其本質上是一個兩點邊值問題,由于三體動力學的復雜性,不存在解析解,必須利用數值方法進行求解。一般的思路是:首先根據邊值約束,采用簡化模型設計初值,之后考慮精確攝動模型,采用微分修正等數值方法求解。
對于地月轉移軌道中途修正問題,可以將實際速度作為初值,若探測器實際速度誤差過大,初值精度變差,修正算法的求解性能降低,因此有必要對初值進行改進。若修正時刻探測器的位置在月球段,探測器主要受月球引力的作用,初值設計本質上是求解二體Lambert問題,文獻[6-8]給出了經典的求解方法,若在地球段,需要考慮月球引力的影響,目前設計方法有圓錐拼接方法[9-15],但圓錐曲線拼接方法精度有限,導致修正算法的求解效率降低。
在使用微分修正求解精確解時,需要計算狀態轉移矩陣,在常用的數值方法中[3,16-18],首先需要推導加速度相對位置速度的偏導數矩陣,然后同時積分動力學方程就可以得到狀態轉移矩陣,但在每次迭代計算中,需要積分42個方程,導致求解時間增加。
針對以上研究現狀,本文首先利用偽狀態理論,結合Vinti預報方法[19],設計高精度的初值,在求解精確解時,提出了一種基于偽狀態理論的狀態轉移矩陣解析算法。該算法通過設計高精度的初值,降低了求解修正速度增量的難度,而且避免了傳統數值方法計算狀態轉移矩陣的復雜性。
由于修正時刻探測器的位置無法改變,只能將速度和轉移時間作為設計變量,本文采用轉移時間固定的修正策略,因此初始修正參數只有速度3個分量,即:
(1)
式中:v01,v02和v03分別為修正時刻探測器在地球慣性系下的速度分量。
對于月球探測任務,通常對到達月球時刻的月心距、軌道傾角以及飛行角有嚴格要求,因此選取目標軌道參數為:
(2)
式中:rp,iL以及γL分別為到達月球時刻探測器的月心距、月球固聯系下的軌道傾角以及相對月球的飛行角。
偽狀態理論是文獻[20]提出的一種地月中心引力作用下的軌道近似計算方法。
令tI時刻探測器相對地球的狀態為ΣI=(RI,VI),利用偽狀態理論計算tk時刻探測器相對地球的狀態Σk=(Rk,Vk)和相對月球的狀態σk=(rk,vk)。定義以tk時刻月球位置為中心的偽狀態轉換球(Pseudostate Transformation Sphere,PTS),半徑為RPTS,在本文算例中,其值取為24個地球半徑。假設RI在PTS外部,Rk在PTS內部,如圖1所示。

圖1 偽狀態理論Fig.1 The pseudostate theory
具體計算步驟如下:


(3)


本節給出了地球慣性系下的中途修正初值的求解方法。令修正時刻探測器位置R0,到達月球時刻tL,轉移時間Δt,迭代變量RS,求解過程如下:
1)計算tL時刻月球相對地球的位置速度分別為RM和VM,令RS=RM。
2)利用Lambert算法,求解R0到RS,轉移時間Δt的地球二體軌道,可得到R0和RS處的速度V0和VS。
3)計算tL時刻探測器相對月球速度在月球固聯系下的經度aL和緯度δL,求解tL時刻探測器在月球固聯系下的升交點經度ΩL,
sin(aL-ΩL)=tanδL/taniL
(4)
根據式(4),ΩL有兩個解,對應到達月球的升降軌方式,則角動量hL在月球固聯系下的表達式為:
hL=[siniLcosΩL-siniLsinΩLcosiL]T
(5)
4)令探測器在PTS邊界時相對月球的位置和速度分別為rS和vS,對于地月轉移軌道,rS在近月點前,大小為RPTS,根據能量守恒,可得近月點速度大小為:
(6)

(7)
令rS到近月點的轉移時間為tS,根據偽狀態理論可得
(8)

對于地月轉移軌道,地球扁率攝動是一項不可忽視的攝動源,因此在精確初值設計時,需要考慮地球扁率的影響。
在求解地球扁率攝動的解析方法中,Vinti預報方法精度高,速度快,因此采用Vinti預報方法來修正地球扁率攝動的影響,求解過程如下:
1)求解R0到RS,轉移時間為Δt的地球二體軌道,可以得到R0地球速度V0k。
2)以R0和V0k為初值,利用Vinti預報方法計算tL時刻的位置RSk。
3)更新V0k,其中Jk可近似為二體狀態轉移矩陣:
(9)
4)若‖RS-Rsk‖大于給定閾值,返回步驟2),否則,輸出V0k作為初值。
為得到精確解,必須在更復雜的攝動環境下通過數值積分計算精確的轉移軌道,然后采用微分修正方法對初值進行修正。在精確的轉移軌道計算中,本文考慮了地球非球形項,日月中心引力的影響,其中,地球非球形引力取8×8階,月球和太陽星歷均采用DE405。
初始修正參數P和目標軌道參數Q存在確定關系,記為Q=g(P),則存在以下關系:
(10)
根據微分修正方法,在每次迭代中初始修正參數的調整量為:
(11)
式中:ΔQ為目標軌道參數的偏差量。令狀態變量X=[R,V]T,R和V分別表示探測器在地球慣性下的位置和速度,則式(10)可以重寫為:
(12)
式中:Xi,Xf分別為狀態變量X的初始狀態和終止狀態。式(12)誤差傳遞矩陣計算中,除狀態轉移矩陣Φfi外均有確定的形式。三體軌道動力學特性導致Φfi沒有解析表達式,通常采用數值方法計算。為提高求解精確解的效率,本節給出了一種基于偽狀態理論的狀態轉移矩陣近似解析計算方法。
根據偽狀態理論,存在以下關系:
(13)
式(13)第一項具體表達式為:
(14)

(15)
式(13)第二項可以分解為以下兩項:
(16)
(17)
式中:
(18)
根據關系式:
(19)
對式(19)求導可得:
(20)
本節給出了一個地月轉移軌道中途修正的求解算例,首先根據本文給出的初值設計方法給出初值,然后利用解析狀態轉移矩陣求解精確解。
假設探測器在2022年1月1日進行修正,此時探測器在地球慣性坐標系下位置為[-4 094.687,4 730.547,2 588.812] km,轉移時間為72 h,近月點高度為600 km,達到月球時軌道傾角25°,到達月球方式為升。
在無奇異情況下,Lambert問題具有長短程兩個解,所以目標軌道參數確定以后,中途修正有兩個解,相應的軌跡如圖2所示,其對應的地心張角不同。經過計算可以得到兩個可行解。其中:短程解為v1=[-8.564 25,-5.888 53,-2.785 82] km/s,長程解為v2=[9.184 99,5.093 68,2.351 90] km/s,v1和v2夾角為174.2°,方向幾乎相反。

圖2 地月轉移軌道Fig.2 Trans-lunar trajectory
通常長程解被忽略,但在某些任務中,如探測器并非直接由運載送入地月轉移軌道,而是在近地軌道上依靠自身推進系統進入地月轉移軌道,長程解可作為地月轉移軌道的選擇,本文算法直接可以設計地月轉移軌道。
下面通過一個算例來說明本文算法的求解性能。選擇一條標準地月轉移軌道,其發射時刻為2022年1月1日,轉移時間120 h。近地點高度為300 km,地球出發時軌道傾角為25°,地球出發方式為降軌,近月點高度為500 km,到達月球時軌道傾角為35°,到達月球方式為升軌,在地球慣性系下,探測器出發時刻的標稱位置速度分別為[ -5 848.290,2 683.276,1 760.631] km和[-5.135 05,-8.845 35,-3.576 42] km/s。假設探測器進入地月轉移軌道時半長軸偏差為10 000 km,在入軌24 h后進行修正。
如果將實際速度作為修正初值,其搜索過程如表1所示,由于誤差過大,導致迭代發散,無法得到精確解。若采用精確的初值,經過8次迭代就可得到精確解,迭代過程如表2所示。通過初值設計,提高了初值精度,保證了修正算法的收斂性和效率。

表1 近月點狀態偏差迭代過程Table 1 Iteration of the error of perilune

表2 利用精確初值的近月點狀態偏差迭代過程Table 2 Iteration of the error of perilune using the initial solution with high accuracy
表3給出了相同計算條件下(Intel Core 2.53G),分別利用解析法計算狀態轉移矩陣(簡稱解析法)和數值法計算狀態轉移矩陣(簡稱數值法)求解精確解的計算結果,盡管兩種方法均可以得到精確解,但若采用解析法求解狀態轉移矩陣,計算時間更少,求解效率更高,而且算法結構簡單,避免數值計算的復雜性。

表3 計算結果Table 3 Simulation results
本文針對地月轉移軌道中途修正問題,提出了一種求解修正速度增量的方法。該算法通過改善初值的精度和簡化狀態轉移矩陣,提高了地月轉移中途修正的求解效率,可有效解決地月轉移中途修正問題,可為月地轉移等深空探測軌道設計以及中途修正提供借鑒和參考。