大气密度模型:Jacchia-Robert

本文介绍轨道力学中常用的大气密度模型 Jacchia-Robert，主要用于近地卫星大气阻力计算。STK 高精度轨道积分器（HPOP）中大气密度模型默认即为 Jacchia-Robert；相关参数包括阻力系数 \(C_D\)、面质比 \(A/m\)、日太阳辐射通量 \(F10.7\)、81 天平均通量 \(\bar{F}10.7\) 以及地磁指数 \(K_p\)。注意：本文不讨论 125 km 以下的大气密度（属飞机、高超滑行飞行器范畴）。

STK中的大气阻力设置

大气阻力公式

在近地轨道（300 km–2000 km），地球大气的阻力不容忽视，尤其对于 500 km 高度以下的卫星，阻力作用明显。大气阻力对卫星产生一个阻力的作用，阻力公式为：

\[ \vec{F}_{drag}=-\frac{1}{2}\rho\left(\frac{C_DA}{m}\right)V_r\vec{V}_r \]

其中：

\(\rho\)：卫星所在位置处的大气密度
\(A/m\)：即卫星的面质比系数，即卫星的迎风面积和质量的比值
\(C_D\)：大气阻力系数。（一般取 2.2）
\(\vec{V}_r\)：卫星在地心惯性系中相对大气的速度

对于大气密度而言，虽然上面仅仅是一个参数，但是却与许多参数有关。

大气密度概述

最简单，也最容易计算的大气密度模型为静止大气模型，即考虑地球为一圆球体，大气密度随着高度呈指数衰减。大气密度仅仅与高度相关，与地固系下具体位置无关，也与太阳辐射和地磁通量无关。

显然，静止大气模型容易计算，但是却与实际情形不符合。对于近地轨道，地球大气密度与太阳辐射和地磁通量密切相关，从而导致地球上空每个地方的密度都不相同，且随着时间的变化而变化，见下图。

大气密度的变化

简单来说，大气密度与大气温度密切相关，而大气温度主要受到太阳辐射和地磁活动影响（称之为空间环境扰动），从而造成地球上空每个地方的温度都不一样！！

空间环境扰动可以引起大气密度的剧烈变化，若大气密度在短时间内快速上升，卫星受到的大气阻力也会突然增加，从而加快卫星轨道的衰减。比如，美国“哥伦比亚”号航天飞机在 1981 年 4 月 12 日飞行时，遇到一次剧烈的空间环境扰动事件，陡增的大气密度导致该航天飞机下降到较低轨道的时间比预期快了 60%。

再比如 1989 年 3 月份的“卡林顿事件”中，空间环境的剧烈扰动引起大气密度急剧增加，使 840 km 高度的大气密度增加了 9 倍。大气密度的剧增引起了美国的太阳峰年卫星（SMM）在整个事件期间的运行高度下降了 5 km，从而提前陨落。

卡林顿事件期间，空间环境出现剧烈扰动，大气密度急剧增加

假定 T 时刻，卫星在地固系下的位置为 \(\vec{r}_{sat}\)，那么卫星所在位置的大气密度主要受到以下几点影响：

太阳辐射对大气的温度影响（\(F10.7\) 参数和 \(\bar{F}10.7\) 参数）
地磁辐射对大气的温度影响（\(K_p\) 指数）
T 时刻太阳的位置（\(\vec{r}_{sun}\)）

因此，计算给定位置 \(\vec{r}_{sat}\) 处的大气密度时，需要的参数为：\(F10.7\)、\(\bar{F}10.7\)、\(K_p\)，以及具体计算大气密度的模型。

下面首先介绍太阳和地磁辐射参数，然后详细介绍 Jacchia-Robert 的大气密度模型。

太阳辐射通量(F10.7）

太阳活动是影响大气的主要因素，常见的太阳活动指数包括太阳黑子数 SSN 及 10.7 厘米射电辐射通量 F10.7。

太阳黑子，是人们最早观测到的太阳活动现象。1843 年，德国天文爱好者施瓦布通过日常观测发现了太阳黑子数量的多少存在 11 年左右的周期。之后，随着观测数据的增加，这一规律不断被证实，并且人们发现黑子数的多少与这个时期的太阳活跃程度相对应。于是，太阳黑子数的这种规律变化成为人们划分太阳活动周期的标志，黑子数量的高峰年称为太阳活动峰年，黑子数最少年称为太阳活动低年，两次低年之间定为一个太阳活动周。

除了太阳黑子数之外，人们还发现了另一种能代表太阳活动周变化的参量—10.7 cm 射电流量（F10.7）。

太阳向宇宙空间发射（太阳高层大气的辐射）的电磁波中，波长从毫米到十米不等，其中，10.7 cm 波长 (2800 MHz) 属于极紫外辐射。地球高层大气吸收太阳极紫外辐射，吸收能量的 20%–30% 用来加热高层大气。

定义 10.7 cm 波长的射电辐射通量的大小为太阳辐射通量指数，也称 F107 指数（常用 F10.7 表示）。F10.7 的单位为“太阳通量单位”(SFU)，具体数值为：

\[ 1\,\mathrm{SFU}=10^{-22}\,\mathrm{W/m^2/Hz} \]

从长期的监测中人们发现，F10.7 和太阳黑子数有很强的相关性，F10.7 值的大小也能很好地代表太阳活动的强弱，并且由于 F10.7 在地面就可以监测获取，长久以来在许多重要的电离层和中高层大气模型中，通常都是以 F10.7 作为输入来表征太阳活动的水平。因此，无论是过去、现在，还是未来，F10.7 监测在太阳活动预报和研究中都将具有举足轻重的地位。

F10.7 指数的范围在 60 到 300 之间，中国气象标准规定了太阳活动水平的分级，按 F107 指数对太阳活动水平进行划分如下：

等级	F10.7数值
很低	\(< 80\)
低	\(80\)–\(100\)
中	\(100\)–\(150\)
高	\(150\)–\(200\)
很高	\(> 200\)

F107 指数分为两种，一种是每天的观测值 \(F10.7\)，另一种是 81 天（3 个太阳自转周期）的平均值 \(\bar{F}10.7\)，见下图。

太阳通量指数F107

地磁指数（kp或Ap）

地磁扰动会引起大气离子化，从而引起大气密度的变化。

全球范围内有 12 个站实时监测地磁的变化，并以相对安静时的磁场作为参考，对地磁观测的三个分量中在水平两分量扰动的平均值，每三小时给出地磁扰动强度的指数，称为三小时磁情指数 \(a_p\)（3-hourly index）。\(a_p\) 的单位为 2 纳特（nT)。

\[ 1\,\mathrm{nT} = 10^{-9}\,\mathrm{Tesla} \]

把一天分为 8 个时间段，每段三小时，每段都有一个指数 \(a_p\)，一天 8 个数值，取其平均值则得到 \(A_p\)，称为日行星振幅 (daily planetary amplitude)。

\(a_p\) 或 \(A_p\) 数值范围为 0–400，数值越大，地磁扰动越强，一般大于 100 就很少见，常见的数值在 10–20 左右。

与 \(a_p\) 数值相对应的是行星指数 \(K_p\)（planetary index），它与 \(a_p\) 为准对数关系，且一一对应，因此 \(K_p\) 也是 3 小时一个数值，数值范围 0–9，越大表示地磁活动越强。

以 \(K_p\) 的数值对地磁暴的等级进行分类如下：

磁暴等级	\(K_p\)数值
平静	\(0\)–\(2\)
扰动	\(2\)–\(3\)
活跃	\(3\)–\(4\)
小磁暴	\(4\)–\(5\)
大磁暴	\(5\)–\(6.5\)
严重磁暴	\(6.5\)–\(9\)

\(a_p\) 和 \(K_p\) 的关系如下表。

地磁指数Kp和ap的关系

两者对应的关系如下图，以对数方式画出。

在这里插入图片描述

地磁指数 \(a_p\) 和 \(K_p\) 的相互转换可以使用 3 次多项式进行插值，此处不再详述。

太阳和地磁指数的数据获取与计算

太阳通量指数 \(F10.7\) 和 \(\bar{F}10.7\)，以及地磁指数 \(a_p\)（\(A_p\)）和 \(K_p\) 在网上可以查询到每天的数据，本文给出几个来源：

其中，Celestrak 网站上给出的数据包含了上述几个指数，包括 3 小时指数和日平均指数，数据格式见下图。

Celestrak网站上的太空天气

计算时，需要根据当前时刻 T 根据给定的数据得到 \(F10.7\)、\(\bar{F}10.7\)、\(K_p\) 等参数。

也可以简化，取常数：

\[ \begin{aligned} F10.7&=\bar{F}10.7=150.0\\ K_p&=3.0 \end{aligned} \]

STK 默认的大气密度设置中，上述三个参数就是这几个常数。

Jacchia-Robert大气密度模型

Jacchia-Roberts 大气密度模型是一种用于描述大气层密度分布的数学模型。该模型由意大利科学家 J. M. Jacchia 和美国科学家 L. G. Roberts 在 20 世纪 60 年代开发。

Jacchia-Roberts 模型旨在提供大气层密度的近似值，以便用于空间科学和工程应用中，例如卫星轨道预测、再入飞行器设计等。该模型基于大量的大气观测数据，并利用统计方法来拟合这些观测数据，以获得对大气密度的预测。

Jacchia-Roberts 模型考虑了多种影响大气密度的因素，包括太阳活动、地球磁场、大气成分等。它将大气层分为多个不同的层次，并使用复杂的数学公式来描述每个层次的密度分布。模型的输入参数包括高度、纬度、经度、时间等。

需要注意的是，Jacchia-Roberts 模型是一种经验模型，它的准确性可能会受到多种因素的影响，例如地理位置、季节变化、太阳活动等。因此，在具体的应用中，可能需要根据实际情况进行调整或采用其他更精确的模型。

总而言之，Jacchia-Roberts 大气密度模型是一种常用的近似模型，用于估计大气层密度分布，特别适用于空间科学和工程领域的应用。

大气密度模型计算的思路为：外层大气温度受太阳辐射和地磁活动的影响并且存在周日变化，首先求出外层大气温度，再利用静态模式求出各大气成分的数密度，对求出的数密度加以地磁、季节纬度项和半年项改正，并最终求出大气密度值，具体过程如下：

基本参数

令 T 时刻，地固系下，卫星的位置 \(\vec{r}_{sat}=(r_I,r_J,r_K)^T\)，太阳的单位位置矢量 \(\vec{r}_{sun}=(r_x,r_y,r_z)^T\)。

那么卫星的大地纬度（近似公式，没有迭代求解）为：

\[ \phi_{gd}=\tan^{-1}\left(\frac{1}{(1-f)^2}\left[\frac{r_K}{\sqrt{r^{2}_{I}+r^{2}_J}}\right]\right) \]

太阳的纬度为：

\[ \delta_s=\tan^{-1}\left(\frac{r_z}{\sqrt{r^{2}_{x}+r^{2}_y}}\right) \]

太阳的相对时角 \(LHA_s\)（此处单位为度）为：

\[ LHA_s=\frac{180^\circ}{\pi}\left(\frac{r_xr_J-r_yr_I}{\left|r_xr_J-r_yr_I\right|}\cos^{-1}\left(\frac{r_xr_I+r_yr_J}{\sqrt{r^{2}_x+r^{2}_y}\sqrt{r^{2}_I+r^{2}_J}}\right)\right) \]

卫星的大地高度 \(h\) 的计算，需要考虑地球椭球体，具体计算公式可以为精确的，也可以近似计算，此处不再给出，注意，在后面的计算中，大地高度 \(h\) 的单位为 km。

后面需要用到的其它常数为：

\[ \begin{aligned} T_0&=183^\circ \\ R_{pole}&=6356.766\,\mathrm{km}\\ g_0 &=9.80665 \\ R &= 8.31432\\ A &=6.022045\times10^{23}\\ \epsilon&=23.439291^\circ=0.4090928 \end{aligned} \]

大气温度

夜间全球外层温度 \(T_c\)：

\[ T_c(K) = 379+3.24\bar{F}_{10.7}+1.3\left[F_{10.7}-\bar{F}_{10.7}\right] \]

由此得到未修正的外层大气温度 \(T_{unc}\)：

\[ T_{unc}=T_c\left(1+0.3\left[\sin^{2.2}(\theta)+\left(\cos^{2.2}(\eta)-\sin^{2.2}(\theta)\right)\cos^3\left(\frac{\tau}{2}\right)\right]\right) \]

上式中：

\[ \begin{aligned} \eta&=\frac{\left|\phi_{gd}-\delta_s\right|}{2}\\ \theta&=\frac{\left|\phi_{gd}+\delta_s\right|}{2}\\ \tau&=LHA_s-37^\circ+6.0^\circ\sin(LHA_s+43.0^\circ) \end{aligned} \]

注意，在计算 \(\tau\) 时，三角函数为度，需要转换为弧度，另外得到 \(\tau\) 后也需要转换为弧度，才能带入 \(T_{unc}\) 中计算。

考虑地磁指数 \(k_p\) 对温度的修正：

\[ \Delta T_{corr}=\begin{cases} 28.0^\circ k_p+0.03e^{k_p}, & h\geq200\,\mathrm{km} \\ 14.0^\circ k_p+0.02e^{k_p}, & h<200\,\mathrm{km} \end{cases} \]

从而得到修正的外层温度 \(T_{corr}\)：

\[ T_{corr}=T_{unc}+\Delta T_{corr} \tag{1} \]

临界点温度 \(T_x\) 为：

\[ T_x=371.6678^\circ+0.0518806T_{corr}-294.3505^\circ e^{-0.00216222T_{corr}} \tag{2} \]

令 \(T_0=183^\circ\)。则 125 km 高度以下的温度为：

\[ T(h)_{0-125}=T_x+\frac{T_x-T_0}{35^4}\sum_{n=0}^4 C_n h^n \tag{3} \]

上式中的系数为：

\[ \begin{aligned} C_0&=-89284375.0;\quad C_1=3542400.0;\quad C_2=-52687.5\\ C_3&=340.5;\quad C_4=-0.8 \end{aligned} \]

125 km 高度以上的温度为（Robert 修正）：

\[ T(h)_{125-}=T_{corr}-(T_{corr}-T_x)e^{-\left(\frac{T_x-T_0}{T_{corr}-T_x}\right)\left(\frac{h-125}{35}\right)\left(\frac{l}{R_{pole}+h}\right)} \tag{4} \]

上式中，系数 \(l\) 为多项式的和：

\[ l=\sum_{j=0}^4 l_j T^{j}_{corr} \tag{5} \]

系数为：

\[ \begin{aligned} l_0&=0.1031445\times10^5\\ l_1&=0.2341230\times10^1\\ l_2&=0.1579202\times10^{-2}\\ l_3&=-0.1252487\times10^{-5}\\ l_4&=0.2462708\times10^{-9} \end{aligned} \]

大气密度

首先考虑几项对大气密度的修正。

地磁修正

高度 200 km 以下（\(h<200\,\mathrm{km}\)）时，考虑地磁的修正为：

\[ (\Delta \log_{10}\rho)_G=0.012k_p+1.2\times10^{-5}e^{k_p} \]

季节纬度变化

首先计算自 1958 年 1 月 1 日零时到当前时刻的累计天数（假设都使用 TAI 时间计算），并可得到累计的世纪数：

\[ T_{1958}=\frac{JD_T-JD_{1958}}{365.2422} \]

\(JD_{1958}\) 对应的儒略日为 2436204.5。

则季节纬度变化的修正为：

\[ (\Delta \log_{10}\rho)_{LT}=0.014(h-90)\sin(2\pi T_{1958}+1.72)\sin(\phi_{gd})\left|\sin(\phi_{gd})\right|e^{-0.0013(h-90)^2} \]

半年变化

令

\[ \tau_{SA}=T_{1958}+0.09544\left(\left[0.5+0.5\sin(2\pi T_{1958}+6.035)\right]^{1.65}-0.5\right) \]

则半年变化的修正为：

\[ \begin{aligned} (\Delta \log_{10}\rho)_{SA}&=(5.876\times10^{-7}h^{2.331}+0.06328)e^{-0.002868h}\\ &\quad\times\left(0.02835+\left[0.3817+0.17829\sin(2\pi \tau_{SA}+4.137)\right]\sin(4\pi \tau_{SA}+4.259)\right) \end{aligned} \]

密度公式

上述三项的修正总和为：

\[ (\Delta \log_{10}\rho)_{corr}=(\Delta \log_{10}\rho)_G+(\Delta \log_{10}\rho)_{LT}+(\Delta \log_{10}\rho)_{SA} \]

Robert (1971) 注意到大气密度在 90 km–125 km 变化比较剧烈，因此他将大气密度分为三个段，不同高度对应的大气密度（单位：\(kg/m^3\)）为：

\[ \rho (h)= \begin{cases} \rho (h)_{90-100}\times1000\times10^{(\Delta \log_{10}\rho)_{corr}}, & 90\,\mathrm{km}<h\leq100\,\mathrm{km}\\ \rho (h)_{100-125}\times1000\times10^{(\Delta \log_{10}\rho)_{corr}}, & 100\,\mathrm{km}< h\leq125\,\mathrm{km}\\ \rho (h)_{125-}\times1000\times10^{(\Delta \log_{10}\rho)_{corr}}, & 125\,\mathrm{km}<h \end{cases} \]

上式中，乘以 1000 是因为原来的公式中大气密度的单位为 \(g/cm^3\)。

由于卫星轨道一般不可能低于 125 km，所以本文只给出 125 km 以上的大气密度，即上式中的 \(\rho (h)_{125-}\)。

125 km 以上密度 ρ(h)₁₂₅₋

大气密度是由 6 种气体分子的密度单独计算的，六种气体分子的相关参数见下表。

大气组成的分子密度

后面公式中，需要使用到上表中的参数 \(M_i\) 和 \(a_i\)。

首先给出各分子在 125 km 高度处的密度 \(\rho_i(125)\)（注意，这里 \(i\) 取值为 1–5，氢气 H 分子暂未包含，后面单独计算）。

\[ \rho_i(125)=\frac{M_i}{A}10^{\sum^6_{j=0} \delta_{ij}T^{j}_{corr}}\tag{6} \]

系数 \(\delta_{ij}\) 由下表给出。

大气各分子密度的拟合系数

并令各分子的 \(\gamma_i\) 系数为：

\[ \gamma_i=\frac{M_i g_0 R^{2}_{pole}}{R l T_{corr}}\left(\frac{T_{corr}-T_x}{T_x-T_0}\right)\left(\frac{35}{6481.766}\right) \]

其中，\(l\) 的取值参见公式 (5)，其余常数参见“基本参数”章节。

则最终的大气密度由 6 种气体分子的密度相加而成，其中前 5 种的和为：

\[ \rho_{1-5} (h)_{125-}=\sum_{i=1}^5\rho_i(125)\left(\frac{T_x}{T(h)}\right)^{1+a_i +\gamma_i}\left(\frac{T_{corr}-T(h)}{T_{corr}-T_x}\right)^{\gamma_i} \tag{7} \]

上式中，对于第 3 种气体 He 的密度 \(\rho_3(125)\) 有修正，修正公式如下：

\[ (\Delta \log_{10}\rho)_{He}=0.65\left|\frac{\delta_s}{\epsilon}\right|\left[\sin^3\left(\frac{\pi}{4}-\frac{\phi_{gd}\delta_s}{2\left|\delta_s\right|}\right)-0.35355\right] \]

则修正后的 \(\rho_3(125)\) 为：

\[ \rho_3(125)=\rho_3(125)\times10^{(\Delta \log_{10}\rho)_{He}} \]

第 6 种气体分子 H 密度 \(\rho_6(h)\) 为（仅在 500 km 以上高度才考虑）：

\[ \rho_6 (h)_{500-}=\rho_H(500)\left(\frac{T_{500}}{T(h)}\right)^{1+\gamma_H}\left(\frac{T_{corr}-T(h)}{T_{corr}-T_{500}}\right)^{\gamma_H} \tag{7} \]

上式中，\(T_{500}\) 表示在 500 km 高度的温度，有 (4) 式给出，另外 \(\rho_H(500)\) 为 500 km 高度的 H 分子的密度，由下式给出：

\[ \rho_H(500)=\frac{M_H}{A}10^{\left[73.13-(39.4-5.5\log_{10}T_{500})(\log_{10}T_{500})\right]} \]

因此，最终 6 种气体分子的密度总和为：

\[ \rho (h)_{125-}=\rho_{1-5} (h)_{125-}+\rho_6 (h)_{500-} \tag{8} \]

因此，最终 Jacchia-Robert 的大气密度模型由上式给出（本文仅给出 125 km 高度以上的公式），注意，当高度大于 2500 km 以上时，可直接令大气密度为 0，即 2500 km 以上高度认为无大气。

算例

时间：2017-01-01 00:00:00 UTC
大地纬度：45 deg
大地经度：0 deg
\(F10.7\): 100
\(\bar{F}\): 100
\(K_p\): 4

则不同大地高度的大气密度如下：

大地高度(km)	大气密度(kg/m³)
125.1	\(1.5899\times10^{-8}\)
300	\(1.3061\times10^{-11}\)
700	\(1.3480\times10^{-14}\)
1500	\(4.0058\times10^{-16}\)

参考

卫星在轨寿命浅谈
空间环境预报的基石——太阳活动监测
地磁活动指数与太阳活动指数
Julia源码: SatelliteToolboxAtmosphericModels.jl
"Fundamentals of Astrodynamics and Applications", 3rd Edition, David A. Vallado; Appendix B
Geosat卫星定轨中的大气阻力摄动
NOAA Kp和ap指数下载