童科偉，程炳琳，汪小衛(wèi)

(1. 中國(guó)運(yùn)載火箭技術(shù)研究院 研究發(fā)展部， 北京 100076; 2. 中國(guó)運(yùn)載火箭技術(shù)研究院， 北京 100076)

雙二體軌道拼接設(shè)計(jì)方法在月球探測(cè)軌道設(shè)計(jì)中具有基礎(chǔ)性作用，國(guó)內(nèi)外學(xué)者以雙二體軌道拼接法為基礎(chǔ)對(duì)月球探測(cè)軌道設(shè)計(jì)方法進(jìn)行了大量的研究[1-8]。

在基于雙二體問(wèn)題求解地月轉(zhuǎn)移軌道時(shí)，Battin、郗曉寧、黃文德、賀波勇等都不約而同地以入口點(diǎn)月心經(jīng)緯度作為中間變量來(lái)拼接地球影響球內(nèi)飛行軌道和月球影響球內(nèi)飛行軌道[1-2,4-5,8]。這種基于月球影響球球坐標(biāo)來(lái)描述入口點(diǎn)位置的方法比較直觀，但是沒(méi)有充分利用轉(zhuǎn)移軌道的物理意義。

對(duì)于著陸地點(diǎn)已知的奔月任務(wù)，賀波勇推導(dǎo)了一種類解析的改進(jìn)雙二體模型，該方法與傳統(tǒng)的從地球計(jì)算到月球不同，在計(jì)算窗口時(shí)針對(duì)已知的月心軌道六根數(shù)逆向計(jì)算地月轉(zhuǎn)移軌道六根數(shù)，計(jì)算量極小[9]。

對(duì)于工程約束不明、具體探測(cè)地點(diǎn)不定的地月轉(zhuǎn)移軌道窗口計(jì)算問(wèn)題，還需要按照傳統(tǒng)的計(jì)算量較大的搜索方法，進(jìn)而得到一些運(yùn)動(dòng)特性。與常規(guī)的以入口點(diǎn)月心經(jīng)緯度作為中間變量的雙二體拼接方法不同，本文提出一種直接基于軌道特征的地月轉(zhuǎn)移飛行軌道的幾何方法，物理意義明確，求解地月軌道更加直觀，通過(guò)兩次降維解耦操作，降低問(wèn)題求解維數(shù)，進(jìn)一步提升計(jì)算效率。

考慮到目前大多數(shù)計(jì)算機(jī)都是多核中央處理器(central processing unit, CPU)，最簡(jiǎn)單的提高算力的途徑是采用并行計(jì)算，充分利用各核的計(jì)算能力，并行計(jì)算技術(shù)在國(guó)內(nèi)外受到高度重視。并行程序設(shè)計(jì)編程語(yǔ)言通常基于消息傳遞接口標(biāo)準(zhǔn)(message passing interface, MPI)、OpenMP等，二者甚至已經(jīng)分別形成了集群并行和單機(jī)多核并行領(lǐng)域的事實(shí)標(biāo)準(zhǔn)[10-13]。

本文以Java相關(guān)的多核并行編程技術(shù)為基礎(chǔ)[13]，設(shè)計(jì)了一種地月轉(zhuǎn)移軌道的多核并行計(jì)算方法，并應(yīng)用于天梯地月轉(zhuǎn)移軌道計(jì)算與特性分析。按照國(guó)際天梯協(xié)會(huì)(international space elevator consortium, ISEC)的定義，天梯是一種把有效載荷從地球表面提升到空間的系統(tǒng)，該系統(tǒng)是一個(gè)長(zhǎng)10萬(wàn)千米的繩系，質(zhì)心在地球靜止軌道(geostationary orbit, GEO)，并固定于地表的某個(gè)錨點(diǎn)。隨著爬升高度的變化，天梯攀爬器在電能作用下抬升高度并獲得地球自轉(zhuǎn)帶來(lái)的圓速度，從而使有效載荷獲得勢(shì)能和動(dòng)能[14]。由于天梯位于赤道上空，軌道傾角受限，與常規(guī)月球轉(zhuǎn)移軌道設(shè)計(jì)有一定區(qū)別。本文以一種幾何切面法雙二體模型來(lái)分析天梯發(fā)射航天器進(jìn)入月球轉(zhuǎn)移軌道的特性。

1 雙二體幾何切面法

雙二體假設(shè)下的地月轉(zhuǎn)移軌道計(jì)算也稱為圓錐曲線拼接法。本質(zhì)上是以月球影響球?yàn)檫吔? 將軌道分成多段來(lái)拼接：進(jìn)入影響球前, 飛行軌跡為地心圓錐曲線(一般為橢圓)；在影響球內(nèi)，飛行軌跡為月心圓錐曲線(一般為雙曲線)。對(duì)于自由返回軌道等任務(wù)，還需考慮飛出影響球后的飛行軌跡，其為地心圓錐曲線(一般也為橢圓)。在進(jìn)入和離開(kāi)影響球的兩個(gè)邊界點(diǎn)(一般稱為入口點(diǎn)和出口點(diǎn))，將這幾段圓錐曲線拼接成完整的運(yùn)動(dòng)軌跡。

圓錐曲線拼接法的設(shè)計(jì)參數(shù)通常以月球影響球入口點(diǎn)/出口點(diǎn)在月心白道坐標(biāo)系下的經(jīng)緯度作為獨(dú)立變量，通過(guò)該經(jīng)緯度來(lái)描述入口點(diǎn)/出口點(diǎn)的空間位置[1-2,4-6,8]。

采用傳統(tǒng)的圓錐曲線拼接法的設(shè)計(jì)參數(shù)計(jì)算轉(zhuǎn)移窗口時(shí)問(wèn)題待求解參數(shù)的維度大，計(jì)算轉(zhuǎn)移窗口時(shí)的計(jì)算時(shí)間過(guò)長(zhǎng)。對(duì)于在天梯系統(tǒng)出發(fā)這一類地月轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題，由于出發(fā)軌道的傾角已知、升交點(diǎn)赤經(jīng)(right ascension of ascending node, RAAN)也可求出，本文借助空間幾何關(guān)系把原始問(wèn)題降維，從而大大降低求解問(wèn)題的規(guī)模。其基本思想是：基于轉(zhuǎn)移軌道面和月球影響球的幾何關(guān)系來(lái)描述空間位置關(guān)系，地月/月地轉(zhuǎn)移軌道面與月球影響球必須相交(否則一定不可能被月球捕獲)。通過(guò)轉(zhuǎn)移軌道面切割月球影響球形成的幾何關(guān)系就可以描述地月/月地轉(zhuǎn)移軌道，結(jié)合Lambert原理就可以求解出所需的轉(zhuǎn)移軌道參數(shù)。這樣就把原來(lái)的三維球面搜索算法降維成了二維平面內(nèi)切割月球影響球形成的圓上的搜索算法。另外本文的描述方式把地月轉(zhuǎn)移軌道求解問(wèn)題解耦成平面內(nèi)轉(zhuǎn)移軌道形狀計(jì)算問(wèn)題以及轉(zhuǎn)移軌道空間定向問(wèn)題，從而大大降低了問(wèn)題求解維度。

1.1 轉(zhuǎn)移三角形形狀計(jì)算問(wèn)題

對(duì)于從天梯出發(fā)的奔月任務(wù)來(lái)說(shuō)，為了擴(kuò)大任務(wù)窗口，可以適當(dāng)允許出發(fā)轉(zhuǎn)移軌道有小傾角，以盡可能節(jié)約推進(jìn)劑，比如要求i∈[-0.5°,0.5°]。由于天梯質(zhì)心位于地球定點(diǎn)位置，相應(yīng)的升交點(diǎn)赤經(jīng)為：

Ω=Ω0+ωEt

其中，Ω0為初始點(diǎn)降交點(diǎn)赤經(jīng)，ωE為地球自轉(zhuǎn)角速度，t為據(jù)初始點(diǎn)時(shí)刻。

如圖1所示，定義地心為E，月心為L(zhǎng)，軌道面與月球影響球切割形成圓的圓心為O。轉(zhuǎn)移軌道在地球近地點(diǎn)出發(fā)，該點(diǎn)定義為P，到達(dá)月球影響球的入口點(diǎn)A。

圖1 軌道面與月球影響球的切割示意圖Fig.1 Demonstration of section of orbital plane and the sphere of influence of moon

轉(zhuǎn)移軌道面與月球影響球切割的交線必然形成一個(gè)圓，該圓上各點(diǎn)到月心的距離等于影響球半徑rs，rs=66 200 km。轉(zhuǎn)移軌道面法方向或軌道動(dòng)量矩方向nOL為：

nOL=[sinisinΩ,-sinicosΩ,cosi]T

(1)

cos∠OLE|nOL·rL|/rL

(2)

其中，rL為地月距離。

由圖(1)所示的幾何關(guān)系可知，順行軌道(i,Ω)和逆行軌道(π-i,Ω+π)都滿足以上幾何關(guān)系(二者軌道面法向相差180°)，相應(yīng)的式(2)取絕對(duì)值，可根據(jù)工程實(shí)際約束來(lái)選取。

l=rLcos∠OLE

(3)

地月轉(zhuǎn)移軌道面與月球影響球必須相交，否則轉(zhuǎn)移軌道與月球影響球沒(méi)有交點(diǎn)，就無(wú)法形成環(huán)月軌道。對(duì)應(yīng)的數(shù)學(xué)約束條件就是：

l≤rs

(4)

發(fā)射窗口搜索時(shí)間區(qū)間出現(xiàn)在滿足式(4)的時(shí)間段，該約束條件可以極大地節(jié)約轉(zhuǎn)移窗口的計(jì)算量。

軌道面與影響球切割而成的圓半徑OA的大小r0為：

(5)

(6)

k-r0≤rEA≤k+r0

(7)

圖2 軌道面與月球影響球的截面平面示意圖Fig.2 Section plane of orbital plane and the sphere of influence of moon

由地月幾何關(guān)系，rEA取值范圍可擴(kuò)大為：

rL-rs≤k-r0≤rEA≤k+r0≤rL+rs

(8)

已知地月平均距離為384 400 km，因此地月距離rL考慮一定余量就可以得到rEA的取值范圍，如取余量Δ=0.1可得：

(rL-rs)(1-Δ)≤rEA≤(rL+rs)(1+Δ)

即

286 380 km≤rEA≤495 660 km

(9)

在給定時(shí)間tPA內(nèi)，從近地出發(fā)點(diǎn)P運(yùn)行到入口點(diǎn)A的地月轉(zhuǎn)移軌道構(gòu)成了一個(gè)Lambert問(wèn)題。根據(jù)Lambert定理，地月轉(zhuǎn)移出發(fā)點(diǎn)P到入口點(diǎn)A的飛行時(shí)間tPA僅與軌道半長(zhǎng)軸、兩點(diǎn)與地心距離之和rPE+rEA以及弦長(zhǎng)rPA有關(guān)。又由余弦定理，給定rPE和rEA，則弦長(zhǎng)rPA和真近點(diǎn)角θPA一一對(duì)應(yīng)。因此地月轉(zhuǎn)移軌道問(wèn)題轉(zhuǎn)化為已知tPA求θPA的問(wèn)題。E、P、A確定了一個(gè)與轉(zhuǎn)移軌道相關(guān)的三角形，稱為轉(zhuǎn)移三角形。

賀波勇推導(dǎo)了tPA和θPA的導(dǎo)數(shù)關(guān)系。將tPA作為設(shè)計(jì)變量，用牛頓迭代法求解θPA[9]。本文的算例測(cè)試發(fā)現(xiàn)，該方法對(duì)初值較敏感，收斂性一般，計(jì)算量較大，且計(jì)算精度也不太高，為此改用Brent法求解。

Brent法在20世紀(jì)60年代，由阿姆斯特丹數(shù)學(xué)中心的Van Wijngaarden和Dekker等研究成功，并在1973年由Brent進(jìn)行了改進(jìn)。該方法特別適合于求解一維非線性方程的根，而無(wú)須提供函數(shù)的導(dǎo)數(shù)。該方法具有超線性收斂特性，又能夠保證二分法的收斂確定性。具體算法參見(jiàn)文獻(xiàn)[15]。

按照Brent算法，選擇tPA作為設(shè)計(jì)變量，求解θPA，二者的求解關(guān)系如下：

地月轉(zhuǎn)移軌道由近地點(diǎn)出發(fā)，存在如下關(guān)系：

rEA=p/(1+ecosθPA)

p=rEP(1+e)

可得到：

(10)

(11)

其中：p為半通徑；e為偏心率；a為半長(zhǎng)軸。

至此由θPA、a、e不難求得轉(zhuǎn)移軌道飛行時(shí)間tPA，可以將求得的tPA與給定的轉(zhuǎn)移軌道飛行時(shí)間tPA0之差作為待求的一維非線性目標(biāo)函數(shù)，對(duì)應(yīng)的未知量θPA可由Brent法求解。

轉(zhuǎn)移三角形形狀求解問(wèn)題在整個(gè)搜索過(guò)程中重復(fù)出現(xiàn)，但由本文模型的描述方式可知，該三角形與軌道傾角、升交點(diǎn)赤經(jīng)和近地點(diǎn)幅角無(wú)關(guān)，因此該求解問(wèn)題與整體窗口搜索算法解耦，可以把轉(zhuǎn)移三角形所有可能的軌道組合提前求解出來(lái)，在接下來(lái)的窗口搜索中就只需要匹配這些已知的轉(zhuǎn)移三角形，從而大大降低問(wèn)題的搜索維度，節(jié)省計(jì)算量。此外轉(zhuǎn)移三角形形狀求解問(wèn)題也從原來(lái)的三維影響球的球面經(jīng)緯度搜索問(wèn)題降維成為平面內(nèi)圓上點(diǎn)的搜索問(wèn)題，進(jìn)一步節(jié)省了計(jì)算量。

1.2 轉(zhuǎn)移三角形的空間定向問(wèn)題

1)對(duì)于如近地軌道空間站出發(fā)的地月轉(zhuǎn)移軌道，i和初始Ω0一般已知，則當(dāng)時(shí)的Ω可由到達(dá)入口點(diǎn)時(shí)刻和飛行時(shí)間來(lái)決定，未知量為近地點(diǎn)幅角ω；

2)對(duì)于天梯奔月任務(wù)來(lái)說(shuō)，天梯位于赤道上方，ω等于0°(順行軌道)或者180°(逆行軌道)，未知量為i。

以上兩種情況均可由Brent算法求解。

由地心慣性系到軌道系的坐標(biāo)旋轉(zhuǎn)矩陣為：

其中：Lx(?)表示繞x軸旋轉(zhuǎn)?角度的坐標(biāo)旋轉(zhuǎn)矩陣；Lz(?)表示繞z軸旋轉(zhuǎn)?角度的坐標(biāo)旋轉(zhuǎn)矩陣；u=ω+θ為緯度幅角，其中θ為真近點(diǎn)角。

地心慣性系下的入口點(diǎn)位置矢量rEA為：

(12)

由圖1可知，rEA與月球位置矢量rL(由星歷計(jì)算)及入口點(diǎn)到月心的矢量rLA滿足如下幾何關(guān)系：

rLA=rEA-rL

而rLA的大小等于月球影響球半徑rs，則：

|rEA-rL|-rs=0

(13)

因此根據(jù)未知的ω或者i，由Brent算法預(yù)先判斷所求區(qū)間是否有解，再求解式(13)對(duì)應(yīng)的一維非線性方程的根，就可以得到滿足條件的未知的ω或者i[15]。即：

1)對(duì)于i給定的轉(zhuǎn)移軌道，ω∈[0°,360°]，由式(12)和式(13)，根據(jù)已知的五個(gè)軌道根數(shù)，由Brent算法求解ω；

2)對(duì)于天梯任務(wù)，i∈[-0.5°,0°]和i∈[0°,0.5°]對(duì)應(yīng)的逆行和順行軌道，同樣由式(12)和式(13)，根據(jù)已知的其他五個(gè)軌道根數(shù)，由Brent算法求解i。

然后根據(jù)式(1)～(6)計(jì)算轉(zhuǎn)移三角形數(shù)據(jù)庫(kù)中待評(píng)估的轉(zhuǎn)移三角形是否滿足式(7)。

以上求解軌道面參數(shù)的過(guò)程稱為轉(zhuǎn)移三角形的空間定向問(wèn)題。

以上求解出了地月轉(zhuǎn)移段軌道六個(gè)軌道根數(shù)，還需要進(jìn)一步判斷該軌道進(jìn)入月球影響球被月球捕獲以后形成的環(huán)月軌道是否滿足給定條件。此時(shí)需把地心慣性系下的入口點(diǎn)處的矢量rLA和vLA轉(zhuǎn)換到月心白道系下[2]：

M=Lz(uL+π)Lx(iL)Lz(ΩL)

其中：M為地心慣性系到月心瞬時(shí)白道系的旋轉(zhuǎn)矩陣。上式中的上標(biāo)L表示月球相關(guān)的參數(shù)。

基于上式求得的月心系下的位置和速度，就可以求得月心段繞月軌道參數(shù)，之后可以判斷是否滿足需要的環(huán)月軌道條件，比如環(huán)月軌道高度約束等。

以上是普通地月轉(zhuǎn)移軌道的求解流程。對(duì)于載人航天任務(wù)常用的自由返回軌道，可基于對(duì)稱性原理，進(jìn)一步求得月地轉(zhuǎn)移段軌道參數(shù)，再判斷是否滿足載人航天任務(wù)常用的自由返回軌道約束條件，如返回地球時(shí)的近地點(diǎn)高度、再入角等。本文只計(jì)算一般的地月轉(zhuǎn)移軌道，自由返回軌道求解的具體流程可參考Battin[1]、黃文德[4-5]、賀波勇[8]等的報(bào)道。

2 天梯奔月解法

2.1 CPU并行計(jì)算技術(shù)簡(jiǎn)介

當(dāng)前CPU主頻的提升已經(jīng)明顯遇到了瓶頸，計(jì)算機(jī)硬件界主要關(guān)注如何發(fā)展多核CPU。目前民用CPU已經(jīng)擁有十多個(gè)核心。針對(duì)CPU核心數(shù)越來(lái)越多的現(xiàn)狀，科學(xué)計(jì)算必須考慮并行/并發(fā)程序設(shè)計(jì)，避免浪費(fèi)計(jì)算資源。

并行計(jì)算一般可分為計(jì)算密集型、數(shù)據(jù)密集型、網(wǎng)絡(luò)密集型[10]，本文主要關(guān)注計(jì)算密集型計(jì)算。1994年誕生的MPI已經(jīng)成為集群并行程序設(shè)計(jì)事實(shí)上的標(biāo)準(zhǔn)。1997年誕生的OpenMP已經(jīng)成為單機(jī)多核并行事實(shí)上的標(biāo)準(zhǔn)[10]。本文的并行程序主要針對(duì)單機(jī)多核并行架構(gòu)，而單機(jī)多核圖形處理器(graphics processing unit, GPU)加速并行架構(gòu)和多機(jī)多核并行架構(gòu)不在本文考慮。常用的并行分解技術(shù)有兩種：遞歸分解和數(shù)據(jù)分解[11-12]。本文著重關(guān)注單機(jī)CPU并行程序設(shè)計(jì)，并主要采用數(shù)據(jù)分解的方式來(lái)設(shè)計(jì)并行程序。

本文綜合考慮到可移植性、編程方便性、程序穩(wěn)健性等，避免采用復(fù)雜的MPI，而是采用Java語(yǔ)言，并行計(jì)算程序設(shè)計(jì)基于Parallel Java 2 Library(PJ2庫(kù))[13]。

PJ2庫(kù)是RIT(Rochester Institute of Technology)的Kaminsky教授開(kāi)發(fā)的Java語(yǔ)言的并行計(jì)算調(diào)度庫(kù)，并應(yīng)用于RIT計(jì)算機(jī)學(xué)院的超級(jí)計(jì)算機(jī)集群上。其主要設(shè)計(jì)思想是簡(jiǎn)化并行計(jì)算程序設(shè)計(jì)的難度，特別是避免使用復(fù)雜的MPI，能夠極好地適用于單機(jī)CPU、GPU并行以及分布式集群異構(gòu)并行計(jì)算。PJ2庫(kù)有五種CPU核心并行調(diào)度機(jī)制：Fixed、Leapfrog、Dynamic、Proportional、Guided。后四種調(diào)度機(jī)制代表了當(dāng)前常見(jiàn)的先進(jìn)并行調(diào)度機(jī)制，具體介紹參見(jiàn)文獻(xiàn)[13]。本文選擇Dynamic調(diào)度機(jī)制作為示例。

2.2 天梯地月轉(zhuǎn)移軌道窗口計(jì)算并行設(shè)計(jì)

天梯系統(tǒng)的基本組成如圖3所示，主要包括天梯繩索、地球表面節(jié)點(diǎn)、GEO節(jié)點(diǎn)、頂點(diǎn)錨、繩索攀爬器、地球表面運(yùn)輸系統(tǒng)和運(yùn)營(yíng)中心[14]。

圖3 天梯系統(tǒng)組成示意圖[15]Fig.3 Composition of space elevator system[15]

天梯繩索是天梯最重要的組成部分，其主要功能為向攀爬器提供依附的途徑，使攀爬器能夠沿著纜繩從地表節(jié)點(diǎn)進(jìn)入空間，同時(shí)其還要平衡整個(gè)天梯的重力與系統(tǒng)繞地球轉(zhuǎn)動(dòng)形成的離心力；攀爬器為運(yùn)輸有效載荷進(jìn)入空間的運(yùn)載器，其沿著繩索爬升，進(jìn)入軌道預(yù)定位置后將有效載荷釋放；地球表面節(jié)點(diǎn)和頂點(diǎn)錨分別位于繩索的兩端，起到固定繩索位置、調(diào)整繩索姿態(tài)的功能。整個(gè)天梯的質(zhì)心設(shè)置在GEO節(jié)點(diǎn)位置，這樣整個(gè)天梯系統(tǒng)就以地球的自轉(zhuǎn)速度圍繞地球旋轉(zhuǎn)，從而保證天梯與地面相對(duì)靜止[14,16]。

由于天梯從地表一直延伸到GEO更遠(yuǎn)端，在執(zhí)行行星際發(fā)射任務(wù)時(shí)，在天梯不同的高度釋放可以獲得不同的能量。比如，在天梯遠(yuǎn)端頂點(diǎn)錨處(距離地心100 000 km)釋放甚至無(wú)須施加速度增量就可以飛往火星[16]。本文設(shè)計(jì)的從天梯飛往月球的轉(zhuǎn)移軌道算法，重點(diǎn)解決在天梯哪一高度釋放才能以最小的代價(jià)飛往月球。

本文假定天梯初始位于升交點(diǎn)赤經(jīng)60°赤道上空。轉(zhuǎn)移軌道約束地月轉(zhuǎn)移軌道真近點(diǎn)角差θPA<π，近月點(diǎn)軌道高度為50 km≤rH≤1 000 km。考慮到2025年3月白道面與赤道面相對(duì)夾角達(dá)到最大值，約為28.72°，本文窗口計(jì)算的時(shí)間段取為2025年3月。

本文采用基于數(shù)據(jù)分解的并行方式，主要設(shè)計(jì)兩級(jí)并行算法：轉(zhuǎn)移三角形形狀并行計(jì)算問(wèn)題和轉(zhuǎn)移三角形的空間定向并行計(jì)算問(wèn)題。具體流程如圖4所示。

圖4 軌道設(shè)計(jì)流程Fig.4 Flow chart of trajectory design

第一個(gè)問(wèn)題是轉(zhuǎn)移三角形形狀計(jì)算問(wèn)題并行化：在初估轉(zhuǎn)移軌道特性時(shí)，rEA根據(jù)式(9)均勻選擇若干個(gè)點(diǎn)作為并行化設(shè)計(jì)參數(shù)，從近地點(diǎn)出發(fā)飛行到入口點(diǎn)的飛行時(shí)間tPA取30～170 h，基于Brent算法求解轉(zhuǎn)移三角形的真近點(diǎn)角，求得的結(jié)果存儲(chǔ)下來(lái)用于下一個(gè)子問(wèn)題的搜索。

第二個(gè)問(wèn)題是轉(zhuǎn)移三角形的空間定向問(wèn)題并行化：天梯地月轉(zhuǎn)移軌道窗口計(jì)算時(shí)，以到達(dá)入口點(diǎn)時(shí)刻tA作為劃分并行顆粒度的設(shè)計(jì)變量，計(jì)算1個(gè)月的時(shí)間周期。以到達(dá)入口點(diǎn)時(shí)刻作為并行顆粒度劃分的主要好處是可以在并行計(jì)算之前計(jì)算好tA時(shí)刻的月球星歷。在計(jì)算空間定向時(shí)，軌道傾角i∈[-0.5°,0.5°]。在前一問(wèn)題中求得的滿足約束的轉(zhuǎn)移三角形數(shù)據(jù)庫(kù)中，求解滿足式(13)的Brent子問(wèn)題，從而求解出對(duì)應(yīng)的轉(zhuǎn)移軌道空間定向參數(shù)軌道傾角。然后判斷所求得的軌道是否滿足地月轉(zhuǎn)移軌道條件，如果滿足則找到一條可行的轉(zhuǎn)移軌道。

2.3 天梯奔月窗口計(jì)算結(jié)果

下面分別計(jì)算GEO點(diǎn)釋放和天梯遠(yuǎn)端頂點(diǎn)錨處釋放奔月的計(jì)算結(jié)果，分別如圖5和圖6所示。本文的算例中，采用一臺(tái)Windows 10系統(tǒng)的電腦(AMD 3700X八核十六線程CPU，主頻3.6 GHz，內(nèi)存8G)，rEA取500個(gè)點(diǎn)，tPA間隔取0.2 h，tA間隔取10 min，并行計(jì)算在天梯GEO點(diǎn)釋放時(shí)1個(gè)月的奔月窗口需要437.57 s，串行計(jì)算需要3 040.32 s，并行加速比為6.95。

由圖5和圖6的計(jì)算結(jié)果表明，從天梯奔月每月內(nèi)存在兩次窗口，每次窗口持續(xù)4～6天。對(duì)于奔月任務(wù)，天梯頂點(diǎn)錨處釋放所需的速度增量反而比在GEO節(jié)點(diǎn)直接奔月所需的速度增量大，這主要是由于在遠(yuǎn)端頂點(diǎn)錨處釋放的地月轉(zhuǎn)移軌道獲得了過(guò)多的能量，此能量甚至能飛行到火星，但不適合用于月球軌道轉(zhuǎn)移。為了詳細(xì)計(jì)算窗口，還需要把rEA、tPA和tA的間隔取得更密集一些，但這里用于了解轉(zhuǎn)移軌道的特性已經(jīng)足夠。

(a) 飛行時(shí)間結(jié)果(a) Results of flight time

(a) 飛行時(shí)間結(jié)果(a) Results of flight time

為尋找到底在何處釋放可以以較小代價(jià)飛往月球，計(jì)算了天梯不同軌道高度釋放的奔月窗口，計(jì)算結(jié)果表明：在天梯距地心51 000 km處釋放進(jìn)入地月轉(zhuǎn)移軌道，所需要的速度增量較小，所需速度增量如圖7所示。在此處釋放飛行器獲得的相對(duì)地心慣性系的線速度為3.719 km/s。

(a) 飛行時(shí)間結(jié)果(a) Results of flight time

為詳細(xì)研究從地心距51 000 km處釋放的地月轉(zhuǎn)移軌道特性，rEA取2 000個(gè)點(diǎn)，tPA間隔取0.1 h，tA間隔取1 min，并行計(jì)算天梯在GEO點(diǎn)釋放時(shí)1個(gè)月的奔月窗口需要30 220.12 s(合503.7 min)。

如圖8所示，在距地心51 000 km處的天梯奔月環(huán)月軌道傾角分布在10°～170°之間。

圖8 月球軌道傾角與升交點(diǎn)赤經(jīng)Fig.8 Inclination and RAAN of Moon orbit

天梯地月轉(zhuǎn)移軌道的月球影響球入口點(diǎn)如圖9所示，在月球影響球上西半球，在南、北緯度均有分布。

圖9 入口點(diǎn)經(jīng)緯度分布Fig.9 Latitude and longitude of entry point

以尋找到的速度增量最小的轉(zhuǎn)移軌道為例，輸入STK進(jìn)行驗(yàn)證，計(jì)算得到的轉(zhuǎn)移軌道根數(shù)如表1所示，到入口點(diǎn)飛行時(shí)間111.9 h，入口點(diǎn)對(duì)應(yīng)歷書時(shí)儒略日(對(duì)應(yīng)STK的JED時(shí)間)為2 460 747.973 611 11，求得的最小速度增量為0.855 m/s，近月點(diǎn)高度為517.040 km。

表1 地月轉(zhuǎn)移軌道入口點(diǎn)根數(shù)

基于商業(yè)軟件STK的Astrogator模塊進(jìn)行驗(yàn)證，STK/Astrogator在地球段和月球段軌道預(yù)報(bào)器分別選地球點(diǎn)質(zhì)量和月球點(diǎn)質(zhì)量模型，STK計(jì)算得到的近月點(diǎn)高度為517.645 km。本文計(jì)算結(jié)果與之接近，顯示了本文算法的正確性。設(shè)計(jì)結(jié)果如圖10所示。

圖10 天梯地月轉(zhuǎn)移軌道Fig.9 Earth-Moon transfer orbit from space elevator

3 結(jié)論

本文提出一種基于飛行軌道面參數(shù)描述地月轉(zhuǎn)移軌道的雙二體模型表達(dá)方式，避免采用傳統(tǒng)的基于月球影響球入口點(diǎn)經(jīng)緯度的描述方式，該表達(dá)方式由軌道定義出發(fā)，通過(guò)地月幾何關(guān)系來(lái)表述，更加明確直觀。通過(guò)表述方式的轉(zhuǎn)換把原來(lái)的三維球面搜索算法降維成為二維平面內(nèi)圓上的搜索算法，結(jié)合一維非線性方程求根算法Brent算法和Lambert原理，可以高效地求解地月轉(zhuǎn)移軌道的形狀參數(shù)。

通過(guò)合理的解耦分解，將原始的軌道窗口多變量搜索問(wèn)題進(jìn)一步分解為兩個(gè)子問(wèn)題：一個(gè)為轉(zhuǎn)移三角形形狀計(jì)算問(wèn)題，通過(guò)轉(zhuǎn)移軌道形成的三角形幾何關(guān)系，結(jié)合Lambert原理求解出轉(zhuǎn)移軌道的形狀參數(shù)(半長(zhǎng)軸、偏心率及真近點(diǎn)角)；另一個(gè)為轉(zhuǎn)移三角形的空間定向問(wèn)題，基于已求得的轉(zhuǎn)移軌道形狀數(shù)據(jù)匹配搜索可行的轉(zhuǎn)移軌道，確定剩余的三個(gè)軌道面相關(guān)參數(shù)。最后為了加速窗口搜索過(guò)程，充分發(fā)揮多核計(jì)算機(jī)算力，形成一種轉(zhuǎn)移三角形形狀計(jì)算問(wèn)題和空間定向問(wèn)題的兩級(jí)并行算法。

基于本文的并行圓錐曲線幾何切面法，計(jì)算了天梯不同高度進(jìn)行奔月任務(wù)的軌道窗口，結(jié)果表明：

1)在天梯遠(yuǎn)端釋放執(zhí)行奔月任務(wù)效果甚至不如在GEO點(diǎn)直接釋放的效果。

2)在距地心51 000 km左右釋放執(zhí)行奔月任務(wù)所需的速度增量極小。

3)天梯出發(fā)的奔月環(huán)月軌道軌道傾角分布在10°～170°之間。入口點(diǎn)在月球影響球上西半球，在南、北緯度均有分布。

4)計(jì)算天梯奔月1個(gè)月的發(fā)射窗口，本文的兩級(jí)并行算法加速比達(dá)6.95。