电离层受到一系列外驱动的影响，包括来自太阳的极紫外辐射、行星际条件、极区电场以及低层的重力波、潮汐、和行星波等，并通过化学和动力学过程与背景中性大气紧密耦合在一起.中性大气参量(温度、成分、风场等)各种尺度的变化，也会在电离层参量里有所体现.要更好地研究电离层日常变化，尤其是进行预报，需要正确描述热层状态.然而，截止到目前热层观测依然十分匮乏，给准确预报电离层带来了困难.

近些年来，众多研究已证实，基于热层电离层耦合数据同化系统，通过同化大量电离层电子密度相关观测，优化热层参量，比如中性质量密度、中性成分以及中性风场等，可以有效提升电离层预报时长和精度(Anderson et al., 2009; Chartier et al., 2016; Chen et al., 2017; Hsu et al., 2018; Lee et al., 2012; Matsuo et al., 2013；Wu et al., 2002).如Hsu等(2014)利用Data Assimilation Research Testbed (DART)同化F3/C(Formosa Satellite 3/Constellation Observing System for Meteorology, Ionosphere, and Climate)电子密度剖面观测数据，表明更新热层成分对电离层预报时长的改善是至关重要的，其改善效果能达到大于24 h，明显优于只更新电离层参量的2~3 h时长.Matsuo等(2013)同样利用DART，证明了同化F3/C电子密度剖面数据能够给热层质量密度带来全球的改善.Chen等(2017)利用DART同化全球2127个地基GPS垂直TEC观测，揭示了纬向中性风场能够在热层电离层耦合同化中被很好地优化，进而能够模拟更符合实际观测的电离层电场日落反转增强现象.

我们将通过同化电离层斜TEC观测，优化中性质量密度以及积分氧氮比，进而提高系统预报能力，并基于全球空地基GNSS、CHAMP及TIMED/GUVI观测构型设计了系列观测系统模拟实验，对系统性能进行了评估.对中性质量密度预报进行了探讨，并对改善预报能力的主导因素进行了分析.

下面从理论模型和同化算法两方面介绍我们开发的数据同化系统.

我们选择NCAR-TIEGCM作为背景模型，因为该理论模式充分考虑了电离层热层物理，代码经过了众多研究检验，具有较高的稳定性和可靠性.关于该模型的详细介绍，可以参考Richmond等(1992).该同化系统中背景模型分辨率是经纬度为5°、高度为半标高、时间为180 s.

集合卡尔曼滤波是一种基于贝叶斯理论的递归数据同化算法，分为更新和预测两个步骤(Evensen，1994; Evensen and Van Leeuwen，1996).更新阶段主要是将观测数据融入模式中，实现当前时刻模型背景值的优化，类似于三维变分算法和卡尔曼滤波算法.该阶段主要包括观测算子计算以及卡尔曼增益的解算.

(1)

(2)

(3)

(4)

公式中x_k^f，y_k表示当前时刻背景值和观测数据，g₁，g₂，…，g_p表示模型网格.P_k^f，ε_k，H_k，K_k，R_k分别表示背景误差协方差、观测误差、观测算子、卡尔曼增益和观测误差协方差.x_k^a代表更新值，也称之为同化后值或分析值.

预测阶段是指模型背景值和误差协方差随时间向前传播，实现下一时刻模型背景状态的估计.不同于卡尔曼滤波和变分同化算法，集合卡尔曼滤波算法采用样本统计的方式对背景误差协方差矩阵进行估计，并假设样本符合正态分布，从而避免对模型的线性近似.

(5)

(6)

在集合卡尔曼滤波数据同化算法中，由于集合体的引入，面临超大矩阵计算和存储，通常依赖超级计算机完成.为了避免大规模计算，我们采用稀疏矩阵算法和线性迭代求解稀疏线性方程组(Saad, 1990)，具体介绍可以参考Yue等(2014).基于该方法，我们能够在一般工作站上开展热层电离层集合卡尔曼滤波数据同化实验(He et al., 2019, 2020).

本研究通过设计观测系统模拟实验来检验同化系统的性能.选择的观测系统包括：图 1a给出了2006年7月12日一天中全球地基IGS TEC观测台站位置和COSMIC掩星观测切点位置；图 1b和1c分别给出了同一天内CHAMP在轨热层中性质量密度观测和TIMED/GUVI积分氧氮比探测位置分布.在观测系统模拟实验中，真实场是NCAR- TIEGCM在太阳活动指数F₁₀₇设置为140 sfu (solar flux unit)、地磁活动指数Kp为1时的模拟结果，并基于上述真实观测构型采样，得到的结果称之为“真实”结果.为了对热层参量优化效果进行评估，我们利用F₁₀₇等于78 sfu以及Kp等于1时驱动模式生成控制场，并利用真实观测构型采样，得到的结果称之为控制结果.同时，我们开展了2006年7月12日00:00 UT至21:00 UT的22个小时集合同化实验以及12日21:00 UT至13日21:00 UT的24个小时预报实验.值得注意的是，在整个同化和预报阶段太阳活动指数F₁₀₇以及地磁活动指数Kp始终为78和1.其他同化参数设置如下：同化时间窗口为1 h；局地化参数在三个空间方向上设置间隔分别是，纬度方向5°、经度方向40°以及高度方向30 km；集合体个数选择为90个；同化更新参量包括电子密度、氧原子离子、中性温度、氧原子和氧气分子百分比(N_e、O⁺、T_N、O、O₂).

图 1

Fig. 1

图 1 (a) 2006年7月12日一天内IGS台站位置(黑色三角形点)和COSMIC掩星观测切点位置(黑色点)分布图；(b) 2006年7月12日一天内CHAMP卫星中性质量密度观测轨道分布图；(c) 2006年7月12日一天内TIMED/GUVI探测的积分氧氮比位置分布图 Fig. 1 (a) The global position distribution of the IGS stations (black triangle dots) and the tangent points of COSMIC occultation (black dots) during the day of 12 July 2006; (b) The position distribution of the neutral mass density observations from CHAMP satellite during the day of 12 July 2006;(c) The position distribution of the Σ[O/N2] observation from TIMED/GUVI during the day of 12 July 2006

通过同化已有的观测参量达到对未观测参量进行优化是集合kalman滤波数据同化算法的精髓所在(Yue et al., 2007)，未观测参量的调整主要是通过包含在背景误差协方差中与已知观测参量之间的相关性来实现.这些被调整的参量在预测阶段经过理论模型的约束，进而实现优化.

中性质量密度(ρ)是关于热层中性温度和中性成分的函数，它是通过下式计算得到的：

(7)

式中，P，T，γo，γo₂分别表示压力(单位：N/m²)，温度(单位：K)，氧原子质量混合比值和氧气分子质量混合比值.R^*是标准气体常数，等于8.314×10³J/K/kmol.为了评估同化电离层TEC观测对这热层参量的优化效果，图 2a给出一个压力面上(pressure level = 20)平均温度值和图 2b给出了全球积分氧氮比均值在整个同化阶段内随时间变化图.

图 2

Fig. 2

图 2 (a) 22个小时集合同化实验内在一个压力面上(pressure level=20)控制实验平均温度(圆形标注线)、“真实”平均温度(三角形标注线)以及被估计的平均温度(正方形标注线)随时间变化图；(b) 22个小时集合同化实验内全球平均的控制实验积分氧氮比值(圆形标注线)、“真实”积分氧氮比值(三角形标注线)以及被估计的平均积分氧氮比值(正方形标注线)随时间变化图 Fig. 2 (a) The time evolution of global averaged neutral temperature at one pressure level (pressure level=20) of the control run, the forecast value of OSSE and the simulated "truth" during the 22 h ensemble assimilation period. The lines with circle, square and triangle show the control, forecast and "truth" results, respectively; (b) The time evolution of global averaged Σ[O/N2] of the control run, the forecast value of OSSE and the simulated "truth" during the 22 h ensemble assimilation period. The lines with circle, square and triangle show the control, forecast and "truth" results, respectively

从图 2a可以看出，在选中的压力面上(pressure level=20)，“真实”温度值大于控制值，这是由于“真实”值相较于控制值是在更大的太阳活动指数驱动下生成的.随着同化的进行，被估计的全球平均温度值逐渐增加，并且更加接近于“真实”值.相较于控制值，在整个集合同化阶段该压力面上被估计的温度值与“真实”值平均偏差为177 K, 这是小于控制结果与“真实”值的平均偏差239 K.同样，根据图 2b，对于被估计的积分氧氮比全球均值，也是随着同化的进行越来越接近于“真实”值.被估计的积分氧氮比的值与“真实”值的平均偏差为0.059，这是明显小于控制结果与“真实”值的平均偏差0.115.这些说明集合卡尔曼滤波的确可以通过同化电子密度优化中性温度和成分.但是，我们可以发现热层参量的改善是有限的，尤其对于热层温度，这主要是由于(1)同化的电子密度观测是有限的；(2)在模式中主导热层温度的参量是太阳活动指数F₁₀₇，而在此次同化模拟实验中我们并没有对F₁₀₇这一驱动项进行改善.

我们利用CHAMP卫星中性质量密度观测和TIMED/GUVI积分氧氮比观测的构型，对同化实验结果进行采样，评估同化电离层TEC对这两种热层参量的优化效果.图 3分别给出了一天内CHAMP卫星轨道位置上的控制值、“真实”值以及被估计的中性质量密度值(图 3a)，以及一天内TIMED/GUVI观测路径上的控制值、“真实”值以及被估计的积分氧氮比值(图 3b).值得注意的是，这里的“真实”中性质量密度值和积分氧氮比值是基于上述两种真实观测构型，通过采样上述NCAR-TIEGCM计算的真实场结果获取的.从图 3a可以看出，在CHAMP卫星观测位置上，“真实”中性质量密度值大于控制值，这是由于“真实”场是在更大的太阳活动指数驱动下生成的.随着同化的进行，被估计的轨道位置处的中性质量密度逐渐变大，并且更加接近于“真实”值，这是由于TEC观测量被逐渐同化.图 3a下图表示的是被估计的中性质量密度相较于控制结果提升的百分比，可以看出被估计的中性质量密度值提升最高能达到接近40%.

图 3

Fig. 3

图 3 (a) 上：CHAMP卫星观测位置处的控制实验中性质量密度值(点线)、“真实”中性质量密度值(实线)以及被估计的中性质量密度值(虚线)在整个同化阶段随时间变化图.下：被估计的中性质量密度相较于控制结果提升的百分比. (b)上：TIMED/GUVI临边探测位置处的控制实验积分氧氮比值(灰色点)、“真实”积分氧氮比值(红色点)以及被估计的积分氧氮比值(青色点)在整个同化阶段随时间变化图.下：被估计的积分氧氮比值相较于控制结果提升的百分比 Fig. 3 (a) The top panel: the time evolution of neutral mass density sampled from CHAMP observations for the control run, the forecast value of OSSE and the simulated "truth" during the whole assimilation period. The dot, dashed and solid lines show the control, forecast and "truth" results, respectively. The bottom panel: the line represents the percentage improvement of the forecast neutral mass densities. (b) The top panel: the time evolution of the value of Σ[O/N2] sampled from TIMED/GUVI limb location for the control run, the forecast value of OSSE and the simulated "truth" during the whole assimilation period. The gray, cyan and red dots show the control, forecast and "truth" results, respectively. The bottom panel: the magenta line represents the percentage improvement of the forecast Σ[O/N2]

图 3b显示，在同化开始阶段，被估计积分氧氮比值逐渐接近于“真实”值，这说明同化TEC有效优化了中性成分.随着同化进行，提升百分比超过了100%，即被估计的积分氧氮比值出现过分改正的情况.为分析原因，图 4给出了在同化阶段各个世界时内“真实”积分氧氮比值与被估计值之间的差值图，图中同时给出了峰值电子密度等值线图.

图 4

Fig. 4

图 4 在整个同化阶段，GUVI临边探测采样位置处的“真实”积分氧氮比值与被估计的积分氧氮比值差值图.图中的青色等值线表示的是赤道异常区位置 Fig. 4 The differential value of Σ[O/N2] sampled from GUVI limb observations between the forecast value of OSSE and the simulated "truth" during the whole assimilation period. The cyan contour lines show the equatorial ionization anomaly region

从图 4可以看出，在同化初期(00:00 UT时刻至07:00 UT时刻)，随着TEC观测数据的不断同化，GUVI探测位置处的积分氧氮比值出现较大幅度的提升.但随着同化的继续进行，被估计的积分氧氮比值与“真实”值开始出现较大偏差，这些差值较大的地方主要发生在白天赤道异常区位置和极区附近.这是由于该区域处电子密度水平梯度变化剧烈，使得已知观测参量与未观测参量之间的相关性可能被错误地统计，进而影响中性成分估计.这与Matsuo等(2013)实验结果是一致的.同时，结合图 2b结果，全球平均积分氧氮比值在整个阶段是提升的，说明在除这些电子密度水平梯度变化剧烈区域以外的其他位置处，中性成分整体上表现为好的改善.

上面已经论证了热层参量(温度、中性成分)能够通过同化电离层TEC观测得到较好的校正，接下来我们尝试对CHAMP卫星观测位置处的中性质量密度预报效果进行评估.这里依然采用观测系统模拟实验，预报时间段选择为7月12日21:00 UT至13日21:00 UT，值得注意的是，在整个24 h预报阶段太阳活动指数F₁₀₇以及地磁活动指数Kp始终为78和1.图 5给出了在选取的预报24 h内CHAMP卫星轨迹位置上的控制值、“真实”值和预报的中性质量密度值，以及预报的中性质量密度相较于控制结果提升的百分比.这里的“真实”中性质量密度值仍然是基于同一天内CHAMP真实观测构型，通过采样上述NCAR-TIEGCM计算的真实场结果获取的.从图中可以看出，随着预报的进行，预报的中性质量密度值逐渐变小，并且更加接近于控制值，这是由于在同化阶段被调整的热层参量随着预报时间的推移慢慢恢复到控制状态值.根据图 5可知，经过了24 h的预报阶段，预报的中性质量密度相较于控制结果提升的百分比仍然可以达到5%左右，这说明中性质量密度预报时间尺度能够达到近24 h.

图 5

Fig. 5

图 5 上：CHAMP卫星观测位置处的控制实验中性质量密度值(点线)、“真实”中性质量密度值(实线)以及预报的中性质量密度值(虚线)在整个预报阶段内随时间变化图.下：预报的中性质量密度相较于控制结果提升的百分比 Fig. 5 The top panel: the time evolution of neutral mass density sampled from CHAMP observations during the forecast stage of the control run, the forecast value of OSSE and the simulated "truth" during the whole forecast period. The dot, dashed and solid lines show the control, forecast and "truth" results, respectively. The bottom panel: the line from the bottom panel represents the percentage improvement of the forecast neutral mass densities

为了研究哪一个热层参量在中性质量密度预报改进中起主导作用，图 6给出了CHAMP卫星观测位置处的“真实”、控制以及预报中性温度预报对比结果.图 7给出了“真实”、控制以及预报的全球平均积分氧氮比值对比结果.从图 6可以看出，随着预报的进行，预报的中性温度逐渐变小，并且很快接近控制值.在接近12个小时预报之后，预报的中性温度相较于控制结果提升的百分比已经低于5%.然而，根据图 7可以看出，在预报的第12个小时，预报的积分氧氮比值相较于控制结果提升的百分比为40.45%，在接近预报结束时刻提升百分比仍保持为28.47%.这说明中性成分可能在提高中性质量密度预报时长上有一定作用.

图 6

Fig. 6

图 6 上：预报阶段CHAMP卫星观测采样位置处的控制实验中性温度值(点线)、“真实”中性温度值(实线)以及预报的中性温度值(虚线)在整个预报阶段内随时间变化图.下：预报的中性温度相较于控制结果提升的百分比 Fig. 6 The top panel: the time evolution of neutral mass density sampled from CHAMP observations for the control run, the forecast value of OSSE and the simulated "truth" during the whole forecast period. The dot, dashed and solid lines show the control, forecast and "truth" results, respectively. The bottom panel: the line represents the percentage improvement of the forecast neutral mass densities

图 7

Fig. 7

图 7 预报阶段全球平均的控制实验积分氧氮比值(圆形标注线)、“真实”平均积分氧氮比值(三角形标注线)以及预报的平均积分氧氮比值(正方形标注线)在整个预报阶段内随时间变化图 Fig. 7 The time evolution of global averaged ∑[O/N2] for the control run, the forecast value of OSSE and the simulated "truth" during the whole forecast period. The lines with circle, square and triangle show the control, forecast and "truth" results, respectively

本文通过观测系统模拟实验，研究发现：基于开发的热层电离层耦合数据同化系统，通过同化电离层TEC观测，可以有效校正热层中性质量密度和积分氧氮比值，进而提高电离层热层的预报能力.鉴于热层观测相对较少，而电离层观测正急剧增加，这也为未来热层参量估算提供了一条新的途径.同时也发现，中性成分在电子密度水平梯度变化剧烈区域效果较差，说明了只同化TEC观测对热层成分估计是不足的，后续为了更好地估计氧氮比参量开展多参量的数据同化是很有必要的.中性质量密度估计存在效果，但更优的密度估计需要对模式驱动——太阳辐射这一参量进行合理的估计.

目前，国内外开发的数据同化现报系统，着眼于区域或全球电离层环境实时监测的目的，可为卫星导航、雷达成像、短波通信等科学研究和工程应用提供及时、准确、有效的电离层参量信息(阿尔察等，2018).本文使用的电离层热层集合卡尔曼滤波同化系统采用稀疏矩阵来进行相关的超大矩阵存储和计算，大大节约了计算量，使得同化实验能够在普通工作站上完成，这为今后实时运行该系统进行电离层环境预报提供了可能性.