CN111125609A - 一种基于双指数驱动的电离层三维电子密度重构方法 - Google Patents

一种基于双指数驱动的电离层三维电子密度重构方法 Download PDF

Info

Publication number
CN111125609A
CN111125609A CN201911328951.0A CN201911328951A CN111125609A CN 111125609 A CN111125609 A CN 111125609A CN 201911328951 A CN201911328951 A CN 201911328951A CN 111125609 A CN111125609 A CN 111125609A
Authority
CN
China
Prior art keywords
layer
ionosphere
height
calculating
electron density
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.)
Granted
Application number
CN201911328951.0A
Other languages
English (en)
Other versions
CN111125609B (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 CN201911328951.0A priority Critical patent/CN111125609B/zh
Publication of CN111125609A publication Critical patent/CN111125609A/zh
Application granted granted Critical
Publication of CN111125609B publication Critical patent/CN111125609B/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
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • 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
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Data Mining & Analysis (AREA)
  • Signal Processing (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Molecular Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种基于双指数驱动的电离层三维电子密度重构方法，包括如下步骤：步骤A，多GNSS数据源电离层垂直总电子含量映射：步骤B，最优化电离层Rz指数确定：步骤C，最优化电离层IG指数确定：步骤D，基于双指数驱动的电离层三维电子密度重构。本发明所公开基于双指数驱动的电离层三维电子密度重构方法，采用国际全球卫星导航系统组织（International GNSS Service，IGS）发布的全球GNSS观测数据和COSMIC数据分析和档案中心（COSMIC Data Analysis and Archive Center，CDAAC）发布的GNSS掩星电子密度剖面数据产品作为电离层三维电子密度重构的数据来源，基于最优化的Rz和IG指数对国际参考电离层模型（International Reference Ionosphere，IRI）进行数据驱动更新，实现高精度电离层三维电子密度重构。

Description

一种基于双指数驱动的电离层三维电子密度重构方法
技术领域
本发明属于空间电离层特征参量现报预报领域，特别涉及该领域中的一种基于双指数驱动的电离层三维电子密度重构方法。
背景技术
电离层作为日地空间环境的重要组成部分，会对穿越其中的无线电波产生折射、反射、散射和吸收等效应，从而影响卫星导航、通信、雷达等诸多无线电信息系统的性能。利用各类技术手段来获取电离层的特征参量，研究电离层中的各种现象，揭示其内含的物理机制，具有重要的科学意义。在诸多电离层的特征参量中，电离层总电子含量(TotalElectron Content,TEC)和电子密度是其中最为关键的探测参数，尤其是电离层电子密度变化。目前，很多无线电信息系统如卫星导航及其增强系统，地基雷达折射修正、天基低频段SAR高精度目标成像的工程应用均需要三维高精度电离层电子密度信息支撑。掌握其时空变化特征和规律，对提升无线电信息系统性能具有重要的实用价值。
根据实际应用场景的不同，实现电子密度重构需要求解的未知系数不尽相同。但针对需要较高分辨率电子密度的应用领域，通常需要求解较大规模的未知系数，这将需要极大的内存和存储空间及非常大的CPU计算量，无法满足实际应用场景对快速的计算时效性、简便的嵌入式应用的需求。为降低计算量，提高运行时效以便于将算法嵌入至系统应用中，降低三维电子密度的重构时需要的待估未知数的数量是非常必要的。
当前，随着全球GPS、GLONASS、北斗、GALILEO等卫星导航系统的发展，以及星载掩星测量技术的进步，综合利用GNSS地基观测数据与GNSS天基掩星测量数据，基于数据驱动技术实现全球高精度电离层三维电子密度重构成为世界各国研究的热点方向。很多研究者开始尝试在各类地基/天基电离层探测的基础上，利用数据驱动技术对经验电离层模型的一些特定输入参数进行驱动更新，一方面降低计算量和存储空间，另一方面也可以较好的提升电离层电子密度重构的精度，以更好的满足工程应用需求。
发明内容
本发明所要解决的技术问题就是提供一种基于双指数驱动的电离层三维电子密度重构方法。
本发明采用如下技术方案：
一种基于双指数驱动的电离层三维电子密度重构方法，其改进之处在于，包括如下步骤：
步骤A，多GNSS数据源电离层垂直总电子含量映射：
步骤A1：根据指定日期，下载IGS发布的GNSS观测文件和GNSS卫星星历文件，对数据进行解压缩，提取观测台站的经度、纬度和高度，分别记为λr,
Figure BDA0002329090390000021
hr，提取GNSS卫星与观测台站间的伪距
Figure BDA0002329090390000022
和载波相位
Figure BDA0002329090390000023
其中i为对应的观测频段i＝1,2，j为对应的卫星编号，r为对应的接收机编号；
步骤A2：对于所有的伪距和载波相位观测值，分别计算无几何距离线性组合，计算方法如下：
Figure BDA0002329090390000024
Figure BDA0002329090390000025
其中：
Figure BDA0002329090390000026
为伪距的无几何距离线性组合，
Figure BDA0002329090390000027
为载波相位的无几何距离线性组合；
步骤A3：在同一观测弧段内，利用载波相位的无几何距离线性组合L4对伪距的无几何距离线性组合P4进行平滑滤波，计算方法表达为：
Figure BDA0002329090390000028
Figure BDA0002329090390000029
其中：
Figure BDA00023290903900000210
表示载波平滑伪距量，是计算
Figure BDA00023290903900000211
时的中间过程量，t表示历元时刻，ωt为t时刻的加权因子；当t＝1时，设定
Figure BDA00023290903900000212
步骤A4：利用GNSS星历文件，计算GNSS卫星的轨道坐标；同时计算卫星和接收机之间的仰角E和方位角A；
步骤A5：假定电离层总电子含量TEC集中在地球上空450km高度上的一个薄层内，计算观测站与GNSS卫星连线在该薄层高度上的穿刺点的地理坐标
Figure BDA00023290903900000213
λ，计算公式如下：
Figure BDA00023290903900000214
其中：
Figure BDA00023290903900000215
λ0为GNSS观测站的纬度和经度坐标，A为方位角；ψ0为观测站与GNSS卫星间的地心夹角，计算公式如下：
Figure BDA0002329090390000031
其中：Re＝6371.2km为地球的平均半径；hI＝450km为电离层薄层高度，E为接收机至卫星间的仰角；
步骤A6：计算倾斜总电子含量与垂直总电子含量的映射函数值，计算方法如下：
M(E)＝[1-(Re cos E/(Re+hI))2]-1/2
步骤A7：利用球谐函数构建电离层垂直总电子含量随空间分布模型，计算方法如下：
Figure BDA0002329090390000032
其中：nmax表示球谐函数展开的最大阶数，
Figure BDA0002329090390000033
为缔合勒让德函数，
Figure BDA0002329090390000034
为穿刺点的地理纬度，s为穿刺点在日固坐标系下的经度，
Figure BDA0002329090390000035
anm,bnm为球谐函数的未知系数；
步骤A8：利用平滑后的伪距的无几何距离线性组合
Figure BDA0002329090390000036
构建总电子含量映射求解方程，计算方法如下：
Figure BDA0002329090390000037
其中：f1为卫星的第一个频率，f2为卫星的第二个频率，DCBr为第r个接收机的硬件延迟，DCBj为第j颗卫星的硬件延迟值，为了对卫星和接收机的硬件延迟进行分离，同时加入以下四个约束条件：
Figure BDA0002329090390000038
其中：
Figure BDA0002329090390000039
表示第j颗GPS卫星的硬件延迟值，NG为GPS卫星的总数量；
Figure BDA00023290903900000310
表示第j颗GLONASS卫星的硬件延迟值，NR为GLONASS卫星的总数量；
Figure BDA00023290903900000311
表示第j颗北斗卫星的硬件延迟值，NC为北斗卫星的总数量；
Figure BDA0002329090390000041
表示第j颗GALILEO卫星的硬件延迟值，NE为GALILEO卫星的总数量；
步骤A9：选择1～2小时为时间窗口，遍历该时间窗口内所有的卫星与接收机，根据步骤A8的方程建立垂直总电子含量的线性矩阵如下：
Y＝AX
式中，Y是所有接收机与卫星间的观测量
Figure BDA0002329090390000042
的集合，A是与球谐函数、接收机和卫星硬件延迟相关的系数矩阵，X是待求解的未知系数矢量，
Figure BDA0002329090390000043
步骤A10：利用最小二乘算法求解拟合方程，解算得到未知系数X，计算方法如下：
XLS＝(ATA)-1ATY
其中：XLS是最小二乘解，其包含待估计的未知系数，包括球谐函数系数anm,bnm，接收机和卫星硬件延迟；
步骤A11：利用解算得到的系数anm,bnm，重新代入到球谐函数模型中，解算得到指定网格点处的电离层垂直TEC值，计算方法如下：
Figure BDA0002329090390000044
其中：
Figure BDA0002329090390000045
λg是网格点的纬度和经度坐标，
Figure BDA0002329090390000046
即为解算得到的垂直总电子含量，m表示累加序号，即m＝0,1,2…n；
步骤B，最优化电离层Rz指数确定：
步骤B1：根据与步骤A1相同的日期，下载CDAAC发布的GNSS掩星电子密度剖面数据，提取掩星碰撞点经度、纬度和高度处对应的电子密度值Ne及掩星发生时刻t，记为
Figure BDA0002329090390000047
步骤B2：对一次完整的掩星事件，提取该掩星事件对应的电离层F2层峰值高度及其时间和位置信息，记为
Figure BDA0002329090390000051
计算方法如下：
Figure BDA0002329090390000052
其中，ho表示观测点处的高度，
Figure BDA0002329090390000053
表示对应的掩星碰撞点的位置和时间信息，R表示一次掩星事件的所有采样点的集合；
步骤B3：对时间窗口内所有掩星事件的电离层F2层峰值高度数据进行质量控制，剔除异常的观测数据，符合要求的掩星峰值高度数据应满足以下公式：
Figure BDA0002329090390000054
其中：DET表示判决条件，
Figure BDA0002329090390000055
Δλ分别表示一次掩星事件期间起始时刻碰撞点纬度、经度与终止时刻碰撞点纬度、经度的绝对差值；
步骤B4：选择国际参考电离层IRI-2016模型为背景电离层模型，设定选择AMTB为IRI的峰值高度模型，根据提取的掩星事件对应的电离层F2层峰值高度及其时间和位置信息，调用电离层经验模型IRI，计算对应时间及地点处的电离层F2层峰值高度，记为
Figure BDA0002329090390000056
步骤B5：设定时间窗口为24小时，获取该时间窗口内所有掩星事件的电离层F2层峰值高度数据序列，计算两者的差值ΔhmF2,i，公式表达为：
Figure BDA0002329090390000057
其中：i表示为当前时间窗口内第i次掩星事件；λo,
Figure BDA0002329090390000058
ho,to表示F2层峰值高度对应的位置和时间信息；
步骤B6：设定Rz指数搜索范围，在该范围内不断调整电离层Rz指数，输出对应Rz指数条件下IRI模型计算的电离层F2层峰值高度值hmF2IRI，直到满足下式为止，即得到最优化Rz参量
Figure BDA0002329090390000059
搜索公式如下：
Figure BDA0002329090390000061
其中，N为设定的时间窗口内所有掩星事件的总数，Rzmin表示Rz搜索时取的最小值，Rzmax表示Rz搜索时取的最大值，设Rzmin＝-40.0，Rzmax＝400.0；
步骤B7：保存电离层
Figure BDA00023290903900000613
指数，以便后续处理过程读取；
步骤C，最优化电离层IG指数确定：
步骤C1：读取步骤A处理得到的电离层垂直TEC值及其对应的网格点坐标，设定
Figure BDA0002329090390000062
λgrid是网格点的纬度和经度坐标，
Figure BDA0002329090390000063
为GNSS观测数据解算得到的电离层垂直TEC；
步骤C2：设定国际参考电离层模型采用国际无线电联合会发布的URSI系数计算F2层临界频率，输入最优化
Figure BDA0002329090390000064
指数，将其作为模型的输入量，计算网格点
Figure BDA0002329090390000065
λgrid对应位置处的电子密度值
Figure BDA0002329090390000066
设定高度h的范围为100-1500km；
步骤C3：对100-1500km高度范围内的电子密度值沿着垂直方向进行积分，得到IRI模型输出的垂直TEC值
Figure BDA0002329090390000067
计算公式如下：
Figure BDA0002329090390000068
步骤C4：计算高度1500km至卫星高度20200km区间的等离子体层的电子密度值，计算方法如下：
Figure BDA0002329090390000069
其中：
Figure BDA00023290903900000610
表示等离子体层的电子密度，h'表示电离层顶部的高度，可近似等于1500km，H为电子密度的标高；
步骤C5：对1500-20200km高度范围内的电子密度值沿着垂直方向进行积分，得到磁层模型输出的垂直TEC值
Figure BDA00023290903900000611
计算公式如下：
Figure BDA00023290903900000612
其中：h1和h2分别表示积分高度的下限和上限；
步骤C6：选择对应的时间窗口，计算该时间窗口内GNSS计算的垂直TEC与IRI模型计算的垂直TEC的差值
Figure BDA0002329090390000071
公式表达为：
Figure BDA0002329090390000072
步骤C7：设定IG指数搜索范围，在该范围内不断调整电离层IG指数，按照步骤C3-C6，计算在输入不同IG指数条件下模型输出的电离层垂直TEC，直到满足下式为止，即得到最优化IG参量
Figure BDA0002329090390000073
搜索公式如下：
Figure BDA0002329090390000074
其中：IGmin表示IG指数搜索时取的最小值，IGmax表示IG指数搜索时取的最大值，设IGmin＝5.0，IGmax＝300.0；
步骤C8：保存所有网格点的电离层
Figure BDA0002329090390000075
指数，以便后续处理过程读取，以实现背景电离层模型IRI的驱动更新；
步骤D，基于双指数驱动的电离层三维电子密度重构：
步骤D1：将步骤B获取的Rz指数
Figure BDA0002329090390000076
输入到IRI模型中，设定IRI模型F2层峰高计算模型为AMTB，计算得到设定网格点
Figure BDA0002329090390000077
λgrid处的电离层峰值高度
Figure BDA0002329090390000078
计算方法如下：
Figure BDA0002329090390000079
其中：
Figure BDA00023290903900000710
表示展开系数，由模型默认输入文件提供，
Figure BDA00023290903900000711
为缔合勒让德函数，M为展开谐波数；
步骤D2：根据网格点的地理纬度和经度，获取网格点的磁倾角I，然后计算网格点所在位置处的修正磁倾角μ，计算公式如下：
Figure BDA00023290903900000712
其中：μ表示网格点所在位置的修正磁倾角，
Figure BDA0002329090390000081
为网格点的地理纬度；
步骤D3：计算地理坐标函数Gk，该函数计算方法：
Figure BDA0002329090390000082
其中：φ表示地理纬度，θ表示地理经度，q(k)为修正磁倾角的阶数，m(k)为经度的阶数，μ为修正磁倾角；
步骤D4：读取URSI系数文件，将步骤C中计算得到的最优化IG指数作为加权系数，按照线性组合的方式计算时间修正系数Ui,k，计算方法如下：
Figure BDA0002329090390000083
其中：
Figure BDA0002329090390000084
为URSI太阳活动低年系数，
Figure BDA0002329090390000085
为太阳活动高年系数，
Figure BDA0002329090390000086
为垂直TEC最优化得到的电离层IG指数；
步骤D5：设定IRI模型F2层峰高计算模型为URSI，计算得到设定网格点
Figure BDA0002329090390000087
λg处的电离层F2层临界频率foF2，计算方法如下：
Figure BDA0002329090390000088
其中：T表示全球时UT(-π≤T≤π)；H表示用于表征foF2日变化的谐波数的最大值，设为常数6；
步骤D6：根据计算得到的电离层峰值高度及临界频率，基于电离层电子密度剖面方程，通过输入不同的时间和地理位置，计算得到三维电子密度剖面并对数据进行输出。
进一步的，所述步骤A4具体包括：
步骤A41：将地基GNSS接收机和GNSS卫星经纬高坐标
Figure BDA0002329090390000089
λ,h转换为空间直角XYZ坐标，分别标记为X0,Y0,Z0和X1,Y1,Z1，转换表达式表示为：
Figure BDA0002329090390000091
其中
Figure BDA0002329090390000092
Re为地球半径，
Figure BDA0002329090390000093
步骤A42：将地基GNSS接收机和GNSS卫星的空间直角坐标转换为站心直角坐标N,E,U，转换表达式表示为：
Figure BDA0002329090390000094
其中：T为旋转矩阵，表示方法如下：
Figure BDA0002329090390000095
其中：λ0,
Figure BDA0002329090390000096
分别对应GNSS接收机的经度和纬度坐标，其单位为弧度；
步骤A43：计算接收机与卫星间的仰角E，计算表达式为：
Figure BDA0002329090390000097
步骤A44：计算接收机与卫星间的方位角A，计算表达式为：
Figure BDA0002329090390000098
进一步的，所述步骤D6具体包括：
步骤D61：将电离层F2层临界频率转换为F2层峰值电子密度值NmF2，转换公式如下：
Figure BDA0002329090390000101
其中：foF2为步骤D5计算得到的电离层F2层临界频率，其单位为Hz；
步骤D62：采用经验函数计算电离层F2层以上高度的电子密度值，计算表达式如下：
Ne1(h)＝NmF2×exp(-Y)
Figure BDA0002329090390000102
Figure BDA0002329090390000103
δ＝[η/(1+Z)-ζ/2]/(η×Z/β×(1+Z)2+ζ/400)
Figure BDA0002329090390000104
其中：Ne1(h)表示F2层以上高度的电离层电子密度，h为高度，ζ、β、η是与地磁纬度、太阳辐射通量及电离层foF2有关的经验参量；
步骤D63：计算电离层F1层与F2层之间高度上的电离层电子密度值Ne2(h)，计算表达式如下：
Figure BDA0002329090390000105
Figure BDA0002329090390000106
其中：B0为F2层的厚度参数，B1为F2层的形状参数；
步骤D64：计算电离层F1层高度区域上的电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000107
Figure BDA0002329090390000111
其中：Ne2(h)为F1层与F2层之间高度上的电离层电子密度值，即步骤D63中的计算公式，C1是F1层的型态参数，hmF1为F1层峰值高度，F1reg为判断F1层是否存在的逻辑参数，该参数由IRI模型给定；
步骤D65：计算电离层中间层高度区域的电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000112
Figure BDA0002329090390000113
T＝D2/(HZ-HEF-D)
其中：中间参数D＝HZ-HST、HZ＝(HST+hmF1)/2、HEF＝HME+HBR；而HST表示为Ne3(h)与E层峰值密度NmE相等时所对应的高度值，Ne3为F1层高度区域上的电离层电子密度值，即步骤D64采用的计算公式；HBR为电离层E层谷区的宽度；
步骤D66：计算电离层E层谷区的电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000114
其中：NmE表示E层峰值密度；Ereg为是否为晚上的逻辑参数，0表示为白天，1代表晚上；参数P由多项式函数计算得到，计算方法如下：
P＝E1x2+E2x3+E3x4+E4x5
x＝h-HmE
式中：HmE白天设定为110km，晚上设定为105km，Ei，i＝1,...,4分别为E层谷区的经验修正系数，由IRI模型给定；
步骤D67：计算电离层D层高度电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000115
其中：NmE为E层峰值密度；NmD为D层峰值密度；hmE为E层峰值高度；HDX为D层临界高度；D1、FP1、FP2、FP3作为经验参数由IRI模型给定；
步骤D68：将步骤D62-D67计算电离层得到的D、E、F1、F2层及F2层以上高度电离层电子密度值进行集合，得到全剖面的电子密度值。
本发明的有益效果是：
本发明所公开基于双指数驱动的电离层三维电子密度重构方法，采用国际全球卫星导航系统组织(International GNSS Service，IGS)发布的全球GNSS观测数据和COSMIC数据分析和档案中心(COSMIC Data Analysis and Archive Center，CDAAC)发布的GNSS掩星电子密度剖面数据产品作为电离层三维电子密度重构的数据来源，基于最优化的Rz和IG指数对国际参考电离层模型(International Reference Ionosphere，IRI)进行数据驱动更新，实现高精度电离层三维电子密度重构。本发明所公开的方法具有数据兼容性好、计算结果稳定可靠、计算速度快、电子密度重构精度高等特点，相比传统的利用气候学电离层模型给出的电离层电子密度结果，本发明所公开的方法能够极大的提高国际参考电离层模型的电离层电子密度计算精度，满足空间科学研究、空间天气态势感知、极端气候条件下(如发生地磁暴、太阳耀斑等状况下)电离层时空变化诊断与分析、卫星通信、卫星导航系统高精度电离层折射误差修正等多个领域对高精度电离层三维电子密度信息的需求。
附图说明
图1是本发明实施例1所公开方法的流程示意图；
图2是本发明实施例1所公开方法中多GNSS数据源电离层垂直总电子含量的映射结果示意图；
图3是本发明实施例1所公开方法中最优化电离层Rz指数的计算结果示意图；
图4是本发明实施例1所公开方法中最优化电离层IG指数的计算结果示意图；
图5是本发明实施例1所公开方法中电离层三维电子密度重构的结果示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
实施例1，本实施例在对实测数据电离层驱动更新技术进行深入研究的基础上，提出了一种基于双指数驱动的电离层三维电子密度重构方法，该方法以国际参考电离层模型IRI为背景模型，利用地基GNSS及天基掩星数据，采用电离层IG指数和Rz指数联合的双指数驱动方法，实现了电离层三维电子密度高精度重构。解决了多GNSS数据源电离层垂直总电子含量映射、最优化电离层Rz指数确定、最优化电离层IG指数确定、基于双指数驱动的电离层三维电子密度重构等关键问题，为实现高精度电离层电子密度信息获取提供了一整套解决方法。
具体的说，如图1所示，本实施例公开了一种基于双指数驱动的电离层三维电子密度重构方法，包括如下步骤：
步骤A，多GNSS(Global Navigation Satellite System)数据源电离层垂直总电子含量映射：
步骤A1：根据指定日期，下载IGS发布的GNSS观测文件和GNSS卫星星历文件，对数据进行解压缩，提取观测台站的经度、纬度和高度，分别记为λr,
Figure BDA0002329090390000131
hr，提取GNSS卫星与观测台站间的伪距
Figure BDA0002329090390000132
和载波相位
Figure BDA0002329090390000133
其中i为对应的观测频段i＝1,2，j为对应的卫星编号，r为对应的接收机编号；
步骤A2：对于所有的伪距和载波相位观测值，分别计算无几何距离线性组合(Geometry-Free Combination)，计算方法如下：
Figure BDA0002329090390000134
Figure BDA0002329090390000135
其中：
Figure BDA0002329090390000136
为伪距的无几何距离线性组合，
Figure BDA0002329090390000137
为载波相位的无几何距离线性组合；
步骤A3：在同一观测弧段内，利用载波相位的无几何距离线性组合L4对伪距的无几何距离线性组合P4进行平滑滤波，计算方法表达为：
Figure BDA0002329090390000138
Figure BDA0002329090390000139
其中：
Figure BDA00023290903900001310
表示载波平滑伪距量，是计算
Figure BDA00023290903900001311
时的中间过程量，t表示历元时刻，ωt为t时刻的加权因子；当t＝1时，设定
Figure BDA00023290903900001312
步骤A4：根据接收机坐标与利用GNSS星历文件计算得到的GNSS卫星轨道坐标，计算卫星和接收机之间的仰角E和方位角A；
步骤A4具体包括：
步骤A41：将地基GNSS接收机和GNSS卫星经纬高坐标
Figure BDA0002329090390000141
λ,h转换为空间直角XYZ坐标，分别标记为X0,Y0,Z0和X1,Y1,Z1，转换表达式表示为：
Figure BDA0002329090390000142
其中
Figure BDA0002329090390000143
Re为地球半径，
Figure BDA0002329090390000144
步骤A42：将地基GNSS接收机和GNSS卫星的空间直角坐标转换为站心直角坐标N,E,U，转换表达式表示为：
Figure BDA0002329090390000145
其中：T为旋转矩阵，表示方法如下：
Figure BDA0002329090390000146
其中：λ0,
Figure BDA0002329090390000147
分别对应GNSS接收机的经度和纬度坐标，其单位为弧度；
步骤A43：计算接收机与卫星间的仰角E，计算表达式为：
Figure BDA0002329090390000148
步骤A44：计算接收机与卫星间的方位角A，计算表达式为：
Figure BDA0002329090390000151
步骤A5：一般假定电离层总电子含量TEC集中在地球上空450km高度上的一个薄层内，计算观测站与GNSS卫星连线在该薄层高度上的穿刺点的地理坐标(
Figure BDA0002329090390000152
λ，计算公式如下：
Figure BDA0002329090390000153
其中：
Figure BDA0002329090390000154
λ0)为GNSS观测站的纬度和经度坐标，A为方位角；ψ0为观测站与GNSS卫星间的地心夹角，计算公式如下：
Figure BDA0002329090390000155
其中：Re＝6371.2km为地球的平均半径；hI＝450km为电离层薄层高度，E为接收机至卫星间的仰角；
步骤A6：计算倾斜总电子含量与垂直总电子含量的映射函数值，计算方法如下：
M(E)＝[1-(Re cos E/(Re+hI))2]-1/2
步骤A7：利用球谐函数构建电离层垂直总电子含量随空间分布模型，计算方法如下：
Figure BDA0002329090390000156
其中：nmax表示球谐函数展开的最大阶数，
Figure BDA0002329090390000157
为缔合勒让德(AssociatedLegendre)函数，
Figure BDA0002329090390000158
为穿刺点的地理纬度，s为穿刺点在日固坐标系下的经度，
Figure BDA0002329090390000159
anm,bnm为球谐函数的未知系数；
步骤A8：利用平滑后的伪距的无几何距离线性组合
Figure BDA0002329090390000161
构建总电子含量映射求解方程，计算方法如下：
Figure BDA0002329090390000162
其中：f1为卫星的第一个频率，f2为卫星的第二个频率，DCBr为第r个接收机的硬件延迟，DCBj为第j颗卫星的硬件延迟值，为了对卫星和接收机的硬件延迟进行分离，同时加入以下四个约束条件：
Figure BDA0002329090390000163
其中：
Figure BDA0002329090390000164
表示第j颗GPS卫星的硬件延迟值，NG为GPS卫星的总数量；
Figure BDA0002329090390000165
表示第j颗GLONASS卫星的硬件延迟值，NR为GLONASS卫星的总数量；
Figure BDA0002329090390000166
表示第j颗北斗卫星的硬件延迟值，NC为北斗卫星的总数量；
Figure BDA0002329090390000167
表示第j颗GALILEO卫星的硬件延迟值，NE为GALILEO卫星的总数量；
步骤A9：选择1～2小时为时间窗口，遍历该时间窗口内所有的卫星与接收机，根据步骤A8的方程建立垂直总电子含量的线性矩阵如下：
Y＝AX
式中，Y是所有接收机与卫星间的观测量
Figure BDA0002329090390000168
的集合，A是与球谐函数、接收机和卫星硬件延迟相关的系数矩阵，X是待求解的未知系数矢量，
Figure BDA0002329090390000169
步骤A10：利用最小二乘算法求解拟合方程，解算得到未知系数X，计算方法如下：
XLS＝(ATA)-1ATY
其中：XLS是最小二乘解，其包含待估计的未知系数，包括球谐函数系数anm,bnm，接收机和卫星硬件延迟等；
步骤A11：利用解算得到的系数anm,bnm，重新代入到球谐函数模型中，解算得到指定网格点处的电离层垂直TEC值，计算方法如下：
Figure BDA0002329090390000171
其中：
Figure BDA0002329090390000172
λg是网格点的纬度和经度坐标，
Figure BDA0002329090390000173
即为解算得到的垂直总电子含量，m表示累加序号，即m＝0,1,2…n；
多GNSS数据源电离层垂直总电子含量的映射结果如图2所示。
步骤B，最优化电离层Rz指数确定：
步骤B1：根据与步骤A1相同的日期，下载CDAAC发布的GNSS掩星电子密度剖面数据，提取掩星碰撞点(经度、纬度和高度)处对应的电子密度值Ne及掩星发生时刻t，记为
Figure BDA0002329090390000174
步骤B2：对一次完整的掩星事件，提取该掩星事件对应的电离层F2层峰值高度及其时间和位置信息，记为
Figure BDA0002329090390000175
计算方法如下：
Figure BDA0002329090390000176
其中，ho表示观测点处的高度，
Figure BDA0002329090390000177
表示对应的掩星碰撞点的位置和时间信息，R表示一次掩星事件的所有采样点的集合；
步骤B3：对时间窗口内所有掩星事件的电离层F2层峰值高度数据进行质量控制，剔除异常的观测数据，符合要求的掩星峰值高度数据应满足以下公式：
Figure BDA0002329090390000178
其中：DET表示判决条件，
Figure BDA0002329090390000179
Δλ分别表示一次掩星事件期间起始时刻碰撞点纬度、经度与终止时刻碰撞点纬度、经度的绝对差值；
步骤B4：选择国际参考电离层IRI-2016模型为背景电离层模型，设定选择AMTB为IRI的峰值高度模型(AMTB为一种电离层hmF2模型，分别代表四个作者的首字母(Altadill,D.,S.Magdaleno,J.M.Torta,and E.Blanch))，根据提取的掩星事件对应的电离层F2层峰值高度及其时间和位置信息，调用电离层经验模型IRI，计算对应时间及地点处的电离层F2层峰值高度，记为
Figure BDA0002329090390000181
步骤B5：设定时间窗口为24小时，获取该时间窗口内所有掩星事件的电离层F2层峰值高度数据序列，计算两者的差值ΔhmF2,i，公式表达为：
Figure BDA0002329090390000182
其中：i表示为当前时间窗口内第i次掩星事件；λo,
Figure BDA0002329090390000183
ho,to表示F2层峰值高度对应的位置和时间信息；
步骤B6：设定Rz指数搜索范围，在该范围内不断调整电离层Rz指数，输出对应Rz指数条件下IRI模型计算的电离层F2层峰值高度值hmF2IRI，直到满足下式为止，即得到最优化Rz参量
Figure BDA0002329090390000184
搜索公式如下：
Figure BDA0002329090390000185
其中，N为设定的时间窗口内所有掩星事件的总数，Rzmin表示Rz搜索时取的最小值，Rzmax表示Rz搜索时取的最大值，设Rzmin＝-40.0，Rzmax＝400.0；
步骤B7：保存电离层
Figure BDA0002329090390000188
指数，以便后续处理过程读取；
最优化电离层Rz指数的计算结果如图3所示。
步骤C，最优化电离层IG指数确定：
步骤C1：读取步骤A处理得到的电离层垂直TEC值及其对应的网格点坐标，设定
Figure BDA0002329090390000186
λgrid是网格点的纬度和经度坐标，
Figure BDA0002329090390000187
为GNSS观测数据解算得到的电离层垂直TEC；
步骤C2：设定国际参考电离层模型采用国际无线电联合会(International Unionof Radio Science,URSI)发布的URSI系数计算F2层临界频率，输入最优化
Figure BDA0002329090390000191
指数，将其作为模型的输入量，计算网格点
Figure BDA0002329090390000192
λgrid对应位置处的电子密度值
Figure BDA0002329090390000193
设定高度h的范围为100-1500km；
步骤C3：对100-1500km高度范围内的电子密度值沿着垂直方向进行积分，得到IRI模型输出的垂直TEC值
Figure BDA0002329090390000194
计算公式如下：
Figure BDA0002329090390000195
步骤C4：计算高度1500km至卫星高度20200km区间的等离子体层的电子密度值，计算方法如下：
Figure BDA0002329090390000196
其中：
Figure BDA0002329090390000197
表示等离子体层的电子密度，h'表示电离层顶部的高度，可近似等于1500km，H为电子密度的标高(Scale Height)；
步骤C5：对1500-20200km高度范围内的电子密度值沿着垂直方向进行积分，得到磁层模型输出的垂直TEC值
Figure BDA0002329090390000198
计算公式如下：
Figure BDA0002329090390000199
其中：h1和h2分别表示积分高度的下限和上限；
步骤C6：选择对应的时间窗口，计算该时间窗口内GNSS计算的垂直TEC与IRI模型计算的垂直TEC的差值
Figure BDA00023290903900001910
公式表达为：
Figure BDA00023290903900001911
步骤C7：设定IG指数搜索范围，在该范围内不断调整电离层IG指数，按照步骤C3-C6，计算在输入不同IG指数条件下模型输出的电离层垂直TEC，直到满足下式为止，即得到最优化IG参量
Figure BDA00023290903900001912
搜索公式如下：
Figure BDA00023290903900001913
其中：IGmin表示IG指数搜索时取的最小值，IGmax表示IG指数搜索时取的最大值，设IGmin＝5.0，IGmax＝300.0；
步骤C8：保存所有网格点的电离层
Figure BDA0002329090390000201
指数，以便后续处理过程读取，以实现背景电离层模型IRI的驱动更新；
最优化电离层IG指数的计算结果如图4所示。
步骤D，基于双指数驱动的电离层三维电子密度重构：
步骤D1：将步骤B获取的Rz指数
Figure BDA0002329090390000202
输入到IRI模型中，设定IRI模型F2层峰高计算模型为AMTB，计算得到设定网格点
Figure BDA0002329090390000203
λgrid处的电离层峰值高度
Figure BDA0002329090390000204
计算方法如下：
Figure BDA0002329090390000205
其中：
Figure BDA0002329090390000206
表示展开系数，由模型默认输入文件提供，
Figure BDA0002329090390000207
为缔合勒让德函数，M为展开谐波数；
步骤D2：根据网格点的地理纬度和经度，获取网格点的磁倾角I，然后计算网格点所在位置处的修正磁倾角μ，计算公式如下：
Figure BDA0002329090390000208
其中：μ表示网格点所在位置的修正磁倾角，
Figure BDA0002329090390000209
为网格点的地理纬度；
步骤D3：计算地理坐标函数Gk，该函数计算方法：
Figure BDA00023290903900002010
其中：φ表示地理纬度，θ表示地理经度，q(k)为修正磁倾角的阶数，m(k)为经度的阶数，μ为修正磁倾角；
步骤D4：读取URSI系数文件，将步骤C中计算得到的最优化IG指数作为加权系数，按照线性组合的方式计算时间修正系数Ui,k，计算方法如下：
Figure BDA0002329090390000211
其中：
Figure BDA0002329090390000212
为URSI太阳活动低年系数，
Figure BDA0002329090390000213
为太阳活动高年系数，
Figure BDA0002329090390000214
为垂直TEC最优化得到的电离层IG指数；
步骤D5：设定IRI模型F2层峰高计算模型为URSI，计算得到设定网格点
Figure BDA0002329090390000215
λg处的电离层F2层临界频率foF2，计算方法如下：
Figure BDA0002329090390000216
其中：T表示全球时UT(-π≤T≤π)；H表示用于表征foF2日变化的谐波数(harmonics)的最大值，设为常数6；
步骤D6：根据计算得到的电离层峰值高度及临界频率，基于电离层电子密度剖面方程，通过输入不同的时间和地理位置，计算得到三维电子密度剖面并对数据进行输出。
步骤D6具体包括：
步骤D61：将电离层F2层临界频率转换为F2层峰值电子密度值NmF2，转换公式如下：
Figure BDA0002329090390000217
其中：foF2为步骤D5计算得到的电离层F2层临界频率，其单位为Hz；
步骤D62：采用经验函数计算电离层F2层以上高度的电子密度值，计算表达式如下：
Ne1(h)＝NmF2×exp(-Y)
Figure BDA0002329090390000218
Figure BDA0002329090390000219
δ＝[η/(1+Z)-ζ/2]/(η×Z/β×(1+Z)2+ζ/400)
Figure BDA0002329090390000221
其中：Ne1(h)表示F2层以上高度的电离层电子密度，h为高度，ζ、β、η是与地磁纬度、太阳辐射通量及电离层foF2有关的经验参量；
步骤D63：计算电离层F1层与F2层之间高度上的电离层电子密度值Ne2(h)，计算表达式如下：
Figure BDA0002329090390000222
Figure BDA0002329090390000223
其中：B0为F2层的厚度参数(bottom-side thickness parameter)，B1为F2层的形状参数(shape parameter)；
步骤D64：计算电离层F1层高度区域上的电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000224
Figure BDA0002329090390000225
其中：Ne2(h)为F1层与F2层之间高度上的电离层电子密度值，即步骤D63中的计算公式，C1是F1层的型态参数(F1 shape parameter)，hmF1为F1层峰值高度，F1reg为判断F1层是否存在的逻辑参数，该参数由IRI模型给定；
步骤D65：计算电离层中间层高度区域的电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000226
Figure BDA0002329090390000227
T＝D2/(HZ-HEF-D)
其中：中间参数D＝HZ-HST、HZ＝(HST+hmF1)/2、HEF＝HME+HBR；而HST表示为Ne3(h)与E层峰值密度NmE相等时所对应的高度值，Ne3为F1层高度区域上的电离层电子密度值，即步骤D64采用的计算公式；HBR为电离层E层谷区(E-valley region)的宽度；
步骤D66：计算电离层E层谷区的电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000231
其中：NmE表示E层峰值密度；Ereg为是否为晚上的逻辑参数，0表示为白天，1代表晚上；参数P由多项式函数计算得到，计算方法如下：
P＝E1x2+E2x3+E3x4+E4x5
x＝h-HmE
式中：HmE白天设定为110km，晚上设定为105km，Ei，i＝1,...,4分别为E层谷区的经验修正系数，由IRI模型给定；
步骤D67：计算电离层D层高度电离层电子密度值，计算表达式如下：
Figure BDA0002329090390000232
其中：NmE为E层峰值密度；NmD为D层峰值密度；hmE为E层峰值高度；HDX为D层临界高度；D1、FP1、FP2、FP3作为经验参数由IRI模型给定；
步骤D68：将步骤D62-D67计算电离层得到的D、E、F1、F2层及F2层以上高度电离层电子密度值进行集合，得到全剖面的电子密度值。
电离层三维电子密度重构的结果如图5所示。

Claims (3)

1.一种基于双指数驱动的电离层三维电子密度重构方法，其特征在于，包括如下步骤：
步骤A，多GNSS数据源电离层垂直总电子含量映射：
步骤A1：根据指定日期，下载IGS发布的GNSS观测文件和GNSS卫星星历文件，对数据进行解压缩，提取观测台站的经度、纬度和高度，分别记为λr,
Figure FDA0002329090380000011
hr，提取GNSS卫星与观测台站间的伪距
Figure FDA0002329090380000012
和载波相位
Figure FDA0002329090380000013
其中i为对应的观测频段i＝1,2，j为对应的卫星编号，r为对应的接收机编号；
步骤A2：对于所有的伪距和载波相位观测值，分别计算无几何距离线性组合，计算方法如下：
Figure FDA0002329090380000014
Figure FDA0002329090380000015
其中：
Figure FDA0002329090380000016
为伪距的无几何距离线性组合，
Figure FDA0002329090380000017
为载波相位的无几何距离线性组合；
步骤A3：在同一观测弧段内，利用载波相位的无几何距离线性组合L4对伪距的无几何距离线性组合P4进行平滑滤波，计算方法表达为：
Figure FDA0002329090380000018
Figure FDA0002329090380000019
其中：
Figure FDA00023290903800000110
表示载波平滑伪距量，是计算
Figure FDA00023290903800000111
时的中间过程量，t表示历元时刻，ωt为t时刻的加权因子；当t＝1时，设定
Figure FDA00023290903800000112
步骤A4：利用GNSS星历文件，计算GNSS卫星的轨道坐标；同时计算卫星和接收机之间的仰角E和方位角A；
步骤A5：假定电离层总电子含量TEC集中在地球上空450km高度上的一个薄层内，计算观测站与GNSS卫星连线在该薄层高度上的穿刺点的地理坐标(
Figure FDA00023290903800000113
λ，计算公式如下：
Figure FDA00023290903800000114
其中：
Figure FDA00023290903800000115
λ0)为GNSS观测站的纬度和经度坐标，A为方位角；ψ0为观测站与GNSS卫星间的地心夹角，计算公式如下：
Figure FDA0002329090380000021
其中：Re＝6371.2km为地球的平均半径；hI＝450km为电离层薄层高度，E为接收机至卫星间的仰角；
步骤A6：计算倾斜总电子含量与垂直总电子含量的映射函数值，计算方法如下：
M(E)＝[1-(RecosE/(Re+hI))2]-1/2
步骤A7：利用球谐函数构建电离层垂直总电子含量随空间分布模型，计算方法如下：
Figure FDA0002329090380000022
其中：nmax表示球谐函数展开的最大阶数，
Figure FDA0002329090380000023
为缔合勒让德函数，
Figure FDA0002329090380000024
为穿刺点的地理纬度，s为穿刺点在日固坐标系下的经度，
Figure FDA0002329090380000025
anm,bnm为球谐函数的未知系数；
步骤A8：利用平滑后的伪距的无几何距离线性组合
Figure FDA0002329090380000026
构建总电子含量映射求解方程，计算方法如下：
Figure FDA0002329090380000027
其中：f1为卫星的第一个频率，f2为卫星的第二个频率，DCBr为第r个接收机的硬件延迟，DCBj为第j颗卫星的硬件延迟值，为了对卫星和接收机的硬件延迟进行分离，同时加入以下四个约束条件：
Figure FDA0002329090380000028
其中：
Figure FDA0002329090380000029
表示第j颗GPS卫星的硬件延迟值，NG为GPS卫星的总数量；
Figure FDA00023290903800000210
表示第j颗GLONASS卫星的硬件延迟值，NR为GLONASS卫星的总数量；
Figure FDA0002329090380000031
表示第j颗北斗卫星的硬件延迟值，NC为北斗卫星的总数量；
Figure FDA0002329090380000032
表示第j颗GALILEO卫星的硬件延迟值，NE为GALILEO卫星的总数量；
步骤A9：选择1～2小时为时间窗口，遍历该时间窗口内所有的卫星与接收机，根据步骤A8的方程建立垂直总电子含量的线性矩阵如下：
Y＝AX
式中，Y是所有接收机与卫星间的观测量
Figure FDA0002329090380000033
的集合，A是与球谐函数、接收机和卫星硬件延迟相关的系数矩阵，X是待求解的未知系数矢量，
Figure FDA0002329090380000034
步骤A10：利用最小二乘算法求解拟合方程，解算得到未知系数X，计算方法如下：
XLS＝(ATA)-1ATY
其中：XLS是最小二乘解，其包含待估计的未知系数，包括球谐函数系数anm,bnm，接收机和卫星硬件延迟；
步骤A11：利用解算得到的系数anm,bnm，重新代入到球谐函数模型中，解算得到指定网格点处的电离层垂直TEC值，计算方法如下：
Figure FDA0002329090380000035
其中：
Figure FDA0002329090380000036
λg是网格点的纬度和经度坐标，
Figure FDA0002329090380000037
即为解算得到的垂直总电子含量，m表示累加序号，即m＝0,1,2…n；
步骤B，最优化电离层Rz指数确定：
步骤B1：根据与步骤A1相同的日期，下载CDAAC发布的GNSS掩星电子密度剖面数据，提取掩星碰撞点经度、纬度和高度处对应的电子密度值Ne及掩星发生时刻t，记为
Figure FDA0002329090380000038
步骤B2：对一次完整的掩星事件，提取该掩星事件对应的电离层F2层峰值高度及其时间和位置信息，记为
Figure FDA0002329090380000041
计算方法如下：
Figure FDA0002329090380000042
其中，ho表示观测点处的高度，
Figure FDA0002329090380000043
表示对应的掩星碰撞点的位置和时间信息，R表示一次掩星事件的所有采样点的集合；
步骤B3：对时间窗口内所有掩星事件的电离层F2层峰值高度数据进行质量控制，剔除异常的观测数据，符合要求的掩星峰值高度数据应满足以下公式：
Figure FDA0002329090380000044
其中：DET表示判决条件，
Figure FDA0002329090380000045
Δλ分别表示一次掩星事件期间起始时刻碰撞点纬度、经度与终止时刻碰撞点纬度、经度的绝对差值；
步骤B4：选择国际参考电离层IRI-2016模型为背景电离层模型，设定选择AMTB为IRI的峰值高度模型，根据提取的掩星事件对应的电离层F2层峰值高度及其时间和位置信息，调用电离层经验模型IRI，计算对应时间及地点处的电离层F2层峰值高度，记为
Figure FDA0002329090380000046
步骤B5：设定时间窗口为24小时，获取该时间窗口内所有掩星事件的电离层F2层峰值高度数据序列，计算两者的差值ΔhmF2,i，公式表达为：
Figure FDA0002329090380000047
其中：i表示为当前时间窗口内第i次掩星事件；λo,
Figure FDA0002329090380000048
ho,to表示F2层峰值高度对应的位置和时间信息；
步骤B6：设定Rz指数搜索范围，在该范围内不断调整电离层Rz指数，输出对应Rz指数条件下IRI模型计算的电离层F2层峰值高度值hmF2IRI，直到满足下式为止，即得到最优化Rz参量
Figure FDA0002329090380000049
搜索公式如下：
Figure FDA0002329090380000051
其中，N为设定的时间窗口内所有掩星事件的总数，Rzmin表示Rz搜索时取的最小值，Rzmax表示Rz搜索时取的最大值，设Rzmin＝-40.0，Rzmax＝400.0；
步骤B7：保存电离层
Figure FDA0002329090380000052
指数，以便后续处理过程读取；
步骤C，最优化电离层IG指数确定：
步骤C1：读取步骤A处理得到的电离层垂直TEC值及其对应的网格点坐标，设定
Figure FDA0002329090380000053
λgrid是网格点的纬度和经度坐标，
Figure FDA0002329090380000054
为GNSS观测数据解算得到的电离层垂直TEC；
步骤C2：设定国际参考电离层模型采用国际无线电联合会发布的URSI系数计算F2层临界频率，输入最优化
Figure FDA0002329090380000055
指数，将其作为模型的输入量，计算网格点
Figure FDA0002329090380000056
λgrid对应位置处的电子密度值
Figure FDA0002329090380000057
设定高度h的范围为100-1500km；
步骤C3：对100-1500km高度范围内的电子密度值沿着垂直方向进行积分，得到IRI模型输出的垂直TEC值
Figure FDA0002329090380000058
计算公式如下：
Figure FDA0002329090380000059
步骤C4：计算高度1500km至卫星高度20200km区间的等离子体层的电子密度值，计算方法如下：
Figure FDA00023290903800000510
其中：
Figure FDA00023290903800000511
表示等离子体层的电子密度，h'表示电离层顶部的高度，可近似等于1500km，H为电子密度的标高；
步骤C5：对1500-20200km高度范围内的电子密度值沿着垂直方向进行积分，得到磁层模型输出的垂直TEC值
Figure FDA00023290903800000512
计算公式如下：
Figure FDA00023290903800000513
其中：h1和h2分别表示积分高度的下限和上限；
步骤C6：选择对应的时间窗口，计算该时间窗口内GNSS计算的垂直TEC与IRI模型计算的垂直TEC的差值
Figure FDA0002329090380000061
公式表达为：
Figure FDA0002329090380000062
步骤C7：设定IG指数搜索范围，在该范围内不断调整电离层IG指数，按照步骤C3-C6，计算在输入不同IG指数条件下模型输出的电离层垂直TEC，直到满足下式为止，即得到最优化IG参量
Figure FDA0002329090380000063
搜索公式如下：
Figure FDA0002329090380000064
其中：IGmin表示IG指数搜索时取的最小值，IGmax表示IG指数搜索时取的最大值，设IGmin＝5.0，IGmax＝300.0；
步骤C8：保存所有网格点的电离层
Figure FDA0002329090380000065
指数，以便后续处理过程读取，以实现背景电离层模型IRI的驱动更新；
步骤D，基于双指数驱动的电离层三维电子密度重构：
步骤D1：将步骤B获取的Rz指数
Figure FDA0002329090380000066
输入到IRI模型中，设定IRI模型F2层峰高计算模型为AMTB，计算得到设定网格点
Figure FDA0002329090380000067
λgrid处的电离层峰值高度
Figure FDA0002329090380000068
计算方法如下：
Figure FDA0002329090380000069
其中：
Figure FDA00023290903800000610
表示展开系数，由模型默认输入文件提供，
Figure FDA00023290903800000611
为缔合勒让德函数，M为展开谐波数；
步骤D2：根据网格点的地理纬度和经度，获取网格点的磁倾角I，然后计算网格点所在位置处的修正磁倾角μ，计算公式如下：
Figure FDA00023290903800000612
其中：μ表示网格点所在位置的修正磁倾角，
Figure FDA0002329090380000071
为网格点的地理纬度；
步骤D3：计算地理坐标函数Gk，该函数计算方法：
Figure FDA0002329090380000072
其中：φ表示地理纬度，θ表示地理经度，q(k)为修正磁倾角的阶数，m(k)为经度的阶数，μ为修正磁倾角；
步骤D4：读取URSI系数文件，将步骤C中计算得到的最优化IG指数作为加权系数，按照线性组合的方式计算时间修正系数Ui,k，计算方法如下：
Figure FDA0002329090380000073
其中：
Figure FDA0002329090380000074
为URSI太阳活动低年系数，
Figure FDA0002329090380000075
为太阳活动高年系数，
Figure FDA0002329090380000076
为垂直TEC最优化得到的电离层IG指数；
步骤D5：设定IRI模型F2层峰高计算模型为URSI，计算得到设定网格点
Figure FDA0002329090380000077
λg处的电离层F2层临界频率foF2，计算方法如下：
Figure FDA0002329090380000078
其中：T表示全球时UT(-π≤T≤π)；H表示用于表征foF2日变化的谐波数的最大值，设为常数6；
步骤D6：根据计算得到的电离层峰值高度及临界频率，基于电离层电子密度剖面方程，通过输入不同的时间和地理位置，计算得到三维电子密度剖面并对数据进行输出。
2.根据权利要求1所述基于双指数驱动的电离层三维电子密度重构方法，其特征在于，所述步骤A4具体包括：
步骤A41：将地基GNSS接收机和GNSS卫星经纬高坐标
Figure FDA0002329090380000079
λ,h转换为空间直角XYZ坐标，分别标记为X0,Y0,Z0和X1,Y1,Z1，转换表达式表示为：
Figure FDA0002329090380000081
其中
Figure FDA0002329090380000082
Re为地球半径，
Figure FDA0002329090380000083
步骤A42：将地基GNSS接收机和GNSS卫星的空间直角坐标转换为站心直角坐标N,E,U，转换表达式表示为：
Figure FDA0002329090380000084
其中：T为旋转矩阵，表示方法如下：
Figure FDA0002329090380000085
其中：λ0,
Figure FDA0002329090380000086
分别对应GNSS接收机的经度和纬度坐标，其单位为弧度；
步骤A43：计算接收机与卫星间的仰角E，计算表达式为：
Figure FDA0002329090380000087
步骤A44：计算接收机与卫星间的方位角A，计算表达式为：
Figure FDA0002329090380000088
3.根据权利要求1所述基于双指数驱动的电离层三维电子密度重构方法，其特征在于，所述步骤D6具体包括：
步骤D61：将电离层F2层临界频率转换为F2层峰值电子密度值NmF2，转换公式如下：
Figure FDA0002329090380000091
其中：foF2为步骤D5计算得到的电离层F2层临界频率，其单位为Hz；
步骤D62：采用经验函数计算电离层F2层以上高度的电子密度值，计算表达式如下：
Ne1(h)＝NmF2×exp(-Y)
Figure FDA0002329090380000092
Figure FDA0002329090380000093
Figure FDA0002329090380000094
δ＝[η/(1+Z)-ζ/2]/(η×Z/β×(1+Z)2+ζ/400)
Figure FDA0002329090380000095
其中：Ne1(h)表示F2层以上高度的电离层电子密度，h为高度，ζ、β、η是与地磁纬度、太阳辐射通量及电离层foF2有关的经验参量；
步骤D63：计算电离层F1层与F2层之间高度上的电离层电子密度值Ne2(h)，计算表达式如下：
Figure FDA0002329090380000096
Figure FDA0002329090380000097
其中：B0为F2层的厚度参数，B1为F2层的形状参数；
步骤D64：计算电离层F1层高度区域上的电离层电子密度值，计算表达式如下：
Figure FDA0002329090380000101
Figure FDA0002329090380000102
其中：Ne2(h)为F1层与F2层之间高度上的电离层电子密度值，即步骤D63中的计算公式，C1是F1层的型态参数，hmF1为F1层峰值高度，F1reg为判断F1层是否存在的逻辑参数，该参数由IRI模型给定；
步骤D65：计算电离层中间层高度区域的电离层电子密度值，计算表达式如下：
Figure FDA0002329090380000103
Figure FDA0002329090380000104
T＝D2/(HZ-HEF-D)
其中：中间参数D＝HZ-HST、HZ＝(HST+hmF1)/2、HEF＝HME+HBR；而HST表示为Ne3(h)与E层峰值密度NmE相等时所对应的高度值，Ne3为F1层高度区域上的电离层电子密度值，即步骤D64采用的计算公式；HBR为电离层E层谷区的宽度；
步骤D66：计算电离层E层谷区的电离层电子密度值，计算表达式如下：
Figure FDA0002329090380000105
其中：NmE表示E层峰值密度；Ereg为是否为晚上的逻辑参数，0表示为白天，1代表晚上；参数P由多项式函数计算得到，计算方法如下：
P＝E1x2+E2x3+E3x4+E4x5
x＝h-HmE
式中：HmE白天设定为110km，晚上设定为105km，Ei，i＝1,...,4分别为E层谷区的经验修正系数，由IRI模型给定；
步骤D67：计算电离层D层高度电离层电子密度值，计算表达式如下：
Figure FDA0002329090380000111
其中：NmE为E层峰值密度；NmD为D层峰值密度；hmE为E层峰值高度；HDX为D层临界高度；D1、FP1、FP2、FP3作为经验参数由IRI模型给定；
步骤D68：将步骤D62-D67计算电离层得到的D、E、F1、F2层及F2层以上高度电离层电子密度值进行集合，得到全剖面的电子密度值。
CN201911328951.0A 2019-12-20 2019-12-20 一种基于双指数驱动的电离层三维电子密度重构方法 Active CN111125609B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911328951.0A CN111125609B (zh) 2019-12-20 2019-12-20 一种基于双指数驱动的电离层三维电子密度重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911328951.0A CN111125609B (zh) 2019-12-20 2019-12-20 一种基于双指数驱动的电离层三维电子密度重构方法

Publications (2)

Publication Number Publication Date
CN111125609A true CN111125609A (zh) 2020-05-08
CN111125609B CN111125609B (zh) 2022-11-29

Family

ID=70501035

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911328951.0A Active CN111125609B (zh) 2019-12-20 2019-12-20 一种基于双指数驱动的电离层三维电子密度重构方法

Country Status (1)

Country Link
CN (1) CN111125609B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111667567A (zh) * 2020-05-27 2020-09-15 北京天工科仪空间技术有限公司 电离层三维展示方法和装置
CN112526617A (zh) * 2020-11-19 2021-03-19 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于多源卫星信号的电离层层析成像系统观测数据模拟方法
CN113031048A (zh) * 2021-03-05 2021-06-25 中国科学院近代物理研究所 一种离子束射程快速质控验证的装置及方法
CN113093189A (zh) * 2021-04-02 2021-07-09 中国空间技术研究院 一种提高迭代初值精度的电离层层析成像方法
CN113433524A (zh) * 2021-06-23 2021-09-24 长安大学 一种联合ig值和sar反演高精度电子密度的方法
CN114384564A (zh) * 2021-12-31 2022-04-22 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于多源数据驱动的电离层层析成像方法
CN115015974A (zh) * 2021-12-31 2022-09-06 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于gnss掩星与三频信标的电离层探测性能仿真评估方法
WO2023288150A1 (en) * 2021-07-15 2023-01-19 Qualcomm Incorporated Ionosphere grid history and compression for gnss positioning
CN115792963A (zh) * 2023-02-13 2023-03-14 天津云遥宇航科技有限公司 基于gim模型的电离层校正方法、电子设备及存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105301601A (zh) * 2015-10-09 2016-02-03 中国科学院光电研究院 一种适用于全球区域的gnss电离层延迟三维建模方法
CN109613565A (zh) * 2019-01-14 2019-04-12 中国人民解放军战略支援部队信息工程大学 基于多星座gnss的电离层层析方法及系统
CN110275186A (zh) * 2019-07-11 2019-09-24 武汉大学 Leo卫星增强的gnss电离层归一化与融合建模方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105301601A (zh) * 2015-10-09 2016-02-03 中国科学院光电研究院 一种适用于全球区域的gnss电离层延迟三维建模方法
CN109613565A (zh) * 2019-01-14 2019-04-12 中国人民解放军战略支援部队信息工程大学 基于多星座gnss的电离层层析方法及系统
CN110275186A (zh) * 2019-07-11 2019-09-24 武汉大学 Leo卫星增强的gnss电离层归一化与融合建模方法

Cited By (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111667567A (zh) * 2020-05-27 2020-09-15 北京天工科仪空间技术有限公司 电离层三维展示方法和装置
CN111667567B (zh) * 2020-05-27 2023-11-24 北京天工科仪空间技术有限公司 电离层三维展示方法和装置
CN112526617A (zh) * 2020-11-19 2021-03-19 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于多源卫星信号的电离层层析成像系统观测数据模拟方法
CN113031048B (zh) * 2021-03-05 2022-11-15 中国科学院近代物理研究所 一种离子束射程快速质控验证的装置及方法
CN113031048A (zh) * 2021-03-05 2021-06-25 中国科学院近代物理研究所 一种离子束射程快速质控验证的装置及方法
CN113093189A (zh) * 2021-04-02 2021-07-09 中国空间技术研究院 一种提高迭代初值精度的电离层层析成像方法
CN113093189B (zh) * 2021-04-02 2022-12-09 中国空间技术研究院 一种提高迭代初值精度的电离层层析成像方法
CN113433524A (zh) * 2021-06-23 2021-09-24 长安大学 一种联合ig值和sar反演高精度电子密度的方法
WO2023288150A1 (en) * 2021-07-15 2023-01-19 Qualcomm Incorporated Ionosphere grid history and compression for gnss positioning
US11686850B2 (en) 2021-07-15 2023-06-27 Qualcomm Incorporated Ionosphere grid history and compression for GNSS positioning
CN115015974A (zh) * 2021-12-31 2022-09-06 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于gnss掩星与三频信标的电离层探测性能仿真评估方法
CN114384564A (zh) * 2021-12-31 2022-04-22 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于多源数据驱动的电离层层析成像方法
CN115015974B (zh) * 2021-12-31 2024-05-10 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于gnss掩星与三频信标的电离层探测性能仿真评估方法
CN114384564B (zh) * 2021-12-31 2024-05-14 中国电波传播研究所（中国电子科技集团公司第二十二研究所） 一种基于多源数据驱动的电离层层析成像方法
CN115792963A (zh) * 2023-02-13 2023-03-14 天津云遥宇航科技有限公司 基于gim模型的电离层校正方法、电子设备及存储介质
CN115792963B (zh) * 2023-02-13 2023-05-16 天津云遥宇航科技有限公司 基于gim模型的电离层校正方法、电子设备及存储介质

Also Published As

Publication number Publication date
CN111125609B (zh) 2022-11-29

Similar Documents

Publication Publication Date Title
CN111125609B (zh) 一种基于双指数驱动的电离层三维电子密度重构方法
CN114518586B (zh) 一种基于球谐展开的gnss精密单点定位方法
CN111273335B (zh) 一种基于垂测数据约束的电离层层析成像方法
Schüler On ground-based GPS tropospheric delay estimation
Rocken et al. Near real‐time GPS sensing of atmospheric water vapor
CN110568459B (zh) 基于igs和cors站的区域电离层tec实时监测方法
WO2023197714A1 (zh) 一种适用于动态载体平台的gnss多路径误差削弱方法
Nesterov et al. Ionospheric perturbation indices based on the low-and high-orbiting satellite radio tomography data
CN109917494A (zh) 降雨预报方法、装置、设备和存储介质
CN111123345B (zh) 一种基于gnss测量的经验电离层模型数据驱动方法
CN114545458A (zh) 一种高精度电离层实时建模方法
CN114384564B (zh) 一种基于多源数据驱动的电离层层析成像方法
Sedeek Ionosphere delay remote sensing during geomagnetic storms over Egypt using GPS phase observations
CN115980751A (zh) 一种幂律模型InSAR对流层延迟改正方法
CN110146904B (zh) 一种适用于区域电离层tec的精确建模方法
Liu et al. Evaluation of HY-2A satellite-borne water vapor radiometer with shipborne GPS and GLONASS observations over the Indian Ocean
Sun et al. An investigation into real-time GPS/GLONASS single-frequency precise point positioning and its atmospheric mitigation strategies
Bahadur et al. Real-time single-frequency multi-GNSS positioning with ultra-rapid products
Li et al. Bottomside ionospheric snapshot modeling using the LEO navigation augmentation signal from the Luojia-1A satellite
Yuan et al. Multipath mitigation in GNSS precise point positioning using multipath hierarchy for changing environments
CN115980317B (zh) 基于修正相位的地基gnss-r数据土壤水分估算方法
Iwabuchi et al. PPP and network true real-time 30 sec estimation of ZTD in dense and giant regional GPS network and the application of ZTD for nowcasting of heavy rainfall
CN112528213B (zh) 基于低地球轨道卫星全球电离层总电子含量多层解析方法
Chen et al. Global ionosphere modeling based on GNSS, satellite altimetry, radio occultation, and DORIS data considering ionospheric variation
CN114048585A (zh) 一种电离层模型事后分析方法及装置

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