基于實(shí)測(cè)數(shù)據(jù)的大氣密度反演方法及應(yīng)用研究

黃美麗，王大軼，馮 昊，田百義

(北京空間飛行器總體設(shè)計(jì)部，北京 100094)

0 引 言

對(duì)于低軌衛(wèi)星，軌道高度受大氣影響將不斷降低，并且軌道高度越低，大氣密度越大，軌道衰減速度越快[1]。因此，衛(wèi)星需要定期消耗推進(jìn)劑進(jìn)行軌道控制以維持軌道高度在一定的范圍內(nèi)。軌道維持的推進(jìn)劑預(yù)算是衛(wèi)星軌道設(shè)計(jì)的一個(gè)重要組成部分。

衛(wèi)星軌道維持的推進(jìn)劑消耗量由大氣密度及衛(wèi)星迎風(fēng)面質(zhì)比決定，衛(wèi)星總體構(gòu)型以及軌道確定后，迎風(fēng)面質(zhì)比也就確定了，軌道維持的推進(jìn)劑消耗量取決于大氣密度[2]。

大氣密度是一個(gè)復(fù)雜的多元變量函數(shù)，因此預(yù)測(cè)難度極大[3]。幾十年來(lái)，大氣模型不斷在改進(jìn)和發(fā)展，但由于高層大氣的變化極其復(fù)雜，模型仍存在約15%～30%的誤差，在空間環(huán)境擾動(dòng)期間甚至可達(dá)100%或更高。因此，進(jìn)一步提高大氣模型預(yù)報(bào)精度成為大氣模型研究需要解決的重要問(wèn)題。目前，大氣模型預(yù)報(bào)精度的改進(jìn)主要有兩種方式：一是建立新模型，增加更多更新的探測(cè)數(shù)據(jù)或者對(duì)模型的控制變量進(jìn)行改進(jìn)；二是在原模型基礎(chǔ)上進(jìn)行結(jié)果修正。無(wú)論是哪種改進(jìn)方法，其基礎(chǔ)都是要獲取大量準(zhǔn)確的大氣密度觀測(cè)數(shù)據(jù)[4]。

衛(wèi)星的兩行根數(shù)(TLE)、我國(guó)在軌衛(wèi)星定軌數(shù)據(jù)，都是實(shí)測(cè)軌道數(shù)據(jù)，將一段時(shí)間內(nèi)的軌道數(shù)據(jù)刨除變軌因素，衛(wèi)星軌道參數(shù)的變化反映了衛(wèi)星所受攝動(dòng)力的情況[5]。把其中與大氣阻力無(wú)關(guān)的攝動(dòng)因素去除之后，得到了衛(wèi)星受大氣阻力攝動(dòng)的軌道變化情況，再結(jié)合衛(wèi)星的外形結(jié)構(gòu)、阻力系數(shù)、姿態(tài)變化等因素，可對(duì)大氣密度進(jìn)行反演，得到近實(shí)時(shí)的大氣密度。

近年來(lái)，國(guó)內(nèi)外在由TLE數(shù)據(jù)進(jìn)行大氣密度反演以及在此基礎(chǔ)上的大氣密度模型修正等領(lǐng)域進(jìn)行了不少研究，但多注重于反演方法的推導(dǎo)及個(gè)例的應(yīng)用，缺乏大量數(shù)據(jù)反演獲得中長(zhǎng)期平均大氣密度的研究，并且尚無(wú)將反演方法應(yīng)用于衛(wèi)星推進(jìn)劑預(yù)算的研究[6-12]。

本文針對(duì)我國(guó)典型遙感衛(wèi)星運(yùn)行軌道開(kāi)展中長(zhǎng)期的平均大氣密度反演研究，并利用反演結(jié)果進(jìn)行大氣密度局部修正以及推進(jìn)劑預(yù)算應(yīng)用研究，達(dá)到提高推進(jìn)劑預(yù)算準(zhǔn)確度、在保障衛(wèi)星安全運(yùn)行基礎(chǔ)上最大限度減少推進(jìn)劑預(yù)算的目的。

1 中長(zhǎng)期大氣密度反演算法研究

1.1 低軌衛(wèi)星軌道高度衰減分析

對(duì)于低軌衛(wèi)星，大氣阻力造成軌道機(jī)械能的長(zhǎng)期衰減。對(duì)于一般的低地球軌道，大氣阻力產(chǎn)生與衛(wèi)星速度相反的阻力，對(duì)應(yīng)的阻力加速度大小為：

(1)

式中：fD為大氣阻力加速度，單位km/d2；CD為阻力系數(shù)，對(duì)于一般衛(wèi)星其取值范圍為2.0～2.3；S為衛(wèi)星有效迎風(fēng)面積，單位km2；m為衛(wèi)星質(zhì)量，單位kg；ρ為衛(wèi)星飛行高度處的大氣密度，單位kg/km3；V為衛(wèi)星速度，單位km/d。

大氣阻力沿速度的反向，將造成軌道面內(nèi)參數(shù)半長(zhǎng)軸a，偏心率e,近地點(diǎn)幅角ω和平近點(diǎn)角M的變化，對(duì)應(yīng)的攝動(dòng)方程如下：

(2)

式中：θ為衛(wèi)星真近點(diǎn)角，單位rad；E為衛(wèi)星偏近點(diǎn)角，單位rad；r為衛(wèi)星地心距，單位km；p=a(1-e2)為軌道半通徑，單位km；n為衛(wèi)星軌道運(yùn)動(dòng)平均角速度，單位rad/d；R,T為大氣阻力加速度的徑向和橫向分量。

(3)

對(duì)于橢圓的低地球軌道，其高度的衰減涉及了半長(zhǎng)軸、遠(yuǎn)地點(diǎn)高度、近地點(diǎn)高度和偏心率的衰減，由于衛(wèi)星在橢圓軌道上不同位置的大氣密度也是變化的，導(dǎo)致軌道衰減方程(2)和近似公式(3)不再適用，橢圓軌道各參數(shù)的衰減情況可參考文獻(xiàn)[12]。

1.2 大氣密度反演算法分析

1)積分反演算法

大氣密度反演基于實(shí)測(cè)的衛(wèi)星軌道平根數(shù)變化，采用帶有大氣攝動(dòng)的軌道動(dòng)力學(xué)模型(見(jiàn)式(2)),對(duì)衛(wèi)星所在高度中長(zhǎng)期的等效大氣密度反演，給出大氣密度反演算法理論和流程。該方法適用于近圓軌道中長(zhǎng)期的大氣密度反演。

衛(wèi)星軌道數(shù)據(jù)包括地心J2000坐標(biāo)系下的軌道瞬根和平根數(shù)，軌道參數(shù)包含軌道半長(zhǎng)軸、偏心率、傾角、升交點(diǎn)赤經(jīng)、近地點(diǎn)幅角和平近點(diǎn)角。平根數(shù)中的軌道半長(zhǎng)軸消去了地球非球形攝動(dòng)的短周期和長(zhǎng)周期攝動(dòng)項(xiàng)，主要受大氣阻力攝動(dòng)影響。主要表現(xiàn)是半長(zhǎng)軸在軌逐漸衰減，可根據(jù)半長(zhǎng)軸衰減率反推衛(wèi)星當(dāng)時(shí)所在高度的等效大氣密度情況?；谲壍烂鎯?nèi)參數(shù)攝動(dòng)方程，大氣密度反演方案可按圖1進(jìn)行。

詳細(xì)算法流程描述如下：

(2)對(duì)攝動(dòng)方程(2)積分，可得到衛(wèi)星運(yùn)行ΔT時(shí)間后的平軌道半長(zhǎng)軸at。

2)解析反演算法

(4)

以495 km軌道為例，采用解析反演算法得到的大氣密度和積分反演算法相比，誤差約0.05%，兩種方法相差極小。且解析反演算法與積分反演算法相比，可極大節(jié)約計(jì)算時(shí)間，因此，本文采用解析算法進(jìn)行大氣密度的反演。

上述方法用于反演衛(wèi)星中長(zhǎng)期的大氣密度情況，需對(duì)長(zhǎng)期實(shí)測(cè)的軌道數(shù)據(jù)進(jìn)行簡(jiǎn)單處理，獲取衛(wèi)星在某段時(shí)期內(nèi)的平均衰減率，以便對(duì)該段時(shí)間的平均大氣密度進(jìn)行反演。

2 典型軌道中長(zhǎng)期大氣密度反演分析

2.1 典型軌道選取

迄今為止，我國(guó)已成功發(fā)射了六發(fā)遙感衛(wèi)星(ZY-2(01)、ZY-2(02)、ZY-2(03)、YAOGAN5、YAOGAN12和YAOGAN21),運(yùn)行于495km軌道上，并且有多發(fā)衛(wèi)星在495km附近的軌道上運(yùn)行。未來(lái)，我國(guó)還將往該軌道附近陸續(xù)發(fā)射多顆遙感衛(wèi)星。該軌道上持續(xù)有衛(wèi)星運(yùn)行已超過(guò)11年的一個(gè)太陽(yáng)活動(dòng)周期，進(jìn)行大氣密度反演對(duì)后續(xù)的軌道維持推進(jìn)劑預(yù)算是非常有參考意義的。

因此，本文選擇典型衛(wèi)星高度495km進(jìn)行大氣密度反演，并與MSIS大氣密度模式計(jì)算的大氣密度進(jìn)行對(duì)比，驗(yàn)證本文的反演模型和算法。大氣密度計(jì)算主要取決于太陽(yáng)活動(dòng)率，太陽(yáng)活動(dòng)率以波長(zhǎng)10.7的射電流量F10.7表示。MSIS00模式計(jì)算大氣密度時(shí)采用實(shí)測(cè)的F10.7數(shù)據(jù)，根據(jù)美國(guó)國(guó)家海洋和大氣管理局(NOAA)發(fā)布的實(shí)測(cè)F10.7數(shù)據(jù)，2000年至2016年，實(shí)測(cè)F10.7數(shù)據(jù)如圖2所示。

大氣密度主要依賴(lài)于太陽(yáng)活動(dòng)水平F10.7，F(xiàn)10.7既有約為11年的長(zhǎng)周期變化，又有周期約為27天的短周期變化。對(duì)于衛(wèi)星推進(jìn)劑預(yù)算，只需考慮長(zhǎng)周期變化即可。

2.2 中長(zhǎng)期大氣密度反演結(jié)果

以ZY-2(01)、ZY-2(02)、ZY-2(03)、YAOGAN5、YAOGAN12和YAOGAN21衛(wèi)星為例，衛(wèi)星平均軌道高度均為495km。

以YAOGAN12和YAOGAN21衛(wèi)星為例，2015年期間兩顆星軌道半長(zhǎng)軸衰減情況如圖3所示。

由圖3可知，兩顆星在2015年期間進(jìn)行了5次正常的軌道維持。雙星軌道高度、整星質(zhì)量和有效迎風(fēng)面積較為一致，因此，雙星軌道半長(zhǎng)軸衰減也較為一致。

大氣密度反演過(guò)程中，ΔT的取值應(yīng)小于衛(wèi)星軌道維持周期。ΔT取值越大，所得大氣密度越接近ΔT時(shí)間內(nèi)的平均密度，對(duì)獲取每天的大氣密度和F10.7的關(guān)系不利；而取值越小，對(duì)抑制軌道半長(zhǎng)軸誤差不利，所得大氣密度受軌道定軌誤差影響較大。綜上考慮，495 km軌道大氣密度反演時(shí)，取ΔT=5天。

基于實(shí)測(cè)F10.7計(jì)算所得的MSIS00大氣模式密度和本文反演算法所得的2010年至2016年反演密度情況見(jiàn)圖4。

本文方法所得大氣密度與MSIS大氣密度模式相比，得到如下結(jié)論：

1)MSIS00大氣密度模式與反演大氣密度的變化趨勢(shì)吻合度較高。

2)在太陽(yáng)活動(dòng)高年，MSIS00模型誤差較大，誤差最大可達(dá)60%。

3)在太陽(yáng)活動(dòng)低年，MSIS00模型誤差較小，模型精度可達(dá)10%。

3 大氣密度模型局部修正方法研究

3.1 大氣密度模型局部修正方法

本節(jié)以MSIS00大氣密度模式為例，進(jìn)行大氣模型局部修正方法研究。

通過(guò)在軌數(shù)據(jù)進(jìn)行大氣密度反演，獲得一段時(shí)間內(nèi)的大氣密度數(shù)據(jù)：

ρs=ρs(t)

(5)

通過(guò)大氣密度模型計(jì)算一段時(shí)間內(nèi)的大氣密度數(shù)據(jù)，計(jì)算大氣密度時(shí)采用實(shí)測(cè)的F10.7數(shù)據(jù)。

ρm=ρm(t,F10.7)

(6)

從而得到同一時(shí)間的實(shí)測(cè)大氣密度與模型密度之間的關(guān)系系數(shù)

(7)

利用衛(wèi)星長(zhǎng)期在軌的數(shù)據(jù)，通過(guò)多項(xiàng)式曲線擬合，得到兩種密度之間的關(guān)系系數(shù)ρx與F10.7的關(guān)系：

ρx(F10.7)=A·f(F10.7)

(8)

式中：A為多項(xiàng)式系數(shù)向量。

通過(guò)上述關(guān)系系數(shù)與F10.7的擬合表達(dá)式，可以對(duì)大氣密度模型進(jìn)行修正，從而再用于工程設(shè)計(jì)中，優(yōu)化衛(wèi)星推進(jìn)劑預(yù)算及總體設(shè)計(jì)。

3.2 典型軌道大氣密度模型局部修正

對(duì)典型軌道高度495 km，通過(guò)式(5)～(8)可獲取大氣密度模型修正系數(shù)同F(xiàn)10.7的關(guān)系如圖5所示。

采用5次多項(xiàng)式擬合得出模型修正系數(shù)同F(xiàn)10.7的關(guān)系見(jiàn)下式：

(9)

其中多項(xiàng)式系數(shù)向量A取值如下：

A= [5.538039×10-11,-5.496297×10-8,

2.098230×10-5,-0.003817,0.322250,

-8.822965]

(10)

模型修正前后和反演的大氣密度情況見(jiàn)圖6。

由圖6可知：

1)在太陽(yáng)活動(dòng)高年，模型修正后的大氣密度平均值同反演的大氣密度平均值相比，最大誤差約20%，與模型修正前的50%誤差相比，模型精度有了較大提高。

2)在太陽(yáng)活動(dòng)低年，修正前后的大氣密度模型同反演的大氣密度均有較高的契合度。

3)通過(guò)系數(shù)修正主要提高了太陽(yáng)活動(dòng)高年的大氣密度模型精度。

4 大氣密度模型局部修正在推進(jìn)劑預(yù)算中的應(yīng)用

4.1 低軌衛(wèi)星推進(jìn)劑預(yù)算

低軌衛(wèi)星推進(jìn)劑預(yù)算主要包括對(duì)初軌調(diào)整、軌道機(jī)動(dòng)、軌道維持以及姿態(tài)控制等消耗的推進(jìn)劑進(jìn)行預(yù)算并留有一定的余量。其中，對(duì)大氣阻力引起的軌道衰減進(jìn)行軌道維持所需的推進(jìn)劑預(yù)算是重點(diǎn)及難點(diǎn)，這主要因?yàn)槠渌麕醉?xiàng)推進(jìn)劑預(yù)算一般情況下均可準(zhǔn)確預(yù)估，而軌道維持所需推進(jìn)劑則與大氣密度預(yù)測(cè)密切相關(guān)，而大氣密度預(yù)測(cè)的不確定性很大，是目前整體推進(jìn)劑預(yù)算的難點(diǎn)。

傳統(tǒng)的軌道保持推進(jìn)劑預(yù)算方法主要有兩種：

1)對(duì)于壽命3年左右的衛(wèi)星，一般都按太陽(yáng)活動(dòng)中高年(F10.7=175)的情況，根據(jù)大氣密度模型預(yù)報(bào)的大氣密度進(jìn)行推進(jìn)劑預(yù)算，推進(jìn)劑余量較大。

2)對(duì)于壽命5年以上的衛(wèi)星，一般按一個(gè)太陽(yáng)活動(dòng)周期中高、中、低年組合的形式，根據(jù)預(yù)計(jì)運(yùn)行年份，以大氣密度模型預(yù)報(bào)的大氣密度進(jìn)行推進(jìn)劑預(yù)算，余量相對(duì)較小。

可以看出，兩種方法都受限于大氣密度模型預(yù)報(bào)精度。對(duì)于典型軌道，在具備長(zhǎng)期(一個(gè)太陽(yáng)活動(dòng)周期11年以上)大氣密度反演數(shù)據(jù)的前提下，可利用反演結(jié)果對(duì)大氣密度模型進(jìn)行局部修正，提高大氣模型預(yù)報(bào)精度，從而提高推進(jìn)劑預(yù)算精度。

4.2 應(yīng)用算例及結(jié)果分析

本文以495 km軌道高度上的5年壽命衛(wèi)星推進(jìn)劑預(yù)算作為應(yīng)用算例進(jìn)行大氣密度模型局部修正方法應(yīng)用研究。

根據(jù)前文分析，對(duì)于5年壽命衛(wèi)星，按太陽(yáng)活動(dòng)中高年的情況，根據(jù)大氣密度模型預(yù)報(bào)的大氣密度進(jìn)行推進(jìn)劑預(yù)算。取衛(wèi)星的迎風(fēng)面質(zhì)比為0.008 m2/kg，分析可得大氣密度模型局部修正前后的軌道維持推進(jìn)劑預(yù)算比對(duì)如表1所示[10]。

由表1可知，軌道保持推進(jìn)劑消耗由修正前的148 kg減少為110 kg，比原設(shè)計(jì)結(jié)果節(jié)省38 kg推進(jìn)劑，節(jié)省了近35%。

表1 大氣密度模型局部修正前后軌道維持推進(jìn)劑預(yù)算比對(duì)Table 1 Propellant budget before modified model and after

5 結(jié)束語(yǔ)

本文基于實(shí)測(cè)軌道數(shù)據(jù)進(jìn)行大氣密度反演方法、大氣模型局部修正方法及工程應(yīng)用研究，研究成果具有創(chuàng)新性以及工程實(shí)用性，可在我國(guó)未來(lái)航天器的高精度軌道預(yù)報(bào)、推進(jìn)劑預(yù)算、碰撞規(guī)避等任務(wù)中得以推廣和應(yīng)用，可以為我國(guó)未來(lái)航天器的總體設(shè)計(jì)等提供技術(shù)支持，促進(jìn)形成合理、可行、優(yōu)化的設(shè)計(jì)方案。

需要指出的是，本文采用的面質(zhì)比計(jì)算方法適用于合作目標(biāo)，對(duì)于非合作目標(biāo)，測(cè)軌參數(shù)里包含的B*參數(shù)是一個(gè)綜合了阻力系數(shù)、面質(zhì)比和大氣密度的綜合氣動(dòng)參數(shù)，后續(xù)如何根據(jù)該參數(shù)獲取大氣密度將進(jìn)行深入研究。