本立言，嚴玲玲，謝祥華，張銳，*，王國際

（1.中國科學院微小衛(wèi)星創(chuàng)新研究院，上海201203； 2.上海微小衛(wèi)星工程中心，上海201203）

自中國探月工程開展以來，在地月轉移軌道設計、中途修正等方面取得了豐富的研究成果，相關技術已成功應用于工程實踐［1-3］。根據探月工程總體規(guī)劃，中國將繼續(xù)開展月球無人采樣返回任務，軌道設計及軌道控制面臨著新的挑戰(zhàn)和難點。月地轉移并不是地月轉移簡單的逆過程，相關約束條件更加復雜。根據月球采樣返回任務的需求，從月球停泊軌道出發(fā)直接再入大氣的月地轉移軌道是月球采樣返回任務的首選方案。

月地轉移軌道設計不同于地月轉移軌道。首先通過調整發(fā)射窗口，可容易獲取共面的地月轉移軌道，但由于月球停泊軌道面的固定，月地轉移的初始狀態(tài)不一定滿足共面轉移的條件，則需要進行軌道面的調整，并且探測器再入大氣返回地球時，受回收場的約束，需要匹配落點的地理位置，避免軌道面的調整，保證任務可靠實施。

目前，國內外針對月地轉移軌道設計展開了深入的研究。Ocampo和Saudemont［4］研究了利用單脈沖和三脈沖機動方式從環(huán)月軌道出發(fā)進入月地轉移軌道的問題，但未考慮再入點約束。高玉東等［5］基于雙二體模型假設，利用二分法搜索從月球停泊軌道到地球停泊軌道的共面月地轉移軌道，但未涉及直接再入大氣的約束。張磊等［6-7］以近月點速度、月球停泊軌道升交點赤經及近月點幅角為設計參數，以再入點地心距、再入角和再入傾角為終端參數，通過兩級修正設計滿足約束的月地轉移軌道，在此基礎上通過分別調整再入傾角和飛行時間完成對落點地理位置匹配，求解過程復雜。鄭愛武等［8-10］以探測器出月球影響球的時刻，位置和速度為中間變量，將月地轉移軌道分為地心段和月心段，通過調整影響球邊界處的位置和速度獲得一條連續(xù)的月地轉移軌道，在此基礎上給出了精確月地轉移軌道設計方法和月地返回窗口搜索策略，但初值選擇依賴于設計者的經驗且未涉及月球停泊軌道面調整。汪中生等［11］分析了多種月地轉移方案，并給出了相應的基本設計方法，但其求解效率尚有改進和優(yōu)化的空間。

由于三體問題不存在解析解，月地轉移軌道設計必須采用數值方法。一般的思路是：根據邊值約束，采用簡化模型設計初值，之后考慮精確攝動模型，采用微分修正等數值方法求解精確轉移軌道。在目前的研究成果中，偽狀態(tài)理論［12］是一種計算三體問題的近似解析形式，由于考慮地月球中心引力共同作用，預報精度相對較高，廣泛應用于三體軌道的初值設計。

針對以上研究現狀，本文基于偽狀態(tài)理論，以月球停泊軌道為出發(fā)點，考慮再入點的邊界約束條件，通過簡單迭代設計單脈沖月地轉移軌道的初值，在此基礎上利用微分修正方法求解精確轉移軌道，并對月地返回窗口進行搜索和分析。

1 動力學模型

在地球慣性系下，建立探測器動力學模型，模型中考慮了地球非球形項、月球及太陽的影響，其中，地球非球形引力取8×8階，月球和太陽星歷均采用DE405。

式中：r和v分別為探測器的位置和速度矢量；μe、μm和μh分別為地球、月球和太陽的引力常數；rh和rm分別為太陽和月球的位置；rhd和rmd分別為探測器到太陽和月球的位置；R為地球引力攝動位函數。

2 設計約束

由于探測器運行在月球停泊軌道上，機動時刻探測器的位置無法改變，只能將速度和轉移時間作為自由變量，故初始軌道參數有速度和轉移時間。則初始設計參數選取為對于直接再入大氣的月地轉移軌道，通常對再入時刻的地心距、再入角、再入傾角有嚴格要求，為匹配落點的地理位置，再入軌道面與回收場必須共面，可采用回收場位置與再入軌道角動量間夾角來描述，則目標終端參數選取為

式中：Re為再入時刻的地心距；γe為再入時刻的再入角；ie為再入時刻的再入傾角；φe為再入時刻回收場位置與再入軌道角動量的夾角。

3 偽狀態(tài)理論

偽狀態(tài)理論是W ilson［12］提出的一種考慮地月中心引力作用下軌道近似計算方法。

考慮地球、月球及探測器三體系統(tǒng)，認為探測器對地球和月球運動沒有影響。給定探測器tI時刻相對地心的位置速度ΣI=（RI，VI），利用偽狀態(tài)理論可以計算tk時刻相對地心的位置速度Σk=（Rk，Vk）和相對月心的位置速度σk=（rk，vk）。定義以tk時刻月球位置為中心的偽狀態(tài)轉換球（Pseudostate Transformation Sphere，PTS），半徑為RPTS。假設RI在PTS外部，Rk在PTS內部，如圖1所示。

圖1 偽狀態(tài)理論示意圖Fig.1 Schemetic of pseudostate theory

具體計算步驟如下：

3）月心二體段。根據tc時刻相對月心的位置速度σxc，按月心二體軌道求解tk時刻相對月心的位置速度σk，即可得到相對地心的位置速度Σk。

4 初值設計

本節(jié)給出了月地轉移軌道的初值設計方法。若機動時刻探測器的位置在近月點前，選擇近月點對應的偽狀態(tài)位置作為迭代變量，轉移軌道分為兩段處理，轉移軌道在近月點前按月球二體軌道處理。若機動時刻探測器的位置在近月點后，選擇機動時刻探測器對應的偽狀態(tài)位置作為迭代變量。

定義機動時刻t0探測器在月球慣性系下位置為r0，近月點前轉移時間Δt1，近月點后轉移時間Δt2，偽狀態(tài)位置為Rs。

步驟1計算t0時刻月球位置速度Rm和Vm，令Rs=Rm，Δt1=0。

步驟2定義回收場的經度為λl，緯度為φl，位置為Rl，he為地心段角動量，ke為地軸的單位方向。根據余弦定理，可得

式中：α的符號對應了升降軌方式，β的符號對應了軌道面的方位，決定了再入航程大?。ㄒ妶D2）。

圖2 月地轉移軌道的幾何約束Fig.2 Geometric constraints of Moon-to-Earth transfer orbit

定義t0時刻Rs在地球固連系下的經度為Gs，再入時刻Rl在地球固連系下的經度為Gl=Gs+α+β，記G=Gl-λl，則地心段的飛行時間為

式中：ωe為地球自轉角速度；N為指定的天數。

則Rs為返回軌道地球段的起始點，探測器到達再入點時地心距為Re，再入角為γe，采用Luo等［13］提出的飛行路徑角約束的Lambert算法可以求解該類邊值問題，這里不再贅述。

再入點速度為

令Rs處的飛行角為γs，速度大小為Vs，Vs由式（9）計算：

步驟3令探測器在PTS邊界處相對月球的速度為vs，位置為rs，其大小為PTS半徑，在月球無窮遠處的速度大小為

軌道半長軸為

定義Δf為初始位置r0和v∞間的夾角，由幾何關系可知：

式中：f∞為月球無窮遠處真近點角；f0為機動時刻的真近點角。

令em為偏心率，存在以下關系：

為完成月地轉移，機動時刻探測器速度應為

若f0≥0，根據開普勒飛行時間公式計算r0到rs的轉移時間ts，否則轉向步驟5。

步驟5此時機動時刻探測器位置在近月點前，根據開普勒飛行時間公式分別計算r0到rp的轉移時間Δt1和rp到rs的轉移時間ts，同時更新Rm為近月點時刻的月球位置。

步驟6更新Rs。根據偽狀態(tài)理論得到

5 精確解求解

初始設計參數P和目標終端參數Q存在確定的函數關系為

則Q和P的誤差之間的關系可線性近似為

為提高求解算法的速度，將軌道積分停止條件設置為到達給定再入角γe，則初始設計參數P中剔除了轉移時間Δt，目標終端參數Q中解除了再入角γe的約束，則月地轉移軌道設計降價為一個三維的兩點邊值問題。

令X=［r，v］T，則式（23）可以重寫為

當積分停止條件設置為到達給定再入角，則存在以下關系：

式（27）給出了消除時間變量后偏差傳遞矩陣形式：

式中：除矩陣?Xf/?Xi外均有確定的形式，需要注意的是，對目標終端參數cosφe偏導數計算時，需要考慮由于再入時間變化引起的返回場慣性系下位置變化。由于三體軌道動力學特性導致矩陣?Xf/?Xi沒有解析表達式，本文采用文獻［2-3，14-15］提出的方法求解。

得到誤差傳遞矩陣后，可利用各種迭代算法搜索精確解，本文采用微分修正方法求解精確解。

6 算例與分析

本節(jié)首先給出了一個月地轉移軌道設計算例。相關算法均采用C++編程，程序運行環(huán)境：CPU為Intel Core 2.53G。

假設探測器在2022年1月1日進行單脈沖機動，此時探測器在月球慣性系下位置為［1 937.4，0，0］km，轉移時間為3～4 d，再入點地心高度為120 km，再入傾角為45°，再入角為-6°，返回場的經度為110°，緯度為40°。

根據第4節(jié)的初值設計方法，存在4條月地轉移軌道滿足給定約束條件，下面以其中一條為例，表1給出了偽狀態(tài)位置誤差的3個分量迭代過程。在初值設計時，迭代6次，耗時1 ms，在精確值求解時，迭代4次，耗時3 s，即可滿足設定精度要求，表2給出了再入點狀態(tài)偏差的迭代過程。比較初值［-337.339，2 218.060，914.856］m/s與精確解［-346.252，2 219.854，914.926］m/s，兩者模值的差為2.925m/s，兩者矢量差的模值為9.092m/s，為初始速度大小的0.37%，在任務分析階段，可用初值代替精確解，為任務分析提供參考。

表1 初值設計時偽狀態(tài)位置誤差迭代過程Table 1 Iteration process of pseudostate position error during initial solution design

表2 精確解求解時再入點狀態(tài)偏差迭代過程Table 2 Iteration p rocess of reentry point state error during exact solution search

表3給出滿足給定約束條件的4條月地轉移軌道的關鍵特征參數，圖3和圖4給出相應的3D和2D示意圖。分析結果可知，4條月地轉移軌道通過調整轉移時間滿足了再入軌道與返回場的共面要求，再入點速度大小基本相同，再入航程因軌道方位不同而出現差別。

下面計算2022年1月1日到2022年1月6日從月球停泊軌道出發(fā)的月地返回窗口，其中，月球停泊軌道在2022年1月1日的具體參數：高度為200 km，偏心率為0，傾角為22.4°，升交點赤經為0°，緯度幅角為0°，目標軌道參數與上述算例相同。其中2022年1月1日滿足共面轉移條件，這是月地轉移能量最省方式之一，以此為基準分析速度增量變化趨勢。由于月球停泊軌道周期較小，角速度較快，速度方向變化較大，在每個軌道周期內變軌速度增量變化劇烈，如圖5所示，在每個軌道周期內，存在一緯度幅角使變軌速度增量最小。

表3 月地轉移軌道特征參數Table 3 Characteristic param eters of M oon-to-Earth transfer orbit

圖3 月地轉移軌道3D示意圖Fig.3 Three-dimensional illustration of Moon-to-Earth transfer orbit

圖4 月地轉移軌道2D示意圖Fig.4 Two-dimensional illustration of Moon-to-Earth transfer orbit

在窗口搜索時選取每個軌道周期內最小變軌速度增量，分析其變化趨勢。由圖6可知，結果有明顯的周期性。在每一天內變軌速度增量逐漸增大，變軌速度增量差別最大為40m/s。隨著探測器出發(fā)時刻延遲，月球停泊軌道與月地轉移軌道面夾角增大，如圖7所示。若實現月地轉移，需要調整軌道平面，變軌速度增量逐漸增大，再入航程同時逐漸增大，若月球出發(fā)時間延遲6 d，變軌速度增量增加130m/s。

圖5 一個軌道周期內變軌速度增量變化趨勢Fig.5 Track change velocity increment versus time in one orbit period

圖6 變軌速度增量變化趨勢Fig.6 Track change velocity increment versus time

圖7 軌道平面改變量變化趨勢Fig.7 Change of orbit plane versus time

通過對月地返回窗口搜索可知，單脈沖轉移方式通過調整軌道面，增加了月地返回機會，保證任務順利進行。但隨著對月球探測的繼續(xù)深入，月球探測任務更加復雜，對探測器中止任務能力的需求相應增加，顯然單脈沖轉移方式在軌道平面調整量較大時，需要消耗更多能量，無法滿足任務的要求，因此對多脈沖月地轉移軌道的研究是后續(xù)工作的重點。

7 結 論

本文給出了從月球停泊軌道出發(fā)直接再入大氣的月地轉移軌道設計方法。

1）初值設計算法精度高，耗時少，在任務分析階段可作為代替精確解，為任務分析提供參考依據。精確解求解算法求解效率高，可以作為設計月地轉移軌道及搜索月地返回窗口的一個有力工具。

2）給出了滿足給定約束的4條月地轉移軌道，并對比了關鍵特征參數，為月地轉移軌道的選擇提供依據。

3）通過月地返回窗口分析可知，由于匹配返回場，每天內變軌速度增量逐漸增大。當月球停泊軌道與月地轉移軌道面夾角逐漸增大時，使用單脈沖機動，則軌道面調整量增大，導致燃料消耗增大，限制了月地轉移機會。