Integrity investigation of global ionospheric TEC maps for high-precision positioning

Aside from the ionospheric total electron content (TEC) information, root-mean-square (RMS) maps are also provided as the standard deviations of the corresponding TEC errors in global ionospheric maps (GIMs). As the RMS maps are commonly used as the accuracy indicator of GIMs to optimize the stochastic model of precise point positioning algorithms, it is of crucial importance to investigate the reliability of RMS maps involved in GIMs of different Ionospheric Associated Analysis Centers (IAACs) of the International GNSS Service (IGS), i.e., the integrity of GIMs. We indirectly analyzed the reliability of RMS maps by comparing the actual error of the differential STEC (dSTEC) with the RMS of the dSTEC derived from the RMS maps. With this method, the integrity of seven rapid IGS GIMs (UQRG, CORG, JPRG, WHRG, EHRG, EMRG, and IGRG) and six final GIMs (UPCG, CODG, JPLG, WHUG, ESAG and IGSG) was examined under the maximum and minimum solar activity conditions as well as the geomagnetic storm period. The results reveal that the reliability of the RMS maps is significantly different for the GIMs from different IAACs. Among these GIMs, the values in the RMS maps of UQRG are large, which can be used as ionospheric protection level, while the RMS values in EHRG and ESAG are significantly lower than the realistic RMS. The rapid and final GIMs from CODE, JPL and WHU provide quite reasonable RMS maps. The bounding performance of RMS maps can be influenced by the location of the stations, while the influence of solar activity and the geomagnetic storm is not obvious.


Introduction
The ionospheric delay is one of the crucial error sources in global navigation satellite system (GNSS)-based precise positioning techniques, such as precise point positioning (PPP). For dual-frequency uncombined PPP, the ionospheric delay can be estimated as unknown parameters and ionospheric constraints derived from regional or global ionospheric models have been introduced in the PPP filter as pseudo-observations (Banville et al. 2014). For singlefrequency PPP, a priori ionospheric information has either been used to correct the ionospheric delay or as pseudoobservations to constrain the ionospheric parameters in PPP models (Shi et al. 2012). In these cases, the accuracy information of the corresponding ionospheric models is of crucial importance to set the optimal weighting strategy and to avoid anomalies in positioning errors, which is also helpful to reduce the convergence time and to enhance the positioning accuracy of PPP (Li et al. 2020a). Hence, it is necessary to provide the accuracy information of the ionospheric model, which are then used as an input in the PPP processing.
The global ionospheric maps (GIMs) generated by different Ionospheric Associated Analysis Centers (IAACs) of the International GNSS Service (IGS) have been widely used for ionospheric error mitigation in single-and dualfrequency PPP. As of early 2020, there are seven IAACs that provide GIM products using different techniques . Those GIMs are typically estimated from a worldwide selected subset of hundreds of permanent receivers and are accessible with latency from less than 5 min (real-time GIMs), 1 day (rapid GIMs) and up to 2-4 weeks (final GIMs) (Hernández-Pajares et al. 2009). In addition to the total electron content (TEC), the accuracy indicator is also considered in the design of GIM products. For the rapid and final GIMs, the root-mean-square (RMS) maps, contained within the files following the IONosphere Exchange (IONEX) format (Schaer et al. 1998), are generally generated as the standard deviation estimates of TEC maps. For the real-time GIMs, the field of the vertical TEC (VTEC) quality indicator in the Radio Technical Commission for Maritime Services (RTCM) standard is designed to reflect the quality of GIMs, which is currently set to zero by different IAACs because the real-time GIMs remain experimental under continuous development and improvement (Li et al. 2020b). There have been some works which apply the RMS maps in rapid and final GIM products to determine the weight of ionospheric pseudo-observations in PPP model to reduce the convergence time (Banville et al. 2014;Xiang et al. 2020). As the reliability of RMS maps has not yet been comprehensively studied, empirical variance methods are commonly adopted to determine the weight of ionospheric error (Cai et al. 2017;Shi et al. 2012;Su et al. 2019). The shortcoming of the empirical variance is that it does not reflect the real ionospheric error, especially during periods of high solar activity or geomagnetic storm. In case a small RMS value is used while the ionospheric delay calculated from GIM product is very large, an unknown bias might be introduced in the positioning result. Additionally, if the accuracy information is not reasonable, the random model constructed according to the RMS might be even worse than that of empirical model. Therefore, it is essential to investigate the reliability of RMS maps provided in the GIMs of IGS IAACs.
Recent researches had focused on the assessment and comparison of the accuracy of different GIMs Roma-Dollase et al. 2018;Li et al. 2020a), which is helpful to improve the performance of combined IGS products (IGSG and IGRG as final and rapid, respectively). However, few studies have been carried out to check the reliability of those RMS maps. Banville et al (2019) analyzed the slant TEC (STEC) errors and their correlation with the precision calculated from EMRG for the period of 2015-2018 for one Canadian station ALGO. Ghoddousi-Fard (2020) analyzed the consistence of the std of VTEC and RMS from IGS, COD, and EMR over Canada, but whether the RMS map can represent the actual error of the TEC is not covered. Considering that the performances of the RMS maps from different GIMs have not been comprehensively studied yet, we presented a method in this paper to investigates the reliability of the RMS maps from different GIMs at different locations and under different levels of ionospheric activity condition. The result provides an overview of what extent the RMS values in the GIMs of different IAACs can cover the actual ionospheric errors when used as accuracy indicators, which can also serve as a reference to choose different weighting strategies in the PPP algorithms when applying different GIMs for ionospheric error mitigation. The remainder of this paper is organized as follows: The rapid and final GIM products provided by the different IAACs and the IGS combined one are reviewed at first, followed by the detailed description of the method employed in this paper. Then, three experimental cases were performed to investigate, characterize, and compare the TEC standard deviations of the seven rapid GIMs and six final GIMs during different levels of solar and geomagnetic conditions. Finally, the conclusion and future work are given.

Review of the GIM products
Global ionospheric mapping is an active field where several techniques coexist for the determination of global ionospheric TECs using the dual-frequency GNSS measurements from a worldwide network of receivers. Currently, there are seven IAACs which developed different techniques to generate respective rapid and final GIM products. Table 1 summarizes the GIM products of different IAACs as well as the corresponding computation methods.
The MSLM and SLM in Table 1 is the mapping function currently used in the corresponding IAAC's GIM generation to convert STEC to VTEC. SLM is the standard single-layer model while MSLM is the modified standard single-layer model (Schaer S, 1999). In addition to the GIMs listed in Table 1, IGSG and IGRG are the final and rapid GIMs which resulting from the rank-weighted mean of the IAAC GIMs from CODE, ESA and JPL (Hernández-Pajares et al. 2009). All these GIMs can be downloaded from ftp://ftp.gipp.org. cn/produ ct/ionex / in IONEX format. Among these GIMs, CARG and CASG only provides the flag of availability (1 and 999 for available and unavailable, respectively) in the RMS map of the IONEX file. The NRCan only provides rapid products from 2015. Hence, we choose seven rapid GIMs (EHRG, EMRG, CORG, JPRG, UQRG, WHRG and IGRG) and six final GIMs (ESAG, CORG, JPRG, UPCG, WHUG, and IGSG) to comparatively investigate their integrity.

Methodology
The RMS maps in GIMs are considered the standard deviations of the provided TEC maps. In PPP processing, the stochastic model is generally based on the zero-mean normal distribution assumption. Ideally, if a GIM product is applied to correct the ionospheric delay or generate ionospheric pseudo-observations, the ionospheric delay error calculated by the TEC maps should follow a zero-mean normal distribution with the standard deviation derived from the corresponding RMS maps. Hence, the RMS maps can be adopted to weight the equations and to estimate the positioning error.
In fact, it is likely that the actual ionospheric error does not follow a zero-mean normal distribution. To ensure the integrity, the distribution of the actual error should be bound by the normal distribution defined by the RMS value, as expressed in Eq. (1), with the notation R for the RMS value, f(0, R) for the probability density function of a zero-mean normal distribution with a standard deviation equal to R. It can be seen that for a given positive threshold n·R, the probability of the absolute value of the actual ionospheric error not exceeding the threshold is not lower than the cumulative probability between − n·R and n·R for the zero-mean normal distribution. If Eq. (1) is satisfied for any positive n, this implies that the actual error is bounded by the normal distribution defined by the RMS value.
The actual ionospheric delay error should first be obtained to see whether Eq. (1) is satisfied. Because the slant TEC (STEC) is usually used in the precise positioning, users concern about the slant TEC. The actual error of STEC should be derived from dual-frequency carrier phase measurements to ensure its precision. However, it is difficult to obtain the actual value of the STEC due to the unknown integer ambiguities. In this paper, we adopt differential STEC (dSTEC), which is defined as the difference between the STEC measured at any epoch and the STEC measured at the epoch with the highest elevation in the same phase-continuous transmitter-receiver arc (Hernández-Pajares et al. 2017). The actual error of dSTEC can obtained by comparing the dSTEC computed by a given GIM to the observed dSTEC by the ionospheric combination of carrier phase measurements. The calculation of observed dSTEC can be expressed as: here, using dS o for the observed dSTEC, S k o (t) for the line-ofsight STEC derived from the carrier phase observations for the satellite k at the epoch t, L k I (t) for the ionospheric combination ( L I = L 1 − L 2 ) of the dual-frequency GPS carrier phases (L 1 and L 2 ), t E max for the epoch at which the highest elevation is achieved in the same phase-continuous arc, and α for the ionospheric delay scaling factor of L I from meters to TEC Unit (TECU). This observed dSTEC was reported to be more accurate than 0.1 TECU (Roma-Dollase et al. 2018), which can be regarded as a reference to obtain the error dSTEC derived from GIMs.
The dSTEC calculated from a GIM corresponding to satellite k between epochs t and t E max can be written as: where M k (t) is the function to convert VTEC to STEC and V k GIM (t) is the VTEC corresponding to satellite k at epoch t. The MSLM mapping function, which is a widely used mapping functions in precise positioning algorithm, was used in this paper. The actual error of the dSTEC derived from GIM is defined as: The RMS of dSTEC can be derived from the RMS map in GIM products. Here we assume that S k GIM (t) and S k GIM (t E max ) are independent when the difference in elevation between the highest elevation and that of epoch t is larger than 20°. Under this assumption, the variance of dSTEC can be appropriated by, are the variances of VTEC of the epochs t and t E max obtained from the RMS map of GIMs. Then, the standard deviation of dSTEC ( dS GIM ) under the assumption of normal distribution can be obtained. Equation (1) can be rewritten as: Equation (6) is applied to assess the bounding performance of dSTEC error distribution. For example, if we set n to 1, 2 and 3, the probability of | (dS)| <= n ⋅ dS GIM should then be no lower than 68.27%, 95.45%, and 99.73%, respectively. This probability can be obtained through the statistics of a large amount of data from IGS stations, which should equal to the percentage of the dSTEC error no larger than the corresponding times of dSTEC RMS. In this way, we can indirectly investigate the reliability of the RMS values as standard deviation estimates of the TEC in GIM products.

Experimental data
To study the bounding characteristic of the RMS of dSTEC derived from the RMS maps, three experiments corresponding to maximum and minimum solar cycle conditions and geomagnetic storm periods were performed with observations from IGS network GNSS receivers during the recent solar cycle, from 2010 to 2019. The Ap index can give a measure of the storminess of the Earth's magnetic field, which is obtained by computing an 8-point running average of successive 3-h ap indices during a geomagnetic storm event. The Ap range 30-49, 50-99, > 100 corresponding to minor storm, major storm, and severe storm (Bowman et al. 2008;Svalgaard et al. 2002). Hence, GNSS data with AP indices smaller than 30 for 2014 and 2019 were chosen to investigate the performance of the proposed method during the solar maximum and minimum conditions, respectively. For the geomagnetic storm period, all days with Ap indices larger than 30 during this solar cycle were selected, as shown in Fig. 1. GPS dual-frequency carrier phase measurements from three permanent GNSS stations, p047, lveg, lexa, which are not used in the generation of any involved GIMs, are selected. The stations are located at north mid-latitudes, north low-latitudes, south mid-latitudes respectively, as shown in Fig. 2. The elevation threshold of the observations is set to 15° to avoid the influence of mapping function errors on lower elevation measurements. In the computation of dSTEC, the minimum difference in elevation between the highest elevation epoch and epoch t is set to 20°, which means a long time and space differences between both involved line-of-sight epochs. Since GIMs are sometimes not available or do not provide available RMS maps, the common days that available for all the involved GIMs were used, Table 2 lists the number of days used in the three experimental periods.

Results and discussion
To provide a clear comparison between the RMS of dSTEC and the actual dSTEC error, we employ a diagram similar to the Stanford diagram (Tossaint et al. 2006) to visualize the results, as shown in Fig. 3, named the Ionospheric Model Integrity Diagram (IMID). In the top subpanel of IMID, the actual dSTEC error in TECU is presented on the horizontal axis, while the estimated RMS of dSTEC in TECU Here, we define this percentage as the RMS bounding percentage (RMSBP) and refer to the percentages corresponding to 1RMS, 2RMS and 3RMS as 1RMSBP, 2RMSBP and 3RMSBP, respectively, for the sake of simplicity, which are also shown in the IMIDs. In addition, the histograms of actual dSTEC error are all provided in the bottom of IMIDs to show and compare the accuracy of different GIMs. In this section, we first analyse and discuss the results of rapid GIMs, then the results of final GIMs are given for comparison.

Rapid GIMs
The IMIDs of different rapid GIMs during the solar minimum year of this solar cycle, 2019, were shown in Figs. 3, 4 and 5 for stations p047, lveg and lexa, respectively. By comparing the results of different GIMs in these figures, it can be seen that (1) for UQRG, almost all the points occur in the upper left part divided by the red dotted line. Compared to other GIMs, there is a quite large minimum RMS of dSTEC, which occurs due to the calibration in the RMS maps generation to guarantee the integrity considering the ionospheric protection levels of the actual VTEC error. (2) In regard to CORG, WHRG, EMRG and JPRG, most of the points are located in the upper left area; (3) The number of points located in the right bottom area is large for EHRG for all these three stations; (4) For IGRG, the number of points located in the right bottom area is quite large for station lveg.
The 1RMSBP, 2RMSBP and 3RMSBP values of the three experiments are listed in Table 3 and the histograms of 1RMSBP and 2RMSBP are shown in Figs. 6 and 7. For EMRG, only the results of minimum solar year are given. We observe that most of the RMSBPs of UQRG are 100% for all these stations, while the 1RMSBP, 2RMSBP and 3RMSBP values are all beyond the reference values for CORG. As for JPRG, WHRG and EMRG, few of the RMSBP values are slightly below the reference values. Regarding EHRG, the RMSBPs are comparatively low, and all the values are much lower than the reference values. For IGRG, most of the RMSBPs for station p047 and lexa are slightly lower than the reference values while the results for station lveg is much lower than the reference values.
The RMS and 90th percentile values of the actual dSTEC error for the selected stations are listed in Table 4, to provide a comparison of the accuracy among the GIMs. It can be found that the differences between the RMS and 90th percentile values of GIMs are within 2.63 and 4.5 TECUs, respectively (Table 5). For comparison, the mode, 10th percentile and 90th percentile values of dSTEC RMS for selected stations are listed in Table 6. The variation of the mode, 10th percentile and 90th percentile values for different GIMs are within 14, 10.75 and 16.35 TECUs, respectively. UQRG exhibits the largest value for the estimated RMS of dSTEC, which means that the generation of the RMS maps is very conservative. For CORG, JPRG, WHRG and EMRG, the estimated RMS values of dSTEC are slightly larger than the actual dSTEC error, which means that the RMS maps are relatively reasonable. Regarding EHRG, the estimated RMS of dSTEC is much smaller than the actual dSTEC error, which implies that the corresponding RMS map generation is too optimistic. For IGRG, the estimated RMS of dSTEC is smaller than the true value for most of the experiments, which means that sometimes the RMS map is quite optimistic.

Final GIMs
The results of final GIMs are presented to compare the performance of rapid and final products from the same IAACs. Only the IMIDs of different final GIMs during the solar minimum year, 2019, for station p047 are shown in Fig. 8. We can see that the distribution of IMIDs of the final GIMs for CODE, JPL, WHU, ESA and IGS are similar with the corresponding rapid GIM products. In contrast, the shape of IMID of UPCG is significant different with UQRG, the RMS values of dSTEC are concentrated in 2-5 TECU, which means the RMS of dSTEC for UPCG is not sensitive to the large error of dSTEC. The reason that the RMS of dSTEC for UPCG is much lower than UQRG lies in UPC use different strategies to calculate RMS maps in the generation of UPCG (RMS as a standard deviation proxy) and UQRG (RMS as a VTEC protection level proxy).
For more detailed analysis and comparison, 1RMSBP, 2RMSBP and 3RMSBP values of the three experiments are listed in Table 6 and the histograms of 1RMSBP and  For the final GIMs, the RMS and 90th percentile values of the actual dSTEC error are listed in Table 7. We can see that the variation of the RMS and the 90th percentile for different GIMs are within 1.07 TECUs and 2 TECUs, respectively. The mode, 10th percentile and 90th percentile values of dSTEC RMS are listed in Table 8. For CODG, JPLG, WHUG and IGSG, the differences between the mode, 10th percentile and 90th percentile values are within 3.5, 3.25 and 4.8 TECUs. For ESAG, most of these values are below 1 TECU.
By comparing the results of these three experiments, it can be seen from Tables 4 and 7 that the dSTEC error in the  Tables 3 and 6, it can be concluded that there is no obvious feature for the RMSBPs for the three conditions. By comparing the results of these three selected stations, we can find from Tables 4 and 7 that the dSTEC error for station lveg is larger than other two stations for all the rapid and final GIMs for all three experiments. The results in Tables 3 and 6 indicate that most of the RMSBP values for station lveg are obviously lower than the other two stations. Based on the above analysis, we can initially conclude that the RMSBP performance can be influenced by the location of the station while the influence of solar activity and the geomagnetic storm is not obvious.

Conclusion
GIM products are frequently employed to mitigate the ionospheric delay or for pseudo-observations to constrain the ionospheric parameters in precise positioning techniques, such as PPP. The RMS maps in GIMs are usually considered as an accuracy indicator of the corresponding TEC estimates. This paper investigates the reliability of the RMS maps in seven rapid IGS GIMs (UQRG, CORG, JPRG, WHRG, EHRG, EMRG and IGRG) and six final GIMs (UPCG, CODG, JPLG, WHUG, ESAG and IGSG) by assessing the bounding relationship between the actual dSTEC error and the dSTEC RMS derived from the RMS maps in the GIM products under the zero-mean normal distribution assumption. An ionospheric integrity diagram is designed and employed to intuitively visualize the results. Three typical experiments, under maximum and minimum solar cycle conditions, and for geomagnetic storm periods, were performed. Observation data from IGS stations not used by any of the involved GIMs were    selected to perform this comparative study. Based on these analyses, the main conclusions are as follows: 1. The reliability of the RMS maps is significantly different for GIMs from different IAACs. 2. The rapid and final GIMs from CODE, JPL and WHU provide quite reasonable RMS maps, and the distribution of the actual error is properly bounded by the normal distribution derived from the RMS map, as well as EMRG. 3. The RMS map of UQRG is the most conservative, because it has been calibrated to a large value to ensure its integrity as a sort of ionospheric protection level. In contrast, the RMS map of UPCG is slightly more optimistic than GIMs from CODE, JPL and WHU.
4. EHRG and ESAG reveal highly optimistic estimated RMS values, which indicates a quite low integrity fulfillment. 5. For IGSG and IGRG as combined products, the RMS bounding performance differ greatly for different stations. The RMSBP values for lveg are comparably lower than the other two stations. 6. The bounding performance of RMS maps can be influenced by the location of the stations while the influence of solar activity and the geomagnetic storm is not obvious.
These conclusions should be considered by users who would like to adopt the RMS maps in GIMs as standard deviation estimation of ionospheric error in their positioning algorithm, different processing strategies should be applied to different GIMs to ensure a high positioning integrity. The results in this study can provide some essential knowledge for the development of proper methods to generate reasonable VTEC precision information that can properly cover/ bound the VTEC errors.
This work is only a rough analysis and comparison of the reliability of the RMS maps for different GIMs, more indepth research is needed in the future for using RMS map as a reliable standard deviation of VTEC. Our future study will be focusing on testing and verifying the PPP performance with employing different weighting strategies for the ionospheric delay calculated from different GIMs.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.