一種地月轉(zhuǎn)移軌道中途修正的制導(dǎo)算法

上海微小衛(wèi)星工程中心，上海 201203

探測器在地月轉(zhuǎn)移過程中，由于各種誤差的存在，實(shí)際軌道偏離標(biāo)稱軌道。為確保探測器順利進(jìn)入目標(biāo)軌道，探測器在飛行期間需要進(jìn)行多次修正。中途修正解決的問題有兩方面：如何確定修正時(shí)機(jī)使修正速度總量最少;終端誤差滿足任務(wù)要求及如何確定單次修正速度增量。

針對第一個(gè)方面的研究相對成熟[1-5]，這里只做簡單的介紹。國外研究中有間距比策略[1]，即盡快進(jìn)行第一次中途修正，以后各次修正時(shí)機(jī)的選擇是按照修正效果與前一次修正效果的比值為常數(shù)這一條件確定的。國內(nèi)研究主要圍繞嫦娥探測器展開，文獻(xiàn)[2-3]給出了典型的兩次中途修正策略，即第一次修正在飛行的第一天實(shí)施,第二次修正在到達(dá)月球的前一天實(shí)施。

本文主要針對第二個(gè)問題展開研究。其本質(zhì)上是一個(gè)兩點(diǎn)邊值問題，由于三體動(dòng)力學(xué)的復(fù)雜性，不存在解析解，必須利用數(shù)值方法進(jìn)行求解。一般的思路是：首先根據(jù)邊值約束，采用簡化模型設(shè)計(jì)初值，之后考慮精確攝動(dòng)模型，采用微分修正等數(shù)值方法求解。

對于地月轉(zhuǎn)移軌道中途修正問題，可以將實(shí)際速度作為初值，若探測器實(shí)際速度誤差過大，初值精度變差，修正算法的求解性能降低，因此有必要對初值進(jìn)行改進(jìn)。若修正時(shí)刻探測器的位置在月球段，探測器主要受月球引力的作用，初值設(shè)計(jì)本質(zhì)上是求解二體Lambert問題，文獻(xiàn)[6-8]給出了經(jīng)典的求解方法，若在地球段，需要考慮月球引力的影響，目前設(shè)計(jì)方法有圓錐拼接方法[9-15]，但圓錐曲線拼接方法精度有限，導(dǎo)致修正算法的求解效率降低。

在使用微分修正求解精確解時(shí)，需要計(jì)算狀態(tài)轉(zhuǎn)移矩陣，在常用的數(shù)值方法中[3,16-18]，首先需要推導(dǎo)加速度相對位置速度的偏導(dǎo)數(shù)矩陣，然后同時(shí)積分動(dòng)力學(xué)方程就可以得到狀態(tài)轉(zhuǎn)移矩陣，但在每次迭代計(jì)算中，需要積分42個(gè)方程，導(dǎo)致求解時(shí)間增加。

針對以上研究現(xiàn)狀，本文首先利用偽狀態(tài)理論，結(jié)合Vinti預(yù)報(bào)方法[19]，設(shè)計(jì)高精度的初值，在求解精確解時(shí)，提出了一種基于偽狀態(tài)理論的狀態(tài)轉(zhuǎn)移矩陣解析算法。該算法通過設(shè)計(jì)高精度的初值，降低了求解修正速度增量的難度，而且避免了傳統(tǒng)數(shù)值方法計(jì)算狀態(tài)轉(zhuǎn)移矩陣的復(fù)雜性。

1 問題描述

由于修正時(shí)刻探測器的位置無法改變，只能將速度和轉(zhuǎn)移時(shí)間作為設(shè)計(jì)變量，本文采用轉(zhuǎn)移時(shí)間固定的修正策略，因此初始修正參數(shù)只有速度3個(gè)分量，即：

(1)

式中：v01，v02和v03分別為修正時(shí)刻探測器在地球慣性系下的速度分量。

對于月球探測任務(wù)，通常對到達(dá)月球時(shí)刻的月心距、軌道傾角以及飛行角有嚴(yán)格要求，因此選取目標(biāo)軌道參數(shù)為：

(2)

式中：rp，iL以及γL分別為到達(dá)月球時(shí)刻探測器的月心距、月球固聯(lián)系下的軌道傾角以及相對月球的飛行角。

2 偽狀態(tài)理論

偽狀態(tài)理論是文獻(xiàn)[20]提出的一種地月中心引力作用下的軌道近似計(jì)算方法。

令tI時(shí)刻探測器相對地球的狀態(tài)為ΣI=(RI,VI)，利用偽狀態(tài)理論計(jì)算tk時(shí)刻探測器相對地球的狀態(tài)Σk=(Rk,Vk)和相對月球的狀態(tài)σk=(rk,vk)。定義以tk時(shí)刻月球位置為中心的偽狀態(tài)轉(zhuǎn)換球(Pseudostate Transformation Sphere,PTS)，半徑為RPTS，在本文算例中，其值取為24個(gè)地球半徑。假設(shè)RI在PTS外部，Rk在PTS內(nèi)部，如圖1所示。

圖1 偽狀態(tài)理論Fig.1 The pseudostate theory

具體計(jì)算步驟如下：

(3)

3 初值設(shè)計(jì)

3.1 初值求解

本節(jié)給出了地球慣性系下的中途修正初值的求解方法。令修正時(shí)刻探測器位置R0，到達(dá)月球時(shí)刻tL，轉(zhuǎn)移時(shí)間Δt，迭代變量RS，求解過程如下：

1)計(jì)算tL時(shí)刻月球相對地球的位置速度分別為RM和VM，令RS=RM。

2)利用Lambert算法，求解R0到RS，轉(zhuǎn)移時(shí)間Δt的地球二體軌道，可得到R0和RS處的速度V0和VS。

3)計(jì)算tL時(shí)刻探測器相對月球速度在月球固聯(lián)系下的經(jīng)度aL和緯度δL，求解tL時(shí)刻探測器在月球固聯(lián)系下的升交點(diǎn)經(jīng)度ΩL,

sin(aL-ΩL)=tanδL/taniL

(4)

根據(jù)式(4)，ΩL有兩個(gè)解，對應(yīng)到達(dá)月球的升降軌方式，則角動(dòng)量hL在月球固聯(lián)系下的表達(dá)式為:

hL=[siniLcosΩL-siniLsinΩLcosiL]T

(5)

4)令探測器在PTS邊界時(shí)相對月球的位置和速度分別為rS和vS，對于地月轉(zhuǎn)移軌道，rS在近月點(diǎn)前，大小為RPTS，根據(jù)能量守恒，可得近月點(diǎn)速度大小為：

(6)

(7)

令rS到近月點(diǎn)的轉(zhuǎn)移時(shí)間為tS，根據(jù)偽狀態(tài)理論可得

(8)

3.2 地球扁率影響的修正

對于地月轉(zhuǎn)移軌道，地球扁率攝動(dòng)是一項(xiàng)不可忽視的攝動(dòng)源，因此在精確初值設(shè)計(jì)時(shí)，需要考慮地球扁率的影響。

在求解地球扁率攝動(dòng)的解析方法中，Vinti預(yù)報(bào)方法精度高，速度快，因此采用Vinti預(yù)報(bào)方法來修正地球扁率攝動(dòng)的影響，求解過程如下：

1)求解R0到RS，轉(zhuǎn)移時(shí)間為Δt的地球二體軌道，可以得到R0地球速度V0k。

2)以R0和V0k為初值，利用Vinti預(yù)報(bào)方法計(jì)算tL時(shí)刻的位置RSk。

3)更新V0k，其中Jk可近似為二體狀態(tài)轉(zhuǎn)移矩陣：

(9)

4)若‖RS-Rsk‖大于給定閾值，返回步驟2)，否則，輸出V0k作為初值。

4 精確解求解

為得到精確解，必須在更復(fù)雜的攝動(dòng)環(huán)境下通過數(shù)值積分計(jì)算精確的轉(zhuǎn)移軌道，然后采用微分修正方法對初值進(jìn)行修正。在精確的轉(zhuǎn)移軌道計(jì)算中，本文考慮了地球非球形項(xiàng)，日月中心引力的影響，其中，地球非球形引力取8×8階，月球和太陽星歷均采用DE405。

初始修正參數(shù)P和目標(biāo)軌道參數(shù)Q存在確定關(guān)系，記為Q=g(P)，則存在以下關(guān)系:

(10)

根據(jù)微分修正方法，在每次迭代中初始修正參數(shù)的調(diào)整量為:

(11)

式中：ΔQ為目標(biāo)軌道參數(shù)的偏差量。令狀態(tài)變量X=[R,V]T，R和V分別表示探測器在地球慣性下的位置和速度，則式(10)可以重寫為：

(12)

式中：Xi，Xf分別為狀態(tài)變量X的初始狀態(tài)和終止?fàn)顟B(tài)。式(12)誤差傳遞矩陣計(jì)算中，除狀態(tài)轉(zhuǎn)移矩陣Φfi外均有確定的形式。三體軌道動(dòng)力學(xué)特性導(dǎo)致Φfi沒有解析表達(dá)式，通常采用數(shù)值方法計(jì)算。為提高求解精確解的效率，本節(jié)給出了一種基于偽狀態(tài)理論的狀態(tài)轉(zhuǎn)移矩陣近似解析計(jì)算方法。

根據(jù)偽狀態(tài)理論，存在以下關(guān)系:

(13)

式(13)第一項(xiàng)具體表達(dá)式為:

(14)

(15)

式(13)第二項(xiàng)可以分解為以下兩項(xiàng)：

(16)

(17)

式中：

(18)

根據(jù)關(guān)系式：

(19)

對式(19)求導(dǎo)可得：

(20)

5 算例與分析

本節(jié)給出了一個(gè)地月轉(zhuǎn)移軌道中途修正的求解算例，首先根據(jù)本文給出的初值設(shè)計(jì)方法給出初值，然后利用解析狀態(tài)轉(zhuǎn)移矩陣求解精確解。

假設(shè)探測器在2022年1月1日進(jìn)行修正，此時(shí)探測器在地球慣性坐標(biāo)系下位置為[-4 094.687，4 730.547，2 588.812] km，轉(zhuǎn)移時(shí)間為72 h，近月點(diǎn)高度為600 km，達(dá)到月球時(shí)軌道傾角25°，到達(dá)月球方式為升。

在無奇異情況下，Lambert問題具有長短程兩個(gè)解，所以目標(biāo)軌道參數(shù)確定以后，中途修正有兩個(gè)解，相應(yīng)的軌跡如圖2所示，其對應(yīng)的地心張角不同。經(jīng)過計(jì)算可以得到兩個(gè)可行解。其中：短程解為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 地月轉(zhuǎn)移軌道Fig.2 Trans-lunar trajectory

通常長程解被忽略，但在某些任務(wù)中，如探測器并非直接由運(yùn)載送入地月轉(zhuǎn)移軌道，而是在近地軌道上依靠自身推進(jìn)系統(tǒng)進(jìn)入地月轉(zhuǎn)移軌道，長程解可作為地月轉(zhuǎn)移軌道的選擇，本文算法直接可以設(shè)計(jì)地月轉(zhuǎn)移軌道。

下面通過一個(gè)算例來說明本文算法的求解性能。選擇一條標(biāo)準(zhǔn)地月轉(zhuǎn)移軌道，其發(fā)射時(shí)刻為2022年1月1日，轉(zhuǎn)移時(shí)間120 h。近地點(diǎn)高度為300 km，地球出發(fā)時(shí)軌道傾角為25°，地球出發(fā)方式為降軌，近月點(diǎn)高度為500 km，到達(dá)月球時(shí)軌道傾角為35°，到達(dá)月球方式為升軌，在地球慣性系下，探測器出發(fā)時(shí)刻的標(biāo)稱位置速度分別為[ -5 848.290，2 683.276，1 760.631] km和[-5.135 05，-8.845 35，-3.576 42] km/s。假設(shè)探測器進(jìn)入地月轉(zhuǎn)移軌道時(shí)半長軸偏差為10 000 km，在入軌24 h后進(jìn)行修正。

如果將實(shí)際速度作為修正初值，其搜索過程如表1所示，由于誤差過大，導(dǎo)致迭代發(fā)散，無法得到精確解。若采用精確的初值，經(jīng)過8次迭代就可得到精確解，迭代過程如表2所示。通過初值設(shè)計(jì)，提高了初值精度，保證了修正算法的收斂性和效率。

表1 近月點(diǎn)狀態(tài)偏差迭代過程Table 1 Iteration of the error of perilune

表2 利用精確初值的近月點(diǎn)狀態(tài)偏差迭代過程Table 2 Iteration of the error of perilune using the initial solution with high accuracy

表3給出了相同計(jì)算條件下(Intel Core 2.53G),分別利用解析法計(jì)算狀態(tài)轉(zhuǎn)移矩陣(簡稱解析法)和數(shù)值法計(jì)算狀態(tài)轉(zhuǎn)移矩陣(簡稱數(shù)值法)求解精確解的計(jì)算結(jié)果，盡管兩種方法均可以得到精確解，但若采用解析法求解狀態(tài)轉(zhuǎn)移矩陣，計(jì)算時(shí)間更少，求解效率更高，而且算法結(jié)構(gòu)簡單，避免數(shù)值計(jì)算的復(fù)雜性。

表3 計(jì)算結(jié)果Table 3 Simulation results

6 結(jié)束語

本文針對地月轉(zhuǎn)移軌道中途修正問題，提出了一種求解修正速度增量的方法。該算法通過改善初值的精度和簡化狀態(tài)轉(zhuǎn)移矩陣，提高了地月轉(zhuǎn)移中途修正的求解效率，可有效解決地月轉(zhuǎn)移中途修正問題，可為月地轉(zhuǎn)移等深空探測軌道設(shè)計(jì)以及中途修正提供借鑒和參考。