熱層大氣密度反演與建模研究進(jìn)展

雷久侯，李若曦，任德馨，翁利斌，阮海炳

1 中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院 中國(guó)科學(xué)院近地空間環(huán)境重點(diǎn)實(shí)驗(yàn)室，合肥 230026

2 國(guó)防科技大學(xué)氣象海洋學(xué)院，長(zhǎng)沙 410015

3 南京信息工程大學(xué)遙感與測(cè)繪工程學(xué)院，南京 210044

0 引 言

熱層大氣受到太陽(yáng)輻射、太陽(yáng)風(fēng)-磁層-電離層相互作用、太陽(yáng)能量粒子注入以及低層大氣波動(dòng)等輸入能量的綜合作用，呈現(xiàn)復(fù)雜的時(shí)空變化特征（Hargreaves, 1992）.太陽(yáng)輻射作為熱層大氣十分重要的能量輸入，對(duì)熱層有著顯著的影響（Knipp et al., 2004）.日夜更替造成熱層接收到的太陽(yáng)輻射周日變化、太陽(yáng)自轉(zhuǎn)導(dǎo)致太陽(yáng)輻射的27 天周期變化、地球公轉(zhuǎn)引起的季節(jié)和年際變化，以及與太陽(yáng)活動(dòng)相關(guān)的11年周期變化等，這些因素都會(huì)顯著影響熱層大氣物理化學(xué)性質(zhì)（Jacchia, 1959a;Kasprzak and Newton, 1976; Paetzold and Zsch?rner,1961; Sharp et al., 1978）.如圖1 所示，以北京地區(qū)所在經(jīng)緯度大氣溫度和密度為例，低熱層溫度隨高度增加而迅速增加，直至150～200 km 處溫度逐漸趨于恒定，稱為逃逸層溫度.日側(cè)逃逸層溫度在太陽(yáng)活動(dòng)低年大約為750 K，在高年則可達(dá)1 400 K，溫度差異接近100%.熱層大氣密度隨高度增加呈指數(shù)衰減，而且受太陽(yáng)活動(dòng)的顯著調(diào)制.在600 km 高度區(qū)域，熱層大氣密度在太陽(yáng)活動(dòng)高年時(shí)比低年時(shí)約有兩個(gè)數(shù)量級(jí)的變化.耀斑、磁暴等太陽(yáng)爆發(fā)事件可對(duì)熱層大氣造成短時(shí)間尺度劇烈的擾動(dòng)（Bruinsma et al., 2006; Lei et al., 2011, 2012; Liu et al.,2010, 2011; Lühr et al., 2004; Sutton et al., 2005）.低層大氣擾動(dòng)也會(huì)通過(guò)波動(dòng)等方式影響熱層大氣（Yi?it et al., 2016）.

圖1 利用MSISE00 熱層大氣經(jīng)驗(yàn)?zāi)Ｐ瞳@得的北京地區(qū)所在經(jīng)緯度（116°E，40°N）大氣分別在太陽(yáng)活動(dòng)高年與低年的溫度（a）和密度（b）隨高度分布情況Fig.1 Altitude profiles of atmospheric temperature (a) and density (b) during solar maximum and minimum, based on the MSISE00 model, over Beijing (116°E, 40°N)

研究熱層有助于進(jìn)一步厘清日地相互作用實(shí)質(zhì)，另一方面，熱層是人類空間活動(dòng)的主要區(qū)域，在航天工程領(lǐng)域具有重要的應(yīng)用價(jià)值.大部分航天器的軌道，包括神舟、天宮以及各類資源遙感衛(wèi)星等在內(nèi)，均位于這個(gè)區(qū)域.熱層大氣密度作為熱層環(huán)境重要的物理參數(shù)，對(duì)衛(wèi)星定軌、航天器控制、太空垃圾監(jiān)測(cè)等工程具有十分重要的影響（李若曦，2022；任德馨，2021；翁利斌，2019）.雖然高層大氣十分稀薄，但是對(duì)于高速運(yùn)行的人造航天器來(lái)說(shuō)，也會(huì)受到稀薄大氣顯著的拖曳作用.在軌航天器所受大氣拖曳力與熱層大氣質(zhì)量密度直接相關(guān).隨著高度增加，大氣密度呈指數(shù)衰減，則大氣拖曳力迅速降低.對(duì)于300 km 高度的航天器，大氣拖曳力可以達(dá)到10-5m/s2，而在大約800 km 高度，大氣拖曳力迅速衰減至約10-10m/s2（Montenbruck and Gill, 2000）.熱層大氣的拖曳作用可造成航天器機(jī)械能減少，軌道高度下降，在無(wú)額外動(dòng)力提升軌道高度情況下使航天器最終隕落.此外，熱層大氣受空間天氣事件影響變化劇烈，航天器軌道也會(huì)隨之發(fā)生顯著變化，從而影響航天器控制和預(yù)報(bào)，增加空間目標(biāo)碰撞預(yù)警的不確定性（Li and Lei,2021a）.因此熱層研究對(duì)于加深人類對(duì)日地空間系統(tǒng)的認(rèn)識(shí)有著重要的科學(xué)意義，在航天工程領(lǐng)域也具有顯著的應(yīng)用價(jià)值.同時(shí)，熱層作為電離層的背景大氣，由于離子-中性相互作用等過(guò)程，熱層大氣的變化能夠引起電離層發(fā)生相應(yīng)的改變，因此，研究熱層大氣變化也是深入認(rèn)識(shí)和理解電離層變化特性的基礎(chǔ).

獲取有效觀測(cè)數(shù)據(jù)是熱層研究的基礎(chǔ).地表附近的大氣參數(shù)可以通過(guò)溫濕度計(jì)、氣壓計(jì)以及風(fēng)速計(jì)等設(shè)備十分便捷地測(cè)量.中間層和低熱層大氣參數(shù)則可以通過(guò)衛(wèi)星遙感、雷達(dá)、探空火箭以及氣象氣球等設(shè)備進(jìn)行探測(cè)（張曉芳和嚴(yán)衛(wèi)，2007）.然而，對(duì)于200 km 以上的大氣，以上觀測(cè)手段都難以有效測(cè)量.因此，對(duì)于熱層大氣的觀測(cè)數(shù)據(jù)十分稀少.目前，熱層大氣密度探測(cè)的主要方式是基于航天器軌道或者星載加速度計(jì)數(shù)據(jù)進(jìn)行密度反演.基于衛(wèi)星軌道數(shù)據(jù)的熱層大氣研究開(kāi)始于上世紀(jì)五六十年代（Priester et al., 1967; Prussing and Conway, 1993; Thomson, 1986; Vallado, 2001）.隨著人造航天器的增多，北美防空聯(lián)合司令部（NORAD）開(kāi)始以兩行數(shù)據(jù)（TLE）的形式發(fā)布空間目標(biāo)的軌道跟蹤數(shù)據(jù).很多研究人員基于TLE 數(shù)據(jù)反演了長(zhǎng)時(shí)期的熱層大氣密度，揭示了熱層大氣的氣候?qū)W演化規(guī)律（Emmert, 2009; Marcos et al., 2005; Picone et al., 2005）.隨著星載GPS 接收設(shè)備和加速度計(jì)的發(fā)展和應(yīng)用，基于衛(wèi)星精密定軌（POD）數(shù)據(jù)和加速度計(jì)數(shù)據(jù)，可以獲取更高時(shí)空分辨率和精度的大氣密度觀測(cè)數(shù)據(jù).近年來(lái)，CHAMP、GRACE、GOCE 以及Swarm 等衛(wèi)星的高質(zhì)量大氣密度觀測(cè)結(jié)果，已被廣泛應(yīng)用于熱層大氣研究中，用來(lái)揭示熱層大氣的精細(xì)變化特征，并優(yōu)化大氣模式（Bruinsma and Forbes, 2007, 2009, 2010; Li and Lei,2021a, 2021b; Lei et al., 2010; Liu and Lühr, 2005;Liu et al., 2007; Lühr et al., 2004; Ren and Lei, 2020,2022; Ruan et al., 2018; Weng et al., 2017, 2019）.

衛(wèi)星只能對(duì)全球熱層大氣進(jìn)行采樣式觀測(cè)，獲取時(shí)空連續(xù)的全球熱層大氣數(shù)據(jù)一般依賴于理論或者經(jīng)驗(yàn)?zāi)Ｊ接?jì)算.理論模式一般基于描述質(zhì)量、動(dòng)量以及能量的方程，通過(guò)數(shù)值計(jì)算的方法求解全球尺度的大氣參數(shù).理論模型可以模擬大氣演化過(guò)程及其物理機(jī)制，以及獲取較為精細(xì)的熱層大氣結(jié)構(gòu)，但其輸出的大氣參數(shù)的絕對(duì)值一般與真實(shí)結(jié)果相差較大，因此理論模式一般用于熱層大氣機(jī)理性研究.經(jīng)驗(yàn)?zāi)Ｊ絼t基于參數(shù)化的經(jīng)驗(yàn)方程，采用觀測(cè)數(shù)據(jù)驅(qū)動(dòng)的方式進(jìn)行參數(shù)擬合，雖然物理原理涉及較少，但是輸出結(jié)果與真實(shí)較為接近.經(jīng)驗(yàn)?zāi)Ｊ揭泊嬖谳^大誤差，在地磁平靜時(shí)期誤差一般約為15%，在磁暴期間則更大（Marcos, 2006）.大氣經(jīng)驗(yàn)?zāi)Ｊ降恼`差對(duì)于衛(wèi)星軌道預(yù)報(bào)、軌道維持、碎片預(yù)警以及隕落預(yù)報(bào)等航天工程具有重要影響（李若曦，2022；翁利斌，2019）.

本文將綜述當(dāng)前基于空間目標(biāo)兩行軌道數(shù)據(jù)、航天器精密軌道數(shù)據(jù)以及星載加速度計(jì)數(shù)據(jù)反演熱層大氣密度的方法，總結(jié)目前常用的大氣經(jīng)驗(yàn)?zāi)Ｊ剿惴ê瓦m用范圍，并針對(duì)當(dāng)前最新的模式優(yōu)化方案進(jìn)行重點(diǎn)闡述.需要指出的是，最近幾年學(xué)界在熱層大氣密度反演與建模研究方面取得諸多進(jìn)展，但這篇文章主要綜合我們近期的三篇博士論文（李若曦，2022；任德馨，2021；翁利斌，2019）的進(jìn)展和部分成果，因此本文內(nèi)容難免存在疏漏之處.

1 熱層大氣密度反演方法

1.1 熱層大氣密度反演基本原理

航天器飛行高度一般高于200 km.200 km 以上熱層區(qū)域大氣十分稀薄，大氣分子平均自由程超過(guò)100 m，遠(yuǎn)大于航天器尺度（劉林，1992）.因此，熱層大氣以自由分子流的形式與航天器相互作用，產(chǎn)生的拖曳力為：

其中A為 航天器跡向等效橫截面積，Msat表示衛(wèi)星質(zhì)量，CD代表大氣阻尼系數(shù)，vr表示航天器相對(duì)大氣運(yùn)動(dòng)速度， ρ為航天器附近的熱層大氣質(zhì)量密度.航天器半長(zhǎng)軸在切向攝動(dòng)力的作用下衰減.對(duì)于軌道高度較低的航天器，切向攝動(dòng)力主要為大氣拖曳阻力.根據(jù)軌道攝動(dòng)理論，航天器半長(zhǎng)軸a在大氣拖曳阻力[式（1）]作用下的變化可以表示為（Vallado, 2001）：

其中e、f和n分別指代航天器軌道偏心率、真近點(diǎn)角和平均運(yùn)動(dòng)等軌道根數(shù).根據(jù)式（1）或式（2），可反演熱層大氣密度.其基本原理為，通過(guò)星載加速度計(jì)直接測(cè)量或者通過(guò)航天器軌道衰減間接計(jì)算大氣拖曳力，然后根據(jù)衛(wèi)星形態(tài)、姿態(tài)以及面板材料等參數(shù)估計(jì)大氣阻尼系數(shù)CD，最后計(jì)算大氣密度.

1.2 計(jì)算有效橫截面積和阻尼系數(shù)

根據(jù)式（1）反演熱層大氣密度，需要準(zhǔn)確計(jì)算衛(wèi)星有效橫截面積A與 大氣阻尼系數(shù)CD.大氣分子與每一個(gè)衛(wèi)星面板發(fā)生作用，產(chǎn)生拖曳效應(yīng).根據(jù)Sentman 模型（Sentman, 1961），可以首先計(jì)算航天器每一個(gè)面板所受大氣拖曳力，再將所有面板所受總的拖曳力疊加.每塊衛(wèi)星面板大氣阻尼系數(shù)與橫截面積的乘積為：

vout表示大氣分子流與衛(wèi)星面板碰撞后的速度，vat為大氣相對(duì)衛(wèi)星速度.

S表示vat與大氣分子熱力學(xué)最可幾分布速度之比：

最后疊加航天器所有面板所得結(jié)果，得到航天器整體有效橫截面積與大氣阻尼系數(shù)之乘積.如圖2 所示，以CHAMP 衛(wèi)星為例，衛(wèi)星由13 個(gè)面板組成，其有效橫截面積與大氣阻尼系數(shù)的乘積呈現(xiàn)與衛(wèi)星環(huán)繞周期相關(guān)的周期性變化特征，最大值約為2.35，最小值約1.95.該周期性變化的主要原因與熱層大氣晝夜溫度和密度差異有關(guān)（李若曦，2017）.

圖2 2009年第32 天CHAMP 衛(wèi)星大氣阻尼系數(shù)與有效面積乘積的變化Fig.2 The products of CHAMP drag coefficients and crosssection areas during day 32, 2009

1.3 基于兩行數(shù)據(jù)反演熱層大氣密度

TLE 是美國(guó)北美防空聯(lián)合司令部發(fā)布的空間目標(biāo)軌道跟蹤測(cè)量數(shù)據(jù).該數(shù)據(jù)庫(kù)囊括了自1957年至今超過(guò)50 000 個(gè)空間目標(biāo)的軌道跟蹤數(shù)據(jù).TLE 數(shù)據(jù)包含空間目標(biāo)編號(hào)、軌道根數(shù)以及歷元時(shí)間等信息.基于TLE 數(shù)據(jù)反演熱層大氣密度的理論最初由King-Hele（1987）給出，Picone 等（2005）進(jìn)一步做了優(yōu)化.根據(jù)式（2），理論上通過(guò)TLE 給出的歷元之間半長(zhǎng)軸差即可以反演密度.但是，TLE給出的是航天器瞬時(shí)半長(zhǎng)軸（或密切半長(zhǎng)軸），且地球非球形引力場(chǎng)以及潮汐攝動(dòng)對(duì)航天器半長(zhǎng)軸的影響遠(yuǎn)大于大氣拖曳力.不同TLE 歷元之間半長(zhǎng)軸的變化，除了大氣拖曳力的作用，還包括地球引力場(chǎng)導(dǎo)致的周期性變化.因此，Picone 等（2005）采用平均運(yùn)動(dòng)（mean motion）的變化反演大氣密度.將式（2）中半長(zhǎng)軸a用平均運(yùn)動(dòng)nM替換得到：

式中B=ACD/Msat，定義為逆彈道系數(shù).F為風(fēng)影響參數(shù)，可表示為航天器位置和速度的函數(shù).相鄰兩個(gè)歷元ti至tk之間的平均大氣密度(ik)為：

式中nM可以通過(guò)TLE 數(shù)據(jù)直接讀取.計(jì)算需要計(jì)算積分步長(zhǎng)上的v(t)和F(t).可以使用SGP4（Simplified General Perturbations 4）軌道積分器（Hoots and Roehrich, 1980）進(jìn)行軌道外推，從初始?xì)v元開(kāi)始計(jì)算每個(gè)積分步長(zhǎng)上的航天器位置和速度，進(jìn)而得到v(t)和F(t).基于TLE 數(shù)據(jù)進(jìn)行大氣密度反演時(shí)，一般將大氣逆彈道系數(shù)近似地設(shè)置為常數(shù).圖3 展示了基于Starshine 衛(wèi)星TLE數(shù)據(jù)的反演密度和MSISE00 結(jié)果之比對(duì).Starshine衛(wèi)星形狀為規(guī)則球形，逆彈道系數(shù)基本保持恒定（Lean et al., 2006）.可以看出TLE 反演密度與MSISE00 結(jié)果一致，較好地描述了衛(wèi)星整個(gè)任務(wù)期間軌道密度的變化特征.由于TLE 數(shù)據(jù)時(shí)間分辨率低，因此一般只能反演沿軌道2～5 天的平均大氣密度.TLE 數(shù)據(jù)反演密度主要用于研究熱層大氣的長(zhǎng)期變化規(guī)律，工程上也被用于大氣模式標(biāo)定工作.

圖3 Starshine-1/2/3 衛(wèi)星軌道高度（a）和基于TLE 數(shù)據(jù)的反演密度（b）.Starshine-1/2/3 分別發(fā)射于1999年5月27日、2001年12月5日和2001年9月29日Fig.3 The altitudes of Starshine-1/2/3 satellites (a) and the TLE-based orbital densities (b).Starshine-1/2/3 satellites were launched on 27 May 1999, 5 December 2001 and 29 September, 2001, respectively

1.4 基于精密軌道數(shù)據(jù)反演熱層大氣密度

隨著星載導(dǎo)航設(shè)備的應(yīng)用，越來(lái)越多的航天器可以被精密定軌.二十一世紀(jì)以來(lái)，諸如CHAMP、GRACE 以及Swarm 等衛(wèi)星均具有較高定軌精度.基于這些高精度的航天器軌道數(shù)據(jù)和精確的力學(xué)模型，可以反演獲得具有較高時(shí)空分辨率的熱層大氣密度觀測(cè)數(shù)據(jù).本節(jié)將介紹幾種基于精密軌道數(shù)據(jù)反演熱層大氣密度的方法.

1.4.1 插值微分法

根據(jù)式（1），可以通過(guò)計(jì)算航天器所受大氣拖曳阻力adrag，反演大氣密度.插值微分法的基本原理是對(duì)航天器位置矢量進(jìn)行二階微分，計(jì)算航天器瞬時(shí)加速度，然后減去大氣拖曳力以外的所有作用力，最終得到大氣拖曳阻力（Calabia et al., 2015;Calabia and Jin, 2017）.航天器瞬時(shí)加速度計(jì)算如下：

其中，r和r˙分別代表航天器位置和速度.t-2到t2代表相鄰的5 個(gè)時(shí)刻.ΔT代表t-1到t1的時(shí)間間隔.為了精確計(jì)算航天器所受加速度，時(shí)間間隔 ΔT需足夠低以降低誤差.經(jīng)過(guò)計(jì)算比較，該時(shí)間間隔被設(shè)置為0.05 s.由于航天器精密星歷時(shí)間間隔一般為幾秒到幾十秒，因此需要在此基礎(chǔ)上進(jìn)行插值，以獲取更短時(shí)間間隔的星歷.經(jīng)過(guò)插值微分所得航天器瞬時(shí)加速度，包括地球重力、潮汐、日月引力、光壓、相對(duì)論效應(yīng)以及大氣拖曳力等所有作用力.除大氣拖曳力，其他作用力使用精確力學(xué)模型計(jì)算，然后從總的加速度中扣除，得到大氣拖曳力，最后進(jìn)行密度反演：

其中aG代表包括地球中心引力、非球形引力以及潮汐在內(nèi)的地球引力，aS和aM分別指代太陽(yáng)和月球中心引力場(chǎng)，asp代表光壓以及地球反照壓，arel代表相對(duì)論效應(yīng).

微分插值法可以獲取較高時(shí)空分辨率的熱層大氣密度，如圖4 所示為基于精密軌道數(shù)據(jù)使用插值微分法獲取的GRACE 軌道密度與加速度計(jì)結(jié)果之比對(duì)，可以看出兩者基本一致.加速度計(jì)數(shù)據(jù)反演密度具有最高的精度，一般作為參考密度.需要說(shuō)明的是，插值的準(zhǔn)確性依賴于較高的星歷數(shù)據(jù)密度.

圖4 基于GRACE 衛(wèi)星精密軌道數(shù)據(jù)以及加速度計(jì)數(shù)據(jù)反演的2011—2016年熱層大氣密度（修改自Calabia and Jin, 2017）Fig.4 POD-based thermospheric mass densities from 2011 to 2016.The densities derived from accelerometer data are plotted as references (modified from Calabia and Jin,2017)

1.4.2 半長(zhǎng)軸衰減法

根據(jù)式（2），使用精密軌道數(shù)據(jù)將方程中e、f和n等開(kāi)普勒軌道根數(shù)由衛(wèi)星位置和速度替換，只需計(jì)算航天器在大氣拖曳力作用下的半長(zhǎng)軸衰減率，則大氣密度可以反演得到.使用半長(zhǎng)軸衰減法反演密度最主要的工作就是計(jì)算半長(zhǎng)軸在大氣拖曳阻力作用下的衰減率.下面將詳細(xì)分析的算法.

對(duì)于兩個(gè)連續(xù)的歷元（t0,r0,v0）和（t1,r1,v1），對(duì)應(yīng)的半長(zhǎng)軸分別為a0和a1.考慮除大氣拖曳阻力以外的所有作用力，從歷元（t0,r0,v0）外推至t1時(shí)刻，得到模擬歷元（t1,r1′,v′1），以及對(duì)應(yīng)的模擬半長(zhǎng)軸a′1.因此在Δt=t1-t0時(shí)間段內(nèi)航天器半長(zhǎng)軸在大氣拖曳力作用下的衰減量為a′1-a1（Sang et al., 2012）.需要說(shuō)明的是，這里采用軌道外推的方法計(jì)算半長(zhǎng)軸衰減率的原因是，不同時(shí)刻的半長(zhǎng)軸變化量是大氣拖曳力與地球非球形攝動(dòng)力共同作用的結(jié)果，并不能反映大氣拖曳力的單獨(dú)作用.因此，航天器半長(zhǎng)軸受大氣拖曳力作用的衰減率為：

將上式帶入式（2）即可計(jì)算熱層大氣密度.圖5 展示了基于精密軌道數(shù)據(jù)使用半長(zhǎng)軸衰減法的反演密度（藍(lán)線），可以看出該結(jié)果與加速度計(jì)數(shù)據(jù)結(jié)果（紅線）具有較好的一致性.

圖5 基于CHAMP 精密軌道數(shù)據(jù)（藍(lán)線）與加速度計(jì)數(shù)據(jù)（紅線）之比對(duì)（修改自Sang et al., 2012）Fig.5 The comparison of CHAMP orbital density derived from POD data and that from accelerometer data (modified from Sang et al., 2012)

1.4.3 能量衰減法

半長(zhǎng)軸衰減法和插值微分法本質(zhì)上都是利用力學(xué)方法，Li 和 Lei（2021b）指出可以基于能量守恒原理反演大氣密度.衛(wèi)星在軌運(yùn)行，受到中心引力和其他攝動(dòng)力的作用.該攝動(dòng)力分為保守?cái)z動(dòng)力和非保守?cái)z動(dòng)力.航天器在保守力作用下機(jī)械能守恒，在非保守力作用下機(jī)械能發(fā)生改變.對(duì)于低軌航天器，非保守力主要為大氣拖曳力、光壓以及地球反照壓.因此，可以通過(guò)計(jì)算衛(wèi)星能量的衰減反演熱層大氣密度.這種方法不需要進(jìn)行軌道外推，能夠大幅降低計(jì)算成本.而且可以根據(jù)星歷精度選擇適當(dāng)時(shí)間步長(zhǎng)，對(duì)于數(shù)據(jù)的要求更低，適用范圍更廣.在軌航天器總機(jī)械能可以表示為：

等式右邊分別代表航天器動(dòng)能、地球引力勢(shì)能、日月引力勢(shì)能以及時(shí)變引力勢(shì)能，時(shí)變引力勢(shì)能主要由地球自轉(zhuǎn)和各種地球潮汐造成.包含時(shí)變引力勢(shì)和地球重力場(chǎng)在內(nèi)的航天器總的引力勢(shì)能變化可以用地球引力對(duì)航天器總的做功來(lái)表示.航天器在所選取弧段內(nèi)的能量變化可以寫(xiě)為：

式中g(shù)eci表示慣性坐標(biāo)系下的地球引力.下標(biāo)1 和2分別代表所選取弧段的初末位置.這里機(jī)械能的衰減等于總的非保守力做功ΔE=vΔt·an，an表示非保守力加速度.圖6 展示了基于精密軌道數(shù)據(jù)使用能量衰減法反演的GRACE 衛(wèi)星非保守力加速度與加速度計(jì)結(jié)果之比對(duì)，兩者結(jié)果較一致.為了得到大氣拖曳阻力還需要從非保守力an中減去光壓asrp和地球反照壓atir等其他非保守力.將所得大氣拖曳力帶入式（1），即可反演所選取弧度上的大氣平均密度.

圖6 2006年1月1日基于GRACE 衛(wèi)星精密軌道數(shù)據(jù)反演得到的非保守力加速度（藍(lán)線）與加速度計(jì)測(cè)量結(jié)果（紅線）之比對(duì)（修改自Li and Lei, 2021b）Fig.6 The comparison of GRACE along-track non-conservative forces derived from POD data on 1 January, 2006,and those from the accelerometer (modified from Li and Lei, 2021b)

1.5 基于加速度計(jì)數(shù)據(jù)反演熱層大氣密度

自上世紀(jì)六七十年代開(kāi)始，加速度計(jì)開(kāi)始作為空間載荷進(jìn)入太空，用于測(cè)量航天器所受到的非保守力（Boudon et al., 1979; Champion and Marcos,1973）.二十一世紀(jì)以來(lái)，一系列低軌衛(wèi)星如CHAMP、GRACE 等攜帶更為精密的加速度計(jì)進(jìn)入太空，當(dāng)前投入使用的星載加速度計(jì)主要由法國(guó)航空航天研究所（ONERA）制造.表1 羅列了部分ONERA 制造的星載加速度計(jì)主要技術(shù)指標(biāo)（Christophe et al., 2015; 馬云, 2019; Touboul et al.,2012）.基于這些高精度加速度計(jì)數(shù)據(jù)，研究人員得以獲取高時(shí)空分辨率、高精度的大氣密度觀測(cè)數(shù)據(jù)，用于研究熱層大氣小尺度變化特征，進(jìn)一步提升大氣模式準(zhǔn)確性.

表1 ONERA 靜電加速度計(jì)主要技術(shù)指標(biāo)Table 1 The technical indicators of some on-board electrostatic accelerometers made by ONERA

加速度計(jì)數(shù)據(jù)一般需要在軌標(biāo)定，以下為基于精密軌道數(shù)據(jù)的線性標(biāo)定方法（Bruinsma et al.,2004）：

式中，aPOD表示使用精密軌道數(shù)據(jù)反演的等效非保守力加速度，aAcc代表加速度計(jì)測(cè)量結(jié)果，Scale和Bias 代表尺度和偏差標(biāo)定因子.根據(jù)不同的衛(wèi)星，選取不同的估計(jì)弧長(zhǎng).以GRACE 為例，選擇一天作為一個(gè)估計(jì)弧段，通過(guò)最小二乘法估計(jì)標(biāo)定因子，然后使用該標(biāo)定因子標(biāo)定加速度計(jì)數(shù)據(jù).加速度計(jì)測(cè)量的是所有非保守力加速度，為了得到大氣拖曳力還需要減去光壓、地球反照壓等其他非保守力.根據(jù)式（1），熱層大氣密度可如下得到：

加速度計(jì)數(shù)據(jù)反演密度具有高精度，高時(shí)空分辨率等優(yōu)勢(shì).圖7 展示了基于GRACE 衛(wèi)星精密軌道數(shù)據(jù)（圖7a、7b）以及加速度計(jì)數(shù)據(jù)（圖7e、7f）反演的2017年9月磁暴期間熱層大氣密度變化.可以看出精密軌道數(shù)據(jù)反演的熱層大氣密度與加速度計(jì)反演密度均能描述磁暴期間熱層大氣的基本變化特征.磁暴期間，極區(qū)熱層大氣密度迅速上升，并激發(fā)了大尺度的熱層大氣波動(dòng)現(xiàn)象，將擾動(dòng)傳播至低緯區(qū)域.加速度計(jì)數(shù)據(jù)反演密度時(shí)空分辨率最高，能夠給出更為精細(xì)的熱層大氣變化特征.如圖7f 所示，加速度計(jì)數(shù)據(jù)反演密度捕捉到了在9月7日平靜時(shí)期熱層大氣存在的大尺度波動(dòng)現(xiàn)象，這可能與當(dāng)天較高的AE 指數(shù)相關(guān).

圖7 基于GRACE 衛(wèi)星精密軌道數(shù)據(jù)反演得到的2017年9月磁暴期間的熱層大氣密度（a, b）與加速度計(jì)反演密度（e,f）.密度單位為10-12 kg/m3.磁暴期間行星際磁場(chǎng)By、Bz 以及AE 指數(shù)在子圖（c, d）和（g, h）中給出（修改自Li and Lei, 2021b）Fig.7 Thermospheric mass densities retrieved from GRACE POD data and those from the accelerometer data during the 2017 September storm (in units of 10-12 kg/m3).The interplanetary magnetic field By, Bz components, and the AE index are shown in the bottom for reference (modified from Li and Lei, 2021b)

2 常用熱層大氣經(jīng)驗(yàn)?zāi)Ｊ?/h2>

在航天器定軌、預(yù)報(bào)以及空間碎片預(yù)警等航天工程中，一般使用經(jīng)驗(yàn)?zāi)Ｊ接?jì)算熱層大氣密度.本節(jié)我們將介紹當(dāng)前主要的大氣經(jīng)驗(yàn)?zāi)Ｊ剑M(jìn)一步介紹最新模式優(yōu)化方法及其結(jié)果.根據(jù)靜力學(xué)平衡假設(shè)，熱層大氣密度 ρ隨高度h的變化為（Montenbruck and Gill, 2000）：

其中 ρ0為參考高度的大氣密度，H表示大氣標(biāo)高.R為氣體常數(shù)，T(h)代表大氣溫度，ma為大氣相對(duì)分子質(zhì)量，g(h)為地球重力加速度.根據(jù)式（16）和（17），熱層不同種類分子隨高度的衰減速率不同，因此不同成分隨高度的分布也不同.熱層大氣經(jīng)驗(yàn)?zāi)Ｊ浇５幕驹硎腔诖罅康臏囟?、密度觀測(cè)數(shù)據(jù)，擬合半經(jīng)驗(yàn)半物理的數(shù)學(xué)公式，得到大氣溫度廓線，再計(jì)算不同分子成分的密度分布.

使用經(jīng)驗(yàn)?zāi)Ｊ綍r(shí)，根據(jù)輸入的太陽(yáng)、地磁活動(dòng)指數(shù)以及時(shí)空信息可以獲得熱層大氣密度結(jié)果.由于其易于使用且能夠較好再現(xiàn)大氣密度變化情況，已被廣泛應(yīng)用于航天工程和科學(xué)研究中.目前常用的經(jīng)驗(yàn)?zāi)Ｊ街饕蠮acchia、DTM、MSIS 等系列模式，經(jīng)過(guò)不斷發(fā)展改進(jìn)更新，這些模式的模擬能力和精度也在不斷提高.

2.1 Jacchia 系列模式

Jacchia 系列模型從上世紀(jì)六十年代開(kāi)始開(kāi)發(fā)，其發(fā)展歷程和主要情況如表2 所示.下面以J71 模型為例簡(jiǎn)要介紹Jacchia 模式算法（Jacchia, 1971;Montenbruck and Gill, 2000）.首先根據(jù)太陽(yáng)活動(dòng)以及地磁活動(dòng)計(jì)算全球逃逸層溫度T∞分布：

表2 Jacchia 系列模式基本情況Table 2 Brief introductions of Jacchia models

其中，T1代表太陽(yáng)活動(dòng)周期、半年、季節(jié)以及周日變化對(duì)逃逸層溫度的影響.ΔT∞代表地磁活動(dòng)對(duì)逃逸層溫度的影響.溫度廓線下邊界設(shè)置為90 km 高度處溫度T0=183 K.在獲取整個(gè)溫度剖面的基礎(chǔ)上，可根據(jù)溫度與密度的函數(shù)關(guān)系計(jì)算密度分布.Gill（1996）建立了以二次多項(xiàng)式表示密度隨高度變化的函數(shù)，通過(guò)插值計(jì)算大氣密度，可以更為高效的計(jì)算全球密度分布.

2.2 MSIS 系列模式

MSIS（Mass Spectrometer Incoherent Scatter）模式基于非相干散射雷達(dá)、探空火箭以及質(zhì)譜儀等觀測(cè)數(shù)據(jù)建立，發(fā)展歷程和主要情況如表3 所示.下面以MSIS86 模式為例簡(jiǎn)要介紹該模式計(jì)算方案（Hedin, 1987）.與Jacchia 模式相似，MSIS 也需要首先計(jì)算熱層大氣溫度分布，然后利用各成分對(duì)溫度的依賴性得到不同高度大氣成分分布，最后將所有成分疊加獲得大氣密度.熱層大氣溫度隨高度變化函數(shù)為：

表3 MSIS 系列模式基本情況Table 3 Brief introductions of MSIS models

式中T∞=T∞[1+G(L)] 表示逃逸層溫度，Tl=Tl[1+G(L)]為120 km 處的大氣溫度.其中，和分別表示逃逸層平均溫度和120 km 處平均溫度，在模式中設(shè)定為常數(shù)，G(L)為經(jīng)驗(yàn)函數(shù)，表征時(shí)間、緯度、太陽(yáng)和磁場(chǎng)活動(dòng)等參數(shù)對(duì)大氣溫度的影響 ，ε (z,zl)=(z-zl)(Rp-zl)/(Rp+zl)，Rp=6 356.77 km.σ=/(T∞-Tl)，其中=[1+G(L)]，為溫度隨高度變化梯度.

熱層大氣的密度分布可以看做是大氣分子擴(kuò)散與混合共同作用的結(jié)果.假設(shè)120 km 高度處不同氣體成分的數(shù)密度為常數(shù)，則成分隨高度的變化可以表示為：

其中nd表示擴(kuò)散函數(shù)，nm為混合函數(shù).參數(shù)C1(z)···Cn(z)表示其他修正參數(shù).最后根據(jù)大氣各成分的相對(duì)分子質(zhì)量Mi，將不同成分質(zhì)量密度疊加得到大氣密度：

其中NA=6.022×1023mol-1，為阿伏加德羅常數(shù).

MSIS90 版本將該模式的計(jì)算高度延展到地面（Hedin, 1991），后來(lái)美國(guó)海軍研究實(shí)驗(yàn)室繼續(xù)開(kāi)發(fā)了NRLMSISE00 版本（Picone et al., 2002），使用的數(shù)據(jù)包括質(zhì)譜儀和地基非相干散射雷達(dá)觀測(cè)數(shù)據(jù)，以及星載加速度計(jì)數(shù)據(jù)、衛(wèi)星軌道數(shù)據(jù)反演密度和衛(wèi)星軌道衰減數(shù)據(jù)庫(kù)等，同時(shí)還引入了異常氧成分信息.

2.3 DTM 系列模式

DTM（Drag Temperature Model）模式所用建模數(shù)據(jù)主要為質(zhì)譜儀、加速度計(jì)以及衛(wèi)星軌道資料反演的大氣密度，發(fā)展歷程和主要情況如表4 所示.本文以DTM-94 模式為例簡(jiǎn)述其基本算法（Berger et al., 1998）.首先該模式計(jì)算了逃逸層溫度：

表4 DTM 系列模式基本情況Table 4 Brief introductions of DTM models

不同成分的數(shù)密度為：

其中，T∞0為常數(shù)，代表逃逸層平均溫度.a0i表示主要成分在120 km 高度的平均濃度.函數(shù)G(L)與MSIS 模式定義相同.fi(z)為基于擴(kuò)散平衡方程得到的各成分隨高度分布曲線.最后疊加所有成分質(zhì)量密度得到總的大氣密度.

2.4 最新版本信息

目前上述三個(gè)常見(jiàn)熱層大氣經(jīng)驗(yàn)?zāi)Ｊ阶钚掳姹痉謩e為MSIS 2.0、JB2008 和DTM2013，具體信息如表5 所示（Emmert, 2015; 翁利斌, 2019）.

表5 熱層大氣密度經(jīng)驗(yàn)?zāi)Ｐ妥钚掳姹拘畔able 5 Summary of the selected empirical thermosphere models

3 熱層大氣模式優(yōu)化研究進(jìn)展

隨著空間技術(shù)的發(fā)展，特別是精密定軌及星載高精度加速度計(jì)的發(fā)展和應(yīng)用，使用低軌衛(wèi)星原位探測(cè)資料獲得了大量高精度大氣密度反演數(shù)據(jù).近年來(lái)基于這些高精度密度探測(cè)數(shù)據(jù)，熱層大氣模式不斷得到優(yōu)化和發(fā)展.本節(jié)主要介紹基于高精度探測(cè)數(shù)據(jù)的大氣經(jīng)驗(yàn)?zāi)Ｊ絻?yōu)化方案.

3.1 熱層大氣逃逸層溫度建模

大部分熱層大氣經(jīng)驗(yàn)?zāi)Ｊ浇⒃陟o力平衡假設(shè)基礎(chǔ)上，一般先計(jì)算得到溫度，然后根據(jù)大氣溫度與成分的關(guān)系獲取不同高度大氣成分的分布，最終綜合所有大氣成分獲得總的大氣密度.在計(jì)算大氣溫度時(shí)空分布時(shí)，逃逸層溫度是一個(gè)重要參量，直接影響著大氣溫度廓線結(jié)果.實(shí)際的逃逸層溫度難以觀測(cè)，Weng 等（2017）利用熱層大氣探測(cè)密度提取出等效的逃逸層溫度，構(gòu)建了逃逸層溫度模型（Exospheric Temperature Model, ETM），利用建模得到的逃逸層溫度替代原有經(jīng)驗(yàn)?zāi)Ｐ椭械奶右輰訙囟?，然后借助MSISE00 計(jì)算框架獲得熱層大氣溫度、成分和密度信息，從而實(shí)現(xiàn)了對(duì)一定高度范圍內(nèi)的熱層大氣密度的修正優(yōu)化.

3.1.1 逃逸層溫度提取

基于熱層大氣密度探測(cè)數(shù)據(jù)提取逃逸層溫度的基本算法是通過(guò)不斷修正經(jīng)驗(yàn)?zāi)Ｐ椭械奶右輰訙囟?，使得模式?jì)算的大氣密度結(jié)果與探測(cè)數(shù)據(jù)一致，則對(duì)應(yīng)的溫度即等效逃逸層溫度，具體提取流程如圖8 所示.圖9 展示了基于CHAMP 衛(wèi)星密度探測(cè)數(shù)據(jù)提取的等效逃逸層溫度和MSIS 模式結(jié)果的對(duì)比，可以看出，經(jīng)驗(yàn)?zāi)Ｊ浇o出的逃逸層溫度明顯偏高.

圖8 利用MSISE00 經(jīng)驗(yàn)?zāi)Ｊ胶蜔釋哟髿饷芏忍綔y(cè)值提取等效逃逸層溫度流程圖Fig.8 The exospheric temperature derived from MSISE00 and thermospheric density

圖9 CHAMP 衛(wèi)星熱層大氣密度（a）和提取的逃逸層溫度（b），以及對(duì)應(yīng)的MSISE00 模式日均值結(jié)果Fig.9 Thermospheric density (a) and exospheric temperature(b) derived from CHAMP observation and MSISE00 empirical model

3.1.2 多項(xiàng)式建模方法

考慮到衛(wèi)星探測(cè)數(shù)據(jù)時(shí)空分布、計(jì)算效率特別是溫度分布特征等因素，Weng 等（2017）在建模時(shí)忽略了經(jīng)度變化的影響，同時(shí)加入MSISE00 經(jīng)驗(yàn)?zāi)Ｊ降挠?jì)算結(jié)果來(lái)約束模式，以保證結(jié)果的合理性，并有效解決了稀疏探測(cè)數(shù)據(jù)帶來(lái)的建模不穩(wěn)定問(wèn)題，計(jì)算公式如下：

上式中，Tm（m=0, 1, 2, 3）表示：

其中，T*為MSISE00 經(jīng)驗(yàn)?zāi)Ｊ接?jì)算的逃逸層溫度結(jié)果.

式（27）中：

式（28）表示熱層大氣密度日變化和半日變化，t表示地方時(shí)；式（29）表示熱層大氣密度年變化和半年變化，td表示天數(shù).

3.1.3 逃逸層溫度模型結(jié)果與分析

最后基于CHAMP 衛(wèi)星探測(cè)數(shù)據(jù)擬合模型參數(shù).為了獨(dú)立檢驗(yàn)ETM 模式的模擬能力，使用GRACE 衛(wèi)星大氣密度及其提取的逃逸層溫度進(jìn)行對(duì)比評(píng)估.圖10 給出了MSISE00 和ETM 模式模擬的逃逸層溫度和熱層大氣密度與探測(cè)值的相對(duì)誤差統(tǒng)計(jì)結(jié)果.從圖10a 可以看出MSISE00 模式的逃逸層溫度結(jié)果整體高于GRACE 衛(wèi)星探測(cè)值，偏差和標(biāo)準(zhǔn)差分別為5.88%和6.96%.相比較而言，ETM模式的誤差結(jié)果基本呈現(xiàn)正態(tài)分布，偏差和標(biāo)準(zhǔn)差分別為1.24%和4.29%.圖10b 為大氣密度與探測(cè)值的誤差統(tǒng)計(jì)結(jié)果，可以看出MSISE00 模式模擬的熱層大氣密度與探測(cè)值相對(duì)偏差為41.52%，而ETM 模式則為10.58%.經(jīng)過(guò)計(jì)算，這兩個(gè)模式模擬的密度與GRACE 衛(wèi)星熱層大氣密度探測(cè)值的標(biāo)準(zhǔn)差分別為45.27%和24.75%.由此可見(jiàn)，與MSISE00 模式相比，使用GRACE 數(shù)據(jù)進(jìn)行獨(dú)立檢驗(yàn)表明ETM 模式計(jì)算得到的逃逸層溫度和熱層大氣密度精度得到了較大幅度提升.

圖10 MSISE00 以及ETM 模式模擬的逃逸層溫度（a）和GRACE 衛(wèi)星所在高度熱層大氣密度（b）與探測(cè)值的相對(duì)誤差統(tǒng)計(jì)結(jié)果Fig.10 The statistics of relative errors of (a) exospheric temperature and (b) density with GRACE measurements for the ETM (blue line) and MSISE00 (red line)

3.2 熱層大氣動(dòng)態(tài)修正方案

熱層大氣動(dòng)態(tài)修正也是模式優(yōu)化的一種有效方案.其基本原理是融合CHAMP、GRACE 等多源衛(wèi)星探測(cè)數(shù)據(jù)，結(jié)合熱層電離層電動(dòng)力學(xué)耦合模型（TIEGCM）模擬結(jié)果，構(gòu)建物理模型+數(shù)據(jù)驅(qū)動(dòng)的熱層大氣密度動(dòng)態(tài)修正模式（Ruan et al., 2018）.

3.2.1 基于衛(wèi)星探測(cè)的熱層大氣密度數(shù)據(jù)反演逃逸層溫度

以靜力學(xué)平衡條件為依據(jù)，熱層大氣密度的垂直變化可以由熱層大氣頂部溫度和底部溫度表征.動(dòng)態(tài)修正的方法是利用經(jīng)驗(yàn)?zāi)Ｊ剑?jì)算不同大氣溫度條件下的熱層大氣密度垂直分布，建立熱層大氣溫度與熱層大氣密度的定量關(guān)系；利用CHAMP、GRACE 等衛(wèi)星探測(cè)的局地?zé)釋哟髿饷芏?，采用迭代方法，獲取對(duì)應(yīng)的逃逸層溫度，并開(kāi)展相對(duì)應(yīng)的分析和研究，以完成對(duì)熱層密度及逃逸層溫度信息的收集、處理及其時(shí)空分布特征的提取和積累.

3.2.2 基于數(shù)值理論模型，提取熱層大氣正交基函數(shù)

熱層大氣變化的外部驅(qū)動(dòng)源豐富、內(nèi)部耦合復(fù)雜，物理自洽是動(dòng)態(tài)修正的重要考慮因素.根據(jù)目前對(duì)熱層大氣的理論研究成果，使用電離層熱層耦合模式模擬獲取熱層大氣頂部溫度（即逃逸層溫度）隨太陽(yáng)和地磁活動(dòng)水平的時(shí)空分布，獲取大量的模擬數(shù)據(jù)源，使用經(jīng)驗(yàn)正交基方法對(duì)數(shù)據(jù)源進(jìn)行建模，獲取熱層大氣頂部溫度不同模態(tài)的時(shí)空分布特征，并基于這些模態(tài)表征真實(shí)大氣的變化規(guī)律.如圖11展示了第一至第五階正交基函數(shù)（PC1-PC5）隨緯度與地方時(shí)變化的分布圖，及其在理論模型數(shù)據(jù)庫(kù)中所占比重.

圖11 （a-e）第一至第五階主成分系數(shù)（PC1-PC5）隨緯度與地方時(shí)變化的分布圖，（f）及其在理論模型數(shù)據(jù)庫(kù)中占所有變化的比重（修改自Ruan et al., 2018）Fig.11 Variations of the basis functions (a-e) PC1-PC5 as a function of local time and latitude and (f) their relative contributions to the total variance (modified from Ruan et al., 2018)

3.2.3 利用稀疏數(shù)據(jù)驅(qū)動(dòng)技術(shù)，構(gòu)建全球熱層大氣變化

利用觀測(cè)數(shù)據(jù)計(jì)算理論模型提取的正交基函數(shù)幅度可以有效地將理論模型和實(shí)際探測(cè)數(shù)據(jù)結(jié)合起來(lái)，并實(shí)現(xiàn)對(duì)熱層大氣的全球重構(gòu).在實(shí)際操作中，采用最小二乘法，通過(guò)多源衛(wèi)星探測(cè)數(shù)據(jù)來(lái)獲取理論模型正交基函數(shù)的變化幅度，進(jìn)而構(gòu)建全球熱層大氣變化.圖12 為全球熱層動(dòng)態(tài)反演示意圖.圖12a 中藍(lán)色實(shí)線代表CHAMP 衛(wèi)星數(shù)據(jù)反演的逃逸層溫度變化，藍(lán)色虛線代表理論模型模擬的相同位置的逃逸層溫度變化，紅色實(shí)線代表重構(gòu)結(jié)果.從圖中可以看到，由于CHAMP 衛(wèi)星運(yùn)行軌道的晝夜差異，其探測(cè)得到的逃逸層溫度也具有顯著的高低變化，而理論模型（TIEGCM）能夠較好地模擬該變化特征，但具有顯著的絕對(duì)值差異，因此，理論模型并不能直接用于預(yù)測(cè)實(shí)際大氣的變化.與直接的理論模型模擬不同，基于探測(cè)數(shù)據(jù)的重構(gòu)結(jié)果（紅線）不僅具有理論模型的變化特征，而且與觀測(cè)結(jié)果具有較好的絕對(duì)值一致性，能夠滿足實(shí)際大氣預(yù)報(bào)需求.結(jié)果表明，結(jié)合理論模型與探測(cè)數(shù)據(jù)，可有效構(gòu)建全球大氣動(dòng)態(tài)模型.

圖12 探測(cè)數(shù)據(jù)驅(qū)動(dòng)技術(shù)示意圖（修改自Ruan et al.,2018）Fig.12 Schematic diagram for the data-driven process (modified from Ruan et al., 2018)

3.2.4 多源衛(wèi)星數(shù)據(jù)融合的參數(shù)化動(dòng)態(tài)建模

融合CHAMP、GRACE 等多源衛(wèi)星探測(cè)數(shù)據(jù)，基于理論模型提取的正交基函數(shù)及衛(wèi)星實(shí)時(shí)探測(cè)數(shù)據(jù)，通過(guò)最小二乘法擬合的方式獲取正交基函數(shù)幅度系數(shù)，并對(duì)逃逸層溫度進(jìn)行動(dòng)態(tài)反演.Ruan 等（2018）將CHAMP 衛(wèi)星繞地球每運(yùn)行一圈所記錄的探測(cè)數(shù)據(jù)作為數(shù)據(jù)驅(qū)動(dòng)源，擬合理論模型構(gòu)建全球熱層大氣變化.

參數(shù)化過(guò)程具體如下：

其中TSG代表太陽(yáng)活動(dòng)和地磁活動(dòng)的作用，TTime代表各正交基系數(shù)的年變化和世界時(shí)（經(jīng)度）變化.具體公式如下：

以上F10.7為提前1 天的太陽(yáng)10.7 cm 射電流量，F(xiàn)10.7A為F10.7 指數(shù)81 天滑動(dòng)平均值，地磁Ap1—Ap4指數(shù)分別代表天平均值（Ap1）、提前3 小時(shí)的數(shù)值（Ap2）、提前6 小時(shí)的數(shù)值（Ap3）以及提前9 小時(shí)的數(shù)值（Ap4），d代表年天數(shù)，t代表世界時(shí)，其余參數(shù)可通過(guò)最小二乘法擬合獲得.

3.3 熱層數(shù)據(jù)同化與預(yù)報(bào)

3.3.1 數(shù)據(jù)同化

近年來(lái)，數(shù)據(jù)同化方法在地球科學(xué)領(lǐng)域的應(yīng)用得到快速發(fā)展.通過(guò)融合探測(cè)與模擬信息，數(shù)據(jù)同化可以提供較觀測(cè)更完整、較模擬更準(zhǔn)確的目標(biāo)參數(shù)狀態(tài).三維/四維變分同化（Bannister, 2017）與集合卡爾曼濾波（Evensen, 2003）等方法自提出以來(lái)獲得了大量關(guān)注，且已被廣泛應(yīng)用于氣象、海洋和陸面同化系統(tǒng).

為了更好地估計(jì)熱層大氣溫度、密度參數(shù)，研究人員發(fā)展了一系列熱層大氣數(shù)據(jù)同化方法（如,Cantrall et al., 2019; Ruan et al., 2018; Sutton, 2018;Weimer et al., 2020）.Picone 等（2005）及Gondelach和Linares（2020）通過(guò)同化基于空間目標(biāo)軌道跟蹤數(shù)據(jù)（TLE）獲得的大氣密度信息估計(jì)了全球熱層大氣密度分布.CHAMP 和GRACE 加速度計(jì)反演密度也廣泛地被應(yīng)用于熱層大氣密度同化工作中（ Ren and Lei, 2020; Ruan et al., 2018; Sutton,2018）.常用的同化算法包括集合卡爾曼濾波算法（Codrescu et al., 2018; Matsuo et al., 2013）、集合最優(yōu)插值方法（Murray et al., 2015）、集合平方根濾波方法（Cantrall et al., 2019）等，用于融合熱層觀測(cè)與經(jīng)驗(yàn)/理論數(shù)值模擬信息.但是，在這些傳統(tǒng)同化方法中，為避免觀測(cè)信息與距離較遠(yuǎn)處信息的“假相關(guān)”，一般通過(guò)“局地化”方法對(duì)同化過(guò)程進(jìn)行限制，使觀測(cè)信息僅可調(diào)整探測(cè)位置附近的大氣狀態(tài).由于熱層探測(cè)極為缺乏，如低軌衛(wèi)星需～1.5 小時(shí)獲得環(huán)繞地球一圈的探測(cè)數(shù)據(jù)，因此傳統(tǒng)方法在高層大氣數(shù)據(jù)同化方面具有較大的局限性.

熱層作為一個(gè)復(fù)雜的非線性耦合系統(tǒng)，其時(shí)空演化依賴于各種參數(shù)的協(xié)同作用，如溫度、風(fēng)、密度以及中性組分.由于熱層大氣環(huán)流模型（GCM）可以較好地模擬高層大氣不同時(shí)空尺度的復(fù)雜變化，因此考慮基于理論模式的數(shù)據(jù)同化可以提供對(duì)熱層更好的預(yù)測(cè)與預(yù)報(bào)能力.Ren 和Lei（2020）利用CHAMP 衛(wèi)星軌道密度探測(cè)，迭代調(diào)整熱層大氣環(huán)流模型中太陽(yáng)與地磁驅(qū)動(dòng)指數(shù)，從而使模擬與探測(cè)的軌道密度一致，構(gòu)建了基于理論模式的熱層數(shù)據(jù)同化系統(tǒng)PIDA.他們進(jìn)一步使用GRACE 衛(wèi)星軌道密度和GUVI 臨邊觀測(cè)對(duì)PIDA 進(jìn)行評(píng)估.比較結(jié)果顯示，PIDA 同化結(jié)果與CHAMP 軌道密度非常一致.然而，通過(guò)與軌道高度較高的GRACE 衛(wèi)星軌道密度以及TIMED-GUVI 臨邊觀測(cè)進(jìn)行對(duì)比，結(jié)果顯示PIDA 可以很好地估計(jì)全球熱層密度的時(shí)空演化，但該同化密度存在高度依賴的系統(tǒng)偏差.

圖13a、13c 展示了GUVI 臨邊觀測(cè)以及TIEGCM 和PIDA 估計(jì)的平均溫度與密度剖面，相應(yīng)的模擬與觀測(cè)比值的統(tǒng)計(jì)平均在圖13b、13d 中給出.可以看出，與TIEGCM 模型對(duì)比，同化后的熱層大氣密度與觀測(cè)結(jié)果更為吻合，而同化后的溫度偏差增大.主要原因是，Kp 與F10.7 的變化主要影響高熱層的狀態(tài)，通過(guò)調(diào)整PIDA 模型外驅(qū)動(dòng)參數(shù)以符合衛(wèi)星軌道密度探測(cè)的方法，無(wú)法對(duì)低熱層（低于150 km）狀態(tài)產(chǎn)生顯著影響.因此，為了與CHAMP 軌道密度相一致，PIDA 模型低估了熱層大氣溫度，使得熱層大氣標(biāo)高相應(yīng)降低.由于溫度的變化對(duì)熱層密度的影響會(huì)在高度上持續(xù)累計(jì)，從而導(dǎo)致PIDA 同化后的熱層密度在低高度比真實(shí)情況偏高，而在更高的高度比探測(cè)結(jié)果偏低.

圖13 2004年第80—280 天平均的GUVI 臨邊觀測(cè)（a）溫度與（c）密度高度剖面（灰色：觀測(cè)；藍(lán)色：TIEGCM；紅色：PIDA）.（b）溫度與（d）密度的模擬與觀測(cè)比值的統(tǒng)計(jì)平均（藍(lán)色：TIEGCM；紅色：PIDA）.統(tǒng)計(jì)標(biāo)準(zhǔn)差由圖（b）和圖（d）中相應(yīng)的彩色陰影表示（修改自Ren and Lei, 2020）Fig.13 The altitude profiles of the temporal and spatial averaged (a) temperature and (c) density from Limb observations (gray dotted lines) and the corresponding simulation results from TIEGCM (blue dotted lines) and PIDA (red dotted lines).The mean ratio of the simulation results (blue line for TIEGCM; red line for PIDA) to observations for (b) temperature and (d) density.The standard deviations were marked by the corresponding colored shading in (b) and (d) (modified from Ren and Lei, 2020)

在流體靜力學(xué)條件下，熱層密度的高度分布與熱層底邊界的密度和溫度顯著相關(guān).結(jié)果顯示，TIEGCM 明顯高估了低熱層的溫度以及密度.盡管F10.7 與Kp 可以顯著影響高熱層狀態(tài)，但是其變化無(wú)法對(duì)低熱層大氣產(chǎn)生顯著影響.因此，理論模式對(duì)低熱層估計(jì)的較大偏差使得PIDA 對(duì)不同高度密度的估計(jì)仍存在較大偏差.因此，基于理論模式的數(shù)據(jù)同化系統(tǒng)應(yīng)同時(shí)考慮低熱層大氣的狀態(tài)估計(jì).此外，研究還指出，熱層數(shù)據(jù)同化需融合不同高度的多種參數(shù)信息，以消除同化系統(tǒng)偏差.

3.3.2 基于智能優(yōu)化算法的熱層預(yù)報(bào)研究

在傳統(tǒng)的高層大氣數(shù)值預(yù)報(bào)應(yīng)用中，一般先通過(guò)特定數(shù)據(jù)同化方法，對(duì)實(shí)時(shí)熱層大氣狀態(tài)（如質(zhì)量密度、溫度和中性成分等）進(jìn)行盡可能準(zhǔn)確地估計(jì).然后，使用預(yù)報(bào)模式對(duì)熱層大氣狀態(tài)進(jìn)行預(yù)報(bào).由于高層大氣理論模式可以重現(xiàn)不同時(shí)空尺度下復(fù)雜的熱層大氣變化特性，其通常被選用為預(yù)報(bào)模式.但是，由于模式中的某些缺陷（如不確定性參數(shù)、缺少物理過(guò)程等）以及太陽(yáng)與地磁活動(dòng)預(yù)測(cè)困難等限制，理論模式預(yù)報(bào)誤差會(huì)隨時(shí)間迅速增大，從而導(dǎo)致有效預(yù)報(bào)時(shí)長(zhǎng)局限在幾天之內(nèi).對(duì)熱層大氣長(zhǎng)期可靠的預(yù)報(bào)仍是科學(xué)界挑戰(zhàn)性的難題之一.

Ren 和Lei（2022）構(gòu)建了一種新型的基于全球探測(cè)的熱層預(yù)報(bào)模式GOFT，以實(shí)現(xiàn)熱層大氣數(shù)十天的穩(wěn)定預(yù)報(bào).GOFT 使用TIEGCM 作為基本理論模式.GOFT 通過(guò)智能優(yōu)化粒子濾波算法，對(duì)TIEGCM 中可能導(dǎo)致熱層狀態(tài)估計(jì)顯著誤差的一系列不確定性參數(shù)進(jìn)行智能調(diào)整，實(shí)現(xiàn)了熱層大氣預(yù)報(bào).他們利用TIEGCM 以及廣泛使用的熱層經(jīng)驗(yàn)?zāi)Ｊ組SISE00 與GOFT 預(yù)報(bào)結(jié)果進(jìn)行對(duì)比.結(jié)果表明，在太陽(yáng)和地磁活動(dòng)給出的情況下，GOFT 可提供可靠的超30 天預(yù)報(bào).此外，基于GOFT 有望實(shí)現(xiàn)復(fù)雜的非線性耦合高層大氣系統(tǒng)的長(zhǎng)期預(yù)報(bào).

基于智能優(yōu)化粒子濾波算法的熱層預(yù)報(bào)模式GOFT 邏輯示意圖如圖14 所示.粒子濾波是非線性非高斯系統(tǒng)狀態(tài)估計(jì)的常用數(shù)據(jù)同化算法（如，Djuric et al., 2003; van Leeuwen, 2009）.GOFT 模式粒子群組中“粒子”為不確定性參數(shù)被隨機(jī)擾動(dòng)后的TIEGCM 模式.在當(dāng)前版本調(diào)整了12 個(gè)對(duì)熱層狀態(tài)存在顯著影響的不確定性參數(shù)（如潮汐、焦耳加熱因子等）.GOFT 首先在每個(gè)時(shí)間步將粒子更新一定的時(shí)長(zhǎng).然后根據(jù)模擬-觀測(cè)誤差獲得每個(gè)粒子權(quán)重.這些權(quán)重粒子便用于熱層預(yù)報(bào).為了獲得更優(yōu)的粒子群，在下一步同化新的觀測(cè)之前，通過(guò)特定的智能優(yōu)化算法對(duì)模式不確定參數(shù)（即粒子群中的粒子）進(jìn)行選擇變異.

圖14 預(yù)報(bào)模式GOFT 示意圖.彩色圓點(diǎn)表示GOFT 中的粒子，其大小與粒子權(quán)重成正比.紅色曲線表示概率密度函數(shù)（修改自Ren and Lei, 2022）Fig.14 Schematic view of the forecast model GOFT.The solid dot symbols indicate the particles in the GOFT with the symbol sizes proportional to the weighting of the particles.The red curve represents the probability density function (modified from Ren and Lei, 2022)

GOFT 模式的關(guān)鍵目標(biāo)是提供較現(xiàn)有模型更為可靠的熱層長(zhǎng)期預(yù)報(bào).圖15 比較了GOFT 和TIEGCM 對(duì)GUVI 臨邊密度的30 天預(yù)報(bào)誤差.如圖15a 所示，在前15 天預(yù)報(bào)時(shí)段內(nèi)，GOFT 對(duì)390 km 高度軌道密度的預(yù)報(bào)誤差低于15%，相較于TIEGCM 的預(yù)報(bào)誤差分別提高了90%.此外，我們還與熱層經(jīng)驗(yàn)?zāi)Ｐ蚆SISE00 預(yù)報(bào)結(jié)果進(jìn)行比較，對(duì)比結(jié)果顯示GOFT 預(yù)報(bào)精度同樣提高了80%.在后續(xù)15 天預(yù)報(bào)時(shí)段，GOFT 對(duì)GUVI 臨邊密度的預(yù)報(bào)誤差仍低于20%，與TIEGCM 相比提高了～80%～90%，而與MSISE00 相比提高了20%～70%.評(píng)估結(jié)果表明，預(yù)報(bào)模式GOFT 對(duì)熱層質(zhì)量密度的30 天預(yù)報(bào)精度比理論模式TIEGCM 提高了70%～90%.該方法為復(fù)雜的高層大氣預(yù)報(bào)系統(tǒng)提供了未來(lái)的改進(jìn)前景.需要指出的是，GOFT 對(duì)熱層的預(yù)報(bào)需要提供F10.7 和Kp 指數(shù)預(yù)報(bào)，這里使用了真實(shí)的Kp 與F10.7 指數(shù)以評(píng)估GOFT 在預(yù)測(cè)熱層方面的性能.

圖15 GOFT（紅色點(diǎn)線）和TIEGCM（藍(lán)色點(diǎn)線）對(duì)（a）CHAMP 和（b）GRACE 軌道平均密度的30 天平均預(yù)報(bào)誤差隨預(yù)報(bào)時(shí)長(zhǎng)的變化.陰影表示統(tǒng)計(jì)標(biāo)準(zhǔn)差.（c）TIEGCM 和（d）GOFT 對(duì)150～600 km軌道平均臨邊密度的30 天預(yù)報(bào)誤差隨預(yù)報(bào)時(shí)長(zhǎng)和高度的變化（修改自Ren and Lei, 2022）Fig.15 Statistical results for the 30-day forecast from the GOFT and TIEGCM.The GOFT (red line) and TIEGCM (blue line) averaged relative forecasting errors with the standard deviations (corresponding colored shading) for the orbital mean mass density from(a) CHAMP and (b) GRACE satellites during the 30-day forecast interval.The statistical average of the relative forecasting errors for limb orbital mean mass density by TIMED-GUVI from the (c) TIEGCM and (d)GOFT (modified from Ren and Lei, 2022)

4 熱層大氣響應(yīng)典型空間事件特征研究進(jìn)展

熱層大氣對(duì)空間事件的響應(yīng)是當(dāng)前熱層研究的熱點(diǎn).在一些突發(fā)空間事件中，例如磁暴、耀斑、日食等，熱層大氣通常表現(xiàn)出短時(shí)間尺度的變化特征.早期關(guān)于磁暴期間熱層大氣的響應(yīng)特征研究主要源于對(duì)空間目標(biāo)軌道跟蹤數(shù)據(jù)的分析，Jacchia（1959b）基于Sputnik 3 軌道數(shù)據(jù)首次證實(shí)了磁暴會(huì)造成熱層大氣密度的顯著增長(zhǎng)（圖16）.近年來(lái)隨著星載加速度計(jì)的應(yīng)用，高時(shí)空分辨率的大氣密度觀測(cè)數(shù)據(jù)揭示了更多熱層大氣在磁暴期間的精細(xì)變化特征（Bruinsma et al., 2006; Lei et al., 2011,2012; Liu et al., 2010, 2011; Lühr et al., 2004; Sutton et al., 2005）.圖17 展示了CHAMP（400 km）和GRACE（490 km）衛(wèi)星觀測(cè)到的2003年11月19—21日磁暴期間熱層大氣在不同高度上的變化特征.磁暴期間，極區(qū)附近熱層大氣密度迅速增加，隨后極區(qū)附近的密度增強(qiáng)區(qū)域開(kāi)始向赤道傳播，大約4 個(gè)小時(shí)后赤道附近大氣密度增加.CHAMP 觀測(cè)表明南北半球大氣響應(yīng)較為對(duì)稱，而GRACE 則觀測(cè)到熱層大氣在南半球的響應(yīng)強(qiáng)于北半球，因此熱層大氣對(duì)磁暴響應(yīng)的半球不對(duì)稱性可能存在高度和地方時(shí)相關(guān)性（Li and Lei, 2021b）.從全球尺度而言，磁暴期間熱層大氣密度增加了300%～800%.磁暴導(dǎo)致的大氣密度增強(qiáng)會(huì)加劇航天器軌道衰減.極端磁暴條件下，航天器軌道衰減率會(huì)增大至平靜時(shí)期的8 倍，中等強(qiáng)度磁暴也會(huì)造成軌道衰減率100%～150%的增加（Li and Lei, 2021a）.另外，共轉(zhuǎn)相互作用區(qū)（CIR）驅(qū)動(dòng)的地磁擾動(dòng)雖然強(qiáng)度弱，但是持續(xù)時(shí)間較長(zhǎng)，因而也會(huì)對(duì)衛(wèi)星軌道造成顯著影響（Chen et al., 2012, 2014）.磁暴結(jié)束后，熱層大氣進(jìn)入恢復(fù)階段，Lei 等（2012）基于CHAMP 和GRACE 衛(wèi)星觀測(cè)資料指出，由于NO恢復(fù)速度較慢，NO 的冷卻效應(yīng)會(huì)導(dǎo)致熱層大氣在磁暴結(jié)束后出現(xiàn)“過(guò)冷卻”效應(yīng)，磁暴結(jié)束后的熱層大氣密度可能低于磁暴之前平靜時(shí)期水平（圖18）.

圖16 1958年7月磁暴期間的Ap 指數(shù)變化（a），以及空間目標(biāo)1958 δ1（SPUTNIK 3 rocket）運(yùn)行周期變化率（b）.航天器軌道周期變化率可以表征大氣密度的變化（修改自Pr?lss, 2011）Fig.16 Ap index (a) and the orbit period change rate of space object 1958 δ1 (b) during geomagnetic storm in July,1958.The change rate of the space object represents the orbital density (modified from Pr?lss, 2011)

圖17 CHAMP（a, b）和GRACE（c, d）衛(wèi)星觀測(cè)到的熱層大氣對(duì)2003年11月20日磁暴的響應(yīng)（修改自Bruinsma et al.,2006）Fig.17 The thermospheric response to 20 November, 2003 geomagnetic storm observed by CHAMP (a, b) and GRACE (c, d) satellites (modified from Bruinsma et al., 2006)

圖18 2003年10月28—31日磁暴期間（a）行星際磁場(chǎng)Bz 分量、地磁（b）Kp、（c）Dst 指數(shù)，（d）歸一化至390 km 高度的CHAMP 衛(wèi)星白天（紅色）和夜晚（綠色）軌道平均密度，（e）基于TIMED/SABER 觀測(cè)的100～200 km 高度上NO 冷卻率的軌道平均，以及（f）125 km 高度上白天（紅色）、夜晚（藍(lán)色）的NO 平均冷卻率.其中密度單位為10-12 kg/m3.圖（d）中虛線表示最平靜時(shí)期的大氣密度；NO 冷卻率單位為107 erg/cm3/s（修改自Lei et al., 2012）Fig.18 Variations of (a) interplanetary magnetic field Bz, geomagnetic (b) Kp and (c) Dst indices, (d) dayside (red)and nightside (green) orbital averaged densities from CHAMP (normalized to 390 km), orbital averaged NO cooling rate from TIMED/SABER (e) between 100 and 200 km, and (f) dayside (red) and nightside (blue) averaged NO cooling rates at 125 km during 28-31 October 2003.Note that the mass densities in (d) are in units of 10-12 kg/m3, and the dashed lines stand for mass densities during the quietest period on October 28; NO cooling rates in (e-f) are in units of 107 erg/cm3/s (modified from Lei et al., 2012)

太陽(yáng)耀斑爆發(fā)，伴隨著太陽(yáng)EUV 輻射的突然增強(qiáng)，導(dǎo)致低熱層加熱效應(yīng)增加，熱層大氣隨即膨脹.CHAMP 衛(wèi)星觀測(cè)顯示，在2003年10月28日X17 級(jí)耀斑事件中，400 km 高度上熱層大氣密度增加了30%～40%，響應(yīng)時(shí)間約為1～3 小時(shí)（Liu et al., 2007; Pawlowski and Ridley, 2008; Sutton et al.,2005）.

日全食發(fā)生時(shí)，月影區(qū)EUV 輻射迅速降低至接近0，導(dǎo)致月影區(qū)熱層溫度下降，密度降低.研究人員通過(guò)理論模式模擬了日食期間熱層大氣的演化特征（Dang et al., 2018; Harding et al.2018; Lei et al., 2018）.Li 等（2021）進(jìn)一步通過(guò)多顆衛(wèi)星聯(lián)合觀測(cè)，證實(shí)了日食不僅會(huì)造成月影區(qū)10%～25%的局地大氣密度擾動(dòng)，還會(huì)以波動(dòng)的形式對(duì)全球熱層大氣造成持續(xù)性影響（圖19）.

圖19 2010年1月15日日食結(jié)束后GOCE 衛(wèi)星軌道密度相對(duì)背景密度的絕對(duì)變化.背景密度定義為日食前8 個(gè)軌道的平均密度.密度單位為10-12 kg/m3.圖（a-c）表示觀測(cè)密度；（d-f）表示模擬結(jié)果.紅色三角形代表衛(wèi)星位置，右上角時(shí)間為當(dāng)前時(shí)間（修改自Li et al., 2021）Fig.19 The absolute changes of neutral densities after the 15 January 15, 2010 solar eclipse with respect to the background density from (a-c) GOCE observation and that (d-f) from TIE-GCM simulation.Note that the densities are in units of 10-12 kg/m3 and the background density is defined as the mean density of eight orbits before the eclipse.The red triangles stand for satellite positions at the specific times shown in the top right-hand corner of each column (modified from Li et al., 2021)

5 總結(jié)

熱層是日地空間系統(tǒng)重要的組成部分，一方面，研究熱層對(duì)于理解整個(gè)日地空間因果鏈上的相互作用具有重要的科學(xué)意義；另一方面，由于人類大部分的空間活動(dòng)都在熱層，掌握熱層大氣變化規(guī)律，構(gòu)建并優(yōu)化熱層大氣模式對(duì)于保障太空資產(chǎn)和空間活動(dòng)安全、確保航天航空任務(wù)的圓滿完成具有顯著的工程應(yīng)用價(jià)值.

本文介紹了當(dāng)前熱層大氣密度反演的方法，比對(duì)了各種方法的優(yōu)缺點(diǎn).TLE 數(shù)據(jù)具有跟蹤目標(biāo)數(shù)量多、觀測(cè)持續(xù)時(shí)間長(zhǎng)等優(yōu)勢(shì)，但觀測(cè)時(shí)間精度相對(duì)較低，因而TLE 反演密度在熱層大氣氣候?qū)W研究以及模式標(biāo)定領(lǐng)域具有一定優(yōu)勢(shì)；基于精密軌道數(shù)據(jù)可以反演獲得較高分辨率和精度的大氣密度觀測(cè)結(jié)果，且在攜帶星載加速度計(jì)的情況下可用于標(biāo)定加速度計(jì)數(shù)據(jù)；加速度計(jì)數(shù)據(jù)反演密度具有最高的時(shí)空分辨率和精度.

隨著觀測(cè)技術(shù)的發(fā)展，基于最新的觀測(cè)數(shù)據(jù)，熱層大氣模式優(yōu)化工作已取得了長(zhǎng)足的進(jìn)展.使用加速度計(jì)數(shù)據(jù)反演密度修正逃逸層溫度時(shí)，采用多項(xiàng)式或者稀疏矩陣擬合等技術(shù)取得了較好的優(yōu)化效果.此外，多源參數(shù)同化以及粒子濾波等技術(shù)也開(kāi)始應(yīng)用于熱層大氣模式優(yōu)化研究中，推動(dòng)了熱層大氣密度中長(zhǎng)期預(yù)報(bào)的研究進(jìn)展.熱層大氣研究雖然取得了諸多成果，但其數(shù)據(jù)的時(shí)空覆蓋依舊貧乏，熱層變化特性尚需進(jìn)一步探索和認(rèn)識(shí).目前的大氣模式精度也有待進(jìn)一步提高，還需要發(fā)展更多的新方法、新技術(shù)應(yīng)用于熱層大氣研究和應(yīng)用，以提升我們對(duì)熱層的認(rèn)識(shí)，并服務(wù)于衛(wèi)星軌道定位與預(yù)報(bào)等.