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

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

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

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

0 引 言

熱層大氣受到太陽輻射、太陽風(fēng)-磁層-電離層相互作用、太陽能量粒子注入以及低層大氣波動等輸入能量的綜合作用，呈現(xiàn)復(fù)雜的時空變化特征（Hargreaves, 1992）.太陽輻射作為熱層大氣十分重要的能量輸入，對熱層有著顯著的影響（Knipp et al., 2004）.日夜更替造成熱層接收到的太陽輻射周日變化、太陽自轉(zhuǎn)導(dǎo)致太陽輻射的27 天周期變化、地球公轉(zhuǎn)引起的季節(jié)和年際變化，以及與太陽活動相關(guān)的11年周期變化等，這些因素都會顯著影響熱層大氣物理化學(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è)逃逸層溫度在太陽活動低年大約為750 K，在高年則可達(dá)1 400 K，溫度差異接近100%.熱層大氣密度隨高度增加呈指數(shù)衰減，而且受太陽活動的顯著調(diào)制.在600 km 高度區(qū)域，熱層大氣密度在太陽活動高年時比低年時約有兩個數(shù)量級的變化.耀斑、磁暴等太陽爆發(fā)事件可對熱層大氣造成短時間尺度劇烈的擾動（Bruinsma et al., 2006; Lei et al., 2011, 2012; Liu et al.,2010, 2011; Lühr et al., 2004; Sutton et al., 2005）.低層大氣擾動也會通過波動等方式影響熱層大氣（Yi?it et al., 2016）.

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

獲取有效觀測數(shù)據(jù)是熱層研究的基礎(chǔ).地表附近的大氣參數(shù)可以通過溫濕度計、氣壓計以及風(fēng)速計等設(shè)備十分便捷地測量.中間層和低熱層大氣參數(shù)則可以通過衛(wèi)星遙感、雷達(dá)、探空火箭以及氣象氣球等設(shè)備進(jìn)行探測（張曉芳和嚴(yán)衛(wèi)，2007）.然而，對于200 km 以上的大氣，以上觀測手段都難以有效測量.因此，對于熱層大氣的觀測數(shù)據(jù)十分稀少.目前，熱層大氣密度探測的主要方式是基于航天器軌道或者星載加速度計數(shù)據(jù)進(jìn)行密度反演.基于衛(wèi)星軌道數(shù)據(jù)的熱層大氣研究開始于上世紀(jì)五六十年代（Priester et al., 1967; Prussing and Conway, 1993; Thomson, 1986; Vallado, 2001）.隨著人造航天器的增多，北美防空聯(lián)合司令部（NORAD）開始以兩行數(shù)據(jù)（TLE）的形式發(fā)布空間目標(biāo)的軌道跟蹤數(shù)據(jù).很多研究人員基于TLE 數(shù)據(jù)反演了長時期的熱層大氣密度，揭示了熱層大氣的氣候?qū)W演化規(guī)律（Emmert, 2009; Marcos et al., 2005; Picone et al., 2005）.隨著星載GPS 接收設(shè)備和加速度計的發(fā)展和應(yīng)用，基于衛(wèi)星精密定軌（POD）數(shù)據(jù)和加速度計數(shù)據(jù)，可以獲取更高時空分辨率和精度的大氣密度觀測數(shù)據(jù).近年來，CHAMP、GRACE、GOCE 以及Swarm 等衛(wèi)星的高質(zhì)量大氣密度觀測結(jié)果，已被廣泛應(yīng)用于熱層大氣研究中，用來揭示熱層大氣的精細(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)星只能對全球熱層大氣進(jìn)行采樣式觀測，獲取時空連續(xù)的全球熱層大氣數(shù)據(jù)一般依賴于理論或者經(jīng)驗?zāi)Ｊ接嬎?理論模式一般基于描述質(zhì)量、動量以及能量的方程，通過數(shù)值計算的方法求解全球尺度的大氣參數(shù).理論模型可以模擬大氣演化過程及其物理機(jī)制，以及獲取較為精細(xì)的熱層大氣結(jié)構(gòu)，但其輸出的大氣參數(shù)的絕對值一般與真實結(jié)果相差較大，因此理論模式一般用于熱層大氣機(jī)理性研究.經(jīng)驗?zāi)Ｊ絼t基于參數(shù)化的經(jīng)驗方程，采用觀測數(shù)據(jù)驅(qū)動的方式進(jìn)行參數(shù)擬合，雖然物理原理涉及較少，但是輸出結(jié)果與真實較為接近.經(jīng)驗?zāi)Ｊ揭泊嬖谳^大誤差，在地磁平靜時期誤差一般約為15%，在磁暴期間則更大（Marcos, 2006）.大氣經(jīng)驗?zāi)Ｊ降恼`差對于衛(wèi)星軌道預(yù)報、軌道維持、碎片預(yù)警以及隕落預(yù)報等航天工程具有重要影響（李若曦，2022；翁利斌，2019）.

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

1 熱層大氣密度反演方法

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

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

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

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

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

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

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

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

最后疊加航天器所有面板所得結(jié)果，得到航天器整體有效橫截面積與大氣阻尼系數(shù)之乘積.如圖2 所示，以CHAMP 衛(wèi)星為例，衛(wèi)星由13 個面板組成，其有效橫截面積與大氣阻尼系數(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 是美國北美防空聯(lián)合司令部發(fā)布的空間目標(biāo)軌道跟蹤測量數(shù)據(jù).該數(shù)據(jù)庫囊括了自1957年至今超過50 000 個空間目標(biāo)的軌道跟蹤數(shù)據(jù).TLE 數(shù)據(jù)包含空間目標(biāo)編號、軌道根數(shù)以及歷元時間等信息.基于TLE 數(shù)據(jù)反演熱層大氣密度的理論最初由King-Hele（1987）給出，Picone 等（2005）進(jìn)一步做了優(yōu)化.根據(jù)式（2），理論上通過TLE 給出的歷元之間半長軸差即可以反演密度.但是，TLE給出的是航天器瞬時半長軸（或密切半長軸），且地球非球形引力場以及潮汐攝動對航天器半長軸的影響遠(yuǎn)大于大氣拖曳力.不同TLE 歷元之間半長軸的變化，除了大氣拖曳力的作用，還包括地球引力場導(dǎo)致的周期性變化.因此，Picone 等（2005）采用平均運動（mean motion）的變化反演大氣密度.將式（2）中半長軸a用平均運動nM替換得到：

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

式中nM可以通過TLE 數(shù)據(jù)直接讀取.計算需要計算積分步長上的v(t)和F(t).可以使用SGP4（Simplified General Perturbations 4）軌道積分器（Hoots and Roehrich, 1980）進(jìn)行軌道外推，從初始?xì)v元開始計算每個積分步長上的航天器位置和速度，進(jìn)而得到v(t)和F(t).基于TLE 數(shù)據(jù)進(jìn)行大氣密度反演時，一般將大氣逆彈道系數(shù)近似地設(shè)置為常數(shù).圖3 展示了基于Starshine 衛(wèi)星TLE數(shù)據(jù)的反演密度和MSISE00 結(jié)果之比對.Starshine衛(wèi)星形狀為規(guī)則球形，逆彈道系數(shù)基本保持恒定（Lean et al., 2006）.可以看出TLE 反演密度與MSISE00 結(jié)果一致，較好地描述了衛(wèi)星整個任務(wù)期間軌道密度的變化特征.由于TLE 數(shù)據(jù)時間分辨率低，因此一般只能反演沿軌道2～5 天的平均大氣密度.TLE 數(shù)據(jù)反演密度主要用于研究熱層大氣的長期變化規(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)用，越來越多的航天器可以被精密定軌.二十一世紀(jì)以來，諸如CHAMP、GRACE 以及Swarm 等衛(wèi)星均具有較高定軌精度.基于這些高精度的航天器軌道數(shù)據(jù)和精確的力學(xué)模型，可以反演獲得具有較高時空分辨率的熱層大氣密度觀測數(shù)據(jù).本節(jié)將介紹幾種基于精密軌道數(shù)據(jù)反演熱層大氣密度的方法.

1.4.1 插值微分法

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

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

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

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

圖4 基于GRACE 衛(wèi)星精密軌道數(shù)據(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 半長軸衰減法

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

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

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

圖5 基于CHAMP 精密軌道數(shù)據(jù)（藍(lán)線）與加速度計數(shù)據(jù)（紅線）之比對（修改自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ì)上都是利用力學(xué)方法，Li 和 Lei（2021b）指出可以基于能量守恒原理反演大氣密度.衛(wèi)星在軌運行，受到中心引力和其他攝動力的作用.該攝動力分為保守攝動力和非保守攝動力.航天器在保守力作用下機(jī)械能守恒，在非保守力作用下機(jī)械能發(fā)生改變.對于低軌航天器，非保守力主要為大氣拖曳力、光壓以及地球反照壓.因此，可以通過計算衛(wèi)星能量的衰減反演熱層大氣密度.這種方法不需要進(jìn)行軌道外推，能夠大幅降低計算成本.而且可以根據(jù)星歷精度選擇適當(dāng)時間步長，對于數(shù)據(jù)的要求更低，適用范圍更廣.在軌航天器總機(jī)械能可以表示為：

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

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

圖6 2006年1月1日基于GRACE 衛(wèi)星精密軌道數(shù)據(jù)反演得到的非保守力加速度（藍(lán)線）與加速度計測量結(jié)果（紅線）之比對（修改自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 基于加速度計數(shù)據(jù)反演熱層大氣密度

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

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

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

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

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

圖7 基于GRACE 衛(wèi)星精密軌道數(shù)據(jù)反演得到的2017年9月磁暴期間的熱層大氣密度（a, b）與加速度計反演密度（e,f）.密度單位為10-12 kg/m3.磁暴期間行星際磁場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)驗?zāi)Ｊ?/h2>

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

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

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

2.1 Jacchia 系列模式

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

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

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

2.2 MSIS 系列模式

MSIS（Mass Spectrometer Incoherent Scatter）模式基于非相干散射雷達(dá)、探空火箭以及質(zhì)譜儀等觀測數(shù)據(jù)建立，發(fā)展歷程和主要情況如表3 所示.下面以MSIS86 模式為例簡要介紹該模式計算方案（Hedin, 1987）.與Jacchia 模式相似，MSIS 也需要首先計算熱層大氣溫度分布，然后利用各成分對溫度的依賴性得到不同高度大氣成分分布，最后將所有成分疊加獲得大氣密度.熱層大氣溫度隨高度變化函數(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)驗函數(shù)，表征時間、緯度、太陽和磁場活動等參數(shù)對大氣溫度的影響 ，ε (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ù)大氣各成分的相對分子質(zhì)量Mi，將不同成分質(zhì)量密度疊加得到大氣密度：

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

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

2.3 DTM 系列模式

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

表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 最新版本信息

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

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

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

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

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

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

3.1.1 逃逸層溫度提取

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

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

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

3.1.2 多項式建模方法

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

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

其中，T*為MSISE00 經(jīng)驗?zāi)Ｊ接嬎愕奶右輰訙囟冉Y(jié)果.

式（27）中：

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

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

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

圖10 MSISE00 以及ETM 模式模擬的逃逸層溫度（a）和GRACE 衛(wèi)星所在高度熱層大氣密度（b）與探測值的相對誤差統(tǒng)計結(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 熱層大氣動態(tài)修正方案

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

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

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

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

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

圖11 （a-e）第一至第五階主成分系數(shù)（PC1-PC5）隨緯度與地方時變化的分布圖，（f）及其在理論模型數(shù)據(jù)庫中占所有變化的比重（修改自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ū)動技術(shù)，構(gòu)建全球熱層大氣變化

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

圖12 探測數(shù)據(jù)驅(qū)動技術(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ù)化動態(tài)建模

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

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

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

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

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

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

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

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

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

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

圖13 2004年第80—280 天平均的GUVI 臨邊觀測（a）溫度與（c）密度高度剖面（灰色：觀測；藍(lán)色：TIEGCM；紅色：PIDA）.（b）溫度與（d）密度的模擬與觀測比值的統(tǒng)計平均（藍(lán)色：TIEGCM；紅色：PIDA）.統(tǒng)計標(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)，但是其變化無法對低熱層大氣產(chǎn)生顯著影響.因此，理論模式對低熱層估計的較大偏差使得PIDA 對不同高度密度的估計仍存在較大偏差.因此，基于理論模式的數(shù)據(jù)同化系統(tǒng)應(yīng)同時考慮低熱層大氣的狀態(tài)估計.此外，研究還指出，熱層數(shù)據(jù)同化需融合不同高度的多種參數(shù)信息，以消除同化系統(tǒng)偏差.

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

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

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

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

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

圖15 GOFT（紅色點線）和TIEGCM（藍(lán)色點線）對（a）CHAMP 和（b）GRACE 軌道平均密度的30 天平均預(yù)報誤差隨預(yù)報時長的變化.陰影表示統(tǒng)計標(biāo)準(zhǔn)差.（c）TIEGCM 和（d）GOFT 對150～600 km軌道平均臨邊密度的30 天預(yù)報誤差隨預(yù)報時長和高度的變化（修改自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)展

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

圖16 1958年7月磁暴期間的Ap 指數(shù)變化（a），以及空間目標(biāo)1958 δ1（SPUTNIK 3 rocket）運行周期變化率（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)星觀測到的熱層大氣對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）行星際磁場Bz 分量、地磁（b）Kp、（c）Dst 指數(shù)，（d）歸一化至390 km 高度的CHAMP 衛(wèi)星白天（紅色）和夜晚（綠色）軌道平均密度，（e）基于TIMED/SABER 觀測的100～200 km 高度上NO 冷卻率的軌道平均，以及（f）125 km 高度上白天（紅色）、夜晚（藍(lán)色）的NO 平均冷卻率.其中密度單位為10-12 kg/m3.圖（d）中虛線表示最平靜時期的大氣密度；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)

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

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

圖19 2010年1月15日日食結(jié)束后GOCE 衛(wèi)星軌道密度相對背景密度的絕對變化.背景密度定義為日食前8 個軌道的平均密度.密度單位為10-12 kg/m3.圖（a-c）表示觀測密度；（d-f）表示模擬結(jié)果.紅色三角形代表衛(wèi)星位置，右上角時間為當(dāng)前時間（修改自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)重要的組成部分，一方面，研究熱層對于理解整個日地空間因果鏈上的相互作用具有重要的科學(xué)意義；另一方面，由于人類大部分的空間活動都在熱層，掌握熱層大氣變化規(guī)律，構(gòu)建并優(yōu)化熱層大氣模式對于保障太空資產(chǎn)和空間活動安全、確保航天航空任務(wù)的圓滿完成具有顯著的工程應(yīng)用價值.

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

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