一种全球临近空间大气密度建模方法及应用

程旋，肖存英，杨钧烽，胡雄，闫召爱，柳丹

（1.中国科学院国家空间科学中心 空间环境态势感知技术重点实验室，北京100190；2.北京师范大学 天文系，北京100875； 3.中国科学院大学，北京100049）

临近空间大气变化非常复杂，真实大气中包含复杂的大气波动信息，这也是真实大气偏离气候平均状态的主要原因。往往大气模型不能很好地表征这部分波动信息。为了避免模型不准确对实际飞行的影响，在工程上通过增加设计裕度的方式来避免。但是，设计裕度的增大是以牺牲其他方面的设计为代价的，如有效载荷质量和尺寸的减小。因此，一些学者开展了大气参数对高超声速飞行器影响的研究。例如，陈闽慷等［1］通过研究85 km北半球几个典型纬度和典型月份大气波动对驻点热流的影响，发现在极端情况下，SABER观测的包含大气波动的大气参数预测的热流比USSA76高40%以上。李健等［2］利用地球扰动大气模型对再入飞行器的弹道进行仿真，结果表明，大气扰动条件下，大气扰动对飞行器的落点、过载和气动热有很大影响。黄华等［3］利用美国Cape Canaveral靶场上空的气象资料建立了大气环境样本集，并利用样本集数据对再入飞行器的弹道进行了仿真，研究发现，大气热力学参量的扰动能够引起落点在射向的散布，风场的扰动对射向和横向落点的散布均有较大的影响。因此，真实大气中的扰动对临近空间飞行器的影响是不可忽略的。此外，大气扰动的分布存在时空差异，飞行过程中，飞行器在不同高度受到大气扰动的影响是不同的，还需要满足不同约束条件，使得狭窄的飞行器气动布局设计对气动误差的容忍能力较弱，提升飞行器气动力和气动热的预测准确性对气动设计精细化是非常重要的［4］。飞行器的气动特性对大气密度非常敏感，获得可靠的大气密度数据是提升气动热环境预测准确性、减少设计冗余的一种有效途径。

大气模型是获取大气密度数据的重要来源。在航空航天领域，国际上常用的大气模型主要有美国标准大气（USSA76）［5］、COSPAR国际参考大气系列模型（CIRA）［6］、质谱仪和非相干散射雷达系 列 模 型（MSIS）［7-8］和 水 平 风 系 列 模 型（HWM）［9］等，但这些模型存在一定的局限性，均不能给出大气的扰动信息。GRAM 系列模型［10］可以给出不同经纬度和月度下的温度、密度、气压和风场等参量的平均值，还可以模拟相对于平均状态的扰动剖面，但该模型代码未公开。相比国外，中国大气模型的发展水平比较落后，先后发布了国军标模型，如GJB 365.1—87、GJB 366.1—87、GJB 544—88和GJB 5601—2006等。除了GJB 5601—2006外，其余都是国际标准的引进。这些模型只能提供大气的气候平均状态，不能提供大气的扰动状态。大气模型未包含偏离平均值的大气扰动是大气模型与真实大气存在偏差的主要原因［11］。

近年来，中国学者开展了一些大气模型修正及建模方法的探索。刘一博等［12］利用酒泉2010—2011年的探空数据，对GRAM模型在酒泉上空的大气密度月均值和标准偏差进行修正，改善了酒泉上空GRAM 模型的精度。肖存英等［13］利用TIMED／SABER统计获得的11年月平均值和标准偏差，将临近空间大气密度表征为气候平均量与扰动量的加和，建立了北纬38°N大气密度扰动模型，并利用激光雷达观测数据验证了建模方法的可行性。肖存英等［13］的建模方法中采用一阶自回归模型对总扰动量的表征，未对不同尺度的扰动进行区分。实际大气中的潮汐波和行星波均具有一定的周期性，这些波动不是随机的，具有一定的规律，采用基于随机过程的一阶自回归模型对整体扰动进行表征还存在一些改进空间。

本文在文献［13］建模方法的基础上，将大气扰动分别用大尺度扰动和小尺度扰动代替，大尺度扰动采用周期性的余弦函数表征，小尺度扰动采用一阶自回归模型表征，初步建立临近空间参考大气模型，并对大气模型开展测试验证。

1 数据来源与处理方法

TIMED卫星于2001年发射，至今已连续积累长达18年的观测数据。卫星以74.1°的轨道倾角在约625 km的太阳同步轨道运行，轨道周期约为97 min。TIMED卫星上搭载4个有效载荷，即SABER、TIDI（TIMED Doppler Interferometer）、GUVI（Global Ultraviolet Imager）和SEE（Solar ultraviolet Experiment）。其中，SABER探测器通过测量大气中15μm和4.3μm的CO2临边红外辐射来反演大气温度和成分等参量［14］。卫星每天大约测量15轨数据，相邻轨道的经向间距约为24°。由于卫星进动缓慢，覆盖全球24 h地方时需要约60天。卫星每60天进行一次机动以防止卫星指向太阳，这就导致在正午无卫星数据覆盖，卫星的纬度覆盖范围从一个半球52°到另一个半球83°每60天交替变换一次。本文采用TIMED／SABER数据的时间范围为2002年1月至2018年12月。卫星数据最新版本为V2.0版本，相比V1.07版本，改进了反演算法，提升了产品精度。本文采用的数据版本为V2.0版本。在中间层及以下高度，SABER的最大观测误差在±2 K［15］。与激光雷达观测结果相比，SABER温度与激光雷达温度的偏差在3 K以内。与COSMIC掩星数据相比，在38～60 km，最大温度偏差约为5 K［16］。SABER卫星数据以其较高的数据精度，在科学研究和工程技术等领域已广泛应用［1，17-20］。

采用与文献［13，21］相同的处理方法，将卫星数据通过信息范围检查、极值检查、垂直一致性检查和统计学检查等质量控制，对错误信息进行订正或野值剔除。利用式（1）和式（2），按月份统计计算每个网格内大气密度的月平均值和标准偏差。网格数据的高度间隔为1 km，纬度间隔为4°，经度间隔为5°。网格内月平均值和标准偏差的计算公式如下：

式中：N为该网格内的数据量；ρi为第i次大气密度观测量。

2 大气密度的变化特征

利用统计获得的大气密度网格化数据，对大气密度的气候平均状态和扰动状态进行分析，获得大气密度的定量变化特征，结果如图1和图2所示。

图1 1月份和7月份大气密度月平均值相对于年平均值的变化量在120°E随纬度和高度的分布Fig.1 Latitude-altitude distribution of variation of atmospheric density relative to annual average density on 120°E in January and July

图1给出了1月份和7月份120°E大气密度多年月平均统计结果相对于年平均密度的变化量，可以看出大气密度在南北半球不同的分布特征。在1月份90 km以下，大气密度呈现出从南半球向北半球逐渐递减的趋势。在南半球，大气密度在高纬地区约82°S 80 km高度处存在大气密度的极大值，最大值高于年平均值的74.8%。在北半球，大气密度在中高纬地区80 km附近存在大气密度的极小值，大气密度低于年平均值的36.9%。在90 km以上，大气密度呈现出从低纬地区向两极地区递减的趋势。7月份的变化规律与1月份的相反。在90 km以下，大气密度从北半球高纬地区逐渐向南半球高纬地区递减。在北半球高纬地区约80 km高度存在大气密度的极大值，最大值高于年平均值的77.9%，在南半球高纬地区存在大气密度的极小值，最小值低于年平均值的54.3%。在90 km以上，大气密度也呈现出从低纬地区向高纬地区递减的趋势。

图2 1月份和7月份大气密度标准偏差相对于年平均值的变化量在120°E随纬度和高度的分布Fig.2 Latitude-altitude distribution of variation of atmospheric density standard deviations relative to annual average density on 120°E in January and July

图2给出了1月份和7月份120°E大气密度多年月平均扰动量相对于年平均密度值的变化。可以看出，在1月份，大气密度的相对扰动量在80 km以下，北半球中高纬地区大气密度的扰动量较大，尤其在高纬地区，大气密度相对扰动量在平流层顶附近（～45 km）存在极大值，最大值为17.9%。在低纬地区和南半球中高纬地区，60 km以下，大气密度相对扰动量较小，基本在4%以内。在60～80 km，大气密度相对扰动量随高度逐渐增大，因为大气潮汐波的贡献逐渐显著。在80 km以上，从南半球到北半球大气密度的扰动均较强，相对扰动的最大值可达19.9%，在相同高度层，南半球的大气密度扰动量高于北半球。在7月份，大气密度扰动量随纬度和高度的变化规律与1月份相反。在80 km以下，南半球中高纬地区大气密度扰动较大，最大值为15.0%。在20～60 km，低纬地区和北半球中高纬地区，大气密度扰动较小。在60～80 km，大气密度扰动量随高度逐渐增强。在80 km以上，大气密度扰动在南北半球均较强。在同一高度层，北半球的大气密度扰动量高于南半球。从1月份和7月份大气密度扰动的变化规律可以看出，在80 km以下，冬季半球的中高纬地区大气密度扰动较大。在80 km以上，夏季半球的大气密度扰动量高于冬季半球。

美国标准大气模型（USSA76）是飞行器设计和仿真常用的模型，USSA76模型与观测数据之间的差异是实际应用中非常关心的问题。利用大气密度网格化数据，对USSA76模型相对SABER大气密度随纬度和高度的偏差进行了评估，结果如图3所示。在1月份，90 km以下，从北半球到南半球USSA76大气密度的相对偏差从正偏差逐渐向负偏差过渡。在北半球，高纬地区80 km附近存在正相对偏差的极大值，约为76.7%。在南半球，高纬地区80 km附近存在负相对偏差的极大值，约为38.5%。在90 km以上，USSA76大气密度的相对偏差均为正偏差，且高纬地区的偏差高于低纬地区。在7月份90 km以下，从南半球到北半球，USSA76大气密度的相对偏差从正偏差逐渐向负偏差转变。在南半球高纬地区存在正相对偏差的极大值，为132.2%，在北半球高纬地区80 km 附近存在负相对偏差的极大值，为36.1%。在90 km以上，USSA76大气密度相对偏差为正偏差。

从以上结果可以看出，在90 km以下，大气密度的相对偏差逐渐从冬季半球正偏差向夏季半球负偏差过渡，且冬季半球正偏差的最大值高于夏季半球负偏差的最大值。在90 km以上为正偏差。

图3 1月份和7月份USSA76相对于TIMED／SABER大气密度的偏差在120°E随纬度和高度的分布Fig.3 Latitude-altitude distribution of relative error between USSA76 and TIMED／SABER atmospheric density on 120°E in January and July

3 大气密度建模方法

将临近空间大气密度表征为气候平均量与大气扰动量之和：

对给定轨迹而言，沿轨迹的气候平均量和扰动量分别由网格化的月平均值和标准偏差通过对数插值获得。根据空间尺度，将大气扰动划分为大尺度扰动和小尺度扰动。大尺度扰动包括大气潮汐波和大气行星波，具有行星尺度的波长。小尺度扰动包含大气重力波和湍流，重力波的波长由几百到几十公里，湍流的空间尺度更小。大气密度扰动量表征为大尺度扰动量ρ′L与小尺度扰动量ρ′S之和：

式中：r为ρ′S（x2）／σS（x2）和ρ′S（x1）／σS（x1）之间的自相关，相关系数随着时间和空间尺度呈指数衰减；δh和δz分别为相邻轨迹点之间位移的水平分量和垂直分量；q为高斯分布随机数；Lh和Lz分别为与小尺度扰动有关的水平和垂直尺度，Lh的取值范围为7.8～840 km，Lz的取值范围为0.53～24.03 km；τ为与背景风速有关的时间尺度［22］。

大尺度扰动量是由于大气潮汐和行星波引起的扰动，具有一定的周期性。因此，将大尺度扰动量表征为余弦函数的形式：

式中：ρ′L为大尺度扰动量；σL为大尺度扰动的标准偏差；A为振幅因子；ω＝2π／T为随机波动周期为T的角速度；m和n分别为纬向波数和经向波数；φ和θ分别为经度和纬度；z为高度；λ为与高度有关的垂直波长；β为0～2π均匀分布的随机相位。

模型的初始扰动值可以是初始位置实际观测数据计算得到的扰动值，或是从全球网格化标准偏差数据中利用插值方法得到的扰动值。

4 模型测试

为检验大气模型的仿真效果，以2019年中国科学院国家空间科学中心敦煌临近空间激光雷达观测试验的数据为基准，利用建立的大气模型仿真敦煌上空20～100 km的大气密度，并将观测结果与仿真结果进行对比。假设给定轨迹仅随高度变化，不随时间和经纬度变化，初始高度为20 km，终点高度为100 km。为方便与激光雷达观测数据对比，将仿真时刻均设置为2019年4月18日21：30，仿真轨迹的高度间隔为1 km，模型仿真结果如图4所示。图中：ρmean为密度月均值。左图给出的是20～100 km高度大气密度100次蒙特卡罗仿真值相对于沿轨迹的大气密度月平均值的变化量，右图给出的是沿轨迹的大气密度月平均值。可以看出，大气密度的仿真结果沿月平均值基本呈对称分布，仿真值基本在3倍标准偏差范围内。仿真的扰动量随高度逐渐增大，仿真结果展现的规律与现有的研究结果相符。这些扰动主要由不同来源的静态行星波、传播性行星波、潮汐波、重力波和湍流等贡献。此外，较高高度上波流相互作用、波与波之间的非线性相互作用等导致了大气密度偏离气候平均态产生复杂的变化，展现了大气的整体扰动水平。

图4 100次蒙特卡罗仿真的大气密度相对于月平均值的变化量随高度的变化Fig.4 Variation of atmospheric density of 100 Monte Carlo simulations relative to monthly mean density with height

图5给出了大气模型在2019年4月18日21：30激光雷达观测结果与单次仿真结果的对比。图中：ρ／ρUS76为大气密度与标准大气的比值。蓝色实线表示模型仿真值，红色实线表示激光雷达观测值，观测的高度范围为38～80 km。可以看出，在激光雷达观测高度范围内，仿真值与观测值吻合程度较好，两者随高度变化具有相同的变化趋势。单次仿真的结果存在偶然性，因此，对100次蒙特卡罗仿真结果进行统计分析。

将100次蒙特卡罗仿真结果与激光雷达观测结果进行相关系数计算，定量地展现模型仿真值与激光雷达观测值之间的吻合程度。图6给出了100次蒙特卡罗仿真值与激光雷达观测值相关性检验的结果。可以看出，100次仿真结果与激光雷达观测结果之间存在显著的线性相关，回归斜率在0.987 33～1.105 9之间，表明每次的仿真结果均与观测结果在高度上具有相同的变化趋势。

图5 一次仿真结果与激光雷达观测值的比较Fig.5 Comparison between one-time simulation result and lidar observation values

图6 100次蒙特卡罗仿真试验大气密度相关性检验Fig.6 Atmospheric density correlation test of 100 Monte Carlo simulation

表1 100次蒙特卡罗仿真结果与激光雷达观测结果之间的平均相对偏差Tab1e 1 Mean re1ative error between 100 Monte Car1o simu1ation and 1idar obser vation resu1ts

对100次蒙特卡罗仿真结果与激光雷达观测结果之间相对偏差进行统计，结果如表1所示。在65 km以下，模型仿真结果为负偏差，模型仿真值相对激光雷达观测值偏小，在65 km处存在负偏差的极大值，相对偏差可达－5.09%。在70 km和75 km模型的相对偏差较小，在±1%之间。在80 km，模型偏差较大，为7.18%。后续需要采用其他观测数据，进一步对40 km以下和80 km以上模型的可靠性进行评估。

5 模型应用

假设飞行器某次再入飞行的轨迹如图7所示，飞行器再入的初始时间为2019年11月24日10：00，初始再入的高度为100 km，经纬度为（116.0°E，40.0°N），终端高度为20 km，经纬度为（90.7°E，36.7°N），轨迹点之间的时间步长为1 s。初始位置大气密度的平均值由统计获得的网格化月平均值经对数插值计算获得。

图8再现了沿飞行轨迹大气密度的可能状态。可以看出，大气密度的仿真值均在3倍的标准偏差范围内。因为随高度大气密度呈指数减小，在飞行器每次跳跃运动的最高点时对应大气密度仿真结果的最低点，在图8中能够清晰地展现飞行器在飞行的不同阶段大气密度的变化特征。

图7 假设的飞行器再入轨迹Fig.7 An assumed trajectory of a reentry vehicle

图9分别给出了飞行轨迹上大气密度不同尺度的扰动量和总扰动量。可以看出，大气密度大、小尺度扰动量和总扰动量的仿真结果基本呈对称分布，且均未超出3倍标准偏差范围，结果表明，大气模型能够很好地表征飞行轨迹上不同尺度大气密度的扰动量。

图8 沿轨迹的大气密度仿真结果Fig.8 Simulation result of atmospheric density along trajectory

图9 沿轨迹大气密度小尺度扰动量、大尺度扰动量和总扰动量的仿真结果Fig.9 Simulation result of small-scale perturbances，large-scale perturbances and total perturbances of atmospheric density along trajectory

图10 大尺度扰动和小尺度扰动对总扰动的贡献Fig.10 Contribution of large-scale and small-scale perturbances to total perturbances

为分析仿真结果中大尺度扰动和小尺度扰动对总大气扰动的贡献程度，分别统计计算了大尺度扰动和小尺度扰动占总扰动的比例，如图10所示。在60 s以内，飞行高度在90 km以上，小尺度扰动占主导地位，60 s以后，随着高度的降低，总扰动逐渐以大尺度扰动为主导。这一规律与文献［21］的研究结果具有较好的一致性，符合大气的物理规律。在中间层顶和低热层高度，小尺度扰动的贡献占主导地位，主要贡献来自湍流、重力波及其他随机性的小尺度扰动。在70 km以上，大气潮汐波的贡献逐渐显著，在70 km以下，大气扰动的贡献主要来自于行星波和重力波。

6 结 论

基于2002—2018年TIMED／SABER卫星观测的大气密度数据，统计获得大气密度多年月平均值和标准偏差。基于网格化的密度数据，对大气密度的变化规律进行了分析，并对大气密度建模方法进行探索。根据大气波动的特征，将波动特征显著的大尺度波动（行星波和潮汐波）以三角函数的形式进行表征，将更具随机性的小尺度扰动（重力波和湍流）用一阶自回归方法进行表征，初步建立了全球临近空间大气密度模型。得到以下结论：

1）通过分析大气密度在不同纬度的垂直分布，结果表明，在90 km以下，大气密度从夏季半球向冬季半球递减；在90 km以上，大气密度从低纬向高纬递减。在80 km以下，冬季半球的中高纬地区大气密度扰动较大。在80 km以上，夏季半球的大气密度扰动量高于冬季半球。

2）通过计算USSA76与SABER大气密度之间的相对偏差，结果表明，在90 km以下，大气密度的相对偏差逐渐从冬季半球正偏差向夏季半球负偏差过渡，且冬季半球正偏差的最大值高于夏季半球负偏差的最大值，最大值可达132.2%。在90 km以上为正偏差。

3）通过对敦煌上空20～100 km大气密度进行100次蒙特卡罗仿真，并与激光雷达观测结果进行对比，结果表明，本文的建模方法是可行的，观测值与仿真值具有较好的吻合度。

4）通过采用蒙特卡罗方法对给定再入轨迹上的大气密度进行仿真，结果表明，模型可以再现飞行轨迹上大气密度可能的状态。

建立的大气密度模型主要优点在于：数据来源采用最近17年的观测数据，数据量大且数据全球覆盖性较好，有利于建立全球模型；建模方法更符合物理规律，不仅可以给出飞行轨迹上大气密度的平均态和总的扰动状态，还可以分别给出大尺度扰动量和小尺度扰动量。该模型可以作为初步工具，为临近空间飞行器气动设计、控制设计和飞行仿真等提供数据输入。后续将在建模方法中以温度、密度和气压满足的物理规律为约束，考虑不同尺度不同参量之间满足的相关性，对多物理参量进行协同建模，提升模型的可靠性和实用性。

致谢 感谢TIMED／SABER工作组提供的探测数据（http：∥saber.gats-inc.com／data＿services.php）。