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

本发明属于空间电离层特征参量现报预报领域，特别涉及该领域中的一种基于双指数驱动的电离层三维电子密度重构方法。

背景技术：

电离层作为日地空间环境的重要组成部分，会对穿越其中的无线电波产生折射、反射、散射和吸收等效应，从而影响卫星导航、通信、雷达等诸多无线电信息系统的性能。利用各类技术手段来获取电离层的特征参量，研究电离层中的各种现象，揭示其内含的物理机制，具有重要的科学意义。在诸多电离层的特征参量中，电离层总电子含量(totalelectroncontent,tec)和电子密度是其中最为关键的探测参数，尤其是电离层电子密度变化。目前，很多无线电信息系统如卫星导航及其增强系统，地基雷达折射修正、天基低频段sar高精度目标成像的工程应用均需要三维高精度电离层电子密度信息支撑。掌握其时空变化特征和规律，对提升无线电信息系统性能具有重要的实用价值。

根据实际应用场景的不同，实现电子密度重构需要求解的未知系数不尽相同。但针对需要较高分辨率电子密度的应用领域，通常需要求解较大规模的未知系数，这将需要极大的内存和存储空间及非常大的cpu计算量，无法满足实际应用场景对快速的计算时效性、简便的嵌入式应用的需求。为降低计算量，提高运行时效以便于将算法嵌入至系统应用中，降低三维电子密度的重构时需要的待估未知数的数量是非常必要的。

当前，随着全球gps、glonass、北斗、galileo等卫星导航系统的发展，以及星载掩星测量技术的进步，综合利用gnss地基观测数据与gnss天基掩星测量数据，基于数据驱动技术实现全球高精度电离层三维电子密度重构成为世界各国研究的热点方向。很多研究者开始尝试在各类地基/天基电离层探测的基础上，利用数据驱动技术对经验电离层模型的一些特定输入参数进行驱动更新，一方面降低计算量和存储空间，另一方面也可以较好的提升电离层电子密度重构的精度，以更好的满足工程应用需求。

技术实现要素：

本发明所要解决的技术问题就是提供一种基于双指数驱动的电离层三维电子密度重构方法。

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

一种基于双指数驱动的电离层三维电子密度重构方法，其改进之处在于，包括如下步骤：

步骤a，多gnss数据源电离层垂直总电子含量映射：

步骤a1：根据指定日期，下载igs发布的gnss观测文件和gnss卫星星历文件，对数据进行解压缩，提取观测台站的经度、纬度和高度，分别记为λr,hr，提取gnss卫星与观测台站间的伪距和载波相位其中i为对应的观测频段i＝1,2，j为对应的卫星编号，r为对应的接收机编号；

步骤a2：对于所有的伪距和载波相位观测值，分别计算无几何距离线性组合，计算方法如下：

其中：为伪距的无几何距离线性组合，为载波相位的无几何距离线性组合；

步骤a3：在同一观测弧段内，利用载波相位的无几何距离线性组合l4对伪距的无几何距离线性组合p4进行平滑滤波，计算方法表达为：

其中：表示载波平滑伪距量，是计算时的中间过程量，t表示历元时刻，ωt为t时刻的加权因子；当t＝1时，设定

步骤a4：利用gnss星历文件，计算gnss卫星的轨道坐标；同时计算卫星和接收机之间的仰角e和方位角a；

步骤a5：假定电离层总电子含量tec集中在地球上空450km高度上的一个薄层内，计算观测站与gnss卫星连线在该薄层高度上的穿刺点的地理坐标λ，计算公式如下：

其中：λ0为gnss观测站的纬度和经度坐标，a为方位角；ψ0为观测站与gnss卫星间的地心夹角，计算公式如下：

其中：re＝6371.2km为地球的平均半径；hi＝450km为电离层薄层高度，e为接收机至卫星间的仰角；

步骤a6：计算倾斜总电子含量与垂直总电子含量的映射函数值，计算方法如下：

m(e)＝[1-(recose/(re+hi))²]^-1/2

步骤a7：利用球谐函数构建电离层垂直总电子含量随空间分布模型，计算方法如下：

其中：nmax表示球谐函数展开的最大阶数，为缔合勒让德函数，为穿刺点的地理纬度，s为穿刺点在日固坐标系下的经度，anm,bnm为球谐函数的未知系数；

步骤a8：利用平滑后的伪距的无几何距离线性组合构建总电子含量映射求解方程，计算方法如下：

其中：f1为卫星的第一个频率，f2为卫星的第二个频率，dcb^r为第r个接收机的硬件延迟，dcbj为第j颗卫星的硬件延迟值，为了对卫星和接收机的硬件延迟进行分离，同时加入以下四个约束条件：

其中：表示第j颗gps卫星的硬件延迟值，ng为gps卫星的总数量；表示第j颗glonass卫星的硬件延迟值，nr为glonass卫星的总数量；表示第j颗北斗卫星的硬件延迟值，nc为北斗卫星的总数量；表示第j颗galileo卫星的硬件延迟值，ne为galileo卫星的总数量；

步骤a9：选择1～2小时为时间窗口，遍历该时间窗口内所有的卫星与接收机，根据步骤a8的方程建立垂直总电子含量的线性矩阵如下：

y＝ax

式中，y是所有接收机与卫星间的观测量的集合，a是与球谐函数、接收机和卫星硬件延迟相关的系数矩阵，x是待求解的未知系数矢量，

步骤a10：利用最小二乘算法求解拟合方程，解算得到未知系数x，计算方法如下：

xls＝(a^ta)^-1a^ty

其中：xls是最小二乘解，其包含待估计的未知系数，包括球谐函数系数anm,bnm，接收机和卫星硬件延迟；

步骤a11：利用解算得到的系数anm,bnm，重新代入到球谐函数模型中，解算得到指定网格点处的电离层垂直tec值，计算方法如下：

其中：λg是网格点的纬度和经度坐标，即为解算得到的垂直总电子含量，m表示累加序号，即m＝0,1,2…n；

步骤b，最优化电离层rz指数确定：

步骤b1：根据与步骤a1相同的日期，下载cdaac发布的gnss掩星电子密度剖面数据，提取掩星碰撞点经度、纬度和高度处对应的电子密度值ne及掩星发生时刻t，记为

步骤b2：对一次完整的掩星事件，提取该掩星事件对应的电离层f2层峰值高度及其时间和位置信息，记为计算方法如下：

其中，ho表示观测点处的高度，表示对应的掩星碰撞点的位置和时间信息，r表示一次掩星事件的所有采样点的集合；

步骤b3：对时间窗口内所有掩星事件的电离层f2层峰值高度数据进行质量控制，剔除异常的观测数据，符合要求的掩星峰值高度数据应满足以下公式：

其中：det表示判决条件，δλ分别表示一次掩星事件期间起始时刻碰撞点纬度、经度与终止时刻碰撞点纬度、经度的绝对差值；

步骤b4：选择国际参考电离层iri-2016模型为背景电离层模型，设定选择amtb为iri的峰值高度模型，根据提取的掩星事件对应的电离层f2层峰值高度及其时间和位置信息，调用电离层经验模型iri，计算对应时间及地点处的电离层f2层峰值高度，记为

步骤b5：设定时间窗口为24小时，获取该时间窗口内所有掩星事件的电离层f2层峰值高度数据序列，计算两者的差值δhmf2,i，公式表达为：

其中：i表示为当前时间窗口内第i次掩星事件；λo,ho,to表示f2层峰值高度对应的位置和时间信息；

步骤b6：设定rz指数搜索范围，在该范围内不断调整电离层rz指数，输出对应rz指数条件下iri模型计算的电离层f2层峰值高度值hmf2iri，直到满足下式为止，即得到最优化rz参量搜索公式如下：

其中，n为设定的时间窗口内所有掩星事件的总数，rzmin表示rz搜索时取的最小值，rzmax表示rz搜索时取的最大值，设rzmin＝-40.0，rzmax＝400.0；

步骤b7：保存电离层指数，以便后续处理过程读取；

步骤c，最优化电离层ig指数确定：

步骤c1：读取步骤a处理得到的电离层垂直tec值及其对应的网格点坐标，设定λgrid是网格点的纬度和经度坐标，为gnss观测数据解算得到的电离层垂直tec；

步骤c2：设定国际参考电离层模型采用国际无线电联合会发布的ursi系数计算f2层临界频率，输入最优化指数，将其作为模型的输入量，计算网格点λgrid对应位置处的电子密度值设定高度h的范围为100-1500km；

步骤c3：对100-1500km高度范围内的电子密度值沿着垂直方向进行积分，得到iri模型输出的垂直tec值计算公式如下：

步骤c4：计算高度1500km至卫星高度20200km区间的等离子体层的电子密度值，计算方法如下：

其中：表示等离子体层的电子密度，h'表示电离层顶部的高度，可近似等于1500km，h为电子密度的标高；

步骤c5：对1500-20200km高度范围内的电子密度值沿着垂直方向进行积分，得到磁层模型输出的垂直tec值计算公式如下：

其中：h1和h2分别表示积分高度的下限和上限；

步骤c6：选择对应的时间窗口，计算该时间窗口内gnss计算的垂直tec与iri模型计算的垂直tec的差值公式表达为：

步骤c7：设定ig指数搜索范围，在该范围内不断调整电离层ig指数，按照步骤c3-c6，计算在输入不同ig指数条件下模型输出的电离层垂直tec，直到满足下式为止，即得到最优化ig参量搜索公式如下：

其中：igmin表示ig指数搜索时取的最小值，igmax表示ig指数搜索时取的最大值，设igmin＝5.0，igmax＝300.0；

步骤c8：保存所有网格点的电离层指数，以便后续处理过程读取，以实现背景电离层模型iri的驱动更新；

步骤d，基于双指数驱动的电离层三维电子密度重构：

步骤d1：将步骤b获取的rz指数输入到iri模型中，设定iri模型f2层峰高计算模型为amtb，计算得到设定网格点λgrid处的电离层峰值高度计算方法如下：

其中：表示展开系数，由模型默认输入文件提供，为缔合勒让德函数，m为展开谐波数；

步骤d2：根据网格点的地理纬度和经度，获取网格点的磁倾角i，然后计算网格点所在位置处的修正磁倾角μ，计算公式如下：

其中：μ表示网格点所在位置的修正磁倾角，为网格点的地理纬度；

步骤d3：计算地理坐标函数gk，该函数计算方法：

其中：φ表示地理纬度，θ表示地理经度，q(k)为修正磁倾角的阶数，m(k)为经度的阶数，μ为修正磁倾角；

步骤d4：读取ursi系数文件，将步骤c中计算得到的最优化ig指数作为加权系数，按照线性组合的方式计算时间修正系数ui,k，计算方法如下：

其中：为ursi太阳活动低年系数，为太阳活动高年系数，为垂直tec最优化得到的电离层ig指数；

步骤d5：设定iri模型f2层峰高计算模型为ursi，计算得到设定网格点λg处的电离层f2层临界频率fof2，计算方法如下：

其中：t表示全球时ut(-π≤t≤π)；h表示用于表征fof2日变化的谐波数的最大值，设为常数6；

步骤d6：根据计算得到的电离层峰值高度及临界频率，基于电离层电子密度剖面方程，通过输入不同的时间和地理位置，计算得到三维电子密度剖面并对数据进行输出。

进一步的，所述步骤a4具体包括：

步骤a41：将地基gnss接收机和gnss卫星经纬高坐标λ,h转换为空间直角xyz坐标，分别标记为x0,y0,z0和x1,y1,z1，转换表达式表示为：

其中re为地球半径，

步骤a42：将地基gnss接收机和gnss卫星的空间直角坐标转换为站心直角坐标n,e,u，转换表达式表示为：

其中：t为旋转矩阵，表示方法如下：

其中：λ0,分别对应gnss接收机的经度和纬度坐标，其单位为弧度；

步骤a43：计算接收机与卫星间的仰角e，计算表达式为：

步骤a44：计算接收机与卫星间的方位角a，计算表达式为：

进一步的，所述步骤d6具体包括：

步骤d61：将电离层f2层临界频率转换为f2层峰值电子密度值nmf2，转换公式如下：

其中：fof2为步骤d5计算得到的电离层f2层临界频率，其单位为hz；

步骤d62：采用经验函数计算电离层f2层以上高度的电子密度值，计算表达式如下：

ne1(h)＝nmf2×exp(-y)

δ＝[η/(1+z)-ζ/2]/(η×z/β×(1+z)²+ζ/400)

其中：ne1(h)表示f2层以上高度的电离层电子密度，h为高度，ζ、β、η是与地磁纬度、太阳辐射通量及电离层fof2有关的经验参量；

步骤d63：计算电离层f1层与f2层之间高度上的电离层电子密度值ne2(h)，计算表达式如下：

其中：b0为f2层的厚度参数，b1为f2层的形状参数；

步骤d64：计算电离层f1层高度区域上的电离层电子密度值，计算表达式如下：

其中：ne2(h)为f1层与f2层之间高度上的电离层电子密度值，即步骤d63中的计算公式，c1是f1层的型态参数，hmf1为f1层峰值高度，f1reg为判断f1层是否存在的逻辑参数，该参数由iri模型给定；

步骤d65：计算电离层中间层高度区域的电离层电子密度值，计算表达式如下：

t＝d²/(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层谷区的电离层电子密度值，计算表达式如下：

其中：nme表示e层峰值密度；ereg为是否为晚上的逻辑参数，0表示为白天，1代表晚上；参数p由多项式函数计算得到，计算方法如下：

p＝e1x²+e2x³+e3x⁴+e4x⁵

x＝h-hme

式中：hme白天设定为110km，晚上设定为105km，ei，i＝1,...,4分别为e层谷区的经验修正系数，由iri模型给定；

步骤d67：计算电离层d层高度电离层电子密度值，计算表达式如下：

其中：nme为e层峰值密度；nmd为d层峰值密度；hme为e层峰值高度；hdx为d层临界高度；d1、fp1、fp2、fp3作为经验参数由iri模型给定；

步骤d68：将步骤d62-d67计算电离层得到的d、e、f1、f2层及f2层以上高度电离层电子密度值进行集合，得到全剖面的电子密度值。

本发明的有益效果是：

本发明所公开基于双指数驱动的电离层三维电子密度重构方法，采用国际全球卫星导航系统组织(internationalgnssservice，igs)发布的全球gnss观测数据和cosmic数据分析和档案中心(cosmicdataanalysisandarchivecenter，cdaac)发布的gnss掩星电子密度剖面数据产品作为电离层三维电子密度重构的数据来源，基于最优化的rz和ig指数对国际参考电离层模型(internationalreferenceionosphere，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(globalnavigationsatellitesystem)数据源电离层垂直总电子含量映射：

步骤a1：根据指定日期，下载igs发布的gnss观测文件和gnss卫星星历文件，对数据进行解压缩，提取观测台站的经度、纬度和高度，分别记为λr,hr，提取gnss卫星与观测台站间的伪距和载波相位其中i为对应的观测频段i＝1,2，j为对应的卫星编号，r为对应的接收机编号；

步骤a2：对于所有的伪距和载波相位观测值，分别计算无几何距离线性组合(geometry-freecombination)，计算方法如下：

其中：为伪距的无几何距离线性组合，为载波相位的无几何距离线性组合；

步骤a3：在同一观测弧段内，利用载波相位的无几何距离线性组合l4对伪距的无几何距离线性组合p4进行平滑滤波，计算方法表达为：

其中：表示载波平滑伪距量，是计算时的中间过程量，t表示历元时刻，ωt为t时刻的加权因子；当t＝1时，设定

步骤a4：根据接收机坐标与利用gnss星历文件计算得到的gnss卫星轨道坐标，计算卫星和接收机之间的仰角e和方位角a；

步骤a4具体包括：

步骤a41：将地基gnss接收机和gnss卫星经纬高坐标λ,h转换为空间直角xyz坐标，分别标记为x0,y0,z0和x1,y1,z1，转换表达式表示为：

其中re为地球半径，

步骤a42：将地基gnss接收机和gnss卫星的空间直角坐标转换为站心直角坐标n,e,u，转换表达式表示为：

其中：t为旋转矩阵，表示方法如下：

其中：λ0,分别对应gnss接收机的经度和纬度坐标，其单位为弧度；

步骤a43：计算接收机与卫星间的仰角e，计算表达式为：

步骤a44：计算接收机与卫星间的方位角a，计算表达式为：

步骤a5：一般假定电离层总电子含量tec集中在地球上空450km高度上的一个薄层内，计算观测站与gnss卫星连线在该薄层高度上的穿刺点的地理坐标(λ，计算公式如下：

其中：λ0)为gnss观测站的纬度和经度坐标，a为方位角；ψ0为观测站与gnss卫星间的地心夹角，计算公式如下：

其中：re＝6371.2km为地球的平均半径；hi＝450km为电离层薄层高度，e为接收机至卫星间的仰角；

步骤a6：计算倾斜总电子含量与垂直总电子含量的映射函数值，计算方法如下：

m(e)＝[1-(recose/(re+hi))²]^-1/2

步骤a7：利用球谐函数构建电离层垂直总电子含量随空间分布模型，计算方法如下：

其中：nmax表示球谐函数展开的最大阶数，为缔合勒让德(associatedlegendre)函数，为穿刺点的地理纬度，s为穿刺点在日固坐标系下的经度，anm,bnm为球谐函数的未知系数；

步骤a8：利用平滑后的伪距的无几何距离线性组合构建总电子含量映射求解方程，计算方法如下：

其中：f1为卫星的第一个频率，f2为卫星的第二个频率，dcb^r为第r个接收机的硬件延迟，dcbj为第j颗卫星的硬件延迟值，为了对卫星和接收机的硬件延迟进行分离，同时加入以下四个约束条件：

其中：表示第j颗gps卫星的硬件延迟值，ng为gps卫星的总数量；表示第j颗glonass卫星的硬件延迟值，nr为glonass卫星的总数量；表示第j颗北斗卫星的硬件延迟值，nc为北斗卫星的总数量；表示第j颗galileo卫星的硬件延迟值，ne为galileo卫星的总数量；

步骤a9：选择1～2小时为时间窗口，遍历该时间窗口内所有的卫星与接收机，根据步骤a8的方程建立垂直总电子含量的线性矩阵如下：

y＝ax

式中，y是所有接收机与卫星间的观测量的集合，a是与球谐函数、接收机和卫星硬件延迟相关的系数矩阵，x是待求解的未知系数矢量，

步骤a10：利用最小二乘算法求解拟合方程，解算得到未知系数x，计算方法如下：

xls＝(a^ta)^-1a^ty

其中：xls是最小二乘解，其包含待估计的未知系数，包括球谐函数系数anm,bnm，接收机和卫星硬件延迟等；

步骤a11：利用解算得到的系数anm,bnm，重新代入到球谐函数模型中，解算得到指定网格点处的电离层垂直tec值，计算方法如下：

其中：λg是网格点的纬度和经度坐标，即为解算得到的垂直总电子含量，m表示累加序号，即m＝0,1,2…n；

多gnss数据源电离层垂直总电子含量的映射结果如图2所示。

步骤b，最优化电离层rz指数确定：

步骤b1：根据与步骤a1相同的日期，下载cdaac发布的gnss掩星电子密度剖面数据，提取掩星碰撞点(经度、纬度和高度)处对应的电子密度值ne及掩星发生时刻t，记为

步骤b2：对一次完整的掩星事件，提取该掩星事件对应的电离层f2层峰值高度及其时间和位置信息，记为计算方法如下：

其中，ho表示观测点处的高度，表示对应的掩星碰撞点的位置和时间信息，r表示一次掩星事件的所有采样点的集合；

步骤b3：对时间窗口内所有掩星事件的电离层f2层峰值高度数据进行质量控制，剔除异常的观测数据，符合要求的掩星峰值高度数据应满足以下公式：

其中：det表示判决条件，δλ分别表示一次掩星事件期间起始时刻碰撞点纬度、经度与终止时刻碰撞点纬度、经度的绝对差值；

步骤b4：选择国际参考电离层iri-2016模型为背景电离层模型，设定选择amtb为iri的峰值高度模型(amtb为一种电离层hmf2模型，分别代表四个作者的首字母(altadill,d.,s.magdaleno,j.m.torta,ande.blanch))，根据提取的掩星事件对应的电离层f2层峰值高度及其时间和位置信息，调用电离层经验模型iri，计算对应时间及地点处的电离层f2层峰值高度，记为

步骤b5：设定时间窗口为24小时，获取该时间窗口内所有掩星事件的电离层f2层峰值高度数据序列，计算两者的差值δhmf2,i，公式表达为：

其中：i表示为当前时间窗口内第i次掩星事件；λo,ho,to表示f2层峰值高度对应的位置和时间信息；

步骤b6：设定rz指数搜索范围，在该范围内不断调整电离层rz指数，输出对应rz指数条件下iri模型计算的电离层f2层峰值高度值hmf2iri，直到满足下式为止，即得到最优化rz参量搜索公式如下：

其中，n为设定的时间窗口内所有掩星事件的总数，rzmin表示rz搜索时取的最小值，rzmax表示rz搜索时取的最大值，设rzmin＝-40.0，rzmax＝400.0；

步骤b7：保存电离层指数，以便后续处理过程读取；

最优化电离层rz指数的计算结果如图3所示。

步骤c，最优化电离层ig指数确定：

步骤c1：读取步骤a处理得到的电离层垂直tec值及其对应的网格点坐标，设定λgrid是网格点的纬度和经度坐标，为gnss观测数据解算得到的电离层垂直tec；

步骤c2：设定国际参考电离层模型采用国际无线电联合会(internationalunionofradioscience,ursi)发布的ursi系数计算f2层临界频率，输入最优化指数，将其作为模型的输入量，计算网格点λgrid对应位置处的电子密度值设定高度h的范围为100-1500km；

步骤c3：对100-1500km高度范围内的电子密度值沿着垂直方向进行积分，得到iri模型输出的垂直tec值计算公式如下：

步骤c4：计算高度1500km至卫星高度20200km区间的等离子体层的电子密度值，计算方法如下：

其中：表示等离子体层的电子密度，h'表示电离层顶部的高度，可近似等于1500km，h为电子密度的标高(scaleheight)；

步骤c5：对1500-20200km高度范围内的电子密度值沿着垂直方向进行积分，得到磁层模型输出的垂直tec值计算公式如下：

其中：h1和h2分别表示积分高度的下限和上限；

步骤c6：选择对应的时间窗口，计算该时间窗口内gnss计算的垂直tec与iri模型计算的垂直tec的差值公式表达为：

步骤c7：设定ig指数搜索范围，在该范围内不断调整电离层ig指数，按照步骤c3-c6，计算在输入不同ig指数条件下模型输出的电离层垂直tec，直到满足下式为止，即得到最优化ig参量搜索公式如下：

其中：igmin表示ig指数搜索时取的最小值，igmax表示ig指数搜索时取的最大值，设igmin＝5.0，igmax＝300.0；

步骤c8：保存所有网格点的电离层指数，以便后续处理过程读取，以实现背景电离层模型iri的驱动更新；

最优化电离层ig指数的计算结果如图4所示。

步骤d，基于双指数驱动的电离层三维电子密度重构：

步骤d1：将步骤b获取的rz指数输入到iri模型中，设定iri模型f2层峰高计算模型为amtb，计算得到设定网格点λgrid处的电离层峰值高度计算方法如下：

其中：表示展开系数，由模型默认输入文件提供，为缔合勒让德函数，m为展开谐波数；

步骤d2：根据网格点的地理纬度和经度，获取网格点的磁倾角i，然后计算网格点所在位置处的修正磁倾角μ，计算公式如下：

其中：μ表示网格点所在位置的修正磁倾角，为网格点的地理纬度；

步骤d3：计算地理坐标函数gk，该函数计算方法：

其中：φ表示地理纬度，θ表示地理经度，q(k)为修正磁倾角的阶数，m(k)为经度的阶数，μ为修正磁倾角；

步骤d4：读取ursi系数文件，将步骤c中计算得到的最优化ig指数作为加权系数，按照线性组合的方式计算时间修正系数ui,k，计算方法如下：

其中：为ursi太阳活动低年系数，为太阳活动高年系数，为垂直tec最优化得到的电离层ig指数；

步骤d5：设定iri模型f2层峰高计算模型为ursi，计算得到设定网格点λg处的电离层f2层临界频率fof2，计算方法如下：

其中：t表示全球时ut(-π≤t≤π)；h表示用于表征fof2日变化的谐波数(harmonics)的最大值，设为常数6；

步骤d6：根据计算得到的电离层峰值高度及临界频率，基于电离层电子密度剖面方程，通过输入不同的时间和地理位置，计算得到三维电子密度剖面并对数据进行输出。

步骤d6具体包括：

步骤d61：将电离层f2层临界频率转换为f2层峰值电子密度值nmf2，转换公式如下：

其中：fof2为步骤d5计算得到的电离层f2层临界频率，其单位为hz；

步骤d62：采用经验函数计算电离层f2层以上高度的电子密度值，计算表达式如下：

ne1(h)＝nmf2×exp(-y)

δ＝[η/(1+z)-ζ/2]/(η×z/β×(1+z)²+ζ/400)

其中：ne1(h)表示f2层以上高度的电离层电子密度，h为高度，ζ、β、η是与地磁纬度、太阳辐射通量及电离层fof2有关的经验参量；

步骤d63：计算电离层f1层与f2层之间高度上的电离层电子密度值ne2(h)，计算表达式如下：

其中：b0为f2层的厚度参数(bottom-sidethicknessparameter)，b1为f2层的形状参数(shapeparameter)；

步骤d64：计算电离层f1层高度区域上的电离层电子密度值，计算表达式如下：

其中：ne2(h)为f1层与f2层之间高度上的电离层电子密度值，即步骤d63中的计算公式，c1是f1层的型态参数(f1shapeparameter)，hmf1为f1层峰值高度，f1reg为判断f1层是否存在的逻辑参数，该参数由iri模型给定；

步骤d65：计算电离层中间层高度区域的电离层电子密度值，计算表达式如下：

t＝d²/(hz-hef-d)

其中：中间参数d＝hz-hst、hz＝(hst+hmf1)/2、hef＝hme+hbr；而hst表示为ne3(h)与e层峰值密度nme相等时所对应的高度值，ne3为f1层高度区域上的电离层电子密度值，即步骤d64采用的计算公式；hbr为电离层e层谷区(e-valleyregion)的宽度；

步骤d66：计算电离层e层谷区的电离层电子密度值，计算表达式如下：

其中：nme表示e层峰值密度；ereg为是否为晚上的逻辑参数，0表示为白天，1代表晚上；参数p由多项式函数计算得到，计算方法如下：

p＝e1x²+e2x³+e3x⁴+e4x⁵

x＝h-hme

式中：hme白天设定为110km，晚上设定为105km，ei，i＝1,...,4分别为e层谷区的经验修正系数，由iri模型给定；

步骤d67：计算电离层d层高度电离层电子密度值，计算表达式如下：

其中：nme为e层峰值密度；nmd为d层峰值密度；hme为e层峰值高度；hdx为d层临界高度；d1、fp1、fp2、fp3作为经验参数由iri模型给定；

步骤d68：将步骤d62-d67计算电离层得到的d、e、f1、f2层及f2层以上高度电离层电子密度值进行集合，得到全剖面的电子密度值。

电离层三维电子密度重构的结果如图5所示。