孫沖, 袁建平, 方群, 崔堯, 汪明曉, 劉易斯

(1．西北工業(yè)大學(xué)航天飛行動(dòng)力學(xué)技術(shù)國家級(jí)重點(diǎn)實(shí)驗(yàn)室，陜西西安710072;2．西北工業(yè)大學(xué)深圳研究生院，深圳518057;3．上海宇航系統(tǒng)工程研究所，上海201109; 4．國家航天局空間碎片監(jiān)測與應(yīng)用中心，北京100012)

連續(xù)小推力軌道設(shè)計(jì)在空間機(jī)動(dòng)軌道設(shè)計(jì)方面有廣泛的應(yīng)用[1-2]。與傳統(tǒng)脈沖推力或大推力相比，連續(xù)小推力軌道機(jī)動(dòng)容易控制，機(jī)動(dòng)能力強(qiáng)；另外其比沖小，實(shí)現(xiàn)軌道機(jī)動(dòng)所需要的燃料很少，因而成為研究熱點(diǎn)。文獻(xiàn)[3]中，Mcinnes首次提出了使用太陽帆來改變引力場的大小，設(shè)計(jì)了懸浮軌道，并推導(dǎo)了實(shí)現(xiàn)懸浮軌道所需推力的解析解。Malcolm采用太陽帆產(chǎn)生切向固定推力，設(shè)計(jì)了從地球軌道到太陽極軌的轉(zhuǎn)移軌道，解決了軌道優(yōu)化問題[4]。然而，由于單獨(dú)采用電推進(jìn)或者太陽帆產(chǎn)生的推力較小，軌道轉(zhuǎn)移常常需要很長時(shí)間。針對這個(gè)問題，Mengali等在文獻(xiàn)[5]中提出了太陽帆和太陽帆/電推進(jìn)結(jié)合的混合推進(jìn)方法，并采用間接優(yōu)化方法進(jìn)行了軌道優(yōu)化。在其文章中，太陽帆產(chǎn)生徑向推力，用于改變引力場的引力常數(shù)，而與此同時(shí)電推進(jìn)產(chǎn)生固定切向推力。在假定太陽帆能夠產(chǎn)生50%的中心引力(μ2=0.5μ1)的條件下，混合推力能夠減少燃料的消耗并且縮短空間機(jī)動(dòng)任務(wù)時(shí)間。然而該方法是作為一種間接優(yōu)化方法，并沒有考慮推力約束問題，很難應(yīng)用于實(shí)際空間機(jī)動(dòng)軌道設(shè)計(jì)中。

本文針對連續(xù)小推力軌道設(shè)計(jì)問題，提出了一種基于虛擬引力場的混合連續(xù)小推力軌道設(shè)計(jì)方法[6-7]。該方法將虛擬引力場與平均算法結(jié)合，通過混合連續(xù)小推力形成虛擬引力場，進(jìn)而采用平均算法求解在虛擬引力場下固定大小切向力作用下的機(jī)動(dòng)軌道解析解。進(jìn)一步，本文采用四階龍格庫塔積分方法分析解析解的精度，并提出了誤差因子參數(shù)，從而保證解析解的精度。本文采用虛擬引力場方法，將航天器機(jī)動(dòng)軌道設(shè)計(jì)優(yōu)化問題轉(zhuǎn)化成參數(shù)優(yōu)化問題。該方法主要分為3步：①設(shè)計(jì)并優(yōu)化滿足任務(wù)約束(包括邊界約束和推力約束)的虛擬引力場，以及固定切向推力的大??；②求解實(shí)現(xiàn)虛擬引力場的推力大?。虎劭紤]太陽帆的動(dòng)力學(xué)模型，求解采用太陽帆來實(shí)現(xiàn)虛擬引力場的控制律。

1 模 型

1.1 虛擬中心引力場的定義

虛擬中心引力場是指引力場中心不在地心上，且其引力常數(shù)不等于地心引力常數(shù)的中心引力場。該引力場可用參數(shù)(r0,μ2)來定義,其中r0為虛擬引力場中心與地心相對位置矢量,稱為虛擬引力場位置參數(shù);μ2為虛擬引力場引力常數(shù)[6-7]。

1.2 虛擬中心引力場下虛擬圓錐曲線軌道

根據(jù)虛擬中心引力場參數(shù)r0維數(shù)的不同,虛擬中心引力場可以分為2種,如圖1所示;

圖1 虛擬中心引力場示意圖

設(shè)定地心引力場引力常數(shù)為μ1,航天器在地心慣性坐標(biāo)系下位置矢量為r1,假定在航天器上施加推力加速度為at,則有:

(1)

假設(shè)地心慣性坐標(biāo)系為OXYZ,虛擬中心引力場坐標(biāo)系為O′X′Y′Z。推力與地心引力場的合力形成的虛擬中心引力場參數(shù)為(r0,μ2),則在虛擬中心引力場坐標(biāo)系下有:

(2)

(3)

由(2)式、(3)式可知,在虛擬中心引力場坐標(biāo)系下,航天器運(yùn)動(dòng)方程可表示為如下形式:

(4)

實(shí)現(xiàn)虛擬中心引力場所需推力為

(5)

將(3)式帶入(5)式, 可得

(6)

由上述分析可知,通過調(diào)整推力,可形成一種虛擬中心引力場。在該虛擬中心引力場中,航天器運(yùn)動(dòng)軌跡為圓錐曲線。若該虛擬圓錐曲線軌道滿足軌道轉(zhuǎn)移約束,則可實(shí)現(xiàn)軌道機(jī)動(dòng)。

1.3 固定切向力作用下軌道攝動(dòng)方程

在切向推力作用下,平面軌道高斯攝動(dòng)方程為

(7)

(8)

徑向推力以及周向推力和切向推力之間的關(guān)系為

(9)

而真近點(diǎn)角和偏近點(diǎn)角的關(guān)系為

(10)

因此,將公式(7)中自變量轉(zhuǎn)換為偏近點(diǎn)角,則軌道半長軸和偏心率的變化[9]為

(11)

式中，E為偏近點(diǎn)角,ft為連續(xù)切向推力。在給定軌道交會(huì)的初始點(diǎn)A和終點(diǎn)B的條件下,采用混合推力方法設(shè)計(jì)機(jī)動(dòng)軌道就是求解引力場參數(shù)(μ2,rvg)以及切向推力ft。若選擇合適參數(shù)X=(μ2,rvg,fsep)時(shí),在虛擬引力場中和切向推力的同時(shí)作用下,航天器在虛擬引力場可從A點(diǎn)轉(zhuǎn)移到目標(biāo)點(diǎn)B點(diǎn),從而滿足軌道交會(huì)的邊界約束。這就是采用基于虛擬引力場方法的混合連續(xù)小推力軌道設(shè)計(jì)的基本思路。

采用虛擬引力場的混合連續(xù)小推力軌道設(shè)計(jì)方法來設(shè)計(jì)軌道,實(shí)際上是設(shè)計(jì)在虛擬引力場中的常值切向推力轉(zhuǎn)移軌道。轉(zhuǎn)移軌道需要滿足兩點(diǎn)邊界值約束,以及推力約束。當(dāng)初始軌道的偏心率很小,推力較小的條件下[9],在虛擬引力場中的切向推力軌道近似解析解為

(12)

當(dāng)初始點(diǎn)和目標(biāo)點(diǎn)已知時(shí),上式中的參數(shù)E已知。邊界約束可以轉(zhuǎn)化為

(13)

由此可見,采用虛擬引力場下的連續(xù)推力軌道設(shè)計(jì)方法,可將所有滿足邊界約束的軌道用一組參數(shù)X=(μvg,rvg,fsep)表示,即轉(zhuǎn)移軌道參數(shù)化。這樣轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化問題就可以轉(zhuǎn)化成為參數(shù)設(shè)計(jì)與優(yōu)化問題。

由上分析可知,采用虛擬引力場下的混合連續(xù)小推力軌道設(shè)計(jì)方法,需要給航天器施加推力,一部分推力用于形成虛擬引力場,另一部分推力用于形成虛擬引力場中的切向力。具體推力分配:將形成虛擬引力場所需要的連續(xù)推力與在切向力的合力分解為徑向力和切向力,然后采用太陽帆和電推進(jìn)分別產(chǎn)生分解的徑向力和切向力。由公式(6)可知,實(shí)現(xiàn)虛擬引力場所需要的推力是有虛擬引力場參數(shù)決定的如下式所示,

Trvg=F(μvg,rvg,fsep)

1.4 太陽帆/電推進(jìn)混合最近最優(yōu)軌道設(shè)計(jì)

1.4.1 太陽帆推進(jìn)模型

由公式(5)可知,實(shí)現(xiàn)虛擬引力場需要的力隨著軌道半徑的變化而變化。因而采用太陽帆來實(shí)現(xiàn)虛擬引力場,其推力大小和方向也都是實(shí)時(shí)變化的。在此采用太陽帆推進(jìn)的理想模型[4]

Fsail=2PAcos2αn

(14)

式中，P=(S0/c)(r0/r)2,為太陽光壓率。

采用理想推進(jìn)模型,太陽帆產(chǎn)生的推力表示如下

(15)

式中，a1,a2,a3均為常數(shù)。而太陽帆產(chǎn)生的力在軌道的徑向和切向的分量為

(16)

由此可知太陽帆推進(jìn)可以實(shí)現(xiàn)部分徑向和周向的正推力,而剩余部分需要用電推進(jìn)來實(shí)現(xiàn)。

1.4.2 太陽帆電推進(jìn)模型

本文選擇文獻(xiàn)[5]中的電推進(jìn)器,其產(chǎn)生的小推力和電推進(jìn)器功率參數(shù)相關(guān),而需要的燃料消耗也由點(diǎn)推進(jìn)器的功率決定。點(diǎn)推進(jìn)器的功率是可變的,如公式(17)所示

(17)

則其產(chǎn)生的推力加速度在徑向周向分量為:

(18)

1.4.3 推力分配

采用虛擬引力場混合小推力軌道設(shè)計(jì)方法設(shè)計(jì)軌道,需要推力來形成虛擬引力場以及在虛擬引力場中的切向力。本文采用太陽帆/電推進(jìn)混合推力來實(shí)現(xiàn)虛擬引力場,并提供在虛擬引力場中的常值推力。由于形成虛擬引力場所需的徑向力和周向力是固定的,而太陽帆產(chǎn)生的推力在徑向和周向是耦合的,因此僅僅采用太陽帆無法直接提供形成虛擬引力場需要的推力。在此,本文在選擇太陽帆的同時(shí),另外選擇安裝2個(gè)電推進(jìn)器,分別用于產(chǎn)生徑向和周向的推力,使其與太陽帆光壓推進(jìn)產(chǎn)生的徑向力和切向力相互匹配,共同作用產(chǎn)生形成引力勢場,并提供常值推力。在此假定形成引力勢場所需要的徑向和周向的推力加速度為Tr,Tc,而太陽帆產(chǎn)生的徑向和周向的推力加速度為assr,assc。電推進(jìn)產(chǎn)生的徑向和周向的推力加速度為asepr,asepc,則有如下公式

(19)

由太陽帆的數(shù)學(xué)模型可知,其產(chǎn)生的徑向力和軸向力與太陽帆的位置和太陽帆姿態(tài)有關(guān)。采用虛擬引力場和常值推力結(jié)合的方法,其太陽帆的在每一時(shí)刻的位置是一定的,因此太陽帆產(chǎn)生的推力,實(shí)際上是由太陽帆的姿態(tài)角來確定的。為了能夠最大程度利用太陽帆產(chǎn)生的推力,從而節(jié)省燃料,需要對其姿態(tài)角進(jìn)行優(yōu)化,在此假定太陽帆的姿態(tài)角α和軌道機(jī)動(dòng)的轉(zhuǎn)移角θ有如下的關(guān)系:

α=a2·θ2+a1·θ+a0

(20)

式中，a0,a1,a2為待優(yōu)化的系數(shù)項(xiàng)。

當(dāng)給定虛擬引力場后,軌道轉(zhuǎn)移過程中需要電推進(jìn)器產(chǎn)生推力的大小為:

(21)

對比本文與文獻(xiàn)[5]可知,本文的推力分配與其不同。文獻(xiàn)[5]中軌道機(jī)動(dòng)的徑向推力完全由太陽帆提供,而切向力由電推進(jìn)提供。本文采用虛擬引力場方法設(shè)計(jì)機(jī)動(dòng)軌道,太陽帆推力一部分提供徑向力一部分提供切向力,這樣可以減少電推進(jìn)所需提供的推力。因而本文的推力分配方法比較靈活,可節(jié)省軌道機(jī)動(dòng)所需燃料。

1.4.4 混合小推力軌道設(shè)計(jì)與優(yōu)化

由上述分析可知,采用基于虛擬引力場的混合連續(xù)小推力軌道設(shè)計(jì)方法,可求出大量的滿足邊界約束的轉(zhuǎn)移軌道,如圖3所示。假定虛擬引力場參數(shù)以及切向力參數(shù)的參數(shù)范圍為

圖2 混合小推力軌道交會(huì)

圖3 混合小推力轉(zhuǎn)移軌道優(yōu)化示意圖

則軌道優(yōu)化問題即可轉(zhuǎn)化為參數(shù)優(yōu)化問題,優(yōu)化的參數(shù)為虛擬引力場參數(shù)和切向推力加速度,以及太陽帆光壓推進(jìn)的推力系數(shù),而優(yōu)化的目標(biāo)函數(shù)可以為軌道轉(zhuǎn)移過程中,電推進(jìn)器消耗燃料質(zhì)量:

J=min(mfuel)

(22)

式中，mfuel為任務(wù)結(jié)束時(shí)航天器的質(zhì)量。

或是轉(zhuǎn)移軌道的任務(wù)時(shí)間:

J=min(Tf)

(23)

或是同時(shí)考慮燃料消耗和飛行時(shí)間的權(quán)重函數(shù)

J=min(k1·mfuel+k2·Tf)

(24)

式中，k1,k2為權(quán)重,k1+k2=1。

結(jié)合上述分析可知,整體軌道設(shè)計(jì)流程包括以下幾個(gè)步驟:

a) 給出初始軌道和目標(biāo)軌道;

b) 給出滿足推力約束的虛擬引力場參數(shù)范圍和切向推力范圍;

c) 在b)的范圍中,選擇虛擬引力場參數(shù)。

d) 將初始軌道狀態(tài)(ri,vi)轉(zhuǎn)化為虛擬引力

場中的軌道狀態(tài)(rvg,vvg);

e) 優(yōu)化切向力參數(shù),求解滿足邊界約束條件的轉(zhuǎn)移軌道;

f) 重復(fù)a)～e),得到大量的滿足約束的轉(zhuǎn)移軌道;

g) 采用Matlab軟件中的Fmincon函數(shù)優(yōu)化轉(zhuǎn)移軌道。

2 誤差分析

由第1節(jié)分析可知,采用虛擬引力場下混合連續(xù)小推力軌道設(shè)計(jì)方法,可求解出滿足推力約束和邊界值約束的解析解。然而在此需強(qiáng)調(diào),公式(10)中的解析解是存在誤差的。在此為了保證該解析方法的精度,本節(jié)采用四階龍格庫塔數(shù)值積分方法求解了軌道積分?jǐn)?shù)值解,求出了解析解的誤差。在此基礎(chǔ)上,分析誤差產(chǎn)生的原因并給出解析解的應(yīng)用范圍。具體的仿真條件見表1。

在此,將解析解的誤差分為位置誤差和速度誤差,并分別用Δεr,Δεv表示

Δεr=ranslytical-rnumerical=f(x)

Δεv=vanalytical-vnumerical=g(x)

對比圖4圖5,以及公式(10)、公式(11)可知,解析解與數(shù)值解的誤差和飛行時(shí)間、引力場參數(shù)、以及軌道參數(shù)有關(guān)。切向推力越大、飛行時(shí)間越長,該解析解的誤差會(huì)越大;初始軌道(虛擬引力場中)的偏心率越大、半長軸越大其誤差越大。這一結(jié)論也與文獻(xiàn)[1]得出的結(jié)論一致。為描述解析解與數(shù)值解的誤差,在此定義一個(gè)誤差因子參數(shù),如(25)式所示

(25)

對比公式(11)及(12)可知,解析解的誤差因子會(huì)隨著切向推力、軌道半長軸、軌道轉(zhuǎn)移時(shí)間(偏近點(diǎn)角)的增大而增大,當(dāng)(e<0.7)時(shí),偏心率越大,誤差參數(shù)越大。如圖5所示,采用解析方法可以得到大量的轉(zhuǎn)移軌道,同時(shí)采用四階龍格庫塔數(shù)值方法,就可以求解并分析這些解析解的速度誤差、位置誤差以及誤差因子之間的關(guān)系。

以地球-火星轉(zhuǎn)移軌道為例,采用虛擬引力場方法可以得到大量的轉(zhuǎn)移軌道,在其中任意選擇n條轉(zhuǎn)移軌道(n=157)。采用數(shù)值積分方法,可以得出這些解析解的速度誤差和位置誤差與誤差因子的關(guān)系,如圖6所示。

圖6 速度、位置誤差隨誤差因子的變化

圖6中,藍(lán)色的線條為位置誤差;紅色的線條為速度誤差;綠色的線條為誤差因子。由圖6可知,速度誤差、位置誤差和誤差因子是同步的,當(dāng)誤差因子很小時(shí),其對應(yīng)的速度誤差和位置誤差也很小;相反,當(dāng)誤差因子比較大時(shí),其速度誤差和位置誤差都很大。由圖6可知，速度誤差和位置誤差都隨著誤差因子的增大而增大。當(dāng)保證誤差因子小于0.002 5時(shí),速度誤差小于0.001%,而位置誤差小于0.002%。由以上分析可知,在采用虛擬引力場常值切向力解析方法設(shè)計(jì)轉(zhuǎn)移軌道時(shí),只需要保證誤差因子k

3 數(shù)值仿真

表1 地球火星軌道轉(zhuǎn)移虛擬引力場參數(shù)和常值推力

表2 太陽帆姿態(tài)角系數(shù)

圖7 虛擬引力場混合推力方法的轉(zhuǎn)移軌道

圖8 切向推力方法的轉(zhuǎn)移軌道(相同時(shí)間)

圖7為采用虛擬引力場混合推力方法設(shè)計(jì)的轉(zhuǎn)移軌道；圖8為在相同轉(zhuǎn)移時(shí)間條件下，采用常值切向力方法設(shè)計(jì)的轉(zhuǎn)移軌道。圖9表述轉(zhuǎn)移軌道所需的周向推力和徑向推力以及太陽帆提供的周向推力和徑向推力。圖10是轉(zhuǎn)移軌道太陽帆姿態(tài)角。圖11是為相同情況下(切向推力大小相等)得到的轉(zhuǎn)移軌道。圖12為實(shí)現(xiàn)虛擬引力場所需要的推力加速度。圖13和圖14表示推力秒流量和推力方向角。由表5可知，相比于常值切向推力，虛擬引力場混合推力方法能夠很大程度上節(jié)省燃料，縮短時(shí)間。假定軌道轉(zhuǎn)移時(shí)間相等，虛擬引力場混合推力方法所需要的燃料僅為常值推力方法的 。采用虛擬引力場混合連續(xù)小推力方法，所需要的電推進(jìn)推力大大減小了，因而虛擬引力場混合推力方法能夠擴(kuò)大小推力軌道設(shè)計(jì)方法的應(yīng)用范圍。另一方面，在同樣大小推力條件下，采用常值推力方法，需要的軌道轉(zhuǎn)移時(shí)間和燃料消耗是相當(dāng)大的。

圖9 需要的推力和太陽帆提供的推力 圖10 太陽帆姿態(tài)角圖11 切向推力轉(zhuǎn)移軌道(與電推進(jìn)相同推力大小)

圖12 電推進(jìn)產(chǎn)生的推力加速度 圖13 電推進(jìn)需要的秒流量圖14 電推進(jìn)方向角

半長軸(AU)偏心率真近點(diǎn)角/rad真近點(diǎn)角/rad 初始值1000 目標(biāo)值1.501.484(第一組)3.108(第二組)

表4 虛擬引力場參數(shù)和切向推力參數(shù)

表5 混合推力方法和切向推力力方法對比

4 結(jié) 論

本文提出了一種基于虛擬引力場的混合小推力設(shè)計(jì)方法。該方法能夠以解析的形式求出大量的滿足軌道邊界約束的轉(zhuǎn)移軌道，以及實(shí)現(xiàn)軌道轉(zhuǎn)移所需的推力。同時(shí)，在給定目標(biāo)函數(shù)的條件下，可對滿足推力約束的轉(zhuǎn)移軌道進(jìn)行優(yōu)化。與其它的小推力軌道設(shè)計(jì)方法相比，該方法有以下的特點(diǎn)：

1) 該軌道設(shè)計(jì)方法是一種解析方法，在保證精度的條件下，能夠減小計(jì)算量，因而該方法可以為數(shù)值優(yōu)化方法(直接法或間接法)提供軌道優(yōu)化初值；

2) 與常值切向推力相比，該方法能夠減小實(shí)現(xiàn)軌道轉(zhuǎn)移所需的推力，擴(kuò)大了小推力軌道應(yīng)用的范圍；

3)與常值切向推力相比，該方法不僅能夠減小軌道轉(zhuǎn)移所消耗的能量，而且能夠減小軌道轉(zhuǎn)移時(shí)間。

綜上所述，相比于傳統(tǒng)的小推力設(shè)計(jì)方法，虛擬引力場混合小推力方法，不僅能夠節(jié)省燃料消耗，而且能減小任務(wù)時(shí)間，并擴(kuò)大了小推力方法的應(yīng)用范圍。為小推力軌道設(shè)計(jì)與優(yōu)化領(lǐng)域提供了新的思路。