耿光有，王玨，侯錫云，余夢(mèng)倫，王建明，張志國

（1.北京航空航天大學(xué) 宇航學(xué)院，北京100083； 2.北京宇航系統(tǒng)工程研究所，北京100076；3.南京大學(xué) 天文與空間科學(xué)學(xué)院，南京210023）

為充分發(fā)揮新一代長征運(yùn)載火箭［1］氫氧入軌級(jí)動(dòng)力系統(tǒng)的高比沖優(yōu)勢(shì)，中國首次探測(cè)火星工程（簡稱探火工程）［2］中，采用了長征運(yùn)載火箭在近地停泊軌道滑行至預(yù)定位置后，二次點(diǎn)火工作將探測(cè)器直接送入地火轉(zhuǎn)移軌道的發(fā)射方式，顯著提高了探測(cè)器有效儀器設(shè)備的裝載量，但同時(shí)也明顯增加了運(yùn)載火箭的發(fā)射難度，如需要大型低溫長征運(yùn)載火箭在發(fā)射日實(shí)現(xiàn)相應(yīng)發(fā)射軌道的零窗口準(zhǔn)時(shí)發(fā)射等。而由于液氫約－253℃的嚴(yán)酷低溫工作條件限制，使得其點(diǎn)火工作條件保障困難，主要包括地面起飛時(shí)預(yù)冷好的點(diǎn)火條件，以及在停泊軌道上較長時(shí)間滑行后滿足二次點(diǎn)火工作的條件等［3-4］。國際上，雖然美國已經(jīng)實(shí)現(xiàn)了低溫入軌級(jí)具備近地滑行一圈后再點(diǎn)火的能力，但即便其成熟火箭Delta IV（采用氫氧推進(jìn)劑），在近年來也多有媒體報(bào)道其地面點(diǎn)火發(fā)射一再推遲。而為了滿足低溫氫氧入軌級(jí)的嚴(yán)酷二次點(diǎn)火工作要求，長征運(yùn)載火箭在當(dāng)前采用了較為保守的設(shè)計(jì)處理措施［3-4］，因此需要嚴(yán)格限制最長滑行時(shí)間不能超過設(shè)計(jì)門限。為了能夠在每26個(gè)月才能有一組較合適奔火任務(wù)的發(fā)射機(jī)會(huì)中實(shí)現(xiàn)及時(shí)發(fā)射，國際上普遍采用了持續(xù)2～3周的發(fā)射日期窗口與每天至少30min的發(fā)射窗口，以及相應(yīng)的發(fā)射軌道設(shè)計(jì)以保證深空探測(cè)任務(wù)的順利實(shí)施［5-6］。新研制的大型長征運(yùn)載火箭在滑行時(shí)間嚴(yán)格受限下，再兼顧測(cè)控布站限制等，完成持續(xù)2～3周發(fā)射日期窗口的探火發(fā)射軌道設(shè)計(jì)尤為艱難，本文即研究如何解決這一問題。

不考慮滑行時(shí)間嚴(yán)格受限和發(fā)射場(chǎng)航落區(qū)嚴(yán)格約束下，地火轉(zhuǎn)移軌道的設(shè)計(jì)已有較成熟的理論［7］，通常根據(jù)Lambert求解逃逸速度并采用圓錐曲線拼接法完成初值設(shè)計(jì)，再借助B平面矢量法［8］與微分修正完成行星際轉(zhuǎn)移軌道設(shè)計(jì)。B平面矢量法由于能夠使計(jì)算初值快速收斂而得以普遍應(yīng)用。國內(nèi)外此方面的文獻(xiàn)較多，如文獻(xiàn)［9］結(jié)合B平面矢量法給出了地球到火星轉(zhuǎn)移軌道的近地發(fā)射段雙曲線軌道的幾何約束分析，但沒有長征運(yùn)載火箭需要關(guān)注的停泊軌道最長滑行時(shí)間等具體約束限制下的深入分析研究等。文獻(xiàn)［10］分析了2009—2024年地球到火星間往返發(fā)射窗口機(jī)會(huì)及發(fā)射能量C3與逃逸赤緯等參數(shù)的數(shù)值列表及扼要分析等，但同樣不考慮運(yùn)載火箭停泊軌道的滑行時(shí)間等約束的限制。文獻(xiàn)［11］主要以運(yùn)載火箭發(fā)射重量及C3為約束，分析了探測(cè)器軌道設(shè)計(jì)，指出火星探測(cè)器全部的雙曲線出發(fā)軌道近地點(diǎn)包絡(luò)為一沿地球逃逸速度矢量VE∞為軸、以r0sinβ為半徑的圓軌跡（r0為近地點(diǎn)高度，β為近地點(diǎn)拱線與雙曲漸近線夾角）。文獻(xiàn)［12］根據(jù)探測(cè)器火星目標(biāo)軌道高度、傾角及天頂距約束分析了地球出發(fā)軌道速度、近地點(diǎn)幅角與赤經(jīng)等參數(shù)設(shè)計(jì)需求，然而發(fā)射場(chǎng)位置、射向及運(yùn)載火箭滑行時(shí)間等約束對(duì)發(fā)射軌道的限制關(guān)系并未涉及。文獻(xiàn)［13］將B平面矢量法推廣至制導(dǎo)策略設(shè)計(jì)中。文獻(xiàn)［14］分析了火星探測(cè)轉(zhuǎn)移軌道的偏差傳播，指出從近地點(diǎn)到近火點(diǎn)的位置誤差放大約100萬倍、速度誤差放大約10萬倍，實(shí)際飛行中途修正不可缺省。此外，一般的近地軌道衛(wèi)星發(fā)射任務(wù)，只要求運(yùn)載火箭的起飛時(shí)刻、星箭分離點(diǎn)的半長軸、偏心率、軌道傾角和近地點(diǎn)幅角等5個(gè)參數(shù)滿足要求即可以完成發(fā)射軌道設(shè)計(jì)，而火星等深空探測(cè)發(fā)射任務(wù)的最終目標(biāo)點(diǎn)是近火點(diǎn)目標(biāo)軌道，因此除要求運(yùn)載火箭以上參數(shù)嚴(yán)格滿足要求外，還需要星箭分離點(diǎn)的升交點(diǎn)經(jīng)度與真近點(diǎn)角等2個(gè)參數(shù)也必須滿足探測(cè)器要求，而這2個(gè)參數(shù)與實(shí)際飛行路徑直接關(guān)聯(lián)，因此如果按照傳統(tǒng)的軌道界面分工只能采用雙方反復(fù)拼接的設(shè)計(jì)方法，計(jì)算效率與精度等顯然難以滿足當(dāng)前設(shè)計(jì)需求。實(shí)際上，復(fù)雜、多約束的運(yùn)載火箭發(fā)射軌道設(shè)計(jì)對(duì)探測(cè)任務(wù)能否順利實(shí)施影響顯著。例如，2003年歐洲航天局采用聯(lián)盟號(hào)及Fregat上面級(jí)火箭發(fā)射火星快車探測(cè)器時(shí)，即遭遇到上述問題［15］：起初僅有10天的預(yù)期發(fā)射窗口及1 120 kg的運(yùn)載能力，發(fā)射風(fēng)險(xiǎn)較高，后來經(jīng)過聯(lián)盟號(hào)火箭優(yōu)化停泊軌道傾角及二次點(diǎn)火條件，以及優(yōu)化調(diào)整了整流罩與火箭三子級(jí)落區(qū)，最終將運(yùn)載能力增加到1 200 kg，發(fā)射日期窗口拓展到了30天，并且優(yōu)化后簡化了控制系統(tǒng)裝訂參數(shù)——采用9套飛行程序滿足了持續(xù)30天的發(fā)射日期窗口，僅需要探測(cè)器在飛行第10天時(shí)深空機(jī)動(dòng)10m／s，這些措施有力保障了探測(cè)工程的順利實(shí)施。美國DeltaⅡ7925型火箭發(fā)射探測(cè)火星等深空任務(wù)時(shí)，通過固定93°和99°兩個(gè)射向保障了航落區(qū)安全，簡化了控制系統(tǒng)裝訂參數(shù)，優(yōu)化確定了連續(xù)多日的發(fā)射日期窗口［5-6］。新一代長征運(yùn)載火箭從海南發(fā)射火星探測(cè)器時(shí)，由于低溫入軌級(jí)2次起動(dòng)間最長滑行時(shí)間嚴(yán)格受限［3-4］，使入軌點(diǎn)位置（即軌道近地點(diǎn)幅角）的調(diào)節(jié)范圍明顯減小，為滿足相應(yīng)發(fā)射年份下2～3周的有效發(fā)射日期窗口，分析發(fā)現(xiàn)不排除結(jié)合探測(cè)器深空速度脈沖機(jī)動(dòng)以完成發(fā)射軌道總體優(yōu)化［15-16］，故當(dāng)前文獻(xiàn)資料中的設(shè)計(jì)算法難以直接滿足長征運(yùn)載火箭發(fā)射軌道深入優(yōu)化的需要。依此背景，本文重點(diǎn)研究了長征運(yùn)載火箭在當(dāng)前滑行時(shí)間、發(fā)射場(chǎng)射向及落區(qū)等多約束限制下，發(fā)射探測(cè)器行星際轉(zhuǎn)移軌道的優(yōu)化設(shè)計(jì)算法。該算法直接采用探測(cè)器近火點(diǎn)目標(biāo)軌道參數(shù)，并結(jié)合近地出發(fā)軌道諸多設(shè)計(jì)約束開展設(shè)計(jì)，在高精度力學(xué)計(jì)算模型基礎(chǔ)上，提出雙向微分修正算法，從而也一并解決了在轉(zhuǎn)移軌道的深空對(duì)接點(diǎn)處增加速度脈沖機(jī)動(dòng)，以進(jìn)一步完成發(fā)射軌道優(yōu)化的問題。

1 長征運(yùn)載火箭海南發(fā)射場(chǎng)約束

工程中從飛行安全角度考慮，一般允許的射向A0范圍約為90°～115°。

注意到，地球自轉(zhuǎn)（赤道上轉(zhuǎn)速為465.1m／s）對(duì)實(shí)現(xiàn)長征運(yùn)載火箭目標(biāo)入軌速度的貢獻(xiàn)Vlaunch為

式中：фL為發(fā)射場(chǎng)緯度。

依據(jù)齊奧爾科夫斯基公式［17］，動(dòng)力比沖Isp，火箭起始質(zhì)量m0，燃燒結(jié)束質(zhì)量mk產(chǎn)生的速度增量ΔV為

對(duì)緯度約20°的發(fā)射場(chǎng)，如果射向從90°增加到115°，按照氫氧發(fā)動(dòng)機(jī)比沖約4 300 m／s，結(jié)合探測(cè)器與火箭末級(jí)入軌質(zhì)量的大致相對(duì)關(guān)系［18］，易知發(fā)射火星探測(cè)器等深空任務(wù)的運(yùn)載能力損失約3%，故可能情況下需要優(yōu)先選擇東射向發(fā)射。

2 火星探測(cè)任務(wù)發(fā)射軌道設(shè)計(jì)

2.1 發(fā)射約束下初始軌道搜索

2.1.1 Lambert問題求解

忽略地球和火星的引力及實(shí)際尺寸，認(rèn)為轉(zhuǎn)移軌道直接從地球（在日心坐標(biāo)系下的位置）出發(fā)，到達(dá)火星（在日心坐標(biāo)系下的位置），如圖1（a）所 示，該 問 題 為 典 型 的 日 心 Lambert問題：

式中：tM和tE分別為抵達(dá)和出發(fā)時(shí)刻；μ為日心天體引力常數(shù)；rM和rE分別為抵達(dá)和出發(fā)時(shí)刻地球與火星距日心距離；c為地球與火星的直線距離；a為轉(zhuǎn)移軌道半長軸。

關(guān)于Lambert問題的求解已有非常成熟的算法，本文不再贅述，推薦采用Izzo的計(jì)算方法［19］，如果軌道近似為拋物線，推薦采用Gooding的算法以避免Izzo算法中的奇點(diǎn)問題［20-21］。圖1（a）中，Vd和Va分別為Lambert轉(zhuǎn)移軌道影響球處出發(fā)和抵達(dá)速度矢量，VE和VM分別為地球和火星的速度矢量，V0和V1分別為Lambert轉(zhuǎn)移軌道相對(duì)于地球與火星影響球的速度矢量。通過V0和V1可以進(jìn)一步設(shè)計(jì)地球影響球（SOI）內(nèi)的雙曲線飛行軌跡及火星影響球內(nèi)的雙曲線飛行軌跡，從而將整個(gè)轉(zhuǎn)移軌道由地球影響球內(nèi)的雙曲線弧段＋日心Lambert弧段＋火星影響球內(nèi)雙曲線弧段構(gòu)成（見圖1（b）），此即經(jīng)典圓錐曲線拼接模型。

圖1 從地球轉(zhuǎn)移至火星的飛行軌道示意圖Fig.1 Schematic diagram of Earth-to-Mars flight trajectory

2.1.2 出入影響球的雙曲線軌道初值設(shè)計(jì)

本文著重介紹地球影響球內(nèi)雙曲線弧段設(shè)計(jì)方法，火星影響球內(nèi)雙曲線弧段設(shè)計(jì)方法類似。有別于傳統(tǒng)的借助B平面矢量法，本文算法直接根據(jù)近地端發(fā)射場(chǎng)的約束及火星端目標(biāo)軌道的約束設(shè)計(jì)雙曲線弧段。通過日心Lambert問題求解，得探測(cè)器相對(duì)于地球影響球邊界速度矢量Vo（地球與火星端分別對(duì)應(yīng)圖1中V0和V1），地球段采用下述計(jì)算單位：

式中：aE為地球赤道半徑；mE為地球質(zhì)量；［］表示括號(hào)內(nèi)變量在不同區(qū)域飛行段是不同參量，此即指無量綱化，如地球段、日心飛行段、火星飛行段等。

設(shè)定近地點(diǎn)高度hpE，從式（6）得偏心率eo：

雙曲線軌道傾角io與發(fā)射場(chǎng)射向AL及發(fā)射場(chǎng)緯度фL關(guān)系為

從式（8）中，得雙曲線軌道近地點(diǎn)幅角ωo和升交點(diǎn)赤經(jīng)Ωo為

式中：p＝cosiosinωosinhEo－ΓcosiocosωocoshEo；q＝ΓsinωocoshEo＋cosωosinhEo；μE為地球引力常數(shù)為速度矢量Vo的分量；Eo為地球影響球半徑ro處偏近點(diǎn)角，即

結(jié)合式（10）得對(duì)應(yīng)的真近點(diǎn)角fo：

式（8b）一般同時(shí)求得2個(gè)近地點(diǎn)幅角ωo，代入式（8a）分別得對(duì)應(yīng)的升交點(diǎn)赤經(jīng)Ωo。根據(jù)式（11）計(jì)算雙曲線軌道緯度幅角uo：

為便于區(qū)別，工程上通常以u(píng)o∈（－π／2，π／2）作為升交出發(fā)，以u(píng)o∈（π／2，3π／2）作為降交出發(fā)。

由式（5）～式（11）可以給出地球至火星轉(zhuǎn)移軌道在脫離地球影響球處的Kepler軌道根數(shù)ao、eo、io、Ωo、ωo、fo，而星箭分離點(diǎn)的地心軌道根數(shù)與該軌道只是真近點(diǎn)角fI不同，fI由火箭發(fā)射段確定。研究表明，不同射向及滑行時(shí)間對(duì)應(yīng)的fI變化不大，可近似認(rèn)為是常值。

綜上，在地球影響球內(nèi)，以地球和探測(cè)器組成二體問題動(dòng)力學(xué)模型，根據(jù)日心Lambert弧段給出的地球影響球邊界處的速度矢量Vo，可得J2000系下雙曲線出發(fā)軌道星箭分離點(diǎn)軌道根數(shù)ao、eo、io、Ωo、ωo、fI，并由之得位置矢量rJ2000和速度矢量再通過式（12）坐標(biāo)轉(zhuǎn)換［22］得星箭分離時(shí)刻與地固CGCS2000系瞬時(shí)重合（之后保持慣性空間不動(dòng)）的地心赤道慣性系（以下標(biāo)G表示）下位置rG和速度

式中：Q（t）為歲差章動(dòng)矩陣；R（t）為地球轉(zhuǎn)角矩陣；W（t）為極移矩陣。由再得相應(yīng)的軌道根數(shù)需求：ao、eo、ioE、ΩoE、ωoE、fI。

抵達(dá)近火點(diǎn)雙曲線軌道計(jì)算在火星J2000平赤道坐標(biāo)系（用下標(biāo)ME表示）下描述，該坐標(biāo)系與J2000火心坐標(biāo)系（用下標(biāo)MJ表示）轉(zhuǎn)換關(guān)系［22］見式（13），IAU年度工作報(bào)告對(duì)火星北極指向α0、δ0的精確數(shù)據(jù)每三年互聯(lián)網(wǎng)發(fā)布一次。

式中：矩陣M下標(biāo)x，y，z分別表示繞相應(yīng)坐標(biāo)軸的轉(zhuǎn)換矩陣［17，22］?；鹦怯绊懬騼?nèi)的雙曲線弧段設(shè)計(jì)方法與地球影響球的雙曲弧段設(shè)計(jì)方法類似但不是本文關(guān)注的重點(diǎn)，所以火星端的雙曲線軌道設(shè)計(jì)這里不再贅述。

2.1.3 長征運(yùn)載火箭滑行時(shí)間及發(fā)射場(chǎng)設(shè)計(jì)約束

長征運(yùn)載火箭發(fā)射軌道星下點(diǎn)軌跡幾何示意見圖2，采用圓停泊軌道，忽略各攝動(dòng)因素，得

圖2 探測(cè)火星地球出發(fā)軌道的幾何圖示Fig.2 Earth-to-Mars launch trajectory geometry illustration

設(shè)θ1為從長征運(yùn)載火箭起飛到進(jìn)入停泊軌道前的工作弧段，θ2為停泊軌道滑行結(jié)束至探測(cè)器分離的弧段（分析表明θ1、θ2可以近似為常值），其他符號(hào)含義同前，據(jù)各參量物理含義，易得所需長征運(yùn)載火箭入軌級(jí)2次起動(dòng)間滑行時(shí)間Th為

圖2中，Υ表示春分點(diǎn)，αo和δo分別為地球至火星轉(zhuǎn)移軌道在地球影響球處的逃逸赤經(jīng)和赤緯，其他符號(hào)含義均在上下文中有明確說明。根據(jù)長征運(yùn)載火箭發(fā)射探測(cè)器入軌時(shí)刻升交點(diǎn)赤經(jīng)Ω及起飛時(shí)刻發(fā)射場(chǎng)恒星時(shí)角αL定義，結(jié)合圖2球面三角關(guān)系易得

式中：λL為發(fā)射場(chǎng)經(jīng)度；Ωe為探測(cè)器分離點(diǎn)時(shí)刻發(fā)射軌道計(jì)算的升交點(diǎn)經(jīng)度；ωe為地球自轉(zhuǎn)角速度；tf0為火箭滑行時(shí)間Th＝0 s時(shí)從起飛到探測(cè)器分離的時(shí)間。對(duì)選定的射向AL，預(yù)先選擇很小的閾值εi＞0（i＝1，2），調(diào)整發(fā)射當(dāng)日出發(fā)時(shí)刻tE使：

易得滿足長征運(yùn)載火箭最長滑行時(shí)間限制的飛行軌道參數(shù)及起飛時(shí)刻初解。

2.1.4 射向?qū)︼w行軌道及航程角的影響

長征運(yùn)載火箭關(guān)注的發(fā)射日出發(fā)軌道及航程角與射向的初步關(guān)系分析如下。結(jié)合圖2，得飛行航程角θ（θ＝ωoE＋fo－d）與地球影響球處逃逸赤緯δo、射向AL及發(fā)射場(chǎng)緯度фL關(guān)系：

設(shè)фL＝20°，可得圖3。最長滑行時(shí)間限定了近地點(diǎn)幅角ωoE與飛行航程角θ的變化范圍，每一組航程角對(duì)應(yīng)一個(gè)起飛時(shí)刻。

結(jié)合圖3，在射向范圍0°～180°情況下，結(jié)合發(fā)射軌道基礎(chǔ)理論，分析易知：

2）當(dāng)δo＝фL，僅有一組射向隨航程角θ減小而增大的解；當(dāng)δo＝－фL，除了一組射向隨航程角θ減小而增大的解外，還有航程角恰好為180°、射向任意的另外一組解。

圖3 射向A L、飛行航程角θ與δo關(guān)系Fig.3 Contourmap ofδo with respect to launch azimuth A L and flight range angleθ

2.1.5 結(jié)合Pork-chop圖完成初始發(fā)射軌道優(yōu)選

采用上述加入滑行時(shí)間限制再精確雙向微分修正的設(shè)計(jì)算法可以較快速完成初步發(fā)射軌道搜索，通過對(duì)近地點(diǎn)與近火點(diǎn)速度增量輔以Porkchop圖［7］分析，即可以初步完成多設(shè)計(jì)約束下，長征運(yùn)載火箭發(fā)射火星探測(cè)軌道的初步運(yùn)載能力確定、連續(xù)發(fā)射日期窗口初選和初始發(fā)射軌道設(shè)計(jì)。限于篇幅，此部分內(nèi)容從略。

2.2 精確設(shè)計(jì)發(fā)射軌道

2.2.1 精確動(dòng)力學(xué)方程

在起飛時(shí)刻地心慣性系下，長征運(yùn)載火箭飛行的質(zhì)心動(dòng)力學(xué)矢量方程［17，23］為

式中：m為火箭質(zhì)量；rG0為慣性系原點(diǎn)到火箭質(zhì)心矢徑；P為火箭推力；R f為作用在火箭上的氣動(dòng)力；g為作用在火箭上的地球引力加速度。通過式（20b）轉(zhuǎn)換得入軌時(shí)刻地心慣性系下位置和速度：

進(jìn)而得相應(yīng)的發(fā)射軌道根數(shù)：ao、eo、ioE、Ωe、ωoE、fI。

從探測(cè)器分離開始，設(shè)定整個(gè)地火轉(zhuǎn)移軌道都是在日心J2000坐標(biāo)系下完成積分，動(dòng)力學(xué)方程為

式中：μS為太陽引力系數(shù)；j∈⊙表示太陽系八大行星及月球等；μj為太陽系下擾動(dòng)天體引力系數(shù)，對(duì)于地火轉(zhuǎn)移軌道設(shè)計(jì)，需要考慮八大行星及月球的影響分別為地球和火星的非球形攝動(dòng)力［24-25］；δE、δM為開關(guān)函數(shù)，當(dāng)探測(cè)器在地球（火星）影響球界內(nèi)取1，其余取0；Fothers指地球輻射壓、高層空間大氣、火星稀薄大氣、小行星引力等；FSRP為太陽輻射壓，設(shè)定探測(cè)器面質(zhì)比為S／m，S為表面積，ms為質(zhì)量，忽略行星陰影區(qū)影響，則

其中：κ為探測(cè)器反射系數(shù)；c0為光速；ρ⊙為距離r⊙處的太陽輻射流。

采用RKF7（8）階龍格庫塔算法，分析計(jì)算表明，探測(cè)器在靠近地球或火星飛行段均采用10×10階非球形引力項(xiàng)，比兩者均采用J2項(xiàng)時(shí)設(shè)計(jì)出的轉(zhuǎn)移軌道在近地端的速度差異小于0.1m／s、位置差異小于30 m，在近火端的速度差異約0.3m／s、位置差異約140m；光壓對(duì)探測(cè)器的攝動(dòng)影響量一般可能超過近地端高階攝動(dòng)影響而小于火星端高階項(xiàng)影響；注意到光壓影響與探測(cè)器具體形態(tài)關(guān)系較密切等因素，故綜合后下面以考慮大行星及月球等中心天體引力、地球與火星飛行段均采用10×10階次非球形引力項(xiàng)給出分析算例。

2.2.2 雙向微分修正算法實(shí)現(xiàn)精確對(duì)接

記星箭分離點(diǎn)時(shí)刻（記為t0）實(shí)際飛行軌道與參考軌道的狀態(tài)量（即位置和速度）偏差為ΔX（t0），小偏差條件下，t時(shí)刻狀態(tài)量偏差ΔX（t）近似滿足：

式中：Φ（t，t0）為參考軌道t0時(shí)刻至t時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣［7］。由于從近地點(diǎn)到近火點(diǎn)轉(zhuǎn)移軌道的誤差傳遞特性急劇放大，因此直接采用近火星目標(biāo)軌道根數(shù)求解近地點(diǎn)出發(fā)軌道極其困難，一般需要轉(zhuǎn)化處理。

若分別從近地點(diǎn)、近火點(diǎn)向地球至火星轉(zhuǎn)移軌道中間的某個(gè)時(shí)間點(diǎn)tm積分，記地球出發(fā)端的狀態(tài)量偏差為ΔX（t0），火星到達(dá)端的狀態(tài)量偏差為ΔX（tf），類似式（21），有

式中：ΔXE（tm）為從地球端正向積分至tm時(shí)刻軌道的狀態(tài)量偏差；ΔXM（tm）為從火星端逆向積分至tm時(shí)刻軌道的狀態(tài)量偏差。對(duì)同一條轉(zhuǎn)移軌道，要求這2個(gè)偏差相同，即

再將式（21）中的t取為tf，并結(jié)合式（22）易得）

即時(shí)刻t0～tf的狀態(tài)轉(zhuǎn)移矩陣Φ（tf，t0）可分解為2個(gè)分弧段的狀態(tài)轉(zhuǎn)移矩陣ΦE（tm，t0）和ΦM（tm，tf）；由于將狀態(tài)轉(zhuǎn)移矩陣Φ（t，t0）一分為二后，兩積分段的誤差傳遞特性急劇降低，即可以較順利地實(shí)現(xiàn)問題求解。實(shí)際上B平面矢量法的本質(zhì)原理亦是如此，只是固定為目標(biāo)天體影響球邊界處的B平面參數(shù)。分析表明，由于日心轉(zhuǎn)移弧段的良好解析性，使中間積分結(jié)合點(diǎn)可以自由移動(dòng)而不影響計(jì)算結(jié)果，所以就解決了后續(xù)在結(jié)合點(diǎn)處增加深空機(jī)動(dòng)以進(jìn)一步優(yōu)化長征火箭發(fā)射軌道的問題。下面對(duì)算法進(jìn)一步介紹。

以rJ2000和.rJ2000代表在地心J2000坐標(biāo)系下地球出發(fā)時(shí)探測(cè)器與長征運(yùn)載火箭分離點(diǎn)的位置和速度，以rMJ和.rMJ代表探測(cè)器火心J2000坐標(biāo)系下抵達(dá)近火點(diǎn)的位置和速度，組成一個(gè)12維的狀態(tài)矢量Y，即

式中：符號(hào)R表示相對(duì)日心系的位置，下標(biāo)d表示地球影響球界出發(fā)點(diǎn)，a表示火星影響球界抵達(dá)點(diǎn)。記下標(biāo)E表示探測(cè)器分離點(diǎn)時(shí)刻地球的日心系位置和速度，下標(biāo)M表示抵達(dá)近火點(diǎn)時(shí)刻火星的日心位置和速度，則

在地球出發(fā)點(diǎn)，存在3個(gè)約束：

在近火雙曲線軌道段，也存在3個(gè)類似約束：

式中：上標(biāo)“＋”表示正向積分軌道，上標(biāo)“－”表示反向積分軌道。基于長征運(yùn)載火箭發(fā)射軌道優(yōu)化中，后續(xù)關(guān)注的深空速度脈沖機(jī)動(dòng)優(yōu)化問題，分析表明，Tm可以從地球影響球邊界至火星影響球邊界之間自由選擇，相應(yīng)根據(jù)深空對(duì)接點(diǎn)處速度脈沖增量為g10～12賦值［16］；對(duì)于無需深空脈沖速度機(jī)動(dòng)的情況，一般可取Tm＝（T1＋T2）／2。

將式（27）～式（29）統(tǒng)一記為G（Y）＝（g1，g2，…，g12）＝0此時(shí)迭代方程為

式中：

可將雅可比矩陣分解為如下形式：

式中：

12項(xiàng)約束對(duì)應(yīng)12個(gè)變量，通過迭代求解即完成火星目標(biāo)軌道的地球出發(fā)軌道設(shè)計(jì)。

2.3 計(jì)算流程

考慮發(fā)射場(chǎng)約束，分析與精確設(shè)計(jì)長征運(yùn)載火箭發(fā)射火星探測(cè)器轉(zhuǎn)移軌道的設(shè)計(jì)流程見圖4。

圖4 長征運(yùn)載火箭設(shè)計(jì)探測(cè)火星軌道示意圖Fig.4 Schematic diagram of Long March launch vehicle for Earth-to-Mars trajectory design

3 分析及算例

計(jì)算采用的火星探測(cè)器目標(biāo)軌道約束為：近火點(diǎn)高度500 km，傾角93°；結(jié)合初步選定的射向，近地約束為：射向AL＝107°，出發(fā)近地點(diǎn)高度200 km，以UTC時(shí)間2020-06-19出發(fā)，轉(zhuǎn)移時(shí)間198天，2021-01-02T23：59升交抵達(dá)火星為算例，計(jì)算得地球影響球邊界處逃逸赤緯約24°，由圖3知有滑行時(shí)間Th相近的2組初解，采用雙向微分修正得表1、表2精確解，表中EJ2000表示地心J2000赤道系，ME表示J2000火心IAU平赤道系［22］，SJ2000代表距近火點(diǎn)時(shí)刻tm＝90天時(shí)日心J2000系參數(shù)，ΔV指對(duì)應(yīng)地球或火星端各自圓停泊軌道的速度增量。

分析表明：①tm可以大幅移動(dòng)而不影響優(yōu)化設(shè)計(jì)結(jié)果；②在結(jié)合點(diǎn)處速度增量可以為0 m／s或增加到超過500m／s，都能找到良好求解。

若限定最長滑行時(shí)間不能超過2 000 s［3-4］，則只有表1中解滿足要求；而如果限定最長滑行時(shí)間不超過1 100 s或更短，則需要采用上述加入滑行時(shí)間限制再精確雙向微分修正的設(shè)計(jì)算法，通過加大射向或結(jié)合探測(cè)器深空機(jī)動(dòng)等完成發(fā)射軌道優(yōu)化；具體的分析是繁瑣復(fù)雜的工程問題，已不在本文討論范圍。

將表1、表2中計(jì)算結(jié)果在同等條件下，與STK仿真分析軟件積分對(duì)算表明，二者計(jì)算精度與所需時(shí)長相當(dāng)，算例中從近地點(diǎn)到深空對(duì)接點(diǎn)、從深空對(duì)接點(diǎn)到近火點(diǎn)，位置差異約為1 km。

表1 發(fā)射軌道T h ＝1 849 s，探測(cè)器05∶25釋放，地球升交出發(fā)Table 1 T rajectory T h ＝1 849 s，p robe released on 05∶25，ascending from the Earth

表2 發(fā)射軌道T h ＝2 312 s，探測(cè)器03∶15釋放，地球降交出發(fā)Table 2 T rajectory T h ＝2 312 s，p robe released on 03∶15，descending from the Ear th

4 結(jié) 論

本文給出了長征運(yùn)載火箭允許滑行時(shí)間受限和發(fā)射場(chǎng)可用射向范圍約束下，長征運(yùn)載火箭發(fā)射火星探測(cè)任務(wù)從地面起飛至近火點(diǎn)軌道設(shè)計(jì)的完整模型算法。

1）解決了適用于工程運(yùn)算所需精度的分析研究模型問題。

2）解決了火星探測(cè)任務(wù)的發(fā)射軌道優(yōu)化設(shè)計(jì)問題，確定了運(yùn)載能力、優(yōu)選了發(fā)射日期窗口，完成了航落區(qū)優(yōu)化選擇等。

雙向微分修正算法具有結(jié)合點(diǎn)可以大幅移動(dòng)，以便于后續(xù)分析中探測(cè)器采用深空機(jī)動(dòng)后，深入開展長征運(yùn)載火箭發(fā)射軌道再優(yōu)化分析的優(yōu)勢(shì)。工程應(yīng)用表明此方法穩(wěn)定、可靠好用。除了可以用于地火轉(zhuǎn)移發(fā)射軌道設(shè)計(jì)外，本文算法還可用于長征運(yùn)載火箭發(fā)射其他深空探測(cè)器的發(fā)射軌道設(shè)計(jì)。

致謝衷心感謝國家航天局、國家留學(xué)基金管理委員會(huì)對(duì)本課題的資助與支持，以及中國科學(xué)院大學(xué)李明濤教授在STK軟件驗(yàn)證中給予的幫助。