Real-Time Global Ionospheric Map and Its Application in Single-Frequency Positioning

The prevalence of real-time, low-cost, single-frequency, decimeter-level positioning has increased with the development of global navigation satellite systems (GNSSs). Ionospheric delay accounts for most errors in real-time single-frequency GNSS positioning. To eliminate ionospheric interference in real-time single-frequency precise point positioning (RT-SF-PPP), global ionospheric vertical total electron content (VTEC) product is designed in the next stage of the International GNSS Service (IGS) real-time service (RTS). In this study, real-time generation of a global ionospheric map (GIM) based on IGS RTS is proposed and assessed. There are three crucial steps in the process of generating a real-time global ionospheric map (RTGIM): estimating station differential code bias (DCB) using the precise point positioning (PPP) method, deriving slant total electron content (STEC) from PPP with raw observations, and modeling global vertical total electron content (VTEC). Experiments were carried out to validate the algorithm’s effectiveness. First, one month’s data from 16 globally distributed IGS stations were used to validate the performance of DCB estimation with the PPP method. Second, 30 IGS stations were used to verify the accuracy of static PPP with raw observations. Third, the modeling of residuals was assessed in high and quiet ionospheric activity periods. Afterwards, the quality of RTGIM products was assessed from two aspects: (1) comparison with the Center for Orbit Determination in Europe (CODE) global ionospheric map (GIM) products and (2) determination of the performance of RT-SF-PPP with the RTGIM. Experimental results show that DCB estimation using the PPP method can realize an average accuracy of 0.2 ns; static PPP with raw observations can achieve an accuracy of 0.7, 1.2, and 2.1 cm in the north, east, and up components, respectively. The average standard deviations (STDs) of the model residuals are 2.07 and 2.17 TEC units (TECU) for moderate and high ionospheric activity periods. Moreover, the average root-mean-square (RMS) error of RTGIM products is 2.4 TECU for the one-month moderate ionospheric period. Nevertheless, for the high ionospheric period, the RMS is greater than the RMS in the moderate period. A sub-meter-level horizontal accuracy and meter-level vertical accuracy can be achieved when the RTGIM is employed in RT-SF-PPP.


Introduction
The use of global navigation satellite systems (GNSSs) has widely increased in fields such as geodesy, navigation, precision agriculture, and hazard monitoring [1][2][3][4]. High-accuracy fields using spherical harmonic functions. However, their ionospheric model is still a regional ionospheric model, which cannot be used by global real-time users.
In this study, we focus on the generation of an RTGIM by using observations from current real-time stations in real-time mode. First, IGS RTS products are introduced in the following section. Then, we introduce the methodology used for estimating differential code bias (DCB), deriving STEC PPP with raw observations, and modeling real-time VTEC, as well as the real-time process. Moreover, we report the experiments that were carried out to verify the methods' effectiveness. Finally, the quality of RTGIM products was assessed by comparing them with the CODE GIM and carrying out RT-SF-PPP experiments.

IGS RTS Products
There are 8 real-time IGS individual ACs: the Bundesamt für Kartographie und Geodäsie (BKG), the Centre National d'Etudes Spatiales (CNES), Deutsches Zentrum für Luft und Raumfahrt (DLR), the European Space Operations Centre (ESOC), the GeoForschungsZentrum (GFZ), the GMV Aerospace and Defense (GMV), the Natural Resources Canada (NRCan), and Wuhan University (WUHAN). Correction streams from the IGS' individual ACs are listed in Table 1. For most ACs, two RTS products with different reference points are available. The reference point is an antenna phase center (APC) or center of mass (CoM). RTS products are generated by ACs and then directly broadcasted to real-time users. Also, as listed in Table 2, ESOC and BKG are two IGS RTS Combination Centers (CCs) that provide IGC01/IGS01, IGS02, and IGS03 products. The combination IGC01/IGS01 is a single epoch combination product that is the result of combining CLK10, CLK16, CLK20, CLK22, CLK53, CLK70, CLK80, and CLK93. IGS02 is a Kalman Filter combination product generated by combining CLK10, CLK16, CLK20, CLK22, CLK53, CLK70, and CLK80. Moreover, another Kalman Filter combination product, IGS03, is the result of combining CLK11, CLK91, CLK20, and CLK80.
Zhang et al. [7] studied the accuracies of 9 RTS streams, namely, IGS01, IGS03, CLK01, CLK15, CLK22, CLK52, CLK70, CLK81, and CLK90. Results showed that the accuracy of orbit products ranges from 3.8 to 7.5 cm for different RTS products, and clock accuracy ranges from 1.9 to 5.6 cm. The individual RTS product CLK90 and combined product IGS03 have been proved to have the best quality. For these overall RTS products, only the products CLK90, CLK91, CLK92, and CLK93, all of which are provided by the CNES, contain a VTEC message. Nie et al. [8] assessed the VTEC products and concluded that sub-meter-level horizontal and meter-level vertical positioning solutions with RT-SF-PPP could be obtained. Thus, in this study, precise orbit and clock products from IGS03 were adopted for the generation of the RTGIM, and the accuracy of the CNES VTEC was used as a reference.

Methodology
An infinitesimal thin shell is adopted for most of the current VTEC ionospheric models. A VTEC model can be defined as spherical harmonic expansions; the model is widely used, similar to the usage frequency of CODE GIM products and CNES real-time VTEC products. To build a global VTEC model, STEC derived from globally distributed stations is required, and it can be obtained from PPP with raw observations method [13]. Moreover, DCBs of stations are necessary for PPP with raw observations method and can be estimated in the post-PPP process [12]. In this section, the derivation of STEC from PPP with raw observations method and the estimation of DCB with the PPP method are presented. Afterward, global ionospheric spherical harmonic modeling is described. Finally, data processing for generating an RTGIM is briefly introduced.

PPP with Raw Observations
The linearized equations for PPP with raw observations can be expressed as [11,13,14] ∆P s r, f = e s r · ∆ x + c∆t r + µ fÎ where the superscript s and subscripts r and f denote specific satellites, receivers, and frequencies, respectively; ∆P s r, f and ∆L s r, f are the observed-minus-computed (O-C) raw code and phase observations with corrections (tidal effects, antenna PCO/PCV, phase wind-up, etc.); e s r is the unit vector from satellite to receiver; ∆ x = [∆x, ∆y, ∆z] is the station increment vector; c denotes the speed of light; ∆t r is the clock offset for the receivers; µ f = (λ s f /λ s 1 ) 2 is the frequency-related factor; I s r is the ionospheric delay on L1 frequency; T r is the site-specific zenith tropospheric delay; m s r is the mapping function; λ s f is the wavelength of the frequency f ;N s r, f is the float phase ambiguity absorbing receiver and satellite phase hardware delays; p and L are the observation noises, multipath effects, and other unmodeled errors for code and phase observations, respectively. Satellite orbit and clock errors are corrected with IGS03 RTS products. The symbols with aˆon top denote reparametrized estimable parameters [11]: where ∆t r is the receiver clock containing code hardware delays;Î s r,1 is the ionospheric delay biased by the satellite and receiver DCBs, where DCB s = b s 1 − b s 2 and DCB r = b r,1 − b r,2 ;N s r, f is the float phase ambiguity containing code and phase hardware delays. More details of the method for deriving STEC by PPP with raw observations can be found in Liu et al. [11], Tu et al. [12], Zhang et al. [13].
As per the raw observation equations (Equation (1)), if DCBs are known, the unknown parameters are the station coordinate increment vectors, the slant ionospheric delay, and the float phase ambiguity.
Satellite DCB can be obtained from the CODE. However, DCBs of some stations are unknown. Thus, DCBs of stations are estimated before applying real-time PPP.

DCB Estimation by PPP Method
Differential Code Biases (DCBs) are the systematic errors or biases between two GNSS code observations at the same or different frequencies. To estimate receiver DCBs, Equation (1) can be rewritten as follows: In Equation (3), the satellite DCB is known; if the ionosphere prior constraint is added, the receiver DCB can be estimated in the PPP process [12]. The CODE GIM was used as the prior constraint.
For the stations without P1 and P2 observations, C1 and C2 were used instead. Thus, satellite P1-C1 and P2-C2 DCBs were added during the PPP process [15].

Spherical Harmonic Function Model
In Equations (1) and (2), the ionospheric delay I S r,1 can be estimated in the PPP process. Based on the infinitesimal thin shell hypothesis, VTEC can be obtained from the ionospheric delay I S r,1 [16]: where MF(z) is the ionospheric mapping function; z is the satellite elevation angle; R is the Earth's radius; H is the attitude of the ionosphere thin shell. H and α can be set by the users. Here, H and α were set to 506.7 km and 0.9782. The spherical harmonic function, which is often used to express global ionospheric models, is expressed as follows [3]: where φ is the geocentric latitude of the ionospheric pierce point (IPP), and λ is the solar-fixed longitude of the IPP. N is the degree of spherical harmonic function and was set to 12 in this study; P m n is the regularization Legendre series; A m n and B m n are the spherical harmonic coefficients for estimation. VTEC is vertical total electron content (TEC) at the IPP.

Real-Time Data Process
The RTGIM generation process is shown in Figure 1. To generate GIM in real time, receiving and decoding the real-time streams of RTS and observations is the first step. An RTCM networked transport via an Internet protocol (NTRIP [17]) client, SGGNtrip, is used. SGGNtrip was developed at the School of Geodesy and Geomatics (SGG), Wuhan University. The client can receive, decode, and store (in a database) RTCM messages, such as observations, broadcast ephemeris, orbit and clock corrections, and other RTCM messages. The client is designed based on a database, which means that once real-time data are stored in the database, the data can be accessed efficiently for different applications, including the PPP process, real-time ionospheric monitoring, or DCB estimation.
Orbit and clock, as well as global observations streams, are necessary for real-time PPP. However, orbit and clock streams are always later than observations because of product latency [7]. In real-time PPP, once an observation of an epoch is received, the epoch is processed. The newest received corrections of orbit and clock are applied to observation time.
DCBs of stations are estimated every day by a post-PPP process. Monthly P1-C1 and P1-P2 DCBs of satellites from the CODE are downloaded every month and used in the post-PPP process. Because DCB is stable for several days, the estimated DCB value is not employed until the difference from the used DCB is too large.
STEC derived from the real-time PPP process is also saved in a database. The modeling process queries the most recent 1-h STEC data from the database every 5 min. Afterward, VTEC is converted by a mapping function from STEC, and the spherical harmonic coefficients are estimated. RTGIM products are generated with a resolution of 5 (longitude) × 2.5 degrees (latitude), which is the same as the resolution of the CODE GIM products. In the global ionospheric modeling process, outlier detection is a crucial factor that affects ionospheric accuracy. There are two outlier detection means in the generation of an RTGIM. The first is applied in the real-time PPP process. Coordinates are estimated in real time, and the reference coordinates of stations are known beforehand. Once the difference is greater than 10 cm on the horizontal component and 25 cm on the vertical component, the epoch is regarded as an invalid solution, and the STEC at this epoch is marked as an outlier. The second is in global ionospheric modeling. If a model residual of STEC is greater than 3 times the STD, it is also marked as an outlier, and twice iteration will be adopted for spherical harmonic coefficient estimation.
Actual computation time and product latency are crucial for RTGIM products. There are three time-consuming procedures in RTGIM generation: real-time observation receiving, PPP computation, and spherical harmonic coefficient estimation. In general, the latency of a real-time observation is about 2-3 s; PPP computation costs less than 3 s for about 80 stations. The most time-consuming procedure is the estimation of spherical harmonic coefficients, which costs about 63 s for a 12-degree spherical harmonic function. Thus, the overall time-cost is about 70 s for each product update.

Experiments and Results
In this section, experiments are divided into two parts. First, experiments of DCB estimation, static PPP with raw observations, and global ionospheric modeling are presented to verify the effectiveness of the RTGIM generation method. Second, the quality of RTGIM products is assessed.

Experiments for Algorithm Validation
Static PPP with raw observations is a crucial step that is needed to obtain STEC. DCBs of satellites and stations are necessary for the PPP process. DCB estimation using the PPP method was first validated. Afterward, to validate the performance of static PPP with raw observations, global distribution station experiments were carried out. Finally, the modeling method was evaluated in different ionospheric activity periods.

Test of DCB Estimation
DCBs are necessary for estimating slant ionospheric delay in real-time PPP. Jin et al. [16] proposed a DCB estimation method that was based on the CCL method and global ionospheric modeling. Experimental results showed that DCBs of station accuracy were better than 0.3 ns for most stations.
To investigate the accuracy of the DCB estimation using the PPP method, experiments using globally distributed stations were carried out. Observation data from 16 global distribution IGS stations were collected from day of year (DOY) 121 (1 May) to DOY 151 (31 May), 2018. The orbit and clock obtained from IGS03 were fixed, and satellites P1-P2, P1-C1, and P2-C2 from May 2018 (provided by the CODE) were adopted in the PPP process. The prior information constraint of the ionospheric delay was also required for DCB estimation. Thus, the CODE GIM from the same day was used. Finally, DCBs, tropospheric parameters, coordinates, and ambiguities were estimated.
The Moreover, the mean bias and standard deviation (STD) of the DCB differences are presented in Figure 3. As we can see, the largest value is 0.6 ns for the station COCO, and the smallest value is 0.02 ns for the station DAV1.

Test of Static PPP
To investigate the accuracy of real-time static PPP as determined with raw observations, 30 stations from DOY 135 (May 15), 2018, were collected for the experiment. The orbit and clock from IGS03 were employed. P1-P2 and P1-C1 DCBs of satellites, as well as the estimated P1-P2 DCBs of stations, were fixed in the PPP process. More detailed static PPP settings are provided in Table 3. The IGS station coordinates from the IGS weekly solution or from SOPAC were used as a reference. The positioning errors of 4 stations are shown in Figure 4. As shown in the figure, the accuracy of the station HARB for the three components converges easily after a short time.

Test of Ionospheric Modeling
To conduct a global ionospheric model, STEC values derived from global stations by the PPP method are required. Data from real-time stations on DOY 135, 2018, were used for global ionospheric modeling. As shown in Figure 6, the real-time station distribution is uneven. The highest-density stations are in Western Europe and Oceania, with more than half of all stations distributed from within these two regions. In the Eastern Asian region, there are more stations than in North America and Northern Africa. Moreover, few stations are distributed in the Middle East, South Asia, and the Arabian Peninsula. The uneven spread of the stations will cause an uneven IPP distribution, which can lead to negative values in blank areas. To avoid negative values, a prior constraint was added to the modeling process. The RTGIM of the previous day for the same time of day was used as the constraint, and the VTEC accuracy at each grid point was set to 10 TEC units (TECU). Spherical harmonic coefficients were estimated every 5 min. One-hour STEC values before the modeling epoch were used. Frequency distribution histograms of model residuals were plotted every 1 h to verify the accuracy of the models. As shown in Figure 8, the distributions of the residuals are within normal expectations. Almost all the mean values are less than 0.02 TECU, and the STD values are less than 2.5 TECU. The minimum STD value is 1.68 TECU at 10:00 UT, and the maximum STD value is 2.47 TECU at 6:00 UT. The mean of the STD values is 2.07 TECU. In contrast, the residual distribution on DOY 238, 2018, is presented in Figure 9. As shown in the figure, the residual distributions are also within normal expectations. The mean of the residuals varies from 0 to 0.08 TECU. In particular, from 13:00 UT to 23:00 UT, the mean of the residuals is greater than that from 1:00 UT to 12:00 UT. The minimum STD value is 1.8 TECU at 20:00 UT, and the maximum STD is 2.7 TECU at 10:00 UT. The mean of the STD values is 2.17 TECU, which is slightly greater than the value on DOY 135.

Quality Assessment of RTGIM Products
The quality assessment of RTGIM contains two aspects: first, the accuracy compared with the reference product; second, the positioning accuracy when the RTGIM product is adopted.

Comparison with CODE GIM Products
To assess the accuracy of RTGIM products, the CODE GIM was set as the reference, and the differences between each grid points relative to the CODE GIM were calculated. The global RMS values for all grid points at each epoch from DOY 121 to DOY 151 are presented in Figure 11. As shown in the figure, the minimum value is 1.3 TECU at 20:00 UT on DOY 140, and the maximum value is 4.9 TECU at 5:00 UT on DOY 123. The mean value is 2.4 TECU, and the STD value is 0.7 TECU. To further verify the accuracy of the RTGIM in different ionospheric activity periods, bias and RMS for different latitudes at 12:00 UT on DOY 135 and DOY 238 are shown in Figures 12 and 13, respectively. As shown in the figures, the bias and RMS of the RTGIM share similar spatial features. For the bias map at 12:00 UT on DOY 135, the bias values vary from −0.87 to 1.11 TECU, and the RMS values vary from 0.45 to 3.06 TECU. The maximum bias is around 30 degrees north latitude, and the maximum RMS is around 5 degrees north latitude. At 12:00 UT on DOY 238, the bias values vary from −1.33 to 3.45 TECU, and the RMS values vary from 0.11 to 6.73 TECU. The maximum bias is around 10 degrees north latitude, and the maximum RMS is around 5 degrees north latitude. The RMS values of the global grid points vary from 1.8 to 3.1 TECU for 12:00 UT on DOY 135 and DOY 238, 2018, respectively. By comparing RMS values between the two days, we can see that the accuracy during the moderate ionospheric activity period is better than that during the high activity period. In the study by Nie et al. [8], the CNES VTEC products were assessed by comparing them with the IGS final GIM products. The ranges of the bias and RMS for the different latitudes are listed in Table 4. As shown in the table, two different degrees of spherical harmonic functions were used for period 1 and period 2. For the 12-degree function in period 2, the bias at different latitudes varies from −2.82 to 0.06 TECU, and the RMS for the CNES product varies from 0.97 to 3.01 TECU, which is similar to the RTGIM in the moderate period. For the high ionospheric activity period, the accuracy is slightly worse than that in the quiet period of the RTGIM.

Single-Frequency PPP Performance in User Domain
A real-time ionospheric model is also necessary for RT-SF-PPP. To investigate positioning accuracy when adopting the RTGIM in the user domain, RT-SF-PPP experiments with 31 globally distributed stations were conducted. Please note that RT-SF-PPP was simulated by processing the GPS data in post mode. Orbit and clock corrections from IGS03 were applied to correct orbit and clock errors, and ionospheric delay errors were corrected by RTGIM. Tropospheric errors were corrected by the Saastamoinen model [18]. The RT-SF-PPP settings are listed in Table 5. The precise coordinates of the stations from the IGS weekly SINEX solution or from SOPAC were used as the reference. Positioning errors of RT-SF-PPP are shown in Figure 14 MKEA  MQZG  GLPS  CHPI  CORD  OHI3  MAW1  STHL  SUTH  ZAMB  SEY2  MBAR  ULAB  DAEJ  IISC  LCK4  ARTU  TASH  YAR2  SYDN  MAS1  FLRS  NICO  ZIM2  LAMA  MDO1  DUBO  WILL  IQAL  GODE Figure 15. According to the global map of RMS on DOY 135, the accuracy is uneven for different regions. The RMS varies from 0.3 to 7.9 TECU. For most regions of the world, the RMS is less than 3 TECU. Nevertheless, for some regions, such as the Western Pacific, Northern Africa, Southern Asia, the Arabian Peninsula, and the Northern Indian Ocean, the RMS values are greater than those in other regions. As presented in Figure 15a, the 3-D RMS values of the global stations vary from 0.5 to 2 m. There are 6 stations with an accuracy worse than 1.4 m, namely, IISC, MBAR, GLPS, SEY2, WILL, and MAS1. As we can see from the two maps, the 6 stations are in low-accuracy regions. Thus, the positioning accuracy is obviously affected by the accuracy of the ionospheric product.

Conclusions
Ionospheric delay correction is necessary for RT-SF-PPP. In this study, the generation of an RTGIM based on IGS RTS products and real-time IGS stations is proposed. Experiments were carried out to validate the algorithm and assess product quality.
Generating a real-time global ionospheric map requires three crucial steps: DCB estimation with the PPP method, STEC estimation using PPP with raw observations method, and global ionospheric modeling. To assess the DCB estimated by the PPP method, one month's data from 16 stations were collected and used to estimate DCBs for each day between DOY 121, 2018, and DOY 151, 2018. Experimental results show that the accuracy of DCBs for the stations is about 0.2 ns, and the average STD value is 0.1 ns; these values represent an accuracy that is similar to those derived using the CCL DCB estimation method [16]. STEC was derived from static PPP with raw observations. To assess the performance of the PPP process, 30 globally distributed stations, as well as IGS03 RTS products, were used for positioning experiments. Results show that the mean accuracies for the stations are 0.7, 1.2, and 2.1 cm in the three components (north, east, and up, respectively). To assess the VTEC modeling accuracy for different ionospheric activity periods, VTEC values derived from real-time stations on DOYs 135 and 238 were used in the modeling tests. Results show that the STDs of the model residuals are 2.07 and 2.17 TECU on average in the moderate and high ionospheric periods, respectively. The modeling accuracy in the high activity period is slightly greater than that in the moderate period.
Product quality was assessed from two aspects. First, the differences from the reference for one month's data (from DOY 121 to DOY 151, 2018) were evaluated. Second, RT-SF-PPP experiments were performed by using 31 globally distributed IGS stations. With the CODE GIM products as references, the RMS values of different grid points at each epoch vary from 1.3 to 4.9 TECU. The mean value is 2.4 TECU, and the STD is 0.7 TECU. To further verify the RTGIM in the high ionospheric activity period, the products on DOY 238 were evaluated. Bias and RMS were calculated for different latitudes. The RMS varies from 0.45 to 3.06 TECU on DOY 135 and from 0.45 to 3.06 TECU on DOY 238. The accuracy in the high ionospheric period is lower than that in the moderate period. The RTGIM shares an approximate accuracy with the CNES VTEC product. The RT-SF-PPP experiments show that the mean RMS values in the north, east, and up components are 0.46, 0.52, and 1.01 m. The performance is also close to that of the CNES VTEC product. The poor-performing stations are in the poor RTGIM accuracy regions, where there is low coverage by the reference stations.
Therefore, it is concluded that it is feasible to generate real-time global ionospheric map products based on IGS RTS products and real-time IGS stations. Accuracy at the sub-meter level and meter level on horizontal and vertical components, respectively, can be achieved for RT-SF-PPP.