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

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

1 引言

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

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

由于密度模式精度對軌道計算，尤其是軌道預報影響顯著，因此模式精度的評估往往是工程應用部門關(guān)注的焦點。Marcos分別對Jacchia70/71/77、MSIS77/79/83/86 和 Jacchia71、DTM2000、NRLMSISE00和JB2006進行了系統(tǒng)分析比較，Vallado從應用層面比較了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)驗模式，也沒有類似于馬歇爾工程熱層模型(Marshall Engineering Thermosphere Model，MET)的工程專用模式，因此，在實際工程應用前如何評估模式精度及計算效率，并據(jù)此選擇合適的經(jīng)驗模式具有重要的理論和現(xiàn)實意義。本文利用CHAMP衛(wèi)星的歷史數(shù)據(jù)和中國自主衛(wèi)星探測數(shù)據(jù)，結(jié)合不同類型航天器的軌道預報效果，綜合評估Jacchia-Robert(JR)、NRLMSISE00(MSISE00)、JB2008和DTM2013 4種不同模式的精度，分析它們在不同應用場景下的優(yōu)勢和特點，為后續(xù)的工程應用以及關(guān)聯(lián)的空間環(huán)境探測提供參考。

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

綜合考慮數(shù)據(jù)的精度和可獲得性，采用2類數(shù)據(jù)作為基準開展模式的精度評估:①使用廣泛且精度得到業(yè)內(nèi)認可的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)星是德國地學研究中心研制的重力測地衛(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)星是中國首顆以大氣密度探測為主要目標的微納衛(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é)果與真值進行比較。定義變量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é)果與實測密度進行比較。由于數(shù)據(jù)量較大，為提高計算效率且不失統(tǒng)計特性，對數(shù)據(jù)進行日平均降采樣處理，其相對誤差的均值和標準偏差如圖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年誤差擴大超過1個數(shù)量級。圖2給出了數(shù)據(jù)時間對應的不同波段太陽輻射流量，太陽輻射的整體水平呈逐年下降趨勢，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ù)的準確性和可靠性，李勰等提出了一種基于動力學的在軌標校方法，標校后數(shù)據(jù)的系統(tǒng)誤差優(yōu)于15%，因此作為評估模式精度的基準數(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個模式相對誤差的標準偏差基本相當，在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)水平相當?shù)?005年CHAMP衛(wèi)星數(shù)據(jù)的比較結(jié)果，可以看出不同模式基于2類數(shù)據(jù)的評估結(jié)論基本一致:JB2008的表現(xiàn)均為最優(yōu)，JR和MSISE00模式精度相當，較JB2008而言，有比較明顯的半年變化的系統(tǒng)偏差。此外，JB2008還呈現(xiàn)除了與其他2個模式不同的誤差特性，即JB2008計算結(jié)果較觀測值整體偏低，而其他2個模式結(jié)果偏高。

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

4 基于軌道的精度評估

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

4.1 評估原理與方法

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

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

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

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

4.2 天宮一號目標飛行器

天宮一號目標飛行器是中國首個自主研制的載人空間試驗平臺，于2011年9月29日發(fā)射入軌。本文選取2012年1月～4月的測量數(shù)據(jù)，利用1天的弧段進行軌道計算，分別采用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)誤差水平。由于天宮一號目標飛行器的構(gòu)型比較復雜，其阻尼系數(shù)的真實值較難確定，因此無法直接給出哪個模式最優(yōu)，只能進行相對比較。如果模式計算結(jié)果高于實際密度，為了彌補模式誤差對軌道的影響，解算阻尼系數(shù)將低于真實值，反之亦然。因此，從不同模式解算結(jié)果來看，DTM2013模式計算密度整體偏高，而JR模式計算密度整體偏低，2個模式平均偏差可達24.2%，而相對偏差最小的2個模式是DTM2013和JB2008，平均偏差為4.58%，這體現(xiàn)了新近開發(fā)模式與早期發(fā)布模式的差異。

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

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

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

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

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

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

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

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

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

將所選14 d CHAMP衛(wèi)星RSO軌道作為觀測量進行軌道解算，不同模式解算軌道位置的均方根誤差如圖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中位置均方根誤差再進行均值和標準偏差統(tǒng)計，結(jié)果見表3。由于數(shù)據(jù)樣本有限，誤差的統(tǒng)計特性不足，但依然反映出JB2008和DTM2013模式的誤差相對較低。

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

5 結(jié)論

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

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

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

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

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

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

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

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