熱層大氣密度動態(tài)修正模型研究

劉 衛(wèi)，羅冰顯，龔建村，王榮蘭，容建剛

（1.中國科學(xué)院 國家空間科學(xué)中心 空間天氣學(xué)國家重點實驗室；2.中國科學(xué)院 國家空間科學(xué)中心 太陽活動與空間天氣重點實驗室：北京 100190；3.中國科學(xué)院微小衛(wèi)星創(chuàng)新研究院，上海 201210；4.中國科學(xué)院大學(xué)，北京 100049）

0 引言

熱層指地球大氣層中距地高度60～1000 km的區(qū)域，是地球低軌道航天器運行密集區(qū)和重要環(huán)境。大氣密度預(yù)報誤差是影響低軌航天器定軌和軌道預(yù)報精度的關(guān)鍵因素，其中地磁暴等空間環(huán)境擾動期問題尤為突出。2022 年2 月初，SpaceX 公司超40 顆“星鏈”衛(wèi)星相繼因電推進器推力無法抵御地磁暴期大氣密度陡增引起的大氣阻力變化，在約300 km 高度的軌道爬升過程中墜毀于大氣層中[1]。

1957 年第一顆人造地球衛(wèi)星進入太空后，通過軌道跟蹤發(fā)現(xiàn)了大氣阻力對航天器軌道運動的影響。隨后，以高層大氣為研究對象的新學(xué)科——高層大氣物理學(xué)誕生[2]。20 世紀(jì)60 年代，基于衛(wèi)星軌道衰變數(shù)據(jù)，Jacchia 等構(gòu)建了第1 代HP、Jacchia經(jīng)驗大氣模型[3]。目前已發(fā)展并形成以JB2008、DTM2020、MSIS2.0 和GOST04 等為代表的熱層大氣密度模型群[4-7]?，F(xiàn)有熱層大氣模型所考慮的因素只能涵蓋已知的大氣密度變化規(guī)律，而高層大氣復(fù)雜多變，即使對于平靜空間環(huán)境，模型標(biāo)準(zhǔn)差達15%～20%，擾動期間甚至可超過100%。不規(guī)則的熱層大氣密度增強也會對衛(wèi)星阻力和空間碎片監(jiān)測預(yù)警產(chǎn)生重大影響[8-10]；現(xiàn)有經(jīng)驗?zāi)Ｐ蜔o法滿足諸如空間態(tài)勢感知等特定場景的高精度需求。

動態(tài)大氣模型校正（DCA）是一種改進或修正現(xiàn)有大氣模型的技術(shù)[11]，一般用于熱層大氣密度高精度短臨預(yù)報，可以提供大氣密度變化的信息和相應(yīng)的統(tǒng)計數(shù)據(jù)。DCA 技術(shù)自20 世紀(jì)80 年代首先由Nazarenko 提出以來，在航天器軌道確定應(yīng)用方面受到持續(xù)關(guān)注。目前DCA 方法有基于加速度計和軌道跟蹤等數(shù)據(jù)，日常提供密度校正，用于200～800 km 熱層中性密度的預(yù)報[12-15]。針對熱層大氣密度高精度預(yù)測的場景，各航天大國均有專屬的解決方案。例如：美國空軍的高精度衛(wèi)星拖曳模型（HASDM）[13,16]，修正后的預(yù)報1 天相對誤差標(biāo)準(zhǔn)差平均下降25%[14]；俄羅斯空間觀測中心基于275～650 km 多顆空間目標(biāo)的軌道跟蹤數(shù)據(jù)，進行直接密度修正（DDCM），應(yīng)用效果良好[17-20]。

本文從大氣模型擴散平衡原理著手，對熱層大氣溫度垂直剖面上的拐點和外逸層2 個邊界溫度計算函數(shù)，按空間經(jīng)緯度球諧系數(shù)展開修正。首先，使用推導(dǎo)的熱層模型密度對修正系數(shù)偏導(dǎo)數(shù)公式計算測量矩陣，結(jié)合密度測量殘差建立動態(tài)修正的理論基礎(chǔ)；并利用泰勒公式對條件方程展開，只保留一階項的線性化處理，通過迭代求解修正系數(shù)。然后，基于MATLAB 軟件構(gòu)建熱層大氣密度動態(tài)修正模型系統(tǒng)，并在磁暴和磁靜期不同環(huán)境下進行應(yīng)用測試。討論動態(tài)修正系數(shù)的現(xiàn)報和預(yù)報應(yīng)用方式；通過分析邊界溫度修正量的序列，總結(jié)變化規(guī)律；結(jié)合國內(nèi)熱層大氣模型探測計劃，對熱層大氣密度修正系統(tǒng)的應(yīng)用前景進行總結(jié)和展望。

1 動態(tài)修正模型

1.1 修正原理

熱層模型的動態(tài)修正首先要選擇待修正的基礎(chǔ)模型。模型要求理論清晰，在表征熱層大氣上具有拓展性，無過多擬合公式。1964 年，Jacchia 發(fā)布了第1 個熱層大氣模型Jacchia-64，1970 年更新為Jacchia1970 (J70)[21]，其理論清晰、拓展性好。近年來，在美國空軍空間司令部推動下，基于半周年密度變化新算法和新的空間環(huán)境指數(shù)，J70 模型又拓展得到JB2006/2008 模型[4]。劉衛(wèi)等[22]利用J77、DTM78、MSIS00、JB2008 等不同熱層模型進行軌道預(yù)報，以精密軌道為參考，評估低地球軌道（LEO）衛(wèi)星的軌道預(yù)報誤差，證明Jacchia 系列模型性能優(yōu)越。綜上，本研究選擇J70 作為待修正的基礎(chǔ)模型。

J70 模型的外逸層溫度T∞為定義參數(shù)，不同地點的外逸層溫度由太陽活動指數(shù)F10.7及其81 d 均值Fˉ10.7，結(jié)合周日變化模型計算獲得。同時基于時間t和高度的修正，包括半周年、季節(jié)變化，以及緯度φ、地方時θ和地磁活動ap變化項。外逸層溫度作為溫度剖面的輸入，形成擴散方程積分的基礎(chǔ)。J70 模型計算熱層密度的第1 步是獲取頂點的溫度，

表征的是地磁指數(shù)為0 時全球夜側(cè)溫度的極小值。

外逸層溫度T∞是通過Tc乘以周日變化函數(shù)D(t,φ,θ)附加地磁活動函數(shù)ΔTG(ap)得到[3]，即

J70 模型溫度剖面通過邊界條件定義，下邊界90 km 處的溫度T=183 K，拐點（125 km 處）溫度Tx的梯度由dT/dh=0 給出。Tx是4 個經(jīng)驗參數(shù)（a,b,c,k）的函數(shù)，即

定義完備的溫度剖面后，通過積分?jǐn)U散方程獲取大氣主要成分N2、Ar、He、O2、O 的數(shù)密度，質(zhì)量密度ρi(z)可通過數(shù)密度獲得[21,23]，其中下標(biāo)i=1, 2,3, 4, 5 分別對應(yīng)成分N2、Ar、He、O2、O。

1998 年，Marcos 等首次通過估算Tc對熱層大氣密度模型進行修正，得到修正量ΔTc，使模型預(yù)報數(shù)據(jù)與衛(wèi)星跟蹤數(shù)據(jù)具有更好的一致性[24]。2005 年，Storz 等[16]對該方法進行了發(fā)展，估算ΔTc和ΔTx的球諧系數(shù)，對2 個全局溫度進行修正，溫度修正函數(shù)[16,25]為

因球諧修正系數(shù)計算過程所使用的測量數(shù)據(jù)構(gòu)成的約束條件數(shù)量有限，無法估算無窮高階的球諧系數(shù)，要對修正函數(shù)進行截斷處理。修正模型中，邊界溫度修正量ΔTx取1 階，ΔTc取2 階，共13 項修正系數(shù)，修正后的兩邊界溫度為：

構(gòu)建的條件方程是非線性的，球諧系數(shù)求解時，首先要對修正函數(shù)進行線性化處理，然后迭代估計。將修正模型密度ρM(X)在X0處作泰勒展開，

其中偏導(dǎo)數(shù)矩陣A根據(jù)鏈導(dǎo)法則分為密度對溫度和溫度對修正系數(shù)兩部分計算，其解析表達式的具體形式見文獻[26]。

1.2 系統(tǒng)框架

以J70 模型為基礎(chǔ)模型，熱層拐點和外逸層溫度補償量的球諧展開系數(shù)為自變量，利用實測密度作為驅(qū)動（或校正）數(shù)據(jù)，構(gòu)建條件方程。利用泰勒展開進行線性化，綜合利用正則變換等數(shù)值方法，對方程的“病態(tài)”進行處理。迭代求解動態(tài)修正系數(shù)，構(gòu)建熱層大氣密度動態(tài)修正系統(tǒng)，數(shù)據(jù)處理的流程如圖1 所示。

圖1 熱層大氣密度動態(tài)修正流程Fig.1 Flow chart of the dynamic correction system for determining thermospheric mass density

1.3 系統(tǒng)功能

構(gòu)建的熱層大氣密度動態(tài)修正系統(tǒng)主界面如圖2 所示，系統(tǒng)主要功能包括：

圖2 熱層大氣模型動態(tài)修正系統(tǒng)運行結(jié)果圖形顯示界面Fig.2 Graphic interface of the computational results of dynamic correction system of thermosphere model

1）校正數(shù)據(jù)覆蓋性分析?？稍趧討B(tài)系數(shù)求解前，評估數(shù)據(jù)的空間覆蓋質(zhì)量。通過大量算例分析，認(rèn)為需要有單軌道面單星的6 h 或雙軌道面3 h及以上空間覆蓋質(zhì)量方能滿足動態(tài)系數(shù)解算的基本需求。但空間覆蓋和動態(tài)系數(shù)解算質(zhì)量間的定量關(guān)系仍需深入研究。

2）按弧長連續(xù)解算動態(tài)修正系數(shù)。按設(shè)定弧長進行批量動態(tài)修正系數(shù)序列解算，可基于長期密度數(shù)據(jù)對熱層邊界溫度建模，進一步構(gòu)建熱層密度模型；也可用于短期（3 d）動態(tài)修正系數(shù)預(yù)測的輸入。

3）動態(tài)修正系數(shù)收斂性分析?？稍谛Ｕ龜?shù)據(jù)覆蓋性差，迭代過程可能發(fā)生動態(tài)修正系數(shù)不收斂的情況下，有效剔除不收斂弧段。

4）修正效果評估?？苫谑潞蟮拇髿饷芏葴y量數(shù)據(jù)，對模型的動態(tài)修正效果進行評估。

5）結(jié)果的圖形展示?？蓪⒂嬎慊蚍治鼋Y(jié)果以圖形方式直觀形象地展示給用戶，便于理解。

以2003 年10 月29 日 0:00—3:00，CHAMP 和GRACE A/B 衛(wèi)星3 h 弧段約700 個測量點的數(shù)據(jù)為校正數(shù)據(jù)，修正系統(tǒng)的運行結(jié)果如圖2 所示。圖2(a)為校正數(shù)據(jù)的空間分布。圖2(b)為解算的球諧系數(shù)。圖2(c)為校正數(shù)據(jù)時段兩邊界溫度的補償量：因校正數(shù)據(jù)來自分布在2 個軌道面的衛(wèi)星，圖中可見邊界溫度分布在2 條幅值相同、具有一定相位差的三角函數(shù)曲線上?？山Y(jié)合邊界溫度的物理范圍，分析修正量的合理性。圖2(d)為校正前/后的模型預(yù)測以及實測密度數(shù)據(jù)，其中紅色、黑色和藍色線條分別表示首次迭代、最終收斂和實測密度值。圖2(e)為校正前/后（圖中以紅色/黑色區(qū)分）的密度殘差分布，可以看到校正后的模型密度殘差顯著下降。

2 應(yīng)用測試

熱層大氣模型動態(tài)修正中，動態(tài)修正系數(shù)的獲取包括現(xiàn)報和預(yù)報2 種方式。定義前者為以相鄰弧段解算的動態(tài)修正系數(shù)直接進行后續(xù)弧段大氣密度模型修正；后者為以前序相鄰弧段解算的動態(tài)修正系數(shù)序列為輸入，利用動態(tài)系數(shù)預(yù)測模型計算相鄰后續(xù)弧段的動態(tài)修正系數(shù)，然后進行大氣密度模型修正。

為量化評估動態(tài)修正系統(tǒng)的性能，可以使用相對誤差ε表示模型預(yù)報值偏離實測值的百分比，ε越小則表示系統(tǒng)的性能越優(yōu)；還可以使用模型預(yù)報與實測密度的比值r，其值越接近1，則表示系統(tǒng)的性能越優(yōu)。這2 個系統(tǒng)性能評估量的具體表達式為：

2.1 動態(tài)修正系數(shù)現(xiàn)報

動態(tài)修正系數(shù)現(xiàn)報測試中，使用CHAMP 和GRACE A/B 這3 顆衛(wèi)星在2001—2010 期間的數(shù)據(jù)，按磁暴和磁靜分類分別選取6 個事件進行評估測試。磁暴事件的篩選標(biāo)準(zhǔn)為當(dāng)天3 h 的ap≥236；磁靜事件篩選標(biāo)準(zhǔn)為連續(xù)10 d 的ap≤48。以3 顆衛(wèi)星當(dāng)天24 h 實測為校正數(shù)據(jù)，進行動態(tài)修正系數(shù)解算。后續(xù)連續(xù)0～72 h 實測作為評估標(biāo)準(zhǔn)，分別對J70 模型、MSIS00 模型和J70 模型校正后（DCA）的預(yù)報密度誤差百分比進行統(tǒng)計，給出各模型預(yù)報0～72 h 相對誤差均值，統(tǒng)計結(jié)果見表1 和表2，表中出現(xiàn)修正前后相對誤差均值符號變化情況，這可能與評測數(shù)據(jù)地方時有關(guān)。

表1 磁暴期動態(tài)修正結(jié)果評測統(tǒng)計Table 1 Statistics of DCA results during magnetic storms

表2 磁靜期動態(tài)修正結(jié)果評測統(tǒng)計Table 2 Statistics of DCA results during magnetic quiet

由表中數(shù)據(jù)可以看到：6 次地磁暴事件中，DCA算法全部有效，修正后模型的密度預(yù)報性能均優(yōu)于MSIS00 模型，預(yù)報相對誤差平均降低33.8 個百分點；6 次地磁平靜事件中，除2003 年1 月3 日事件外，DCA 算法全部有效，較MSIS00 模型預(yù)報相對誤差平均降低9.1 個百分點?？梢姡? d 內(nèi)動態(tài)修正后的模型預(yù)報殘差較未修正模型顯著下降，這對軌道預(yù)報性能的提升非常關(guān)鍵。由修正失敗的個例可見，將計算所得修正系數(shù)直接用于后續(xù)的動態(tài)修正預(yù)報時，并不能確保大氣密度預(yù)報準(zhǔn)確性的提升，這可能與計算修正系數(shù)時的空間環(huán)境和后續(xù)預(yù)報時的空間環(huán)境之間是否存在突變等因素相關(guān)。

2.2 動態(tài)修正系數(shù)預(yù)報

動態(tài)修正系數(shù)預(yù)報測試中，以3 h 弧長分8 個弧段計算2003 年10 月29 日的動態(tài)修正系數(shù)，計算的拐點溫度和外逸層溫度修正量分別如圖3 和圖4 所示。

圖3 以3 h 弧長計算的拐點溫度修正量ΔTx（北極俯視圖）Fig.3 Corrected inflection point temperature ΔTx by using 3 h arc length (Arctic top view)

圖4 以3 h 弧長計算的外逸層溫度修正量ΔTc（北極俯視圖）Fig.4 Corrected exosphere temperature ΔTc by using 3 h arc length (Arctic top view)

由圖可見：拐點溫度修正量ΔTx具有呈扇形的兩極值區(qū)，相位差180°，隨時間推移兩極值中軸線順時針旋轉(zhuǎn)（北極俯視），這與地球自轉(zhuǎn)反向；在旋轉(zhuǎn)過程中兩極值表現(xiàn)也不同，角寬度有展寬、有收縮。對于外逸層溫度修正量ΔTc也存在兩扇形極值區(qū)，但兩極值區(qū)的相位差不是180°，且隨時間變化。即，兩邊界溫度都存在相位旋轉(zhuǎn)和幅值變化特性，結(jié)合空間環(huán)境指數(shù)可進行動態(tài)修正系數(shù)的短期變化建模。

3 小結(jié)及建議

討論動態(tài)修正的效果離不開具有一定時空覆蓋的熱層密度實時測量數(shù)據(jù)，分析動態(tài)修正系統(tǒng)的應(yīng)用前景必須結(jié)合熱層實測數(shù)據(jù)的來源。目前在軌空間目標(biāo)數(shù)量已超過20 000 顆，其中多數(shù)為運行在熱層大氣之中的低地球軌道空間目標(biāo)。選擇一定數(shù)量、合理分布的空間目標(biāo)，基于其軌道衰減來反演大氣密度，無疑是獲取實時大氣密度數(shù)據(jù)的絕佳途徑。美國空軍挑選氣動構(gòu)型穩(wěn)定、軌道面均勻分布的約80 顆空間目標(biāo)的軌道精密跟蹤數(shù)據(jù)，根據(jù)其軌道衰減反演全球覆蓋性的熱層大氣密度數(shù)據(jù)，并以此數(shù)據(jù)驅(qū)動其高精度衛(wèi)星拖曳模型（HASDM）。該模式已實現(xiàn)業(yè)務(wù)化運行，且能夠滿足軍事領(lǐng)域高精度需求。

但HASDM 的運行模式是建立在美空軍地基監(jiān)測網(wǎng)設(shè)備類型完備和全球布站的基礎(chǔ)上的，我們不能照搬，需結(jié)合具體情況來運營我們的動態(tài)修正系統(tǒng)。目前，國內(nèi)有2 項進行中的熱層大氣密度探測計劃，可結(jié)合熱層大氣模型動態(tài)修正系統(tǒng)進行業(yè)務(wù)化運行。具體建議如下：

1）氣動力標(biāo)校微小衛(wèi)星項目。該項目旨在研制并發(fā)射低成本、氣動穩(wěn)定、攜帶GNSS、定位精準(zhǔn)的大氣反演和彈道系數(shù)標(biāo)校微小衛(wèi)星（如圖5 所示），既可以基于自身精密星歷反演大氣密度，又可以標(biāo)校相近軌道空間目標(biāo)彈道信息，再與我國空間目標(biāo)監(jiān)視網(wǎng)軌道跟蹤數(shù)據(jù)聯(lián)合進行大氣密度的解算。這樣，只需少量標(biāo)校星進行最小規(guī)模的組網(wǎng)，并根據(jù)軌道衰減情況適時補充，即可滿足數(shù)據(jù)驅(qū)動需求，構(gòu)建高精度大氣密度動態(tài)修正系統(tǒng)。

圖5 氣動力標(biāo)校微小衛(wèi)星Fig.5 Aerodynamically calibrated microsatellite

2）熱層大氣探測衛(wèi)星星座?？梢栽撔亲l(wèi)星實時探測的熱層密度作為熱層大氣密度動態(tài)修正系統(tǒng)的數(shù)據(jù)驅(qū)動，動態(tài)實時進行熱層模型修正，以適應(yīng)高精度的熱層大氣密度需求場景。

4 結(jié)束語

本文以提高熱層大氣密度模型的分析精度為需求，對溫度剖面中逃逸層及拐點兩邊界溫度進行理論分析，通過CHAMP 和GRACE 衛(wèi)星觀測的大氣密度數(shù)據(jù)進行擬合，對現(xiàn)有大氣經(jīng)驗?zāi)Ｐ瓦M行了優(yōu)化，修正后的大氣模型比以往模型性能提升，預(yù)報相對誤差顯著下降。通過動態(tài)修正實例和修正模型的軌道預(yù)報應(yīng)用，有如下結(jié)論：

1）動態(tài)修正系數(shù)解算對校正數(shù)據(jù)的時空覆蓋具有一定需求；

2）動態(tài)修正后可顯著降低模型相對誤差均值，對提高軌道預(yù)報精度非常有效；

3）解算的動態(tài)修正系數(shù)直接用于模型密度預(yù)報時并不保證奏效，需要建立動態(tài)修正系數(shù)序列的短期預(yù)報模型；

4）兩邊界溫度修正量均具有兩扇形極值結(jié)構(gòu)，時變上具有相位旋轉(zhuǎn)和幅值變化特性。

以上結(jié)論可對熱層大氣動態(tài)修正模型構(gòu)建的工程實踐提供理論支撐和應(yīng)用參考。要說明的是，本文中的空間環(huán)境參數(shù)為實測值，意在算法開發(fā)階段避免因空間環(huán)境參數(shù)預(yù)報的不準(zhǔn)確影響算法參數(shù)估計的穩(wěn)定性問題。待算法開發(fā)完畢，后續(xù)實際業(yè)務(wù)應(yīng)用中，空間環(huán)境參數(shù)須采用預(yù)報值，并在動態(tài)系數(shù)預(yù)測、建模等方面開展進一步的分析工作。