周姜濱 葉 松 張華明

北京航天自動(dòng)控制研究所，北京 100854

空間飛行器因?yàn)樾羌蛛x或長時(shí)間飛行等原因會(huì)產(chǎn)生一定的軌道偏差，為了能夠順利完成后續(xù)任務(wù)，需要進(jìn)行軌道修正來消除這一偏差[1]。經(jīng)典的Hohmann轉(zhuǎn)移、雙橢圓轉(zhuǎn)移及蘭伯特變軌等軌道轉(zhuǎn)移與修正理論和方法[2-4]都要求滿足二體假設(shè)條件，并假設(shè)發(fā)動(dòng)機(jī)推力充分大，在瞬間就能獲得所需要的速度增量，因此難以直接應(yīng)用于工程實(shí)踐之中。田野[5]等提出了基于軌道方程線性化的中段飛行軌道修正方法，采用對(duì)線性化方程增加J2修正項(xiàng)的方法，減小軌道攝動(dòng)引起的終端偏差。任金磊[6]等提出了基于速度增益制導(dǎo)的大橢圓轉(zhuǎn)移中段制導(dǎo)方法，通過優(yōu)化修正時(shí)間設(shè)計(jì)了終端約束小、燃料消耗少的中段修正優(yōu)化軌道。

提出了一種考慮攝動(dòng)影響的連續(xù)推力軌道修正制導(dǎo)方法，考慮地球非球形引力及大氣阻力、日月引力和太陽光壓等攝動(dòng)因素，將軌道修正問題轉(zhuǎn)換為二體假設(shè)條件下的蘭伯特變軌問題，通過迭代計(jì)算給出連續(xù)推力軌道修正制導(dǎo)指令，由地面站將制導(dǎo)指令上行至空間飛行器，空間飛行器接收并執(zhí)行該指令，實(shí)現(xiàn)軌道修正，并開展數(shù)學(xué)仿真驗(yàn)證。

1 蘭伯特變軌問題

蘭伯特變軌可以描述為：給定空間飛行器初始時(shí)刻的位置矢量rc0和速度矢量vc0，目標(biāo)點(diǎn)位置矢量rtf和速度矢量vtf，以及軌道轉(zhuǎn)移時(shí)間tf，確定作用在空間飛行器上的2次速度增量Δv1和Δv2。求解蘭伯特變軌問題其實(shí)就是求解Gauss問題。Gauss問題的定義如下：給定2個(gè)位置矢量r1和r2，以及空間飛行器從r1運(yùn)行到r2的飛行時(shí)間t和初始運(yùn)動(dòng)方向，求解空間飛行器在r1和r2處的速度矢量v1和v2。顯然，求解了Gauss問題，就解決了蘭伯特變軌問題[4]。

蘭伯特變軌過程如圖1所示，空間飛行器在初始時(shí)刻應(yīng)具有的速度為

圖1 蘭伯特變軌

到達(dá)目標(biāo)點(diǎn)時(shí)刻的速度為

式中，a為橢圓轉(zhuǎn)移軌道的半長軸；θ為rc0至rtf的夾角。

橢圓轉(zhuǎn)移軌道的半長軸a可表示為

其中，

式中，c為空間飛行器初始時(shí)刻的位置至目標(biāo)點(diǎn)的距離；λ是蘭伯特參數(shù)，為迭代變量，其取值范圍為-π<λ<π。

蘭伯特變軌所需要的時(shí)間為

其中，

蘭伯特變軌所需要的2個(gè)速度增量分別為

Δv1=v0-vc0
Δv2=vtf-vf

在計(jì)算中需先對(duì)蘭伯特參數(shù)λ迭代求解，得到最省燃料消耗時(shí)的λopt，然后計(jì)算出橢圓轉(zhuǎn)移軌道的半長軸a。

2 高精度軌道預(yù)報(bào)

空間飛行器軌道預(yù)報(bào)通常由解析法(如擬平均根數(shù)法)和數(shù)值積分法(如Cowell法)，隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值積分法得到廣泛應(yīng)用。用J2000.0地心慣性坐標(biāo)系中的位置和速度描述空間飛行器狀態(tài)，則其動(dòng)力學(xué)模型的狀態(tài)向量微分方程可寫成[5]

其中，ax,ay,az為空間飛行器的加速度，有

式中，aC,aE,aD,aS,aL,aR分別表示控制力、地球引力、大氣阻力、日月引力和太陽光壓產(chǎn)生的加速度；aO為其他攝動(dòng)加速度，此處將其忽略不計(jì)?？刂屏τ芍茖?dǎo)系統(tǒng)給出，其余可建立準(zhǔn)確的數(shù)學(xué)模型，利用四階龍格-庫塔積分方法對(duì)上述微分方程進(jìn)行數(shù)值積分計(jì)算，得到空間飛行器的位置速度，完成軌道預(yù)報(bào)。

下面分別對(duì)地球引力、大氣阻力、日月引力和太陽光壓加速度進(jìn)行數(shù)學(xué)建模。

1)地球引力加速度

用完全正?；那蛑C函數(shù)表示的地球引力場(chǎng)位函數(shù)V可寫成

根據(jù)地球引力場(chǎng)位函數(shù)V即可得到地球引力加速度

2)大氣阻力加速度

大氣阻力加速度計(jì)算公式為

其中，CD為阻力系數(shù)，通常取2.2；SD為參考面積，m為質(zhì)量；ρ為大氣密度；v為相對(duì)于大氣的速度。

3)太陽引力加速度

太陽引力計(jì)算公式為

其中，GMS為太陽引力常數(shù)。

4)月球引力加速度

同太陽引力加速度，月球引力計(jì)算公式為

其中，GML為月球引力常數(shù)。

5)太陽光壓加速度

太陽光壓計(jì)算公式為

其中，CR為表面反射系數(shù)；SR/m為太陽光壓面質(zhì)比；ρSR為作用在距離太陽一個(gè)天文單位處單位面積黑體上的光壓力，一般取4.560×10-6N/m2。

3 考慮攝動(dòng)影響的連續(xù)推力軌道修正制導(dǎo)方法

蘭伯特變軌利用2個(gè)速度脈沖實(shí)現(xiàn)了在二體假設(shè)條件下的軌道修正，而在實(shí)際工程應(yīng)用中用于軌道修正的推力往往較小，需要連續(xù)作用，工作時(shí)間較長，不能近似為速度脈沖，另一方面，空間飛行器還會(huì)受到包括地球非球形引力、大氣阻力、日月引力和太陽光壓在內(nèi)的各種攝動(dòng)力的作用，并不滿足理想的二體假設(shè)。為此，提出考慮攝動(dòng)影響的連續(xù)推力軌道修正制導(dǎo)方法，具體流程如圖2所示，相應(yīng)描述如下：

1)將軌道修正問題轉(zhuǎn)換為二體假設(shè)條件下的蘭伯特變軌問題，根據(jù)初始位置矢量、(虛擬)目標(biāo)點(diǎn)位置矢量和飛行時(shí)間求解所需要的速度脈沖；

2)將該速度脈沖等效為連續(xù)推力；

3)考慮地球非球形引力、大氣阻力、日月引力和太陽光壓等攝動(dòng)因素，計(jì)算在連續(xù)推力作用下空間飛行器的實(shí)際位置；

4)利用實(shí)際位置與目標(biāo)點(diǎn)位置的偏差構(gòu)造新虛擬目標(biāo)點(diǎn)，即虛擬目標(biāo)點(diǎn)位置矢量減實(shí)際位置與目標(biāo)點(diǎn)位置的偏差得到新虛擬目標(biāo)點(diǎn)的位置矢量；

5)重復(fù)步驟1)～4)，計(jì)算至實(shí)際位置與目標(biāo)點(diǎn)位置的偏差在允許范圍內(nèi)，給出連續(xù)推力軌道修正制導(dǎo)指令。

圖2 考慮攝動(dòng)影響的連續(xù)推力軌道修正制導(dǎo)流程

4 數(shù)學(xué)仿真與分析

為驗(yàn)證提出的連續(xù)推力軌道修正制導(dǎo)方法的有效性，開展了相應(yīng)的數(shù)學(xué)仿真驗(yàn)證，圖3～7給出了某算例500次蒙特卡洛打靶數(shù)學(xué)仿真的計(jì)算結(jié)果。從仿真結(jié)果可以看出該制導(dǎo)算法具有較高的制導(dǎo)精度，半長軸偏差的平均值為-4.529m，偏心率偏差的平均值為2.905×10-7，升交點(diǎn)赤經(jīng)偏差的平均值為-1.312×10-5(°)，軌道傾角偏差的平均值為-1.460×10-8(°)，近地點(diǎn)幅角偏差的平均值為-1.671×10-2(°)。

圖3 半長軸偏差

圖4 偏心率偏差

圖5 升交點(diǎn)赤經(jīng)偏差

圖6 軌道傾角偏差

圖7 近地點(diǎn)幅角偏差

5 結(jié)論

為了消除在軌空間飛行器軌道偏差，本文從工程應(yīng)用角度著眼提出了一種考慮攝動(dòng)影響的連續(xù)推力軌道修正制導(dǎo)方法，該方法易于實(shí)現(xiàn)，仿真表明制導(dǎo)精度較高，能夠以上行制導(dǎo)指令的方式實(shí)現(xiàn)連續(xù)推力軌道修正，具有一定的實(shí)用價(jià)值。