基于X射線能譜擬合的地球中性大氣數密度反演模擬及誤差分析

余道淳, 李保權, 劉亞寧, 李海濤*

1 中國科學院國家空間科學中心， 北京 100190 2 中國科學院大學， 北京 100049

0 引言

大氣密度的精確預報對低軌道衛(wèi)星的機動規(guī)劃、精密定軌、衛(wèi)星壽命預測以及再入飛行器的返回控制等十分重要(苗娟等,2016;Kalafatoglu Eyigule et al., 2019; 任志鵬, 2020). 隨著人們對地球大氣密度需求的不斷增加，尤其是地球中高層大氣密度，發(fā)展了多種半經驗大氣模型，比如Jacchia參考大氣(Jacchia, 1970)，DTM模型(Bruinsma et al., 2003)，MSIS模型(Picone et al., 2002; Emmert et al., 2021)等，但由于地球中高層大氣的復雜變化，半經驗模型往往存在約30%甚至更多的標準偏差(Doornbos et al., 2008; 陳旭杏等, 2013; 萬田等, 2015). 因此，需要發(fā)展更有效的方法來進一步獲取更精確的大氣密度. 比如大氣密度的原位探測方法，這是利用搭載在探空氣球或探空火箭上的載荷直接探測大氣密度的方法，但探空氣球在上升到距離地面35 km的高空時，會因內部氣壓過大而發(fā)生爆炸，所以探空氣球探測高度集中在35 km以下(黃炯等, 2014)，探空火箭的探測高度多集中在距離地面約80～100 km的高度范圍(Katsuda et al., 2021)，而且探空氣球和探空火箭作為探測大氣密度的原位測量方法，存在探測持續(xù)時間短、空間覆蓋性差等弊端，同時探測結果會因探空氣球和探空火箭空氣動力學擾動的存在而造成較大誤差；另外通過大氣遙感探測的手段也可以獲取地球大氣密度，NASA研制的大氣遙感衛(wèi)星儀器SABER(Sounding of the Atmosphere using Broadband Emission Radiometry)通過臨邊探測CO2的紅外輻射進行地球大氣層150 km以下高度范圍內中性大氣密度的反演(Rezac et al., 2015; 介陽陽等, 2018)，但大氣密度的數據質量并不樂觀，尤其是100 km以上的大氣密度(Emmert et al., 2021)，除了反演得到地球中性大氣密度，SABER通過紫外-紅外波段的臨邊探測的方法也可以獲得60～180 km海拔高度范圍內的NO2、O3、CO2、H2O等大氣單一組分的密度廓線，除了衛(wèi)星大氣遙感手段之外，地基大氣遙感探測技術也被廣泛用以大氣密度反演，比如VHF雷達、瑞利激光雷達等，它們可測量30～90 km范圍內的地球中性大氣密度，且空間分辨率高，探測精度高，但它也難以實現全球范圍內大氣中性密度的長期監(jiān)測，且受惡劣天氣的影響大(王英鑒, 1997). 星光掩星技術作為一種遙感探測技術在20世紀70年代逐步發(fā)展起來，并用于反演大氣密度等大氣參數(孫明晨等，2020a).目前，基于紫外(UV)(Meyer et al., 2005; Lumpe et al., 2007)、可見光(Bauer et al., 2012)、紅外(IR)(齊瑾等, 2007; No?l et al., 2010)、無線電(胡雄等, 2005; 宮曉艷等, 2007; Lei et al., 2007; Chou et al., 2017； Wang et al., 2018)等波段的星光掩星技術反演大氣密度已有大量的研究，但是基于紫外、可見光、紅外等波段的星光掩星技術受大氣化學、大氣熱力學的影響，模擬吸收過程十分復雜，而且基于紫外、可見光、紅外波段的星光掩星技術根據大氣成分在光譜上表現出的不同的吸收特性，主要用以測量大氣單一組分或痕量氣體的密度廓線，難以測量總的中性大氣密度，比如：紫外波段掩星技術可實現NO2、O2、O3、H2等的測量，可見光波段掩星技術可實現NO2、NO3、O2等的測量，紅外波段掩星技術可實現氣溶膠、H2O、CO2、CH4等的測量(孫明晨等，2020b).無線電GPS掩星技術根據地球大氣對無線電信號的折射作用可反演得到地球總的中性大氣密度，反演結果精度高，反演的地球中性大氣密度的海拔高度范圍集中在0～60 km之間(曾楨等, 2004).通過以上分析，發(fā)現目前的探測手段難以長時間探測地球低熱層的中性大氣密度，尤其是100～200 km高度范圍內的中性大氣密度，因此該區(qū)域成為了地球大氣的“忽略層”(Oberheide et al., 2011).

X射線掩星探測技術是一門新型的交叉學科(Katsuda et al., 2021)，而且利用X射線掩星技術探測地球大氣密度有很多優(yōu)點， X射線光子直接與原子(包括分子中的原子)的K層或L層的電子發(fā)生作用，因此X射線掩星技術不受大氣化學，大氣熱力學等因素的影響，從而降低了模擬吸收過程的數學復雜性(Determan et al., 2007). 基于X射線掩星技術對行星大氣也開展了一系列研究(Determan et al., 2007; Rahmati et al., 2020; Katsuda et al., 2021)，發(fā)現基于X射線掩星技術可以反演得到地球高中層和低熱層總的中性大氣密度，這是其他探測手段難以長時間探測到的區(qū)域(Baron et al., 2020; Zeitler et al., 2021). Determan等(2007)利用RXTE/PCA觀測蟹狀星云的X射線掩星數據以及ARGOS/USA觀測Cygnus X-2的X射線掩星數據分別反演得到了100～120 km以及70～90 km海拔高度范圍內地球總的中性大氣密度，這是基于X射線掩星技術開展地球大氣參數反演的首次研究；Katsuda等(2021)通過分析Suzaku和Hitomi觀測蟹狀星云的219次掩星數據，反演得到了70～200 km海拔高度范圍內低緯度地區(qū)的平均中性密度.本文中，我們提出了一種基于X射線能譜擬合反演大氣中性密度的新方法，這不同于之前基于X射線反演大氣密度的方法，Determan等(2007)通過擬合X射線掩星過程中的光變曲線反演得到地球中性大氣密度，Katsuda等(2021)首先通過X射線掩星過程中的能譜反演出某個海拔范圍內的大氣柱密度，然后反推出地球大氣數密度，而本文中的方法直接通過擬合X射線掩星過程中的能譜，給出NRLMSISE-00的修正值，從而實現地球中性大氣密度的反演.為了驗證該方法的可行性并計算該方法反演結果的系統(tǒng)誤差，具體如下：通過模擬NICER觀測蟹狀星云的地球大氣掩星過程，得到不同高度范圍內的能譜模型和仿真數據，利用非線性最小二乘擬合方法反演得到大氣密度，對最佳擬合模型和仿真數據的擬合優(yōu)度進行統(tǒng)計分析，并計算大氣密度反演結果和真值的測量誤差. 作為NICER的主要科學儀器，X射線定時儀(XTI)擁有前所未有的吞吐量和背景噪聲的低敏感性(Gendreau et al., 2016)，可以模擬生成高信噪比的能譜數據. 通過以上仿真計算，發(fā)現基于能量范圍為1～10 keV的能譜數據，可以反演得到100～200 km高度范圍內的地球大氣中性密度，反演結果精度較高，這為基于實測能譜數據開展大氣密度的反演工作提供了理論支持.

1 前向模型與仿真數據

1.1 前向模型的構建

地球大氣的X射線掩星探測是指利用X射線衛(wèi)星上的X射線探測器臨邊接收無窮遠處的X射線源發(fā)射的X射線光子信號，當視線方向上的X射線穿過地球大氣層時，由于X射線光子與原子之間的相互作用，X射線強度會被地球大氣衰減. X射線掩星的觀測幾何如圖1所示.

圖1 地球大氣的X射線掩星的觀測幾何Fig.1 Observation geometry of Earth′s atmosphere occultation of X-ray source

圖1中，地心為O，碰撞參數為a，它表示視線方向上距離地球表面最近的一點，該點到地球表面的距離為切點高度，用h表示，RE為地球半徑，l0為X射線衛(wèi)星的位置.

比爾-朗伯定律(Beer-Lambert Law)是光吸收的基本定律，它用來描述X射線在地球大氣中的衰減:

I(E,h)=I0(E)e-τ(E,h),

(1)

其中，I(E,h)是經過地球大氣衰減的X射線能譜，I0(E)是未衰減的X射線能譜，τ(E,h)為光學厚度，它可以表示為

(2)

其中，s表示大氣的組分，σs(E)為大氣中各組分的X射線截面，ns(h)為視線方向上大氣各組分在切點高度處的數密度.

通過模擬NICER觀測蟹狀星云(Crab Nebula)的地球大氣掩星過程，獲得掩星過程中能譜和光變曲線的前向模型(式(3)). NICER是一臺安裝在國際空間站上用來研究中子星的望遠鏡，X射線定時儀(XTI)是NICER的主要科學儀器，由56個X射線光子探測器組成，具有高吞吐量和背景噪聲的低敏感性(Gendreau et al., 2016)，可以探測軟X射線能段范圍內的光子計數. 我們選擇蟹狀星云(Crab Nenula)作為模擬觀測的目標源，因為蟹狀星云作為天空中最亮的X射線源之一，其能譜是一個簡單的冪律譜而且具有非常高的穩(wěn)定性，蟹狀星云也常作為標準X射線源進行X射線天文衛(wèi)星的標定工作(Meyer et al., 2010).

FM(E,h)=RI0(E)e-τ(E,h)+B,

(3)

式中，FM(E,h)表示掩星過程中光變曲線和能譜的前向模型，它是能量和切點高度的函數，在能量維度對FM積分，得到光變曲線前向模型，在切點高度維度對FM積分，得到能譜前向模型.R為探測器響應矩陣，它的每一行表示某一入射能量范圍內的光子在各個電子學通道內被探測到的概率密度分布，B為背景噪聲，I0(E)和τ(E,h)的具體描述見公式(1).

通過數據分析軟件HEAsoft 6.29c(https:∥heasarc.gsfc.nasa.gov/lheasoft/download.html)和標定數據庫CALDB xti20170707(https:∥heasarc.gsfc.nasa.gov/docs/heasarc/caldb/nicer/)對表1中的NICER觀測數據進行處理，得到探測器響應矩陣R，如圖2所示.

表1 選取的觀測數據與觀測源的信息Table 1 The information of the selected observations and target source

圖2 基于觀測數據的NICER的探測器響應矩陣Fig.2 NICER detector response matrix based on observed data

通過能譜擬合的標準軟件包Xspec中的wabs*powerlaw模型(Yan et al., 2018)模擬得到蟹狀星云未衰減的原始能譜(I0(E))，該能譜表示蟹狀星云發(fā)射的未穿過地球大氣的X射線能譜，也稱為參考能譜，按照以下數值固定模型中的參數：等效氫柱密度NH=3.6×1021cm-2，冪率的光子指數PhoIndex=2.11，冪率系數N=8.76keV-1·cm-2·s-1(Li et al., 2020)，得到蟹狀星云在1～10 keV的參考能譜，如圖3所示.

圖3 蟹狀星云在1～10 keV的未衰減的能譜Fig.3 Unattenuated energy spectrum of the Crab Nebula in the energy range of 1～10 keV

將50～400 km海拔范圍內的大氣層分層，將相鄰兩層之間的距離固定為500 m.基于NRLMSISE-00(Picone et al., 2002)給出50～400 km海拔范圍內大氣各組分的離散數密度，NRLMSISE-00模型大氣密度作為真值用以生成仿真數據，在本文中，我們只考慮地球大氣中氮(N，N2)和氧(O，O2)兩種主要的元素成分，基于剝洋蔥法利用Abel積分求出地球大氣每一層視線方向上的大氣柱密度(張斯敏等，2022；Sun et al. , 2022)，然后通過XCOM數據庫計算得到地球大氣各組分的X射線截面(σs(E))，從而根據公式(2)求出光學厚度(τ(E,h))，然后通過公式(3)求出掩星過程中XTI/NICER接收的X射線能譜和光變曲線的前向模型，為了簡化計算，公式(3)中的背景噪聲B固定為0.

圖4 1～10 keV能量范圍內光變曲線模型與仿真數據Fig.4 Model lightcurve and simulation data in the energy range of 1～10 keV

1.2 仿真數據的生成

前向模型(FM(E,h))是一個關于能量和切點高度的二維矩陣，將1～10 keV的能道合并，就可以得到能量范圍在1～10 keV的模型光變曲線，在該模型上添加均值為0，標準偏差為模型計數的開方的高斯噪聲，得到光變曲線的仿真數據，如圖4所示，發(fā)現在切點高度為200 km處，X射線光子計數開始衰減，在切點高度為100 km處，X射線光子計數完全衰減，因此通過仿真計算發(fā)現，在1～10 keV能量范圍內，蟹狀星云的X射線掩星發(fā)生在地球大氣層約100～200 km的海拔高度范圍內，因此通過能譜擬合反演得到地球總的中性大氣密度分布在100～200 km的海拔高度范圍內.

將NRLMSISE-00模型的大氣密度作為真值，按照切點高度維度對前向模型(FM(E,h))并道，得到模型能譜，在每個模型能譜上添加均值為0，標準偏差為模型光子計數開方的高斯噪聲，得到不同切點高度范圍內的能譜仿真數據，如圖5所示，通過能譜的仿真數據發(fā)現，在掩星過程中，隨著海拔高度的降低X射線能譜的光子計數減小，這和海拔高度降低大氣密度的增加有關.

圖5 不同切點高度范圍內的能譜仿真數據Fig.5 Simulation energy spectrum data in the different altitude ranges

2 密度反演

本文中，采用非線性最小二乘擬合的方法對不同高度的能譜進行擬合從而得到地球中性大氣密度，選擇NRLMSIS 2.0(Emmert et al., 2021)模型的大氣密度作為模型能譜的初值(ns(h))，在初值前乘以一個修正因子γ，對公式(2)作以下變形:

(4)

式中，對于各變量的描述參看公式(2).將公式(4)和公式(3)結合可以求出模型能譜，其中修正因子γ作為自由參數.

將NRLMSIS 2.0大氣密度作為初值的模型能譜和NRLMSISE-00大氣密度作為真值的能譜仿真數據進行非線性最小二乘擬合，不同海拔高度范圍內的模型能譜和能譜仿真數據的擬合結果如圖6a—j所示，每個子圖的上面板中，紅色實線為最佳擬合模型，帶有誤差棒的藍點為能譜仿真數據，每個子圖的下面板展示了不同海拔高度范圍內最佳擬合模型和能譜仿真數據之間的標準化殘差.不同海拔高度范圍內，模型能譜中修正因子γ的最佳擬合值如表2所示，通過計算不同海拔高度范圍內的最佳擬合能譜和能譜仿真數據之間的2/dof(degree of freedom)和p值(p-value)對擬合結果進行統(tǒng)計分析，發(fā)現不同海拔高度范圍內的最佳擬合模型和能譜仿真數據之間2/dof都接近1，而且p值的計算結果不太小，如表2所示，因此不同海拔高度范圍內的模型能譜和能譜仿真數據的擬合結果很好.

表2 不同海拔高度范圍內修正因子γ的擬合結果 以及最佳擬合模型和仿真數據之間的擬合優(yōu)度Table 2 The fitting results of the correction factor γ in different altitude ranges and the goodness of fit between the best fitting model and simulation data

根據不同海拔高度范圍內的模型能譜和能譜仿真數據之間的擬合結果，反演得到100～200 km海拔高度范圍內的地球大氣數密度，如圖7所示，紅色實線為NRLMSISE-00的大氣密度值，本文中它作為真值，藍色實線為NRLMSIS 2.0的大氣密度值，本文中它作為模型能譜的初值，綠色實線為基于能譜擬合得到的大氣密度反演結果，為了更好的展示三條大氣密度曲線的相對關系，三條曲線都通過除以真值進行歸一化，如圖7右面板所示，其中大氣密度反演結果的誤差(±1σ、±2σ、±3σ)也被標記在圖中.發(fā)現在115～180 km海拔高度范圍內，大氣密度的反演結果和真值吻合較好，測量誤差介于-3.67%～1.92%之間. 在100～115 km海拔高度范圍內，大氣密度的反演結果較真值偏大，測量誤差介于-0.067%～10.22%之間，這和我們選取的模型能譜和能譜仿真數據的空間尺度有關，為了獲得較大信噪比(SNR)的能譜數據，在100～200 km海拔高度范圍內每隔10 km進行能譜的提取，所以擬合結果受初值形狀的影響較大，因此尋求空間尺度和能譜數據信噪比之間的平衡十分重要，可以通過選取更高吞吐量的X射線探測儀器來獲取更小空間尺度上的能譜數據，從而減弱甚至消除初值形狀帶來的影響. 為了支持我們的觀點，我們在100～115 km海拔高度范圍內按照每隔3 km的空間尺度提取能譜仿真數據，并反演得到100～115 km范圍內的大氣密度，如圖8所示為提取出的能譜數據和最佳擬合模型的對比，表3為3 km空間尺度下不同海拔高度范圍內，修正因子γ的最佳擬合值以及能譜仿真數據與最佳擬合模型之間的擬合優(yōu)度，發(fā)現不同海拔高度范圍內能譜仿真數據和最佳擬合模型之間的擬合優(yōu)度很好，圖9為不同空間尺度下，反演結果與真值、初值的對比，其中綠色實線表示10 km空間尺度下的反演結果，青色實線表示3 km空間尺度下的反演結果，并且計算得到3 km空間尺度下，100～115 km海拔高度范圍內反演結果的測量誤差介于-1.78%～4.99%之間，發(fā)現3 km空間尺度下反演結果與真值之間的測量誤差小于10 km空間尺度下反演結果與真值之間的測量誤差.在180～200 km海拔高度范圍內，基于能譜擬合的地球大氣反演結果較真值偏小，測量誤差介于-8.03%～-0.305%之間，這和該海拔高度范圍內X射線消光不顯著有關，在該高度范圍內大部分X射線光子穿透大氣層，相對于低海拔高度范圍內X射線光子的衰減情況，180～200 km海拔高度范圍內有更少的X射線光子被吸收或散射，這使得反演結果存在較大的不確定性(Yu et al., 2022).通過以上分析，發(fā)現利用X射線能譜擬合的方法可以反演得到高精度的地球中性大氣密度，該仿真計算也為基于實測數據的能譜擬合反演地球中性大氣密度提供了理論支持.

圖6 (a)—(j)表示不同海拔高度范圍內模型能譜和仿真數據的擬合結果Fig.6 (a)—(j) represents the fitting results of model energy spectrum and simulation data in different altitudes

圖7 大氣密度反演結果與真值、初值的比較Fig.7 Comparison of atmospheric density retrieved results with truth value and initial values

表3 3 km空間尺度下100～115 km海拔高度范圍內 修正因子γ的擬合結果以及最佳擬合模型和 仿真數據之間的擬合優(yōu)度Table 3 The fitting results of the correction factor γ and the goodness of fit between the best fitting model and simulation data in the range of 100～115 km at 3 km spatial scale

圖8 3 km空間尺度下100～115 km范圍內最佳擬合能譜和能譜仿真數據的對比Fig.8 Comparison of the best fit energy spectrum and energy spectrum simulation data in the range of 100～115 km at 3 km spatial scale

圖9 不同空間尺度下大氣密度反演結果與真值、初值的比較(青色實線：3 km空間尺度反演結果，綠色實線： 10 km空間尺度反演結果)Fig.9 Comparison of atmospheric density inversion results with true values and initial values at different spatial scales (cyan solid line: retrieved result under 3 km spatial scale, green solid line: retrieved result under 10 km spatial scale)

3 結論

通過對X射線掩星過程中的能譜進行建模，生成能譜仿真數據，利用非線性最小二乘擬合法實現模型能譜和仿真數據的擬合，并利用統(tǒng)計學方法評價模型和仿真數據之間的擬合優(yōu)度，最終將地球大氣密度反演結果和真值進行比較，得到以下結論：

(1) 提出了一種X射線掩星過程中的能譜建模方法，通過模擬生成能譜仿真數據，利用非線性最小二乘方法擬合模型能譜和仿真數據反演得到地球中性大氣密度，為基于實測數據反演地球中性大氣密度提供了理論支持.

(2) 利用X射線掩星過程中地球對X射線的衰減性質，發(fā)現基于不同海拔高度范圍內1～10 keV的X射線能譜，可以通過能譜擬合的方法反演得到地球低熱層(100～200 km)總的中性大氣密度，這是其他手段難以探測到的區(qū)域，另外未來可以基于大量來自以前、現在和未來的X射線衛(wèi)星的掩星數據，通過X射線能譜擬合的方法分析該區(qū)域內地球大氣密度的時空變化特性.