基于航天器探測密度與軌道數(shù)據(jù)的密度模式精度評估

李 勰，滿海鈞，劉舒蒔，陳光明

(1.中國科學(xué)院國家空間科學(xué)中心，北京100190；2.北京航天飛行控制中心，北京100094；3.航天飛行動力學(xué)技術(shù)重點實驗室，北京100094；4.中國科學(xué)院大學(xué)，北京100049)

1 引言

熱層大氣密度的不確定性是影響低軌航天器定軌預(yù)報精度的最關(guān)鍵因素，而描述熱層大氣密度變化的經(jīng)驗?zāi)Ｊ绞芙?shù)據(jù)以及對于密度變化自身物理機制認(rèn)識的局限，雖歷經(jīng)半個世紀(jì)的發(fā)展，依然無法突破15%的精度瓶頸。經(jīng)驗密度模式依據(jù)建立時采用的數(shù)據(jù)源，主要分為兩大類:①基于衛(wèi)星阻力數(shù)據(jù)，典型代表是美國史密松天文臺發(fā)布的Jacchia系列和法國宇航局發(fā)布的DTM系列；②基于地面非相干散射雷達(dá)和衛(wèi)星原位探測數(shù)據(jù)，典型代表是美國的MSIS系列。上述3個系列的經(jīng)驗密度模式在國內(nèi)外航天工程任務(wù)中得到了廣泛應(yīng)用。

進(jìn)入21世紀(jì)以來，隨著各類工程、科學(xué)和軍事應(yīng)用對航天器軌道精度的要求越來越高，傳統(tǒng)的經(jīng)驗密度模式成為制約精度提升的一個關(guān)鍵因素。為此在原有模式基礎(chǔ)上，通過引入新的觀測數(shù)據(jù)和輻射指數(shù)，對部分方程進(jìn)行修改，美國和歐洲相繼對上述不同模式進(jìn)行了升級。Picone等在MSIS90基礎(chǔ)上融合新的觀測數(shù)據(jù)，并考慮500 km以上不規(guī)則原子氧和離子氧的影響，開發(fā)了NRLMSISE00模式，成為當(dāng)前科學(xué)和工程應(yīng)用領(lǐng)域使用最廣泛的模式；Bowman等在Jacchia71基礎(chǔ)上，引入新的太陽輻射指數(shù)和地磁指數(shù)，建立了Jacchia-Bowman模式，作為北美防空司令部軌道計算的標(biāo)準(zhǔn)模式；法國宇航局在歐盟第七研究計劃支持下，通過融合GRACE、GOCE衛(wèi)星數(shù)據(jù)，推出了DTM2013模式，并宣稱該模式是目前誤差最小的模式。

由于密度模式精度對軌道計算，尤其是軌道預(yù)報影響顯著，因此模式精度的評估往往是工程應(yīng)用部門關(guān)注的焦點。Marcos分別對Jacchia70/71/77、MSIS77/79/83/86 和 Jacchia71、DTM2000、NRLMSISE00和JB2006進(jìn)行了系統(tǒng)分析比較，Vallado從應(yīng)用層面比較了Jacchia60/70/71、MSIS86/90和NRLMSISE00等模式。由于驅(qū)動JB2008模式輸入?yún)?shù)的發(fā)布存在約3個月延遲，且最新的DTM2013模式不對外公開發(fā)布，因此國內(nèi)工程領(lǐng)域?qū)嵱玫拿芏饶Ｊ饺允菙?shù)十年前的產(chǎn)品，除汪宏波等利用CAHMP衛(wèi)星加速度計數(shù)據(jù)分析了不同輻射指數(shù)對Jacchia71、DTM94、NRLMSISE00和JB2006模式精度的影響外，國內(nèi)對于新近發(fā)布模式的精度及使用建議關(guān)注較少。

中國在熱層大氣建模領(lǐng)域尚沒有自主的經(jīng)驗?zāi)Ｊ剑矝]有類似于馬歇爾工程熱層模型(Marshall Engineering Thermosphere Model，MET)的工程專用模式，因此，在實際工程應(yīng)用前如何評估模式精度及計算效率，并據(jù)此選擇合適的經(jīng)驗?zāi)Ｊ骄哂兄匾睦碚摵同F(xiàn)實意義。本文利用CHAMP衛(wèi)星的歷史數(shù)據(jù)和中國自主衛(wèi)星探測數(shù)據(jù)，結(jié)合不同類型航天器的軌道預(yù)報效果，綜合評估Jacchia-Robert(JR)、NRLMSISE00(MSISE00)、JB2008和DTM2013 4種不同模式的精度，分析它們在不同應(yīng)用場景下的優(yōu)勢和特點，為后續(xù)的工程應(yīng)用以及關(guān)聯(lián)的空間環(huán)境探測提供參考。

2 評估數(shù)據(jù)描述

綜合考慮數(shù)據(jù)的精度和可獲得性，采用2類數(shù)據(jù)作為基準(zhǔn)開展模式的精度評估:①使用廣泛且精度得到業(yè)內(nèi)認(rèn)可的CHAMP衛(wèi)星數(shù)據(jù)；②中國自主熱層大氣密度探測與精密定軌試驗衛(wèi)星(Atmospheric density detecting and Precise Orbit Determination，APOD)的觀測數(shù)據(jù)。上述數(shù)據(jù)覆蓋了400～500 km高度的大部分區(qū)域，也包含了不同太陽輻射水平和地磁條件，具有較好的代表性。

2.1 CHAMP衛(wèi)星反演數(shù)據(jù)

CHAMP衛(wèi)星是德國地學(xué)研究中心研制的重力測地衛(wèi)星。由于CHAMP衛(wèi)星的軌道特性，5個月左右就能覆蓋全部緯度范圍內(nèi)的所有地方，有利于考察模式誤差的全球分布特征。利用其搭載的高精度加速度計，Bruinsma等和Sutton等均開展了密度反演工作。本文采用Sutton反演的密度數(shù)據(jù)開展工作(密度下載地址為http://sisko.colorado.edu/suttion/data.htm l)。數(shù)據(jù)時間區(qū)間為2002～2008年，分辨率為45 s，涵蓋第23太陽活動周的大部分時段。

2.2 APOD衛(wèi)星原位探測數(shù)據(jù)

APOD衛(wèi)星是中國首顆以大氣密度探測為主要目標(biāo)的微納衛(wèi)星，在國內(nèi)首次獲取了500 km高度以下的長期連續(xù)觀測數(shù)據(jù)。由于晨昏軌道特性，APOD衛(wèi)星不同緯度的地方時基本不變，數(shù)據(jù)主要集中在地方時早晚6點30分左右，這也是APOD衛(wèi)星有別于其他大氣密度探測衛(wèi)星的一個顯著特點。本文采用自2015年12月1日至2016年12月31日共計397天的觀測數(shù)據(jù)，數(shù)據(jù)采樣頻率為1 s。

3 基于實測密度的精度評估

3.1 評估方法

將上述實測結(jié)果作為密度真值，對4種不同模式計算，結(jié)果與真值進(jìn)行比較。定義變量R為模式計算值相對于觀測真值的誤差，如式(1)所示:

其中，N為數(shù)據(jù)個數(shù)。而模式相對真實密度數(shù)據(jù)偏差的彌散程度(隨機誤差)σ如式(3)所示:

3.2 CHAMP加速度計反演密度評估結(jié)果

按照CHAMP衛(wèi)星的軌跡和實測空間環(huán)境參數(shù)，采用JR、MSISE00、JB2008和DTM2013 4種不同模式計算，結(jié)果與實測密度進(jìn)行比較。由于數(shù)據(jù)量較大，為提高計算效率且不失統(tǒng)計特性，對數(shù)據(jù)進(jìn)行日平均降采樣處理，其相對誤差的均值和標(biāo)準(zhǔn)偏差如圖1和表1所示。

圖1 不同模式計算結(jié)果與CHAMP衛(wèi)星數(shù)據(jù)的比較Fig.1 Com parisons between differentmodels and CHAMP observations

表1 基于CHAMP數(shù)據(jù)的不同模式相對誤差統(tǒng)計結(jié)果Table 1 Statistical results of relative error for different models based on CHAMP

可以看出:JR與MSISE00模式存在比較明顯的非線性系統(tǒng)偏差，隨著時間推移，系統(tǒng)偏差明顯增加，2008年相對于2002年誤差擴(kuò)大超過1個數(shù)量級。圖2給出了數(shù)據(jù)時間對應(yīng)的不同波段太陽輻射流量，太陽輻射的整體水平呈逐年下降趨勢，2007和2008年處于一個比較低的水平。綜上，JR和MSIS00模式計算結(jié)果隨著太陽輻射的減弱而逐漸偏高；JR和MSIS00模式呈現(xiàn)出了非常明顯的周年/半年變化特征，且變化的幅度較大，成為顯著的誤差來源；而JB2008和DTM2013模式與觀測結(jié)果符合相對較好，相對誤差均值分別為0.4%和5.2%，除2007、2008以外，基本沒有明顯的系統(tǒng)偏差，其隨機偏差也均在10%左右，較JR和MSISE00模式降低約50%；JB2008模式由于引入了新的半年變化，其模式誤差基本沒有呈現(xiàn)周年變化的現(xiàn)象。

圖2 F10.7 cm太陽輻射流量指數(shù)Fig.2 F10.7 Solar radiation flux

3.3 原位探測密度評估結(jié)果

中國自主原位探測載荷主要是質(zhì)譜儀和真空計，相較于加速度計反演結(jié)果，原位探測數(shù)據(jù)精度略低。為了確保探測數(shù)據(jù)的準(zhǔn)確性和可靠性，李勰等提出了一種基于動力學(xué)的在軌標(biāo)校方法，標(biāo)校后數(shù)據(jù)的系統(tǒng)誤差優(yōu)于15%，因此作為評估模式精度的基準(zhǔn)數(shù)據(jù)是合理可行的。需要說明的是由于DTM2013模式?jīng)]有正式對外公開發(fā)布，驅(qū)動模式運行的數(shù)據(jù)更新截止到2014年，本文只給出了JR、JB2008和MSISE00的比較結(jié)果。圖3為3種模式計算日均值與觀測日均值的比較，表2為數(shù)據(jù)統(tǒng)計結(jié)果?？梢钥闯?JR和JB2008模式的系誤差均值分別為2.2%和-1.1%，系統(tǒng)偏差較小。而MSISE00模式誤差均值為12.4%，系統(tǒng)偏差略大。3個模式相對誤差的標(biāo)準(zhǔn)偏差基本相當(dāng)，在15%～22%之間。相比日均值結(jié)果可以發(fā)現(xiàn)，通過平均也可以有效降低數(shù)據(jù)的隨機誤差。

圖3 3種模式相對APOD觀測數(shù)據(jù)日均值誤差Fig.3 Daily average error of 3 different models vs APOD observations

表2 基于APOD數(shù)據(jù)的模式相對誤差統(tǒng)計結(jié)果Tab le 2 Statistical results of relative error for different models based on APOD observations

由于APOD數(shù)據(jù)覆蓋時間有限，選取與APOD數(shù)據(jù)時間(2016年)范圍內(nèi)太陽輻射流量(70～125 sfu)水平相當(dāng)?shù)?005年CHAMP衛(wèi)星數(shù)據(jù)的比較結(jié)果，可以看出不同模式基于2類數(shù)據(jù)的評估結(jié)論基本一致:JB2008的表現(xiàn)均為最優(yōu)，JR和MSISE00模式精度相當(dāng)，較JB2008而言，有比較明顯的半年變化的系統(tǒng)偏差。此外，JB2008還呈現(xiàn)除了與其他2個模式不同的誤差特性，即JB2008計算結(jié)果較觀測值整體偏低，而其他2個模式結(jié)果偏高。

針對APOD衛(wèi)星晨昏軌道特性，圖4給出了3種模式相對誤差的緯度、地方時分布，其中緯度按照5°，地方時按照1.5 h間隔進(jìn)行數(shù)據(jù)平均。從圖中可以看出:3種模式的主要誤差均集中在南北半球的中高緯地區(qū)，且北半球多表現(xiàn)為正向誤差，南半球表現(xiàn)為負(fù)向誤差。此外，晨側(cè)和昏側(cè)的誤差特性也截然相反，晨側(cè)模式結(jié)果低于觀測結(jié)果，而昏側(cè)模式結(jié)果普遍高于觀測結(jié)果，MSISE00模式在昏側(cè)的誤差尤為顯著。從統(tǒng)計數(shù)據(jù)來看，JB2008模式的誤差均值和標(biāo)準(zhǔn)偏差均最小，分別為-0.4%和18.7%(圖4)。

4 基于軌道的精度評估

將不同模式應(yīng)用于軌道計算，通過軌道確定和軌道預(yù)報的效果進(jìn)行分析比較是精度評估的另一種重要和有效手段，可以為工程應(yīng)用提供技術(shù)依據(jù)，具有較強的現(xiàn)實意義。

4.1 評估原理與方法

大氣阻尼加速度計算公式如式(4)所示:

圖4 3種模式相對誤差的緯度地方時分布Fig.4 Relative error of latitude vs local time for 3 differentmodels

軌道預(yù)報精度主要取決于動力學(xué)模型誤差。目前常用的地球重力場、大行星等保守力攝動模型本身精度較高，且不會產(chǎn)生長期變化；而作為非保守力的大氣阻力攝動，其誤差傳播隨時間平方發(fā)散，是預(yù)報誤差的主要來源。因此，可以通過比較軌道預(yù)報的誤差來評估所使用的密度模式的精度。

以航天器對外發(fā)布的精密軌道為基準(zhǔn)，采用不同模式進(jìn)行軌道計算，將計算軌道與基準(zhǔn)軌道進(jìn)行比較。對于數(shù)據(jù)弧段內(nèi)，通過比較解算阻尼系數(shù)評估不同模式的定軌精度差異；對于預(yù)報弧段，統(tǒng)計預(yù)報軌道與基準(zhǔn)軌道誤差的均值和最大值來評估各個模式的預(yù)報精度。評估采用的軌道數(shù)據(jù)分別為CHAMP衛(wèi)星2003年1月～6月的快速科學(xué)軌道(Rapid Scientific Orbit，RSO)數(shù)據(jù)和天宮一號2012年1月～4月的工程軌道數(shù)據(jù)，其標(biāo)稱精度分別為10 cm和10 m。定軌預(yù)報過程中，為了降低其他誤差因素的影響，空間環(huán)境參數(shù)均采用實測值，面質(zhì)比為固定值，阻尼系數(shù)采用定軌解算結(jié)果。

4.2 天宮一號目標(biāo)飛行器

天宮一號目標(biāo)飛行器是中國首個自主研制的載人空間試驗平臺，于2011年9月29日發(fā)射入軌。本文選取2012年1月～4月的測量數(shù)據(jù)，利用1天的弧段進(jìn)行軌道計算，分別采用JR、MSISE00、JB2008和DTM2013 4種模式求解大氣阻尼系數(shù)，結(jié)果如圖5所示。由于模式存在誤差，在軌道解算過程中，阻尼系數(shù)會吸收模式的誤差，確保解算軌道與觀測符合最好。從圖中可以看出，JR、MSISE00、JB2008和DTM2013 4種不同模式解算阻尼系數(shù)的均值分別為2.576、2.308、2.074和2.169，反映不同模式之間的系統(tǒng)誤差水平。由于天宮一號目標(biāo)飛行器的構(gòu)型比較復(fù)雜，其阻尼系數(shù)的真實值較難確定，因此無法直接給出哪個模式最優(yōu)，只能進(jìn)行相對比較。如果模式計算結(jié)果高于實際密度，為了彌補模式誤差對軌道的影響，解算阻尼系數(shù)將低于真實值，反之亦然。因此，從不同模式解算結(jié)果來看，DTM2013模式計算密度整體偏高，而JR模式計算密度整體偏低，2個模式平均偏差可達(dá)24.2%，而相對偏差最小的2個模式是DTM2013和JB2008，平均偏差為4.58%，這體現(xiàn)了新近開發(fā)模式與早期發(fā)布模式的差異。

以不同模式解算阻尼系數(shù)的平均值為基準(zhǔn)，其相對均值的相對偏差可以反映不同模式的隨機偏差，JR、MSISE00、JB2008和DTM2013 4種模式結(jié)果分別為23.7%、22.7%、21.1%和26.7%?？梢钥闯?，4種模式的隨機偏差大致相當(dāng)，在21%～26%之間，隨機偏差最小的是DTM2013模式。

進(jìn)一步利用解算軌道和阻尼系數(shù)進(jìn)行24 h軌道預(yù)報，將預(yù)報軌道與基準(zhǔn)軌道進(jìn)行比較，用以統(tǒng)計最大位置和速度誤差。

圖5 不同模式定軌解算大氣阻尼系數(shù)結(jié)果Fig.5 Estimation of drag coefficient using different models

為了更直觀反映預(yù)報效果，將位置誤差按照500 m間隔進(jìn)行頻度統(tǒng)計，總共121 d數(shù)據(jù)樣本，結(jié)果如圖6所示。統(tǒng)計結(jié)果表明:JR、MSISE00、JB2008和DTM2013 4種模式預(yù)報24 h位置偏差小于500 m的天數(shù)分別為22 d、38 d、36 d和22 d，MSISE00和DTM2013模式明顯優(yōu)于Jacchia系列；預(yù)報誤差小于1 km的天數(shù)分別為51 d、64 d、63 d和59 d，仍然是MSISE00和DTM2013模式較優(yōu)；而預(yù)報誤差大于3 km的天數(shù)分別為15 d、20 d、16 d和13 d，MSISE00模式效果略差。從天宮一號軌道預(yù)報效果來看，基于探測數(shù)據(jù)評估中表現(xiàn)不佳的MSISE00模式反而表現(xiàn)最好，說明實測密度數(shù)據(jù)反映的模式誤差不能完全體現(xiàn)在軌道預(yù)報誤差中。另一方面，由于模式誤差隨預(yù)報時間的平方發(fā)散，如果預(yù)報超過24 h，模式誤差的影響可能會更顯著增加。

4.3 基于CHAMP衛(wèi)星軌道的精度評估

利用CHAMP衛(wèi)星2003年1月至6月(共計181 d)的快速科學(xué)軌道進(jìn)行定軌，將解算軌道預(yù)報24 h與RSO軌道進(jìn)行比較，仍然按照500m的間隔統(tǒng)計最大位置誤差的頻度分布，結(jié)果見圖7。結(jié)果表明:絕大多數(shù)時間的軌道預(yù)報精度均優(yōu)于1.5 km，主要原因是CHAMP衛(wèi)星軌道高度較高，受大氣阻力的影響比天宮一號顯著減小，因此受密度模式誤差的影響也相應(yīng)減小。JR、MSISE00、DTM2013和JB2008模式預(yù)報24 h位置誤差優(yōu)于1 km的天數(shù)分別為139 d、149 d、148 d和138 d，仍然是MSISE00和DTM2013表現(xiàn)優(yōu)于其他2個模式，這與天宮一號軌道預(yù)報評估結(jié)果完全一致。

圖6 不同模式24 h軌道預(yù)報誤差頻率直方圖Fig.6 Frequency distribution of 24 h orbit p rediction error using differentmodels

圖7 不同模式24 h軌道預(yù)報誤差頻率直方圖(CHAMP)Fig.7 Frequency distribution of 24 h orbit prediction error using differentmodels(CHAMP)

2003年為太陽活動高年，為了考察不同模式在太陽地磁活動較強情況下的響應(yīng)，選取全年地磁指數(shù)Ap大于50 nT的時段(共計14 d)進(jìn)行軌道解算結(jié)果比較，其中超過100 nT的有6 d，2003年10月29至31日太陽輻射指數(shù)和地磁指數(shù)均達(dá)到了較強活動水平。

將所選14 d CHAMP衛(wèi)星RSO軌道作為觀測量進(jìn)行軌道解算，不同模式解算軌道位置的均方根誤差如圖8所示。從圖中可以看出:在地磁擾動最強烈的時段(圖8中橫軸第9～11個算例)，所有模式解算軌道的位置均方根誤差均呈現(xiàn)增長趨勢，但區(qū)別不大，而11月20日的地磁擾動沒有10月底劇烈，但4個模式的表現(xiàn)卻出現(xiàn)了比較顯著的差別，JB2008和DTM2013顯著優(yōu)于另外2個模式。

圖8 地磁擾動情況下不同模式解算位置均方根誤差Fig.8 RM S of estimated position error using differentmodels during geomagnetic disturbance

對圖8中位置均方根誤差再進(jìn)行均值和標(biāo)準(zhǔn)偏差統(tǒng)計，結(jié)果見表3。由于數(shù)據(jù)樣本有限，誤差的統(tǒng)計特性不足，但依然反映出JB2008和DTM2013模式的誤差相對較低。

表3 不同模式定軌位置均方根誤差的均值和標(biāo)準(zhǔn)偏差Table 3 M ean and standard deviation of position RMS error for differentmodels

5 結(jié)論

本文從應(yīng)用支持決策的角度，對當(dāng)前主流的JR、MSISE00、JB2008和DTM2013經(jīng)驗密度模式進(jìn)行了全面的精度評估。

1)利用國外權(quán)威發(fā)布和國內(nèi)自主探測的密度數(shù)據(jù)，計算不同模式相對觀測誤差。結(jié)果表明:DTM2013和JB2008的系統(tǒng)偏差和隨機誤差較早期發(fā)布的JR和MSISE00模式有明顯的改善，整體表現(xiàn)更為優(yōu)異。

2)從應(yīng)用角度出發(fā)，以天宮一號目標(biāo)飛行器和CHAMP衛(wèi)星精密軌道數(shù)據(jù)為基準(zhǔn)，計算4種模式軌道預(yù)報相對精密軌道的偏差。結(jié)果表明:MSISE00和DTM2013的效果略優(yōu)于JR和JB2008模式，主要原因之一是24 h的預(yù)報時間對模式誤差傳遞的敏感度有限；在地磁擾動情況下，JB2008和DTM2013的軌道解算精度更高。

3)新發(fā)布模式由于同化了新近獲取的觀測數(shù)據(jù)，在模式精度改善方面取得了明顯進(jìn)步。盡管在文中軌道預(yù)報的應(yīng)用效果不明顯，但通過實測數(shù)據(jù)的比對，仍然證明了DTM2013和JB2008等新模式的精度優(yōu)勢。

綜合后續(xù)的應(yīng)用需求，給出以下建議:

1)高精度觀測數(shù)據(jù)對于改善現(xiàn)有模式精度具有顯著作用。在當(dāng)前微納衛(wèi)星平臺日益發(fā)展的條件下，利用微納平臺開展分布式探測，不斷提升大氣密度觀測能力和水平，將有助于開展模式修正工作。

2)為了改善對太陽輻射和地磁活動響應(yīng)的不足，新近發(fā)布的模式采用了不同的指數(shù)驅(qū)動模式運行，包括很多衛(wèi)星觀測數(shù)據(jù)。而目前中國在這方面的能力還十分有限，提高對太陽的多譜段觀測能力和地磁活動的觀測能力對于改善現(xiàn)有密度模式和提高預(yù)報能力具有重要意義。

3)與新發(fā)布的模式比較，目前國內(nèi)工程中廣泛應(yīng)用的MSISE00模式已經(jīng)沒有明顯的優(yōu)勢。針對后續(xù)應(yīng)用對精度需求的不斷提高，應(yīng)著手開展對現(xiàn)有密度模式的改進(jìn)，并建立專門服務(wù)航天應(yīng)用的工程模型。