CN112649899A - 一种全球电离层数据同化和预报方法 - Google Patents

Abstract

本发明公开了一种全球电离层数据同化和预报方法，包括获取全球分布的地基GNSS接收机、LEO卫星掩星、垂测仪、卫星信标接收机和小型光度计探测数据等步骤。本发明基于数据最优化估计理论，提供了一种全球电离层数据同化和预报方法，该方法采用限带卡尔曼滤波（Band‑Limited Kalman Filter）和高斯—马尔科夫过程（Gauss‑Markov Process），结合稀疏矩阵存储和快速处理算法，有效降低了全球电离层数据同化大规模矩阵运算对计算机内存的需求，同时大大减少数据同化过程所需的运行时间，可实现全球地基和天基多源探测数据的实时同化和短期预报。对于实现业务化的全球电离层环境参量现报和预报，提升无线电信息系统电离层环境信息保障能力具有重要意义。

Description

一种全球电离层数据同化和预报方法

技术领域

本发明属于空间环境现报预报技术应用领域，特别涉及该领域中的一种全球电离层数据同化和预报方法。

背景技术

电离层作为日地空间环境的重要组成部分，会对穿越其中的无线电波产生折射、反射、散射和吸收等效应，从而影响卫星导航、通信、雷达等诸多无线电信息系统的性能。利用各类技术手段来获取电离层的特征参量，研究电离层中的各种现象，从而揭示其背后的物理机制，具有重要的科学乃至政治、军事和经济意义。

随着地基和天基电离层探测网络的不断发展和日益完善，各类电离层探测数据开始逐渐增多，如何有效利用这些数据满足电离层精确现报和预报的要求，从而构建业务化电离层现报预报系统，已成为当前电离层研究的热点方向。数据同化技术作为在现代数值天气预报中广泛应用的一种技术，它能够对各种天地基观测数据进行综合利用，从而把各种时间和空间上既不规则而且零散分布的观测资料“吸收”到背景模式中，达到观测数据与背景模式之间的深度融合。数据同化技术在电离层研究领域开始受到世界各国研究者极大关注。目前，电离层数据同化大多采用三维/四维变分和Kalman滤波算法。变分同化算法既需要进行背景模型的后向预报，也需要计算模型切线性的正向积分和伴随模式的后向积分，实现难度非常大；由于全球电离层数据同化待求解的状态变量的维度非常高(通常大于20万)，卡尔曼滤波同化中同样面临增益矩阵求逆、协方差矩阵传播过程中不可避免的超大规模稀疏矩阵存储和处理的难题，构建业务化的电离层实时数据同化和预报系统尚有不少关键问题待进一步完善和解决。

发明内容

本发明所要解决的技术问题就是提供一种全球电离层数据同化和预报方法，可为构建业务化的电离层数据同化和预报系统提供重要的技术支撑。

本发明采用如下技术方案：

一种全球电离层数据同化和预报方法，其改进之处在于，包括如下步骤：

步骤1，获取全球分布的地基GNSS接收机、LEO卫星掩星、垂测仪、卫星信标接收机和小型光度计探测数据，分别提取电离层总电子含量TEC和电子密度的观测数据，按照设定的时间窗口将数据存储在一个指定的文件内，该文件包括：数据类型编号、观测时刻、测站经纬高坐标、卫星测站经纬高坐标、TEC或电子密度测量值；

步骤2，获取电离层数据同化和预报需要的日地空间环境参量数据，包括太阳辐射通量和地磁指数，将数据同化所对应的太阳辐射通量、地磁指数及模型运行所需要的配置参量，包括运行模式设定、输入的观测数据存储路径、输出数据的存储和网格分辨率设定，一并存储到指定的配置文件内；

步骤3，构建电离层数据同化观测算子，将电离层TEC和电子密度数据按照“线”型和“点”型观测数据进行分类；

对于“点”型电子密度观测数据，插值方法如下：

其中：Ne表示电子密度，K表示级数展开的级数，r表示地理位置，h_k(r)为基函数，a_k为插值系数，基函数计算方法如下：

对于“线”型TEC观测数据，根据接收机与卫星间的几何位置关系，TEC和同化估计的状态参量之间的关系为：

其中：STEC表示电离层总电子含量，N_e表示电离层电子密度值，积分上下限分别为卫星和接收机位置，因此几何转换矩阵计算方法如下：

其中ΔR_in表示第i个TEC观测值在第n个网格内的射线截距，i为TEC观测值的序号，N为电子密度网格总数，K为展开的总级数；

遍历同化窗口内所有电离层电子密度和TEC数据，根据式(1)和式(3)同化观测方程可构建为：

y＝Hx (4)

其中：y为电离层电子密度和TEC数据组成的矢量，矩阵H为观测数据与同化状态参量“电子密度”对应的几何转换矩阵，x为待估计的各网格点的电子密度；

将矩阵H按照行优先稀疏矩阵CSR格式存储，即只存储H矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤4，构建观测误差协方差矩阵，根据不同观测数据间观测误差的不同，计算数据同化的观测误差协方差矩阵R，将观测误差协方差矩阵R用对角矩阵表示：

其中i，j表示网格点的位置，d为观测的TEC或电子密度值，R_ij为观测点间的误差协方差值，α是比例系数；

将矩阵R按照CSR格式存储，即只存储R矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤5，构建背景场误差协方差矩阵，确定电离层相关距离；

电离层垂直方向相关距离的计算方法如下：

分别表示i和j点在垂直方向的电离层相关长度；

电离层水平相关距离计算方法如下：

其中，

和L_λ分别表示纬向和经向的相关长度，θ表示两个网格点之间的方位角，

其中，

表示磁纬度，z表示高度，γ是地方时LT的函数，表示为：

在存储数据同化背景场误差协方差矩阵P时，设定的限带方法是任意两点间的距离超过1200km，其背景协方差为零值，根据(6)和(7)式，数据同化背景场误差协方差矩阵P具体表示为：

其中：

分别代表背景模型在i点和j点的背景值，e为自然指数，z_ij代表第i点和第j点在高度上的距离，r_ij代表第i点和第j点的直线距离，L_z分别表示电离层在高度方向的相关距离，L_H分别表示电离层在水平方向的相关距离，β表示模式误差与模式值之间的比例系数，g_ij代表第i点和第j点在水平方向上的大圆距离，计算方法如下：

式中，

是第i点的纬度和经度，

是第j点的纬度和经度；

将矩阵P按照CSR格式存储，即只存储P矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤6，读取步骤2设定的配置文件参量，读取太阳辐射通量和地磁指数，根据设定的数据同化网格分辨率，运行电离层背景模型，计算背景场电离层电子密度分布x_b；

步骤7，采用限带卡尔曼滤波算法进行地基和天基TEC或电子密度观测数据的同化，模型使用的卡尔曼滤波方程如下：

x_a＝x_b+PH^T(R+HPH^T)^-1(y-Hx_b) (14)

其中x_a表示同化后电子密度分析场；

对矩阵作以下变换：

(R+HPH^T)^-1(y-Hx_b)＝T (15)

(R+HPH^T)T＝(y-Hx_b) (16)

采用广义最小残差法求解方程，使用不完全下三角矩阵—上三角矩阵LU分解算法进行预处理，具体方法如下：

C^-1[(R+HPH^T)T-(y-Hx_b)]＝0 (17)

其中：C是预条件矩阵，同化给出电离层TEC和电子密度现报结果后，将相关现报结果存储到指定的输出路径下；

步骤8，利用高斯—马尔科夫预报方法对同化的电子密度进行预报，具体公式为：

其中，t表示预报时间，K为增益矩阵，L代表转换矩阵，其表现形式为对角矩阵，具体计算方法为：

其中：ΔT为提前预报的时间间隔，τ表示电离层的时间相关尺度，取5h，磁暴期间取值缩短为2h，同化给出电离层TEC和电子密度预报结果后，将相关现报结果存储到指定的输出路径下。

进一步的，步骤(2)中的运行模式为区域或全球。

本发明的有益效果是：

本发明基于数据最优化估计理论，提供了一种全球电离层数据同化和预报方法，该方法采用限带卡尔曼滤波(Band-Limited Kalman Filter)和高斯—马尔科夫过程(Gauss-Markov Process)，结合稀疏矩阵存储和快速处理算法，有效降低了全球电离层数据同化大规模矩阵运算对计算机内存的需求，同时大大减少数据同化过程所需的运行时间，可实现全球地基和天基多源探测数据的实时同化和短期预报。对于实现业务化的全球电离层环境参量现报和预报，提升无线电信息系统电离层环境信息保障能力具有重要意义。

本发明方法能够灵活配置同化区域范围和数据产品的空间分辨率，实现全球中等分辨率，区域高分辨率的电离层现报和预报。

本发明方法能够兼容目前主流的电离层探测数据，包括地基GNSS接收机、LEO卫星掩星、垂测仪、卫星信标接收机、小型光度计等，从而保证同化结果的精度和可靠性。

本发明方法无需超级计算机，仅利用个人工作站笔记本(Intel i7处理器，64G内存，2T硬盘)，即可实现全球电离层TEC和电子密度参量的现报和预报，一次同化过程仅需3-5分钟，从而满足电离层数据同化和预报的实时性要求，且方法具备非常高的运行稳定性。

附图说明

图1是本发明方法的流程示意图；

图2是电离层数据同化给出的电离层TEC对比结果图；

图3是电离层数据同化给出的电离层NmF2对比结果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例1，本实施例公开了一种全球电离层数据同化和预报方法，利用全球分布的地基GNSS接收机、LEO卫星掩星、垂测仪、卫星信标接收机、小型光度计等探测数据，基于限带卡尔曼滤波数据同化和预报方法，实现区域或全球电子密度、电子总含量的现报和预报，为构建业务化的电离层数据同化和预报系统奠定基础。如图1所示，具体包括如下步骤：

步骤1，获取全球分布的地基GNSS接收机、LEO卫星掩星、垂测仪、卫星信标接收机和小型光度计等探测数据，分别提取电离层总电子含量TEC和电子密度的观测数据，为满足数据同化观测算子构建的要求，分别按照设定的时间窗口将数据存储在一个指定的文件内，文件应包括以下要素：数据类型编号、观测时刻、测站经纬高坐标、卫星测站经纬高坐标、TEC或电子密度测量值；

步骤2，获取电离层数据同化和预报需要的日地空间环境参量数据，包括太阳辐射通量和地磁指数等，将数据同化所对应的太阳辐射通量、地磁指数及模型运行所需要的配置参量，包括运行模式设定(区域或全球)、输入的观测数据存储路径、输出数据的存储和网格分辨率设定等，一并存储到指定的配置文件内，便于电离层数据同化和预报模型读取；

步骤3，构建电离层数据同化观测算子，将电离层TEC和电子密度数据按照“线”型和“点”型观测数据进行分类；

对于“点”型电子密度观测数据，由于观测参量与同化待估计的状态参量相同，因此只需要利用简单的插值即可建立几何转换矩阵，插值方法如下：

其中：Ne表示电子密度，K表示级数展开的级数，r表示地理位置，h_k(r)为基函数，a_k为插值系数，基函数计算方法如下：

对于“线”型TEC观测数据，根据接收机与卫星间的几何位置关系，TEC和同化估计的状态参量(电子密度)之间的关系为：

其中：STEC表示电离层总电子含量，N_e表示电离层电子密度值，积分上下限分别为卫星和接收机位置，因此几何转换矩阵计算方法如下：

其中ΔR_in表示第i个TEC观测值在第n个网格内的射线截距，i为TEC观测值的序号，N为电子密度网格总数，K为展开的总级数；

遍历同化窗口内所有电离层电子密度和TEC数据，根据式(1)和式(3)同化观测方程可构建为：

y＝Hx (4)

其中：y为电离层电子密度和TEC数据组成的矢量，矩阵H为观测数据与同化状态参量“电子密度”对应的几何转换矩阵，x为待估计的各网格点的电子密度；

将矩阵H按照行优先稀疏矩阵CSR格式存储，即只存储H矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身，以此减少矩阵存储所需要的内存消耗；

步骤4，构建观测误差协方差矩阵，根据不同观测数据间观测误差的不同，计算数据同化的观测误差协方差矩阵R，将观测数据间视为相互独立，因此观测误差协方差矩阵R可以用对角矩阵表示：

其中i，j表示网格点的位置，d为观测的TEC或电子密度值，R_ij为观测点间的误差协方差值，α是比例系数；

同样地，将矩阵R按照CSR格式存储，即只存储R矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤5，构建背景场误差协方差矩阵，确定电离层水平和垂直方向相关距离，同时根据限带卡尔曼滤波的需要，截断超过相关距离的协方差，计算数据同化的背景场误差协方差矩阵：

电离层垂直方向的相关距离与两个点所在的高度有关，计算方法如下：

分别表示i和j点在垂直方向的电离层相关长度；

电离层水平相关距离可以表示为经度、纬度相关距离及两者间方位角间的函数，计算方法如下：

其中，

和L_λ分别表示纬向和经向的相关长度，θ表示两个网格点之间的方位角，

L_λ和L_z均是与地磁纬度、高度及地方时LT相关的函数，计算方法如下：

其中，

表示磁纬度(deg)，z表示高度(km)，γ是地方时LT的函数，表示为：

在存储数据同化背景场误差协方差矩阵P时，设定的限带方法是任意两点间的距离超过1200km，其背景协方差为零值，根据(6)和(7)式，数据同化背景场误差协方差矩阵P具体表示为：

其中：

分别代表背景模型在i点和j点的背景值，e为自然指数，z_ij代表第i点和第j点在高度上的距离，r_ij代表第i点和第j点的直线距离，L_z分别表示电离层在高度方向的相关距离，L_H分别表示电离层在水平方向的相关距离，β表示模式误差与模式值之间的比例系数，g_ij代表第i点和第j点在水平方向上的大圆距离，计算方法如下：

式中，

是第i点的纬度和经度，

是第j点的纬度和经度；

同样地，将矩阵P按照CSR格式存储，即只存储P矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤6，读取步骤2设定的配置文件参量，读取太阳辐射通量和地磁指数，根据设定的数据同化网格分辨率，运行电离层背景模型，计算背景场电离层电子密度分布x_b；

步骤7，采用限带卡尔曼滤波算法进行地基和天基TEC或电子密度观测数据的同化(利用卡尔曼滤波算法进行全球电离层数据同化现报，给出电离层TEC和电子密度现报结果)，模型使用的卡尔曼滤波方程如下：

x_a＝x_b+PH^T(R+HPH^T)^-1(y-Hx_b) (14)

其中x_b和x_a分别表示同化前的电离层电子密度背景场和同化后电子密度分析场，P和R分别代表背景场误差协方差和观测数据的误差协方差，H和y分别是观测算子几何转换矩阵和观测向量；

由于数据同化矩阵的维数很大，卡尔曼滤波增益矩阵是具有大型未知变量的稀疏线性方程，为了提升计算稳定性，对矩阵作以下变换：

(R+HPH^T)^-1(y-Hx_b)＝T (15)

(R+HPH^T)T＝(y-Hx_b) (16)

采用广义最小残差法求解方程，为了在迭代过程中保证解收敛，使用不完全下三角矩阵—上三角矩阵LU分解算法进行预处理，具体方法如下：

C^-1[(R+HPH^T)T-(y-Hx_b)]＝0 (17)

其中：C是预条件矩阵(Pre-Conditioner Matrix)，同化给出电离层TEC和电子密度现报结果后，将相关现报结果存储到指定的输出路径下；

步骤8，利用高斯—马尔科夫预报方法对同化的电子密度进行预报(进行全球电离层预报)，具体公式为：

其中，t表示预报时间，K为增益矩阵，L代表转换矩阵(transition matrix)，其表现形式为对角矩阵，具体计算方法为：

其中：ΔT为提前预报的时间间隔，τ表示电离层的时间相关尺度，一般情况下取5h，磁暴期间取值可缩短为2h，同化给出电离层TEC和电子密度预报结果后，将相关现报结果存储到指定的输出路径下。

图2是电离层数据同化给出的电离层TEC对比结果图；图3是电离层数据同化给出的电离层NmF2对比结果图。

综上所述，本发明基于数据最优化估计理论，提供了一种全球电离层数据同化和预报方法，该方法采用限带卡尔曼滤波(Band-Limited Kalman Filter)和高斯—马尔科夫过程(Gauss-Markov Process)，结合稀疏矩阵存储和矩阵快速处理算法，有效降低了全球电离层数据同化大规模矩阵运算对计算机内存的需求，同时大大减少数据同化过程所需的运行时间，可实现全球地基和天基多源探测数据的实时同化和短期预报。本发明可为构建业务化的电离层数据同化和预报系统提供重要的技术支撑。

Claims (2)

1.一种全球电离层数据同化和预报方法，其特征在于，包括如下步骤：

步骤1，获取全球分布的地基GNSS接收机、LEO卫星掩星、垂测仪、卫星信标接收机和小型光度计探测数据，分别提取电离层总电子含量TEC和电子密度的观测数据，按照设定的时间窗口将数据存储在一个指定的文件内，该文件包括：数据类型编号、观测时刻、测站经纬高坐标、卫星测站经纬高坐标、TEC或电子密度测量值；

步骤2，获取电离层数据同化和预报需要的日地空间环境参量数据，包括太阳辐射通量和地磁指数，将数据同化所对应的太阳辐射通量、地磁指数及模型运行所需要的配置参量，包括运行模式设定、输入的观测数据存储路径、输出数据的存储和网格分辨率设定，一并存储到指定的配置文件内；

步骤3，构建电离层数据同化观测算子，将电离层TEC和电子密度数据按照“线”型和“点”型观测数据进行分类；

对于“点”型电子密度观测数据，插值方法如下：

其中：Ne表示电子密度，K表示级数展开的级数，r表示地理位置，h_k(r)为基函数，a_k为插值系数，基函数计算方法如下：

对于“线”型TEC观测数据，根据接收机与卫星间的几何位置关系，TEC和同化估计的状态参量之间的关系为：

其中：STEC表示电离层总电子含量，N_e表示电离层电子密度值，积分上下限分别为卫星和接收机位置，因此几何转换矩阵计算方法如下：

其中ΔR_in表示第i个TEC观测值在第n个网格内的射线截距，i为TEC观测值的序号，N为电子密度网格总数，K为展开的总级数；

遍历同化窗口内所有电离层电子密度和TEC数据，根据式(1)和式(3)同化观测方程可构建为：

y＝Hx (4)

其中：y为电离层电子密度和TEC数据组成的矢量，矩阵H为观测数据与同化状态参量“电子密度”对应的几何转换矩阵，x为待估计的各网格点的电子密度；

将矩阵H按照行优先稀疏矩阵CSR格式存储，即只存储H矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤4，构建观测误差协方差矩阵，根据不同观测数据间观测误差的不同，计算数据同化的观测误差协方差矩阵R，将观测误差协方差矩阵R用对角矩阵表示：

其中i，j表示网格点的位置，d为观测的TEC或电子密度值，R_ij为观测点间的误差协方差值，α是比例系数；

将矩阵R按照CSR格式存储，即只存储R矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤5，构建背景场误差协方差矩阵，确定电离层相关距离；

电离层垂直方向相关距离的计算方法如下：

分别表示i和j点在垂直方向的电离层相关长度；

电离层水平相关距离计算方法如下：

其中，

和L_λ分别表示纬向和经向的相关长度，θ表示两个网格点之间的方位角，

其中，

表示磁纬度，z表示高度，γ是地方时LT的函数，表示为：

在存储数据同化背景场误差协方差矩阵P时，设定的限带方法是任意两点间的距离超过1200km，其背景协方差为零值，根据(6)和(7)式，数据同化背景场误差协方差矩阵P具体表示为：

其中：

分别代表背景模型在i点和j点的背景值，e为自然指数，z_ij代表第i点和第j点在高度上的距离，r_ij代表第i点和第j点的直线距离，L_z分别表示电离层在高度方向的相关距离，L_H分别表示电离层在水平方向的相关距离，β表示模式误差与模式值之间的比例系数，g_ij代表第i点和第j点在水平方向上的大圆距离，计算方法如下：

式中，

是第i点的纬度和经度，

是第j点的纬度和经度；

将矩阵P按照CSR格式存储，即只存储P矩阵中非零元素所在的起始行的序号、非零元素所在的列编号及非零元素本身；

步骤6，读取步骤2设定的配置文件参量，读取太阳辐射通量和地磁指数，根据设定的数据同化网格分辨率，运行电离层背景模型，计算背景场电离层电子密度分布x_b；

步骤7，采用限带卡尔曼滤波算法进行地基和天基TEC或电子密度观测数据的同化，模型使用的卡尔曼滤波方程如下：

x_a＝x_b+PH^T(R+HPH^T)^-1(y-Hx_b) (14)

其中x_a表示同化后电子密度分析场；

对矩阵作以下变换：

(R+HPH^T)^-1(y-Hx_b)＝T (15)

(R+HPH^T)T＝(y-Hx_b) (16)

采用广义最小残差法求解方程，使用不完全下三角矩阵—上三角矩阵LU分解算法进行预处理，具体方法如下：

C^-1[(R+HPH^T)T-(y-Hx_b)]＝0 (17)

其中：C是预条件矩阵，同化给出电离层TEC和电子密度现报结果后，将相关现报结果存储到指定的输出路径下；

步骤8，利用高斯—马尔科夫预报方法对同化的电子密度进行预报，具体公式为：

其中，t表示预报时间，K为增益矩阵，L代表转换矩阵，其表现形式为对角矩阵，具体计算方法为：

其中：ΔT为提前预报的时间间隔，τ表示电离层的时间相关尺度，取5h，磁暴期间取值缩短为2h，同化给出电离层TEC和电子密度预报结果后，将相关现报结果存储到指定的输出路径下。

2.根据权利要求1所述的全球电离层数据同化和预报方法，其特征在于：步骤(2)中的运行模式为区域或全球。

Images