CN114417580B - 一种观测系统对全球电离层数据同化性能的影响评估方法 - Google Patents

一种观测系统对全球电离层数据同化性能的影响评估方法 Download PDF

Info

Publication number
CN114417580B
CN114417580B CN202111675121.2A CN202111675121A CN114417580B CN 114417580 B CN114417580 B CN 114417580B CN 202111675121 A CN202111675121 A CN 202111675121A CN 114417580 B CN114417580 B CN 114417580B
Authority
CN
China
Prior art keywords
observation
data
assimilation
satellite
gnss
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111675121.2A
Other languages
English (en)
Other versions
CN114417580A (zh
Inventor
欧明
陈龙江
杨升高
许娜
王妍
吴家燕
陈亮
熊雯
冯健
甄卫民
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Institute of Radio Wave Propagation CETC 22 Research Institute
Original Assignee
China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Institute of Radio Wave Propagation CETC 22 Research Institute filed Critical China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority to CN202111675121.2A priority Critical patent/CN114417580B/zh
Publication of CN114417580A publication Critical patent/CN114417580A/zh
Application granted granted Critical
Publication of CN114417580B publication Critical patent/CN114417580B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/10—Devices for predicting weather conditions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

本发明公开了一种观测系统对全球电离层数据同化性能的影响评估方法，包括如下步骤：步骤A，地基和天基观测系统基本参数输入；步骤B，GNSS卫星和LEO掩星卫星的轨道坐标计算；步骤C，地基和天基电离层观测数据仿真；步骤D，基于Kalman滤波的全球电离层数据同化；步骤E，电离层TEC和电子密度同化性能评估。本发明所公开的方法，对观测系统采用仿真试验手段（OSSE），利用NeQuick模型模拟生成不同观测手段“真实”的观测数据，同时基于Kalman滤波同化算法对各种观测系统的观测数据进行同化处理，并采用总电子含量（TEC）和电子密度（Ne）的技术分（Skill Score）评估方法，对各类观测系统及其组合对全球电离层数据同化性能的影响进行定量评估。

Description

一种观测系统对全球电离层数据同化性能的影响评估方法
技术领域
本发明属于空间电离层环境数值预报领域，特别涉及该领域中的一种观测系统对全球电离层数据同化性能的影响评估方法。
背景技术
作为一种在现代气象数值天气预报中广泛应用的技术，数据同化是实现空间环境的现报和预报从气候学研究转向天气学应用的热点方向。数据同化能够实现对各种来源的观测数据的综合利用，把各种时空上不规则的零散分布的观测资料“吸收”到背景模式中，从而实现观测数据与背景模式的互补融合。随着人类对电离层现报和预报要求的不断提高，数据同化方法开始在电离层研究领域蓬勃发展。
电离层垂直探测仪(简称垂测)是最早也是目前最为常用的电离层探测手段之一，多年以来，国内外先后构建了一系列的电离层垂测网，累计了大量的电离层观测数据，这些为电离层数据同化提供了重要的数据支撑；地基GNSS探测是目前应用最为广泛的电离层探测手段，自上世纪90年代以来，随着全球卫星导航系统的发展，利用地基广泛分布的GNSS接收机监测电离层的时空变化获得了巨大的成功。当前，国际GNSS服务组织(IGS)已经公开发布数据的地基GNSS台站已经超过800个。此外，掩星探测在空间天气领域开始崭露头角，COSMIC、GRACE、SWARM等卫星均搭载了掩星接收机进行电离层探测；美国地球光学公司Planet IQ计划部署18颗卫星组成“持续地球无线电掩星群计划”(CICERO)星座，美国Spire卫星公司通过风险投资的形式，计划完成超过100颗小卫星组成的掩星探测星座。随着LEO掩星数量的逐步增加，基于小卫星星座的GNSS掩星观测数据的空间和时间分辨率正在快速提升，基于天基掩星数据的全球电离层数据同化将成为电离层精确估计的重要数据来源。
观测系统仿真试验(OSSE)作为一种数值试验，通常用于定量评估气象学界在气候分析和天气预报方面的潜在改进，OSSE既可以检验同化方法的有效性，也可以检验观测系统的有效性，是一种常用的试验手段。OSSE的框架通常包括观测系统数据模拟和数据同化两部分。在OSSE中，观测系统数据模拟是通过对观测系统应用前向观测算子模拟生成“真实”的观测数据；而数据同化部分则对这些观测数据进行同化处理和性能评估。OSSE是目前评估并优化观测系统设计最为重要的技术手段之一。
发明内容
本发明所要解决的技术问题就是针对全球电离层观测站网的设计和优化需求，提供一种观测系统对全球电离层数据同化性能的影响评估方法。
本发明采用如下技术方案：
一种观测系统对全球电离层数据同化性能的影响评估方法，其改进之处在于，包括如下步骤：
步骤A，地基和天基观测系统基本参数输入：
步骤A1，读取地基观测系统的配置参数，获取地基垂测、GNSS接收机的地理经纬高坐标和GNSS接收机的观测截止角；
步骤A2，读取天基观测系统的配置参数，获取GNSS卫星和LEO掩星卫星的两行轨道TLE星历、掩星接收机可视角FOV；
步骤B，GNSS卫星和LEO掩星卫星的轨道坐标计算：
步骤B1，解析步骤A2中得到的GNSS卫星和LEO掩星卫星的TLE星历，提取卫星编号、历元时刻、轨道倾角、历元时刻升交点赤经、轨道偏心率、历元时刻近地点幅角、历元时刻的平近点角、平均转速；
步骤B2，利用SGP4模型解算GNSS和LEO掩星卫星的地心坐标系ECI下的经度、纬度和高度坐标；
步骤B3，将卫星地心坐标系下的经度、纬度和高度坐标转换为地心地固坐标系ECEF下的经度、纬度和高度坐标；
步骤C，地基和天基电离层观测数据仿真：
步骤C1，设定太阳辐射指数F10.7、月份和UT时刻作为电离层模型NeQuick的输入；
步骤C2，将地基垂测的地理经纬度输入到NeQuick模型中，计算对应输入条件下，地基垂测站100千米至F2层峰值高度以下高度上的电子密度剖面Neion(r)；
步骤C3，输入地基GNSS台站地心地固坐标系下的纬度、经度、高度
Figure BDA0003450928620000021
和GNSS卫星的纬度、经度、高度
Figure BDA0003450928620000022
计算出两者间的仰角E，计算方法如下：
Δλ＝λd-λs (1)
rd＝Re+hd (2)
rs＝Re+hs (3)
Figure BDA0003450928620000023
Figure BDA0003450928620000024
其中，Re表示地球半径，Δλ表示GNSS台站与GNSS卫星二者间经度的差值，L表示二者之间的大圆距离，rd表示台站与地心之间的距离，rs表示卫星到地心的距离；
步骤C4，判断仰角是否大于GNSS接收机的观测截止角Eobs，若E大于Eobs，则将地基GNSS台站的经、纬、高度作为起始点，GNSS卫星轨道的经、纬、高度作为结束点坐标输入到NeQuick模型中，计算出两点间的倾斜总电子含量STECgnss；
步骤C5，将LEO掩星卫星的纬度、经度、高度
Figure BDA0003450928620000031
和GNSS卫星轨道的纬度、经度、高度
Figure BDA0003450928620000032
分别转化为地心地固坐标系ECEF下的坐标(Xt,Yt,Zt)和(Xs,Ys,Zs)，转换表达式表示为：
Figure BDA0003450928620000033
其中，
Figure BDA0003450928620000034
Re表示地球半径，e2＝0.00669437999013；
步骤C6，计算LEO卫星与LEO卫星-GNSS卫星射线之间的夹角θ和碰撞点高度，判定是否为掩星观测事件；
步骤C7，将满足掩星观测事件的LEO掩星卫星的经、纬、高度作为起始点，GNSS卫星轨道的经、纬、高度作为结束点坐标输入到NeQuick模型中，计算出两点间的倾斜总电子含量STECgro；
步骤C8，在观测数据中加入随机扰动δ，以模拟观测噪声的影响，计算方法如下：
Neobs(r)＝Neion(r)+δ(μion,σion) (7)
STECgobs＝STECgnss+δ(μgnss,σgnss) (8)
STECrobs＝STECgro+δ(μgro,σgro) (9)
其中，Neobs(r)表示位置r处垂测的模拟观测数据，STECgobs表示GNSS模拟观测数据，STECrobs表示掩星接收机模拟观测数据；(μion,σion)为地基垂测的电子密度的平均观测误差和标准差，(μgnss,σgnss)为地基GNSS接收机的TEC的平均观测误差和标准差，(μgro,σgro)为掩星观测TEC的平均观测误差和标准差；
步骤D，基于Kalman滤波的全球电离层数据同化：
步骤D1，将地基GNSS接收机、掩星接收机的电离层TEC和垂测仪的电子密度测量数据分别按照“点”型和“线”型观测数据进行分类；
步骤D2，对于“点”型电子密度观测数据，插值方法如下：
Figure BDA0003450928620000041
其中，Ne表示电离层电子密度，r表示电子密度所在的位置，K表示展开所取的级数；hk(r)为基函数，ak为插值系数，基函数计算方法如下：
Figure BDA0003450928620000042
步骤D3，对于“线”型TEC观测数据，根据接收机与卫星间的几何位置关系，TEC和同化状态参量之间的关系为：
Figure BDA0003450928620000043
其中，STEC表示电离层总电子含量，dr表示沿传播路径积分，积分上下限sat和rec分别为卫星和接收机位置，同化观测矩阵计算方法如下：
Figure BDA0003450928620000044
其中，ΔRin表示第i个TEC观测在第n个网格内的射线截距；
步骤D4，遍历同化窗口内所有电离层电子密度和TEC数据，根据式(11)和式(13)同化观测方程可构建为：
y＝Hx (14)
其中，y为垂测仪电离层电子密度、地基GNSS接收机TEC以及掩星接收机TEC数据序列组成的矢量，矩阵H为同化观测矩阵，x为待估计的各离散网格点的电子密度；
步骤D5，采用卡尔曼滤波同化算法进行地基和天基TEC或电子密度观测数据的同化，同化方程为：
xa＝xb+PHT(R+HPHT)-1(y-Hxb)＝xb+G(y-Hxb) (15)
其中，xb和xa分别表示同化前的电子密度背景场和同化后的电子密度分析场，P和R分别代表背景场误差协方差和观测数据的误差协方差，G为增益矩阵，H和y分别是同化观测矩阵和观测向量；
步骤D6，利用高斯—马尔科夫预报方法进行全球电离层预报，给出电离层TEC和电子密度预报结果，同时存储到指定的输出路径下：
利用高斯—马尔科夫方法对同化的电子密度进行预报，具体公式为：
Figure BDA0003450928620000051
其中，t表示预报时间，K为增益矩阵，L代表转换矩阵，为对角矩阵，具体计算方法为：
Figure BDA0003450928620000052
其中，ΔT为提前预报的时间间隔，τ表示电离层的时间相关尺度；
步骤E，电离层TEC和电子密度同化性能评估：
步骤E1，根据NeQuick模型，采用与步骤C1相同的输入，包括太阳辐射指数F10.7、月份和UT时刻，计算得到真实的三维电子密度和TEC地图值，标记为ParTruth；
步骤E2，根据国际参考电离层模型IRI，采用与步骤C1相同的输入，包括太阳辐射指数F10.7、月份和UT时刻，计算得到参考模型的三维电子密度和TEC地图值，标记为Parmod；
步骤E3，分别提取数据同化给出三维电子密度结果，利用式(12)积分计算垂直TEC地图结果，标记为ParDA；
步骤E4：计算数据同化的技术性能分，计算公式为：
SKS＝1-RMSEDA/RMSEmod (18)
RMSEDA＝RMSE(ParDA-ParTruth) (19)
RMSEmod＝RMSE(Parmod-ParTruth) (20)
其中，SKS表示数据同化的技术性能分，RMSE(﹒)表示计算样本的均方根误差；RMSEDA表示数据同化的均方根误差，RMSEmod表示背景模型的均方根误差；
步骤E5，按照电离层数据同化现报及提取预报结果，按照步骤E1—E4进行反复计算，以评估数据同化模型现报和预报的技术性能分；
步骤E6，采用不同的输入条件，重复步骤C至步骤D，生成不同的观测模拟数据并进行数据同化，以评估不同条件下观测系统的数据同化性能；
步骤E7，定量观测系统对数据同化的影响和贡献：当SKS＝1，则表明数据同化的结果与真实结果完全一致，表明此时的观测系统是理论上最佳的；当SKS＝0～1时，则表明数据同化的效果优于IRI模型的结果；SKS越大，则改进效果越好，观测系统的影响贡献是正面的；当SKS<0，则表明同化结果不如IRI模型，观测系统的影响贡献是负面的。
进一步的，所述步骤C6具体包括：
步骤C61，计算LEO卫星和GNSS卫星传播路径上的矢量夹角θ，计算方法如下：
Figure BDA0003450928620000061
其中，
Figure BDA0003450928620000062
表示LEO卫星的位置矢量，
Figure BDA0003450928620000063
Figure BDA0003450928620000064
表示GNSS卫星的位置矢量，
Figure BDA0003450928620000065
步骤C62，计算
Figure BDA0003450928620000066
在LEO卫星和GNSS卫星传播路径方向的投影截距
Figure BDA0003450928620000067
计算方法如下：
Figure BDA0003450928620000068
步骤C63，计算中间矢量，方法如下：
Figure BDA0003450928620000069
步骤C64，计算LEO卫星和GNSS卫星射线路径碰撞点坐标
Figure BDA00034509286200000610
Figure BDA00034509286200000611
步骤C65，根据碰撞点矢量信息，通过坐标转换即可计算得到碰撞点的高度，在真空近似下，电离层掩星的判断条件为：掩星碰撞点大于60km，小于LEO卫星高度，矢量夹角θ必须在掩星接收机可视角FOV范围内。
本发明的有益效果是：
本发明所公开的方法，对观测系统采用仿真试验手段(OSSE)，利用NeQuick模型模拟生成不同观测手段“真实”的观测数据，同时基于Kalman滤波同化算法对各种观测系统的观测数据进行同化处理，并采用总电子含量(TEC)和电子密度(Ne)的技术分(Skill Score)评估方法，对各类观测系统及其组合对全球电离层数据同化性能的影响进行定量评估。本发明方法可为全球电离层观测站网的设计和优化提供技术支撑。
附图说明
图1是本发明方法的流程示意图；
图2是LEO卫星—GNSS卫星射线传播路径的几何示意图；
图3是不同观测系统条件下的全球电离层数据同化结果对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
随着地基垂测仪、GNSS台站及天基掩星观测数量的日渐增多，如何选择合适的地基和天基观测作为观测资料来源，从而实现高精度的全球电离层数据同化和预报，是世界各国电离层观测系统设计和优化的一个重要研究方向。
实施例1，如图1所示，本实施例公开了一种观测系统对全球电离层数据同化性能的影响评估方法，包括如下步骤：
步骤A，地基和天基观测系统基本参数输入：
步骤A1，读取地基观测系统的配置参数，获取地基垂测、GNSS接收机的地理经纬高坐标和GNSS接收机的观测截止角等信息；
步骤A2，读取天基观测系统的配置参数，获取GNSS卫星和LEO掩星卫星的两行轨道TLE(Two Line Elements)星历、掩星接收机可视角FOV；
步骤B，GNSS卫星和LEO掩星卫星的轨道坐标计算：
步骤B1，解析步骤A2中得到的GNSS卫星和LEO掩星卫星的TLE星历，提取卫星编号、历元时刻、轨道倾角、历元时刻升交点赤经、轨道偏心率、历元时刻近地点幅角、历元时刻的平近点角、平均转速等参量；
步骤B2，利用SGP4(Simplified General Perturbation)模型解算GNSS和LEO掩星卫星的地心坐标系ECI下的经度、纬度和高度坐标；
步骤B3，将卫星地心坐标系下的经度、纬度和高度坐标转换为地心地固坐标系ECEF下的经度、纬度和高度坐标；
步骤C，地基和天基电离层观测数据仿真：
步骤C1，设定太阳辐射指数F10.7、月份和UT时刻作为电离层模型NeQuick的输入；
步骤C2，将地基垂测的地理经纬度输入到NeQuick模型中，计算对应输入条件下，地基垂测站100千米至F2层峰值高度以下高度上的电子密度剖面Neion(r)；
步骤C3，输入地基GNSS台站地心地固坐标系下的纬度、经度、高度
Figure BDA0003450928620000081
和GNSS卫星的纬度、经度、高度
Figure BDA0003450928620000082
计算出两者间的仰角E，计算方法如下：
Δλ＝λd-λs (1)
rd＝Re+hd (2)
rs＝Re+hs (3)
Figure BDA0003450928620000083
Figure BDA0003450928620000084
其中，Re表示地球半径，一般取值6371.2km；Δλ表示GNSS台站与GNSS卫星二者间经度的差值，L表示二者之间的大圆距离(弧度)，rd表示台站与地心之间的距离，rs表示卫星到地心的距离；
步骤C4，判断仰角是否大于GNSS接收机的观测截止角Eobs，若E大于Eobs，则将地基GNSS台站的经、纬、高度作为起始点，GNSS卫星轨道的经、纬、高度作为结束点坐标输入到NeQuick模型中，计算出两点间的倾斜总电子含量STECgnss；
步骤C5，将LEO掩星卫星的纬度、经度、高度
Figure BDA0003450928620000085
和GNSS卫星轨道的纬度、经度、高度
Figure BDA0003450928620000086
分别转化为地心地固坐标系ECEF下的坐标(Xt,Yt,Zt)和(Xs,Ys,Zs)，转换表达式表示为：
Figure BDA0003450928620000087
其中，
Figure BDA0003450928620000088
Re表示地球半径，e2＝0.00669437999013；
步骤C6，计算LEO卫星与LEO卫星—GNSS卫星射线之间的夹角θ和碰撞点高度，判定是否为掩星观测事件；
所述步骤C6具体包括：
步骤C61，计算LEO卫星和GNSS卫星传播路径上的矢量夹角θ，计算方法如下：
Figure BDA0003450928620000091
其中，
Figure BDA0003450928620000092
表示LEO卫星的位置矢量，
Figure BDA0003450928620000093
Figure BDA0003450928620000094
表示GNSS卫星的位置矢量，
Figure BDA0003450928620000095
其它位置矢量的定义如图2所示；
步骤C62，计算
Figure BDA0003450928620000096
在LEO卫星和GNSS卫星传播路径方向的投影截距
Figure BDA0003450928620000097
计算方法如下：
Figure BDA0003450928620000098
步骤C63，计算中间矢量，方法如下：
Figure BDA0003450928620000099
步骤C64，计算LEO卫星和GNSS卫星射线路径碰撞点坐标
Figure BDA00034509286200000910
Figure BDA00034509286200000911
步骤C65，根据碰撞点矢量信息，通过坐标转换即可计算得到碰撞点的高度，在真空近似下，电离层掩星的判断条件为：掩星碰撞点大于60km，小于LEO卫星高度，矢量夹角θ必须在掩星接收机可视角FOV范围内。
步骤C7，将满足掩星观测事件的LEO掩星卫星的经、纬、高度作为起始点，GNSS卫星轨道的经、纬、高度作为结束点坐标输入到NeQuick模型中，计算出两点间的倾斜总电子含量STECgro；
步骤C8，在观测数据中加入随机扰动δ，以模拟观测噪声的影响，计算方法如下：
Neobs(r)＝Neion(r)+δ(μion,σion) (7)
STECgobs＝STECgnss+δ(μgnss,σgnss) (8)
STECrobs＝STECgro+δ(μgro,σgro) (9)
其中，Neobs(r)表示位置r处垂测的模拟观测数据，STECgobs表示GNSS模拟观测数据，STECrobs表示掩星接收机模拟观测数据；(μion,σion)为地基垂测的电子密度的平均观测误差和标准差，(μgnss,σgnss)为地基GNSS接收机的TEC的平均观测误差和标准差，(μgro,σgro)为掩星观测TEC的平均观测误差和标准差；
步骤D，基于Kalman滤波的全球电离层数据同化：
步骤D1，将地基GNSS接收机、掩星接收机的电离层TEC和垂测仪的电子密度测量数据分别按照“点”型和“线”型观测数据进行分类；
步骤D2，对于“点”型电子密度观测数据，由于观测参量与同化待估计的状态参量相同，因此只需要利用简单的插值即可建立同化观测矩阵，插值方法如下：
Figure BDA0003450928620000101
其中，Ne表示电离层电子密度，r表示电子密度所在的位置，K表示展开所取的级数；hk(r)为基函数，ak为插值系数，该系数与插值算法的选取有关。基函数计算方法如下：
Figure BDA0003450928620000102
步骤D3，对于“线”型TEC观测数据，根据接收机与卫星间的几何位置关系，TEC和同化状态参量(电子密度)之间的关系为：
Figure BDA0003450928620000103
其中，STEC表示电离层总电子含量，dr表示沿传播路径积分，积分上下限sat和rec分别为卫星和接收机位置，同化观测矩阵计算方法如下：
Figure BDA0003450928620000104
其中，ΔRin表示第i个TEC观测在第n个网格内的射线截距；
步骤D4，遍历同化窗口内所有电离层电子密度和TEC数据，根据式(11)和式(13)同化观测方程可构建为：
y＝Hx (14)
其中，y为垂测仪电离层电子密度、地基GNSS接收机TEC以及掩星接收机TEC数据序列组成的矢量，矩阵H为同化观测矩阵，x为待估计的各离散网格点的电子密度；
步骤D5，采用卡尔曼滤波同化算法进行地基和天基TEC或电子密度观测数据的同化，同化方程为：
xa＝xb+PHT(R+HPHT)-1(y-Hxb)＝xb+G(y-Hxb) (15)
其中，xb和xa分别表示同化前的电子密度背景场和同化后的电子密度分析场，P和R分别代表背景场误差协方差和观测数据的误差协方差，G为增益矩阵，H和y分别是同化观测矩阵和观测向量；
步骤D6，利用高斯—马尔科夫(Gauss-Markov)预报方法进行全球电离层预报，给出电离层TEC和电子密度预报结果，同时存储到指定的输出路径下：
利用高斯—马尔科夫方法对同化的电子密度进行预报，具体公式为：
Figure BDA0003450928620000111
其中，t表示预报时间，K为增益矩阵，L代表转换矩阵(transition matrix)，其表现形式为对角矩阵，具体计算方法为：
Figure BDA0003450928620000112
其中，ΔT为提前预报的时间间隔，τ表示电离层的时间相关尺度，一般情况下取5h，磁暴期间取值可缩短为2h；
步骤E，电离层TEC和电子密度同化性能评估：
步骤E1，根据NeQuick模型，采用与步骤C1相同的输入，包括太阳辐射指数F10.7、月份和UT时刻，计算得到真实的三维电子密度和TEC地图值，标记为ParTruth；
步骤E2，根据国际参考电离层模型IRI，采用与步骤C1相同的输入，包括太阳辐射指数F10.7、月份和UT时刻，计算得到参考模型的三维电子密度和TEC地图值，标记为Parmod；
步骤E3，分别提取数据同化给出三维电子密度结果，利用式(12)积分计算垂直TEC地图结果，标记为ParDA；
步骤E4：计算数据同化的技术性能分(Skill Score)，计算公式为：
SKS＝1-RMSEDA/RMSEmod (18)
RMSEDA＝RMSE(ParDA-ParTruth) (19)
RMSEmod＝RMSE(Parmod-ParTruth) (20)
其中，SKS表示数据同化的技术性能分，RMSE(﹒)表示计算样本的均方根误差；RMSEDA表示数据同化的均方根误差，RMSEmod表示背景模型的均方根误差，电离层参量是TEC或电子密度；
步骤E5，按照电离层数据同化现报及提取预报结果，按照步骤E1—E4进行反复计算，以评估数据同化模型现报和预报的技术性能分；
步骤E6，采用不同的输入条件，重复步骤C至步骤D，生成不同的观测模拟数据并进行数据同化，以评估不同条件下观测系统的数据同化性能；
步骤E7，定量观测系统对数据同化的影响和贡献：当SKS＝1，则表明数据同化的结果与真实结果完全一致，表明此时的观测系统是理论上最佳的；当SKS＝0～1时，则表明数据同化的效果优于IRI模型的结果；SKS越大，则改进效果越好，观测系统的影响贡献是正面的；当SKS<0，则表明同化结果不如IRI模型，观测系统的影响贡献是负面的。
图3是一例不同观测系统条件下的全球电离层数据同化结果对比图。其中(a)为NeQuick模型模拟的真实电离层TEC分布；(b)为IRI模型给出的背景电离层TEC分布；(c)为869个地基GNSS接收机数据同化结果；(d)为12个COSMIC卫星数据同化结果；(e)为12颗COSMIC卫星与869个地基GNSS卫星联合观测条件下的数据同化结果；(f)为105颗掩星卫星观测条件下的数据同化结果。

Claims (2)

1.一种观测系统对全球电离层数据同化性能的影响评估方法，其特征在于，包括如下步骤：
步骤A，地基和天基观测系统基本参数输入：
步骤A1，读取地基观测系统的配置参数，获取地基垂测、GNSS接收机的地理经纬高坐标和GNSS接收机的观测截止角；
步骤A2，读取天基观测系统的配置参数，获取GNSS卫星和LEO掩星卫星的两行轨道TLE星历、掩星接收机可视角FOV；
步骤B，GNSS卫星和LEO掩星卫星的轨道坐标计算：
步骤B1，解析步骤A2中得到的GNSS卫星和LEO掩星卫星的TLE星历，提取卫星编号、历元时刻、轨道倾角、历元时刻升交点赤经、轨道偏心率、历元时刻近地点幅角、历元时刻的平近点角、平均转速；
步骤B2，利用SGP4模型解算GNSS和LEO掩星卫星的地心坐标系ECI下的经度、纬度和高度坐标；
步骤B3，将卫星地心坐标系下的经度、纬度和高度坐标转换为地心地固坐标系ECEF下的经度、纬度和高度坐标；
步骤C，地基和天基电离层观测数据仿真：
步骤C1，设定太阳辐射指数F10.7、月份和UT时刻作为电离层模型NeQuick的输入；
步骤C2，将地基垂测的地理经纬度输入到NeQuick模型中，计算对应输入条件下，地基垂测站100千米至F2层峰值高度以下高度上的电子密度剖面Neion(r)；
步骤C3，输入地基GNSS台站地心地固坐标系下的纬度、经度、高度
Figure FDA0003450928610000011
和GNSS卫星的纬度、经度、高度
Figure FDA0003450928610000012
计算出两者间的仰角E，计算方法如下：
Δλ＝λd-λs (1)
rd＝Re+hd (2)
rs＝Re+hs (3)
Figure FDA0003450928610000013
Figure FDA0003450928610000014
其中，Re表示地球半径，Δλ表示GNSS台站与GNSS卫星二者间经度的差值，L表示二者之间的大圆距离，rd表示台站与地心之间的距离，rs表示卫星到地心的距离；
步骤C4，判断仰角是否大于GNSS接收机的观测截止角Eobs，若E大于Eobs，则将地基GNSS台站的经、纬、高度作为起始点，GNSS卫星轨道的经、纬、高度作为结束点坐标输入到NeQuick模型中，计算出两点间的倾斜总电子含量STECgnss；
步骤C5，将LEO掩星卫星的纬度、经度、高度
Figure FDA0003450928610000021
和GNSS卫星轨道的纬度、经度、高度
Figure FDA0003450928610000022
分别转化为地心地固坐标系ECEF下的坐标(Xt,Yt,Zt)和(Xs,Ys,Zs)，转换表达式表示为：
Figure FDA0003450928610000023
其中，
Figure FDA0003450928610000024
Re表示地球半径，e2＝0.00669437999013；
步骤C6，计算LEO卫星与LEO卫星—GNSS卫星射线之间的夹角θ和碰撞点高度，判定是否为掩星观测事件；
步骤C7，将满足掩星观测事件的LEO掩星卫星的经、纬、高度作为起始点，GNSS卫星轨道的经、纬、高度作为结束点坐标输入到NeQuick模型中，计算出两点间的倾斜总电子含量STECgro；
步骤C8，在观测数据中加入随机扰动δ，以模拟观测噪声的影响，计算方法如下：
Neobs(r)＝Neion(r)+δ(μion,σion) (7)
STECgobs＝STECgnss+δ(μgnss,σgnss) (8)
STECrobs＝STECgro+δ(μgro,σgro) (9)
其中，Neobs(r)表示位置r处垂测的模拟观测数据，STECgobs表示GNSS模拟观测数据，STECrobs表示掩星接收机模拟观测数据；(μion,σion)为地基垂测的电子密度的平均观测误差和标准差，(μgnss,σgnss)为地基GNSS接收机的TEC的平均观测误差和标准差，(μgro,σgro)为掩星观测TEC的平均观测误差和标准差；
步骤D，基于Kalman滤波的全球电离层数据同化：
步骤D1，将地基GNSS接收机、掩星接收机的电离层TEC和垂测仪的电子密度测量数据分别按照“点”型和“线”型观测数据进行分类；
步骤D2，对于“点”型电子密度观测数据，插值方法如下：
Figure FDA0003450928610000031
其中，Ne表示电离层电子密度，r表示电子密度所在的位置，K表示展开所取的级数；hk(r)为基函数，ak为插值系数，基函数计算方法如下：
Figure FDA0003450928610000032
步骤D3，对于“线”型TEC观测数据，根据接收机与卫星间的几何位置关系，TEC和同化状态参量之间的关系为：
Figure FDA0003450928610000033
其中，STEC表示电离层总电子含量，dr表示沿传播路径积分，积分上下限sat和rec分别为卫星和接收机位置，同化观测矩阵计算方法如下：
Figure FDA0003450928610000034
其中，ΔRin表示第i个TEC观测在第n个网格内的射线截距；
步骤D4，遍历同化窗口内所有电离层电子密度和TEC数据，根据式(11)和式(13)同化观测方程可构建为：
y＝Hx (14)
其中，y为垂测仪电离层电子密度、地基GNSS接收机TEC以及掩星接收机TEC数据序列组成的矢量，矩阵H为同化观测矩阵，x为待估计的各离散网格点的电子密度；
步骤D5，采用卡尔曼滤波同化算法进行地基和天基TEC或电子密度观测数据的同化，同化方程为：
xa＝xb+PHT(R+HPHT)-1(y-Hxb)＝xb+G(y-Hxb) (15)
其中，xb和xa分别表示同化前的电子密度背景场和同化后的电子密度分析场，P和R分别代表背景场误差协方差和观测数据的误差协方差，G为增益矩阵，H和y分别是同化观测矩阵和观测向量；
步骤D6，利用高斯—马尔科夫预报方法进行全球电离层预报，给出电离层TEC和电子密度预报结果，同时存储到指定的输出路径下：
利用高斯—马尔科夫方法对同化的电子密度进行预报，具体公式为：
Figure FDA0003450928610000041
其中，t表示预报时间，L代表转换矩阵，为对角矩阵，具体计算方法为：
Figure FDA0003450928610000042
其中，ΔT为提前预报的时间间隔，τ表示电离层的时间相关尺度；
步骤E，电离层TEC和电子密度同化性能评估：
步骤E1，根据NeQuick模型，采用与步骤C1相同的输入，包括太阳辐射指数F10.7、月份和UT时刻，计算得到真实的三维电子密度和TEC地图值，标记为ParTruth；
步骤E2，根据国际参考电离层模型IRI，采用与步骤C1相同的输入，包括太阳辐射指数F10.7、月份和UT时刻，计算得到参考模型的三维电子密度和TEC地图值，标记为Parmod；
步骤E3，分别提取数据同化给出三维电子密度结果，利用式(12)积分计算垂直TEC地图结果，标记为ParDA；
步骤E4：计算数据同化的技术性能分，计算公式为：
SKS＝1-RMSEDA/RMSEmod (18)
RMSEDA＝RMSE(ParDA-ParTruth) (19)
RMSEmod＝RMSE(Parmod-ParTruth) (20)
其中，SKS表示数据同化的技术性能分，RMSE(﹒)表示计算样本的均方根误差；RMSEDA表示数据同化的均方根误差，RMSEmod表示背景模型的均方根误差；
步骤E5，按照电离层数据同化现报及提取预报结果，按照步骤E1—E4进行反复计算，以评估数据同化模型现报和预报的技术性能分；
步骤E6，采用不同的输入条件，重复步骤C至步骤D，生成不同的观测模拟数据并进行数据同化，以评估不同条件下观测系统的数据同化性能；
步骤E7，定量观测系统对数据同化的影响和贡献：当SKS＝1，则表明数据同化的结果与真实结果完全一致，表明此时的观测系统是理论上最佳的；当SKS＝0～1时，则表明数据同化的效果优于IRI模型的结果；SKS越大，则改进效果越好，观测系统的影响贡献是正面的；当SKS<0，则表明同化结果不如IRI模型，观测系统的影响贡献是负面的。
2.根据权利要求1所述观测系统对全球电离层数据同化性能的影响评估方法，其特征在于，所述步骤C6具体包括：
步骤C61，计算LEO卫星和GNSS卫星传播路径上的矢量夹角θ，计算方法如下：
Figure FDA0003450928610000051
其中，
Figure FDA0003450928610000052
表示LEO卫星的位置矢量，
Figure FDA0003450928610000053
Figure FDA0003450928610000054
表示GNSS卫星的位置矢量，
Figure FDA0003450928610000055
步骤C62，计算
Figure FDA0003450928610000056
在LEO卫星和GNSS卫星传播路径方向的投影截距
Figure FDA0003450928610000057
计算方法如下：
Figure FDA0003450928610000058
步骤C63，计算中间矢量，方法如下：
Figure FDA0003450928610000059
步骤C64，计算LEO卫星和GNSS卫星射线路径碰撞点坐标
Figure FDA00034509286100000510
Figure FDA00034509286100000511
步骤C65，根据碰撞点矢量信息，通过坐标转换即可计算得到碰撞点的高度，在真空近似下，电离层掩星的判断条件为：掩星碰撞点大于60km，小于LEO卫星高度，矢量夹角θ必须在掩星接收机可视角FOV范围内。
CN202111675121.2A 2021-12-31 2021-12-31 一种观测系统对全球电离层数据同化性能的影响评估方法 Active CN114417580B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111675121.2A CN114417580B (zh) 2021-12-31 2021-12-31 一种观测系统对全球电离层数据同化性能的影响评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111675121.2A CN114417580B (zh) 2021-12-31 2021-12-31 一种观测系统对全球电离层数据同化性能的影响评估方法

Publications (2)

Publication Number Publication Date
CN114417580A CN114417580A (zh) 2022-04-29
CN114417580B true CN114417580B (zh) 2022-12-02

Family

ID=81271473

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111675121.2A Active CN114417580B (zh) 2021-12-31 2021-12-31 一种观测系统对全球电离层数据同化性能的影响评估方法

Country Status (1)

Country Link
CN (1) CN114417580B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116170067B (zh) * 2023-04-25 2023-07-21 中国电子科技集团公司第五十四研究所 一种多星跟踪策略生成方法
CN117008154B (zh) * 2023-08-03 2024-03-26 北京航空航天大学 一种基于松弛因子逆时衰减函数的快速电离层层析方法
CN116956470B (zh) * 2023-09-11 2023-11-28 中国空气动力研究与发展中心超高速空气动力研究所 一种基于动态纵横比的leo航天器大气阻力算法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110275186B (zh) * 2019-07-11 2020-04-03 武汉大学 Leo卫星增强的gnss电离层归一化与融合建模方法
CN111123300B (zh) * 2020-01-13 2022-04-01 武汉大学 近实时大范围高精度电离层电子密度三维监测方法及装置
CN112649899B (zh) * 2020-11-19 2023-01-24 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种全球电离层数据同化和预报方法

Also Published As

Publication number Publication date
CN114417580A (zh) 2022-04-29

Similar Documents

Publication Publication Date Title
CN114417580B (zh) 一种观测系统对全球电离层数据同化性能的影响评估方法
CN107390233B (zh) 一种低轨卫星导航增强电离层延迟改正参数方法
Champollion et al. GPS water vapour tomography: preliminary results from the ESCOMPTE field experiment
Feltens et al. Comparative testing of four ionospheric models driven with GPS measurements
CN110020462B (zh) 一种对气象数据进行融合处理并生成数值天气预报的方法
Scherliess et al. Development of a physics-based reduced state Kalman filter for the ionosphere
RU2318222C2 (ru) Способ и система навигации в реальном масштабе времени, использующие три несущих радиосигнала, передаваемых спутником, и ионосферные коррекции
Isaksen et al. ERS scatterometer wind data impact on ECMWF's tropical cyclone forecasts
Kang et al. Development of an observation processing package for data assimilation in KIAPS
Mascitelli et al. Data assimilation of GPS-ZTD into the RAMS model through 3D-Var: preliminary results at the regional scale
CN112526617A (zh) 一种基于多源卫星信号的电离层层析成像系统观测数据模拟方法
CN112632473B (zh) 一种融合地基和空基gnss大气可降水量的计算方法
Seko et al. The meso-γ scale water vapor distribution associated with a thunderstorm calculated from a dense network of GPS receivers
CN115616637B (zh) 一种基于三维格网多径建模的城市复杂环境导航定位方法
Adavi et al. Analyzing different parameterization methods in GNSS tomography using the COST benchmark dataset
Norman et al. Simulating the impact of refractive transverse gradients resulting from a severe troposphere weather event on GPS signal propagation
CN112444825A (zh) 一种基于北斗geo的电离层tec同化建模方法
Izanlou et al. Enhanced Troposphere Tomography: Integration of GNSS and Remote Sensing Data with Optimal Vertical Constraints
CN114037901A (zh) 光伏发电预测导向的实时卫星近红外图像推算方法
Yu et al. Optimizing Global Navigation Satellite Systems network real-time kinematic infrastructure for homogeneous positioning performance from the perspective of tropospheric effects
Fu et al. An Introduction of GNSS Reflectometer Remote Sensing Mission From Yunyao Aerospace Technology Co., Ltd.
Ganguly et al. Real-time characterization of the ionosphere using diverse data and models
Ugur Modeling the neutral atmosphere in continuously operating GNSS networks using OPUS-Projects
Vankadara et al. Performance Analysis of Various Ionospheric Delay Corrections in Single-frequency GPS Positioning solution at a low latitude Indian location
Nilsson et al. GPS tomography using phase observations

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant