CN116203591B - 一种基于多站点联合估计中国区域高精度电离层vtec方法 - Google Patents

Abstract

本发明公开了一种基于多站点联合估计中国区域高精度电离层VTEC方法，包括如下步骤：步骤1，获取多站点GNSS接收机与卫星链路STEC结果：步骤2，计算接收机与卫星间仰角、方位角和地心夹角以及IPP穿刺点处经度、纬度：步骤3，获得自适应网格及差分方程：步骤4，求解各接收机和卫星硬件延迟：步骤5，获取IPP穿刺点处VTEC_real结果，采用改进Kriging插值求取区域电离层VTEC变化。本发明所公开的方法，充分利用了多站点观测数据，相邻IPP穿刺点处VTEC不会存在突变，可用于获取特定区域及周边范围内高精度电离层TEC变化，在电离层空间环境监测以及电波信号延迟效应评估等方面具有重要应用价值。

Description

一种基于多站点联合估计中国区域高精度电离层VTEC方法

技术领域

本发明属于空间环境遥感研究领域，特别涉及该领域中的一种基于多站点联合估计中国区域高精度电离层VTEC方法。

背景技术

电离层作为日地空间环境的重要组成部分，与人类活动有着极其密切的联系。特别是近些年我国航天事业蓬勃发展，电离层作为日地空间中航天器、卫星分布最广的区域，其异常活动会对航天器飞行安全、卫星导航定位等造成巨大影响。电离层TEC即电离层总电子含量是表征电离层活动的重要参数之一，通过反演区域电离层垂直电离层总电子含量(VTEC)变化对于研究空间天气环境具有极其重要的意义。

当前区域电离层VTEC反演多采用基于实测数据驱动电离层模型的数据同化方法，该方法受限于实测数据质量，所采用的电离层模型只能反映电离层平静期变化。部分学者利用GNSS监测站观测数据，采用球谐函数、多项式等硬件延迟(DCB)估计方法反演区域电离层VTEC变化，但未考虑当多站点解算的IPP穿刺点位置相近时VTEC约束问题，具有区域电离层VTEC反演精度欠佳等缺点。

发明内容

本发明所要解决的技术问题就是提供一种基于多站点联合估计中国区域高精度电离层VTEC方法，利用中国区域范围内多站点GNSS监测站观测数据，根据站点分布和卫星位置自适应划分网格，通过假定同一网格内不同站点IPP穿刺点VTEa相同，开展各站点接收机硬件延迟和卫星硬件延迟联合估计，利用获得的多站点IPP穿刺点VTEC，采用改进Kriging插值方法，实现中国区域电离层VTEC高精度反演。

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

一种基于多站点联合估计中国区域高精度电离层VTEC方法，采用多站点联合自适应网格法精确获取IPP穿刺点VTEC，利用改进Kriging插值算法实现区域电离层VTEC高精度反演。涉及的主要技术包括：采用载波相位平滑伪距法获取GNSS接收机与卫星链路高精度电离层总电子含量(STEC)，根据链路分布自适应划分网格数量解算各站点接收机硬件延迟和卫星硬件延迟，进而获取IPP穿刺点处高精度VTEC，采用改进Kriging插值方法，实现区域电离层VTEC高精度反演。其改进之处在于，包括如下步骤：

步骤1，获取多站点GNSS接收机与卫星链路STEC结果：

步骤11，计算码伪距、相伪距STEC：

上式中，STEC_P、STEC_L分别为采用码伪距和相伪距计算得到的STEC，f₁、f₂为GNSS卫星信号频率，λ₁、λ₂为GNSS卫星信号波长，P₁、P₂为码伪距，L₁、L₂为相位测量值(含有整周模糊度)，STEC_P为绝对值精度较低，STEC_L为相对值精度高；

步骤12，通过检测周跳利用STEC_P、STEC_L将整周模糊度进行补偿，计算得到含有硬件延迟的高精度STEC_obs：

STEC_obs＝STEC_real+B_r+B_s

上式中，STEC_real为剔除接收机和卫星硬件延迟后电离层总电子含量，B_r、B_s分别为接收机硬件延迟和卫星硬件延迟，为待求参量；

步骤2，计算接收机与卫星间仰角、方位角和地心夹角以及IPP穿刺点处经度、纬度：

步骤21，将GNSS接收机经度、纬度、高度记为lat_r、lon_r、height_r，将GNSS卫星经度、纬度、高度记为lat_s、lon_s、height_s，通过下式将经纬高坐标(lat、lon、hight)转为地心地固坐标(X、Y、Z)：

上式中，为基准椭球体的曲率半径，R_e为地球半径，为椭球偏心率，a为地球椭球长半轴，c为地球椭球短半轴；

步骤22，通过下式将地心地固坐标(X、Y、Z)转换为东北天坐标(E、N、U)：

上式中，E、N、U为东北天坐标(站心直角坐标)，S为坐标变换矩阵，X_s、Y_s、Z_s为GNSS卫星地心地固坐标，X_r、Y_r、Z_r为GNSS接收机地心地固坐标；

步骤23，计算GNSS接收机与GNSS卫星仰角Ele、方位角Az：

上式中，atan2为4象限反正切函数，每30s计算一次接收机与卫星仰角，每15min对计算结果取均值，剔除仰角低于15°的数据；

步骤24，计算GNSS接收机与GNSS卫星地心夹角ψ：

上式中，h为电离层等效高度；

步骤25，计算GNSS接收机与GNSS卫星IPP穿刺点处经度lat、纬度lon：

lat＝asin(sin(lat_r)cosψ+cos(lat_r)sinψcos(Az))

每30s计算一次IPP穿刺点处经度、纬度，每15min对计算结果取均值；

步骤3，获得自适应网格及差分方程：

步骤31，计算投影函数SF：

步骤32，计算IPP穿刺点处剔除GNSS接收机与GNSS卫星硬件延迟VTEC_real：

VTEC_real＝(STEC_obs-B_r-B_s)·SF

步骤33，根据IPP穿刺点位置及分布，按照经度、纬度自适应划分网格，假定同一网格内不同IPP穿刺点处VTEC_real相同，每15min建立差分方程：

VTEC_{real_IPPi}-VTEC_{real_IPPj}＝0

上式中，VTEC_{real_IPPi}为IPP穿刺点i的VTEC_real值，VTEC_{real_IPPj}为IPP穿刺点j的VTEC_real值，且穿刺点i和j在同一网格内；

步骤4，求解各接收机和卫星硬件延迟：

步骤41，累积一天观测数据，得到超定线性方程组：

Ax＝b

上式中，A为N_obs×(n_r+n_s)的矩阵，n_r为接收机数量，n_s为卫星数量，N_obs为一天累计有效观测量，x为(n_r+n_s)×1的向量，为待求解的未知系数，由n_r个接收机硬件延迟B_r和n_s颗GNSS卫星硬件延迟B_s组成，b为N_obs×1的向量；

步骤42，采用最小二乘法求解未知数x，得到各接收机B_r和卫星B_s：

x＝(A^TA)^-1A^Tb

上式中，上标T为矩阵转置，上标-1为矩阵求逆；

步骤5，获取IPP穿刺点处VTEC_real结果，采用改进Kriging插值求取区域电离层VTEC变化。

本发明的有益效果是：

本发明所公开的方法，充分利用了多站点观测数据，相邻IPP穿刺点处VTEC不会存在突变，可用于获取特定区域及周边范围内高精度电离层TEC变化，在电离层空间环境监测以及电波信号延迟效应评估等方面具有重要应用价值。

附图说明

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

图2给出了某时刻IPP穿刺点分布情况；

图3给出了采用改进Kriging插值引起的误差；

图4给出了区域高精度电离层VTEC变化情况。

具体实施方式

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

实施例1，本实施例公开了一种基于多站点联合估计中国区域高精度电离层VTEC方法，主要利用中国区域范围内多站点GNSS监测站观测数据，通过采用载波相位平滑伪距法获取高精度站星链路STEC，根据站点分布和卫星位置自适应划分网格，并假定同一网格内不同站点IPP穿刺点VTEC相同，解算各站点接收机硬件延迟和卫星硬件延迟，获取IPP穿刺点处高精度绝对VTEC结果，采用改进Kriging插值方法实现中国区域电离层VTEC高精度反演。如图1所示，包括如下步骤：

步骤1，以第i个GNSS接收机和第j颗GNSS卫星链路为例，说明获取GNSS接收机与卫星链路STEC结果方法：

步骤11，计算码伪距、相伪距STEC：

码伪距法计算第i个GNSS接收机和第j颗GNSS卫星链路STEC：

相伪距法计算第i个GNSS接收机和第j颗GNSS卫星链路STEC：

上式中，f₁、f₂为GNSS卫星信号频率，λ₁、λ₂为GNSS卫星信号波长，P₁、P₂为码伪距，L₁、L₂为相位测量值；对于GPS卫星而言f₁＝1575.42MHz、f₂＝1227.60MHz、λ₁＝19.04cm、λ₂＝24.44cm，对于北斗卫星而言f₁＝1575.42MHz、f₂＝1207.14MHz、λ₁＝19.04cm、λ₂＝24.85cm。

步骤12，通过检测周跳利用将整周模糊度进行补偿，计算得到含有硬件延迟的高精度/>

上式中，为剔除接收机和卫星硬件延迟后电离层总电子含量，/>分别为接收机硬件延迟和卫星硬件延迟，为待求参量；

步骤2，计算第i个GNSS接收机与第j颗GNSS卫星间仰角、方位角和地心夹角以及IPP穿刺点处经度、纬度：

步骤21，对于第i个GNSS接收机，其经度、纬度、高度记为将经纬高坐标转为地心地固坐标：

对于第j颗GNSS卫星，其经度、纬度、高度记为将经纬高坐标转为地心地固坐标：

上式中，为基准椭球体的曲率半径，R_e为地球半径，R_e＝6371km，/>为椭球偏心率，a为地球椭球长半轴，a＝6378137m，c为地球椭球短半轴，c＝6356752.3141m；

步骤22，通过下式将第i个接收机地心地固坐标与第j颗卫星地心地固坐标/>转换为东北天站心直角坐标E^ij、N^ij、U^ij：

上式中，S为坐标变换矩阵：

步骤23，计算第i个GNSS接收机与第j颗GNSS卫星仰角Ele^ij和方位角Az^ij：

上式中，atan2为4象限反正切函数，每30s计算一次接收机与卫星仰角，每15min对计算结果取均值，剔除仰角低于15°的数据；

步骤24，计算第i个GNSS接收机与第j颗GNSS卫星地心夹角ψ^ij：

上式中，电离层等效高度h＝450km；

步骤25，计算第i个GNSS接收机与第j颗GNSS卫星IPP穿刺点处经度lat^ij、纬度lon^ij：

每30s计算一次IPP穿刺点处经度、纬度，每15min对计算结果取均值；

步骤3，以第i个GNSS接收机第j颗GNSS卫星链路和第m个GNSS接收机第n颗GNSS卫星链路为例，获得自适应网格及差分方程：

步骤31，计算第i个GNSS接收机与第j颗GNSS卫星链路投影函数SF^ij：

计算第m个GNSS接收机与第n颗GNSS卫星链路投影函数SF^mn：

步骤32，计算第i个GNSS接收机与第j颗GNSS卫星链路IPP穿刺点处

计算第m个GNSS接收机与第n颗GNSS卫星链路IPP穿刺点处

步骤33，按照经度、纬度2.5°×2.5°网格，根据IPP穿刺点位置及分布自适应划分监测网格，假定同一网格内不同IPP穿刺点处相同，每15min某一网格内建立如下差分方程：

上式中，n_r为站点接收机数量，n_s为GNSS卫星数量，k为同一网格内IPP穿刺点个数，根据每15min接收机与卫星链路分布情况，自适应划分网格数量和网格内IPP穿刺点个数。

步骤4，求解各接收机和卫星硬件延迟：

步骤41，累积一天观测数据，得到超定线性方程组：

Ax＝b

上式中，A为N_obs×(n_r+n_s)的矩阵，n_r为接收机数量，n_s为卫星数量，N_obs为一天累计有效观测量，x为(n_r+n_s)×1的向量，为待求解的未知系数，由n_r个接收机硬件延迟B_r和n_s颗GNSS卫星硬件延迟B_s组成，b为N_obs×1的向量；

步骤42，采用最小二乘法求解未知数x，得到各接收机B_r和卫星B_s：

x＝(A^TA)^-1A^Tb

上式中，上标T为矩阵转置，上标-1为矩阵求逆；

步骤5，获取IPP穿刺点处VTEC_real结果，采用改进Kriging插值求取区域电离层VTEC变化。其中经度权重为1，纬度权重为3。图2给出了某时刻IPP穿刺点分布情况；图3给出了采用改进Kriging插值引起的误差；图4给出了区域高精度电离层VTEC变化情况。

Claims (1)

1.一种基于多站点联合估计中国区域高精度电离层VTEC方法，其特征在于，包括如下步骤：

步骤1，获取多站点GNSS接收机与卫星链路STEC结果：

步骤11，计算码伪距、相伪距STEC：

上式中，STEC_P、STEC_L分别为采用码伪距和相伪距计算得到的STEC，f₁、f₂为GNSS卫星信号频率，λ₁、λ₂为GNSS卫星信号波长，P₁、P₂为码伪距，L₁、L₂为相位测量值；

步骤12，通过检测周跳利用STEC_P、STEC_L将整周模糊度进行补偿，计算得到含有硬件延迟的高精度STEC_obs：

STEC_obs＝STEC_real+B_r+B_s

上式中，STEC_real为剔除接收机和卫星硬件延迟后电离层总电子含量，B_r、B_s分别为接收机硬件延迟和卫星硬件延迟，为待求参量；

步骤2，计算接收机与卫星间仰角、方位角和地心夹角以及IPP穿刺点处经度、纬度：

步骤21，将GNSS接收机经度、纬度、高度记为lat_r、lon_r、height_r，将GNSS卫星经度、纬度、高度记为lat_s、lon_s、height_s，通过下式将经纬高坐标(lat、lon、hight)转为地心地固坐标(X、Y、Z)：

上式中，为基准椭球体的曲率半径，R_e为地球半径，为椭球偏心率，a为地球椭球长半轴，c为地球椭球短半轴；

步骤22，通过下式将地心地固坐标(X、Y、Z)转换为东北天坐标(E、N、U)：

上式中，E、N、U为东北天坐标，S为坐标变换矩阵，X_s、Y_s、Z_s为GNSS卫星地心地固坐标，X_r、Y_r、Z_r为GNSS接收机地心地固坐标；

步骤23，计算GNSS接收机与GNSS卫星仰角Ele、方位角Az：

上式中，atan2为4象限反正切函数，每30s计算一次接收机与卫星仰角，每15min对计算结果取均值，剔除仰角低于15°的数据；

步骤24，计算GNSS接收机与GNSS卫星地心夹角ψ：

上式中，h为电离层等效高度；

步骤25，计算GNSS接收机与GNSS卫星IPP穿刺点处经度lat、纬度lon：

lat＝asin(sin(lat_r)cosψ+cos(lat_r)sinψcos(Az))

每30s计算一次IPP穿刺点处经度、纬度，每15min对计算结果取均值；

步骤3，获得自适应网格及差分方程：

步骤31，计算投影函数SF：

步骤32，计算IPP穿刺点处剔除GNSS接收机与GNSS卫星硬件延迟VTEC_real：

VTEC_real＝(STEC_obs-B_r-B_s)·SF

步骤33，根据IPP穿刺点位置及分布，按照经度、纬度自适应划分网格，假定同一网格内不同IPP穿刺点处VTEC_real相同，每15min建立差分方程：

VTEC_{real_IPPi}-VTEC_{real_IPPj}＝0

上式中，VTEC_{real_IPPi}为IPP穿刺点i的VTEC_real值，VTEC_{real_IPPj}为IPP穿刺点j的VTEC_real值，且穿刺点i和j在同一网格内；

步骤4，求解各接收机和卫星硬件延迟：

步骤41，累积一天观测数据，得到超定线性方程组：

Ax＝b

上式中，A为N_obs×(n_r+n_s)的矩阵，n_r为接收机数量，n_s为卫星数量，N_obs为一天累计有效观测量，x为(n_r+n_s)×1的向量，为待求解的未知系数，由n_r个接收机硬件延迟B_r和n_s颗GNSS卫星硬件延迟B_s组成，b为N_obs×1的向量；

步骤42，采用最小二乘法求解未知数x，得到各接收机B_r和卫星B_s：

x＝(A^TA)^-1A^Tb

上式中，上标T为矩阵转置，上标-1为矩阵求逆；

步骤5，获取IPP穿刺点处VTEC_real结果，采用改进Kriging插值求取区域电离层VTEC变化。