基于X射线能谱拟合的地球中性大气数密度反演模拟及误差分析

余道淳, 李保权, 刘亚宁, 李海涛*

1 中国科学院国家空间科学中心， 北京 100190 2 中国科学院大学， 北京 100049

0 引言

大气密度的精确预报对低轨道卫星的机动规划、精密定轨、卫星寿命预测以及再入飞行器的返回控制等十分重要(苗娟等,2016;Kalafatoglu Eyigule et al., 2019; 任志鹏, 2020). 随着人们对地球大气密度需求的不断增加，尤其是地球中高层大气密度，发展了多种半经验大气模型，比如Jacchia参考大气(Jacchia, 1970)，DTM模型(Bruinsma et al., 2003)，MSIS模型(Picone et al., 2002; Emmert et al., 2021)等，但由于地球中高层大气的复杂变化，半经验模型往往存在约30%甚至更多的标准偏差(Doornbos et al., 2008; 陈旭杏等, 2013; 万田等, 2015). 因此，需要发展更有效的方法来进一步获取更精确的大气密度. 比如大气密度的原位探测方法，这是利用搭载在探空气球或探空火箭上的载荷直接探测大气密度的方法，但探空气球在上升到距离地面35 km的高空时，会因内部气压过大而发生爆炸，所以探空气球探测高度集中在35 km以下(黄炯等, 2014)，探空火箭的探测高度多集中在距离地面约80～100 km的高度范围(Katsuda et al., 2021)，而且探空气球和探空火箭作为探测大气密度的原位测量方法，存在探测持续时间短、空间覆盖性差等弊端，同时探测结果会因探空气球和探空火箭空气动力学扰动的存在而造成较大误差；另外通过大气遥感探测的手段也可以获取地球大气密度，NASA研制的大气遥感卫星仪器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等大气单一组分的密度廓线，除了卫星大气遥感手段之外，地基大气遥感探测技术也被广泛用以大气密度反演，比如VHF雷达、瑞利激光雷达等，它们可测量30～90 km范围内的地球中性大气密度，且空间分辨率高，探测精度高，但它也难以实现全球范围内大气中性密度的长期监测，且受恶劣天气的影响大(王英鉴, 1997). 星光掩星技术作为一种遥感探测技术在20世纪70年代逐步发展起来，并用于反演大气密度等大气参数(孙明晨等，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).通过以上分析，发现目前的探测手段难以长时间探测地球低热层的中性大气密度，尤其是100～200 km高度范围内的中性大气密度，因此该区域成为了地球大气的“忽略层”(Oberheide et al., 2011).

X射线掩星探测技术是一门新型的交叉学科(Katsuda et al., 2021)，而且利用X射线掩星技术探测地球大气密度有很多优点， X射线光子直接与原子(包括分子中的原子)的K层或L层的电子发生作用，因此X射线掩星技术不受大气化学，大气热力学等因素的影响，从而降低了模拟吸收过程的数学复杂性(Determan et al., 2007). 基于X射线掩星技术对行星大气也开展了一系列研究(Determan et al., 2007; Rahmati et al., 2020; Katsuda et al., 2021)，发现基于X射线掩星技术可以反演得到地球高中层和低热层总的中性大气密度，这是其他探测手段难以长时间探测到的区域(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海拔高度范围内低纬度地区的平均中性密度.本文中，我们提出了一种基于X射线能谱拟合反演大气中性密度的新方法，这不同于之前基于X射线反演大气密度的方法，Determan等(2007)通过拟合X射线掩星过程中的光变曲线反演得到地球中性大气密度，Katsuda等(2021)首先通过X射线掩星过程中的能谱反演出某个海拔范围内的大气柱密度，然后反推出地球大气数密度，而本文中的方法直接通过拟合X射线掩星过程中的能谱，给出NRLMSISE-00的修正值，从而实现地球中性大气密度的反演.为了验证该方法的可行性并计算该方法反演结果的系统误差，具体如下：通过模拟NICER观测蟹状星云的地球大气掩星过程，得到不同高度范围内的能谱模型和仿真数据，利用非线性最小二乘拟合方法反演得到大气密度，对最佳拟合模型和仿真数据的拟合优度进行统计分析，并计算大气密度反演结果和真值的测量误差. 作为NICER的主要科学仪器，X射线定时仪(XTI)拥有前所未有的吞吐量和背景噪声的低敏感性(Gendreau et al., 2016)，可以模拟生成高信噪比的能谱数据. 通过以上仿真计算，发现基于能量范围为1～10 keV的能谱数据，可以反演得到100～200 km高度范围内的地球大气中性密度，反演结果精度较高，这为基于实测能谱数据开展大气密度的反演工作提供了理论支持.

1 前向模型与仿真数据

1.1 前向模型的构建

地球大气的X射线掩星探测是指利用X射线卫星上的X射线探测器临边接收无穷远处的X射线源发射的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射线卫星的位置.

比尔-朗伯定律(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射线源之一，其能谱是一个简单的幂律谱而且具有非常高的稳定性，蟹状星云也常作为标准X射线源进行X射线天文卫星的标定工作(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))，该能谱表示蟹状星云发射的未穿过地球大气的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所示，发现在切点高度为200 km处，X射线光子计数开始衰减，在切点高度为100 km处，X射线光子计数完全衰减，因此通过仿真计算发现，在1～10 keV能量范围内，蟹状星云的X射线掩星发生在地球大气层约100～200 km的海拔高度范围内，因此通过能谱拟合反演得到地球总的中性大气密度分布在100～200 km的海拔高度范围内.

将NRLMSISE-00模型的大气密度作为真值，按照切点高度维度对前向模型(FM(E,h))并道，得到模型能谱，在每个模型能谱上添加均值为0，标准偏差为模型光子计数开方的高斯噪声，得到不同切点高度范围内的能谱仿真数据，如图5所示，通过能谱的仿真数据发现，在掩星过程中，随着海拔高度的降低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)对拟合结果进行统计分析，发现不同海拔高度范围内的最佳拟合模型和能谱仿真数据之间2/dof都接近1，而且p值的计算结果不太小，如表2所示，因此不同海拔高度范围内的模型能谱和能谱仿真数据的拟合结果很好.

表2 不同海拔高度范围内修正因子γ的拟合结果 以及最佳拟合模型和仿真数据之间的拟合优度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σ)也被标记在图中.发现在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空间尺度下不同海拔高度范围内，修正因子γ的最佳拟合值以及能谱仿真数据与最佳拟合模型之间的拟合优度，发现不同海拔高度范围内能谱仿真数据和最佳拟合模型之间的拟合优度很好，图9为不同空间尺度下，反演结果与真值、初值的对比，其中绿色实线表示10 km空间尺度下的反演结果，青色实线表示3 km空间尺度下的反演结果，并且计算得到3 km空间尺度下，100～115 km海拔高度范围内反演结果的测量误差介于-1.78%～4.99%之间，发现3 km空间尺度下反演结果与真值之间的测量误差小于10 km空间尺度下反演结果与真值之间的测量误差.在180～200 km海拔高度范围内，基于能谱拟合的地球大气反演结果较真值偏小，测量误差介于-8.03%～-0.305%之间，这和该海拔高度范围内X射线消光不显著有关，在该高度范围内大部分X射线光子穿透大气层，相对于低海拔高度范围内X射线光子的衰减情况，180～200 km海拔高度范围内有更少的X射线光子被吸收或散射，这使得反演结果存在较大的不确定性(Yu et al., 2022).通过以上分析，发现利用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海拔高度范围内 修正因子γ的拟合结果以及最佳拟合模型和 仿真数据之间的拟合优度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射线掩星过程中的能谱进行建模，生成能谱仿真数据，利用非线性最小二乘拟合法实现模型能谱和仿真数据的拟合，并利用统计学方法评价模型和仿真数据之间的拟合优度，最终将地球大气密度反演结果和真值进行比较，得到以下结论：

(1) 提出了一种X射线掩星过程中的能谱建模方法，通过模拟生成能谱仿真数据，利用非线性最小二乘方法拟合模型能谱和仿真数据反演得到地球中性大气密度，为基于实测数据反演地球中性大气密度提供了理论支持.

(2) 利用X射线掩星过程中地球对X射线的衰减性质，发现基于不同海拔高度范围内1～10 keV的X射线能谱，可以通过能谱拟合的方法反演得到地球低热层(100～200 km)总的中性大气密度，这是其他手段难以探测到的区域，另外未来可以基于大量来自以前、现在和未来的X射线卫星的掩星数据，通过X射线能谱拟合的方法分析该区域内地球大气密度的时空变化特性.