热层大气密度动态修正模型研究

刘 卫，罗冰显，龚建村，王荣兰，容建刚

（1.中国科学院 国家空间科学中心 空间天气学国家重点实验室；2.中国科学院 国家空间科学中心 太阳活动与空间天气重点实验室：北京 100190；3.中国科学院微小卫星创新研究院，上海 201210；4.中国科学院大学，北京 100049）

0 引言

热层指地球大气层中距地高度60～1000 km的区域，是地球低轨道航天器运行密集区和重要环境。大气密度预报误差是影响低轨航天器定轨和轨道预报精度的关键因素，其中地磁暴等空间环境扰动期问题尤为突出。2022 年2 月初，SpaceX 公司超40 颗“星链”卫星相继因电推进器推力无法抵御地磁暴期大气密度陡增引起的大气阻力变化，在约300 km 高度的轨道爬升过程中坠毁于大气层中[1]。

1957 年第一颗人造地球卫星进入太空后，通过轨道跟踪发现了大气阻力对航天器轨道运动的影响。随后，以高层大气为研究对象的新学科——高层大气物理学诞生[2]。20 世纪60 年代，基于卫星轨道衰变数据，Jacchia 等构建了第1 代HP、Jacchia经验大气模型[3]。目前已发展并形成以JB2008、DTM2020、MSIS2.0 和GOST04 等为代表的热层大气密度模型群[4-7]。现有热层大气模型所考虑的因素只能涵盖已知的大气密度变化规律，而高层大气复杂多变，即使对于平静空间环境，模型标准差达15%～20%，扰动期间甚至可超过100%。不规则的热层大气密度增强也会对卫星阻力和空间碎片监测预警产生重大影响[8-10]；现有经验模型无法满足诸如空间态势感知等特定场景的高精度需求。

动态大气模型校正（DCA）是一种改进或修正现有大气模型的技术[11]，一般用于热层大气密度高精度短临预报，可以提供大气密度变化的信息和相应的统计数据。DCA 技术自20 世纪80 年代首先由Nazarenko 提出以来，在航天器轨道确定应用方面受到持续关注。目前DCA 方法有基于加速度计和轨道跟踪等数据，日常提供密度校正，用于200～800 km 热层中性密度的预报[12-15]。针对热层大气密度高精度预测的场景，各航天大国均有专属的解决方案。例如：美国空军的高精度卫星拖曳模型（HASDM）[13,16]，修正后的预报1 天相对误差标准差平均下降25%[14]；俄罗斯空间观测中心基于275～650 km 多颗空间目标的轨道跟踪数据，进行直接密度修正（DDCM），应用效果良好[17-20]。

本文从大气模型扩散平衡原理着手，对热层大气温度垂直剖面上的拐点和外逸层2 个边界温度计算函数，按空间经纬度球谐系数展开修正。首先，使用推导的热层模型密度对修正系数偏导数公式计算测量矩阵，结合密度测量残差建立动态修正的理论基础；并利用泰勒公式对条件方程展开，只保留一阶项的线性化处理，通过迭代求解修正系数。然后，基于MATLAB 软件构建热层大气密度动态修正模型系统，并在磁暴和磁静期不同环境下进行应用测试。讨论动态修正系数的现报和预报应用方式；通过分析边界温度修正量的序列，总结变化规律；结合国内热层大气模型探测计划，对热层大气密度修正系统的应用前景进行总结和展望。

1 动态修正模型

1.1 修正原理

热层模型的动态修正首先要选择待修正的基础模型。模型要求理论清晰，在表征热层大气上具有拓展性，无过多拟合公式。1964 年，Jacchia 发布了第1 个热层大气模型Jacchia-64，1970 年更新为Jacchia1970 (J70)[21]，其理论清晰、拓展性好。近年来，在美国空军空间司令部推动下，基于半周年密度变化新算法和新的空间环境指数，J70 模型又拓展得到JB2006/2008 模型[4]。刘卫等[22]利用J77、DTM78、MSIS00、JB2008 等不同热层模型进行轨道预报，以精密轨道为参考，评估低地球轨道（LEO）卫星的轨道预报误差，证明Jacchia 系列模型性能优越。综上，本研究选择J70 作为待修正的基础模型。

J70 模型的外逸层温度T∞为定义参数，不同地点的外逸层温度由太阳活动指数F10.7及其81 d 均值F¯10.7，结合周日变化模型计算获得。同时基于时间t和高度的修正，包括半周年、季节变化，以及纬度φ、地方时θ和地磁活动ap变化项。外逸层温度作为温度剖面的输入，形成扩散方程积分的基础。J70 模型计算热层密度的第1 步是获取顶点的温度，

表征的是地磁指数为0 时全球夜侧温度的极小值。

外逸层温度T∞是通过Tc乘以周日变化函数D(t,φ,θ)附加地磁活动函数ΔTG(ap)得到[3]，即

J70 模型温度剖面通过边界条件定义，下边界90 km 处的温度T=183 K，拐点（125 km 处）温度Tx的梯度由dT/dh=0 给出。Tx是4 个经验参数（a,b,c,k）的函数，即

定义完备的温度剖面后，通过积分扩散方程获取大气主要成分N2、Ar、He、O2、O 的数密度，质量密度ρi(z)可通过数密度获得[21,23]，其中下标i=1, 2,3, 4, 5 分别对应成分N2、Ar、He、O2、O。

1998 年，Marcos 等首次通过估算Tc对热层大气密度模型进行修正，得到修正量ΔTc，使模型预报数据与卫星跟踪数据具有更好的一致性[24]。2005 年，Storz 等[16]对该方法进行了发展，估算ΔTc和ΔTx的球谐系数，对2 个全局温度进行修正，温度修正函数[16,25]为

因球谐修正系数计算过程所使用的测量数据构成的约束条件数量有限，无法估算无穷高阶的球谐系数，要对修正函数进行截断处理。修正模型中，边界温度修正量ΔTx取1 阶，ΔTc取2 阶，共13 项修正系数，修正后的两边界温度为：

构建的条件方程是非线性的，球谐系数求解时，首先要对修正函数进行线性化处理，然后迭代估计。将修正模型密度ρM(X)在X0处作泰勒展开，

其中偏导数矩阵A根据链导法则分为密度对温度和温度对修正系数两部分计算，其解析表达式的具体形式见文献[26]。

1.2 系统框架

以J70 模型为基础模型，热层拐点和外逸层温度补偿量的球谐展开系数为自变量，利用实测密度作为驱动（或校正）数据，构建条件方程。利用泰勒展开进行线性化，综合利用正则变换等数值方法，对方程的“病态”进行处理。迭代求解动态修正系数，构建热层大气密度动态修正系统，数据处理的流程如图1 所示。

图1 热层大气密度动态修正流程Fig.1 Flow chart of the dynamic correction system for determining thermospheric mass density

1.3 系统功能

构建的热层大气密度动态修正系统主界面如图2 所示，系统主要功能包括：

图2 热层大气模型动态修正系统运行结果图形显示界面Fig.2 Graphic interface of the computational results of dynamic correction system of thermosphere model

1）校正数据覆盖性分析。可在动态系数求解前，评估数据的空间覆盖质量。通过大量算例分析，认为需要有单轨道面单星的6 h 或双轨道面3 h及以上空间覆盖质量方能满足动态系数解算的基本需求。但空间覆盖和动态系数解算质量间的定量关系仍需深入研究。

2）按弧长连续解算动态修正系数。按设定弧长进行批量动态修正系数序列解算，可基于长期密度数据对热层边界温度建模，进一步构建热层密度模型；也可用于短期（3 d）动态修正系数预测的输入。

3）动态修正系数收敛性分析。可在校正数据覆盖性差，迭代过程可能发生动态修正系数不收敛的情况下，有效剔除不收敛弧段。

4）修正效果评估。可基于事后的大气密度测量数据，对模型的动态修正效果进行评估。

5）结果的图形展示。可将计算或分析结果以图形方式直观形象地展示给用户，便于理解。

以2003 年10 月29 日 0:00—3:00，CHAMP 和GRACE A/B 卫星3 h 弧段约700 个测量点的数据为校正数据，修正系统的运行结果如图2 所示。图2(a)为校正数据的空间分布。图2(b)为解算的球谐系数。图2(c)为校正数据时段两边界温度的补偿量：因校正数据来自分布在2 个轨道面的卫星，图中可见边界温度分布在2 条幅值相同、具有一定相位差的三角函数曲线上。可结合边界温度的物理范围，分析修正量的合理性。图2(d)为校正前/后的模型预测以及实测密度数据，其中红色、黑色和蓝色线条分别表示首次迭代、最终收敛和实测密度值。图2(e)为校正前/后（图中以红色/黑色区分）的密度残差分布，可以看到校正后的模型密度残差显著下降。

2 应用测试

热层大气模型动态修正中，动态修正系数的获取包括现报和预报2 种方式。定义前者为以相邻弧段解算的动态修正系数直接进行后续弧段大气密度模型修正；后者为以前序相邻弧段解算的动态修正系数序列为输入，利用动态系数预测模型计算相邻后续弧段的动态修正系数，然后进行大气密度模型修正。

为量化评估动态修正系统的性能，可以使用相对误差ε表示模型预报值偏离实测值的百分比，ε越小则表示系统的性能越优；还可以使用模型预报与实测密度的比值r，其值越接近1，则表示系统的性能越优。这2 个系统性能评估量的具体表达式为：

2.1 动态修正系数现报

动态修正系数现报测试中，使用CHAMP 和GRACE A/B 这3 颗卫星在2001—2010 期间的数据，按磁暴和磁静分类分别选取6 个事件进行评估测试。磁暴事件的筛选标准为当天3 h 的ap≥236；磁静事件筛选标准为连续10 d 的ap≤48。以3 颗卫星当天24 h 实测为校正数据，进行动态修正系数解算。后续连续0～72 h 实测作为评估标准，分别对J70 模型、MSIS00 模型和J70 模型校正后（DCA）的预报密度误差百分比进行统计，给出各模型预报0～72 h 相对误差均值，统计结果见表1 和表2，表中出现修正前后相对误差均值符号变化情况，这可能与评测数据地方时有关。

表1 磁暴期动态修正结果评测统计Table 1 Statistics of DCA results during magnetic storms

表2 磁静期动态修正结果评测统计Table 2 Statistics of DCA results during magnetic quiet

由表中数据可以看到：6 次地磁暴事件中，DCA算法全部有效，修正后模型的密度预报性能均优于MSIS00 模型，预报相对误差平均降低33.8 个百分点；6 次地磁平静事件中，除2003 年1 月3 日事件外，DCA 算法全部有效，较MSIS00 模型预报相对误差平均降低9.1 个百分点。可见，在3 d 内动态修正后的模型预报残差较未修正模型显著下降，这对轨道预报性能的提升非常关键。由修正失败的个例可见，将计算所得修正系数直接用于后续的动态修正预报时，并不能确保大气密度预报准确性的提升，这可能与计算修正系数时的空间环境和后续预报时的空间环境之间是否存在突变等因素相关。

2.2 动态修正系数预报

动态修正系数预报测试中，以3 h 弧长分8 个弧段计算2003 年10 月29 日的动态修正系数，计算的拐点温度和外逸层温度修正量分别如图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具有呈扇形的两极值区，相位差180°，随时间推移两极值中轴线顺时针旋转（北极俯视），这与地球自转反向；在旋转过程中两极值表现也不同，角宽度有展宽、有收缩。对于外逸层温度修正量ΔTc也存在两扇形极值区，但两极值区的相位差不是180°，且随时间变化。即，两边界温度都存在相位旋转和幅值变化特性，结合空间环境指数可进行动态修正系数的短期变化建模。

3 小结及建议

讨论动态修正的效果离不开具有一定时空覆盖的热层密度实时测量数据，分析动态修正系统的应用前景必须结合热层实测数据的来源。目前在轨空间目标数量已超过20 000 颗，其中多数为运行在热层大气之中的低地球轨道空间目标。选择一定数量、合理分布的空间目标，基于其轨道衰减来反演大气密度，无疑是获取实时大气密度数据的绝佳途径。美国空军挑选气动构型稳定、轨道面均匀分布的约80 颗空间目标的轨道精密跟踪数据，根据其轨道衰减反演全球覆盖性的热层大气密度数据，并以此数据驱动其高精度卫星拖曳模型（HASDM）。该模式已实现业务化运行，且能够满足军事领域高精度需求。

但HASDM 的运行模式是建立在美空军地基监测网设备类型完备和全球布站的基础上的，我们不能照搬，需结合具体情况来运营我们的动态修正系统。目前，国内有2 项进行中的热层大气密度探测计划，可结合热层大气模型动态修正系统进行业务化运行。具体建议如下：

1）气动力标校微小卫星项目。该项目旨在研制并发射低成本、气动稳定、携带GNSS、定位精准的大气反演和弹道系数标校微小卫星（如图5 所示），既可以基于自身精密星历反演大气密度，又可以标校相近轨道空间目标弹道信息，再与我国空间目标监视网轨道跟踪数据联合进行大气密度的解算。这样，只需少量标校星进行最小规模的组网，并根据轨道衰减情况适时补充，即可满足数据驱动需求，构建高精度大气密度动态修正系统。

图5 气动力标校微小卫星Fig.5 Aerodynamically calibrated microsatellite

2）热层大气探测卫星星座。可以该星座卫星实时探测的热层密度作为热层大气密度动态修正系统的数据驱动，动态实时进行热层模型修正，以适应高精度的热层大气密度需求场景。

4 结束语

本文以提高热层大气密度模型的分析精度为需求，对温度剖面中逃逸层及拐点两边界温度进行理论分析，通过CHAMP 和GRACE 卫星观测的大气密度数据进行拟合，对现有大气经验模型进行了优化，修正后的大气模型比以往模型性能提升，预报相对误差显著下降。通过动态修正实例和修正模型的轨道预报应用，有如下结论：

1）动态修正系数解算对校正数据的时空覆盖具有一定需求；

2）动态修正后可显著降低模型相对误差均值，对提高轨道预报精度非常有效；

3）解算的动态修正系数直接用于模型密度预报时并不保证奏效，需要建立动态修正系数序列的短期预报模型；

4）两边界温度修正量均具有两扇形极值结构，时变上具有相位旋转和幅值变化特性。

以上结论可对热层大气动态修正模型构建的工程实践提供理论支撑和应用参考。要说明的是，本文中的空间环境参数为实测值，意在算法开发阶段避免因空间环境参数预报的不准确影响算法参数估计的稳定性问题。待算法开发完毕，后续实际业务应用中，空间环境参数须采用预报值，并在动态系数预测、建模等方面开展进一步的分析工作。