High-precision VTEC derivation with GEONET

This paper proposes a new technique, namely Phase bias-based Small Grid Model (PSGM), to derive absolute ionospheric vertical total electron content (VTEC) with observations of Global Navigation Satellite System Earth Observation Network of Japan (GEONET). The proposed technique deals with the phase observations alone without handling the pseudoranges, which reduces the noise in VTEC estimation. A new parameter, the arc bias (Barc), is introduced to combine the phase ambiguities and differential phase biases. To solve Barc, equations are constructed under the assumption that the VTEC is identical in the same 0.1° × 0.1° grid. The performance of PSGM is evaluated with the observations in solar maximum year 2014. The root mean square error (RMSE) of PSGM is 0.40 TECU in average, the maximum RMSE is 0.73 TECU and the minimum RMSE is 0.26 TECU. The fitting accuracy of the VTEC results is improved compared with most of the existing methods.


Introduction
In the past 30 years, dual-frequency Global Positioning System (GPS) has become the most important approach for absolute vertical total electron content (VTEC) measurement. The differential code bias (DCB) of GPS hardware device is an important factor affecting the VTEC estimation. To get absolute VTEC, scientists have proposed many methods to remove DCB (e.g., Mannucci et al. 1998;Sarma et al. 2008;Zhang et al. 2009). These methods are usually based on the assumption that the DCBs of each satellite/receiver remain unchanged for a period of time. The DCBs were assumed to remain unchanged over several hours (Li et al. 2015), over an entire day (e.g., Arikan et al. 2008;Ma and Maruyama 2003;Li et al. 2018), over several days (Otsuka et al. 2002) or even over several months (Han et al. 2018).
Another basic assumption of the existing VTEC and DCB fitting models is that VTEC is smooth over a short period and a small region. The polynomial model was proposed in early times, which is suitable for local region of single station (Lanyi and Roth 1988;Coco et al. 1991). The VTEC was assumed to conform to a polynomial function of longitude and latitude. Later under the assumption that VTEC is smooth over time, a Kalman filter model was proposed (Sardon et al. 1994). To fit the global ionosphere model, the VTEC was mapped to a spherical surface harmonic function of longitude and latitude with global GPS network (Wilson et al. 1995). A well trained multilayer neural network, which is equivalent to a nonlinear, smooth and interpolated function, was applied to VTEC and DCB estimation (Ma et al. 2005;Perez 2019). Assuming that VTEC is smoothed in a small area, DCBs were estimated by constructing convolution equation set with over 1300 Global Navigation Satellite System Earth Observation Network of Japan (GEONET) stations (Li et al. 2018).
The precision of absolute VTEC estimating is not only very important for scientific research, but also affects the accuracy of space applications such as satellite navigation and communications, and remote sensing. The efforts to obtain absolute VTEC with higher precision, smaller spatial-scale and shorter temporal-scale have never ceased. Using around 300 GPS stations from GEONET, and assuming that the VTEC is identical at any point within a 2° × 2° grid, the VTEC and DCBs were determined by grid method Open Access *Correspondence: liqi@nao.cas.cn 1 National Astronomical Observatories, Chinese Academy of Sciences, Beijing, China Full list of author information is available at the end of the article (Ma and Maruyama 2003) with the standard deviation of satellite DCBs from 0.22 TECU to 1.89 TECU and the standard deviation of receiver DCBs smaller than 5.7 TECU. Zhang et al. (2010) estimated GPS DCBs with 0.5° × 0.1 h grids (0.1 h is equivalent to 1.5° in longitude as earth rotation), and found the RMS of estimated DCBs from the data observed in the low latitude region is larger than middle latitude, the RMS value varies from approximately 1 to 3 TECU, but all RMS values are less than 4 TECU. DCBs obtained by single-station method and multi-station methods were compared by Choi et al. (2013), it is found that the accuracy of DCB fitted by multi-station method is better than that of single-station method. Receiver DCB was usually believed to be affected by temperature, but even after the temperature dependence is removed, a noise level of 1-3 TECU still remains in the DCB estimation (Coster et al. 2013). With Multi-GNSS Observations and Global Ionosphere Maps, DCB's day-to-day variation was estimated for satellites and receivers: DCB's day-to-day variation in satellite was usually smaller than 0.86 TECU and in receiver smaller than 5.7 TECU (Montenbruck et al. 2014).
In this paper, in order to obtain more accurate absolute VTEC, a phase bias-based Small Grid Model (PSGM) is proposed for VTEC derivation with GEONET, which is the world's densest GNSS observation network. New assumptions are introduced as the basis of the model. Pseudorange measurements will not be used at all, but only GPS phase measurements are used for calculation as carrier phase measurement has higher precision. A new parameter B arc is introduced to represent hardware delay biases and integer ambiguities. Finally, the accuracy of PSGM are validated with the observation in solar maximum year 2014, and the estimated VTEC result is compared with the International Reference Ionosphere 2016 (IRI) model (Bilitza et al. 2017) and the VTEC products of University of Bern (UNIBE) (Schaer et al. 1996).

Algorithms of PSGM
Assuming that all the electrons in the ionosphere are concentrated on an infinitely thin spherical shell, then for a pair of receiver and satellite, the basic formulas for VTEC estimation with GPS phase observations are: where χ is the zenith of the optical path (Ma and Maruyama 2003). λ 1 , λ 2 and L 1 , L 2 represent wave lengths (unit: m) and phases (unit: cycle) observations of the corresponding carriers. (1) VTEC = STEC * cos χ, (2) STEC = n 0 (L 1 1 − L 2 2 ) + B arc , f 1 = 1.57542 GHz and f 2 = 1.22760 GHz are the frequencies of the two GPS carriers and k = 80.62 m 3 /s 2 is a constant related to the charge and mass of electron (Ma and Maruyama 2003). STEC represents the optical path integration of electron density. In Formula (2), B arc contains contribution from integer ambiguities and differential phase biases of both carriers based on the hypothesis that all these parameters are invariant within the same phase connected arc (Blewitt 1990;Hernández-Pajares et al. 1999).
The ionospheric shell concerned is divided into small grids according to 0.1° by 0.1° in latitude and longitude, and the VTEC in the same grid at the same time is assumed identical. If the ionospheric piercing points IPP i and IPP j of two arcs i, j are located in the same grid at time t, then there is: By taking Formulas (1) and (2) into (3), an equation containing two unknowns, B arci and B arcj , can be obtained: Using all eligible pairs of IPPs in a day to build equations, these equations can form an overdetermined equation set: where B is solved by least square method: Then, the fitting residuals are: The root mean square error (RMSE) of the fitting residuals is used to measure the accuracy of the model: where L represents the number of equations and Δ n represents the residual of the nth equation.
The observation data are from over 1300 GEONET stations with 1/30 Hz sampling frequency of 364 days in 2014. With 30° cutoff angle of satellite elevation, and after the arcs with cycle slip are removed, more than 40,000 arcs a day are generated. The number of equations per day is about 3 × 10 6 . So A is a sparse matrix of about 3,000,000 × 40,000; the giant matrix Eq. (5) is solved by the high performance computing cluster of the Information and Computing Center of National Astronomical Observatories of China (NAOC). The cluster has one management node, one login node and 60 computing nodes; each node has 1.28 TB DDR4-2133/2400 Ecc memory and two E5-2600v3/v4 series processors; software is Matlab R2012a glnxa64; hundreds of tasks can be processed in parallel. It takes about 16 h to process the observation data for a whole day with single task, from reading the GPS data to calculating Formula (8).
Compared with the DCB based grid method, the B arcbased PSGM depends on a prerequisite of intensive observation data. Only if an arc has at least one neighbor IPP to form Eq. (3), its B arc can appear in the matrix Eq. (5), otherwise its B arc may not be able to be solved out. In the GEONET observation data of 2014, the number of B arc s that cannot be solved accounts for about 8% of total available arcs.
As long as B arc is solved, the VTEC of the whole arc can be obtained by Formulas (1) and (2). If there is only one IPP in a 0.1° × 0.1° grid, the VTEC of the IPP is considered as the VTEC of the grid. If there are multiple IPPs, the median VTEC of these IPPs is taken as the VTEC of the grid.
The observation data cannot guarantee to cover all the grids all the time as the grids are quite small. Taking the grid centered at 135°E 35°N (near Kobe) as an example, the VTEC available ratio is about 60%, that is to say, around 1700 (60% of 2880) epochs in a day have VTEC measurements. The available ratio of grids 135°E 40°N and 135°E 30°N (above the ocean) decrease rapidly to around 15%. This shows that the available ratio decreases rapidly due to observation density reduction.

VTEC result
In order to simplify the problem, the ionospheric shell height parameter is assumed to be a constant h = 400 km (Ma and Maruyama 2003). The GEONET observations of 2014 were fitted with PSGM. For comparison and analysis, we also calculated the VTEC based on the 1° × 1° grid size, which is referred as Phase bias-based Big Grid Model (PBGM) VTEC. The algorithms of PBGM and PSGM are exactly the same except for the grid size. The IRI-VTEC and UNIBE-VTEC were introduced in some figures. Need to note that IRI-VTEC and UNIBE-VTEC can be considered as external reference and validation here rather than the object of criticism and evaluation because they are based on different backgrounds and environments. Figure 1 displays the VTEC variation on 16 magnetically quiet (Kp < 4) days of four seasons of 2014. The diurnal variation in PSGM-VTEC is consistent with UNIBE and IRI in general. The smoothness and continuity of PSGM are very good. The difference between PSGM and PBGM will not exceed 2 TECU in most cases, which shows the self-consistency and stability of the algorithm. The most obvious difference between the two is in Fig. 1b. The diurnal peaks of the PBGM result around noon are slightly higher than those of PSGM, and the fluctuation is more intense as well. This may indicate that PSGM is better than PBGM in accuracy. March is more than 30 TECU higher than that of 15 March and 17 March. This TEC enhancement gradually weakened with increasing latitude. The peak VTEC decreased rapidly from about 100 TECU at 30°N to about 65 TECU at 35°N, that is to say, VTEC had a spatial gradient of up to 7 TECU/°. So in this case, the assumption of constant VTEC in a 1° × 1° grid will undoubtedly lead to larger errors. Reflected in Fig. 2c, the VTEC curve of PBSM contains some burrs. Figure 3 displays the VTEC seasonal variation at the grid point of [135°E, 35°N] for daytime (04:00 UT, 13:00 JST) and nighttime (16:00 UT, 01:00 JST) of 364 days in 2014. In daytime, VTEC is the largest in spring, followed by autumn and winter, and smallest in summer. At nighttime, VTEC is the largest in summer, followed by spring, and then in autumn and winter. The trends of PSGM-VTEC, PBGM-VTEC, IRI-VTEC and UNIBE-VTEC are similar.
The average difference between PSGM and PBGM is 5.04% (1.7 TECU) in daytime and 7.76% (0.88 TECU) in nighttime. In Fig. 3b, the nighttime UNIBE-VTEC is generally larger than PSGM-VTEC, which may be caused by these reasons: (1) the shell height parameter is 450 km in UNIBE, VTEC may be slightly overestimated than 400 km (Wang et al. 2016); (2) the UNIBE-VTEC is provided every 2.5° of latitudes, lower latitude coverage (33.75°N to 35°N) may enlarge the VTEC estimation at 35°N; (3) the assumption that DCB in one day is unchanged may overestimate nighttime VTEC compared with the assumption that phase bias of an arc is unchanged.
The curve of PSGM-VTEC is quite consistent with that of F10.7 from April to August from Fig. 3. The correlation coefficient of the PSGM-VTEC and the F10.7 is 0.63 in nighttime during April to August while it is only − 0.22 for the whole year 2014. In daytime, the correlation coefficient is 0.52 during April to August and 0.30 for the whole year. It shows that solar radiation is an important  Figure 4 displays the two-dimensional VTEC maps at each hour from 10:00 UT to 19:00 UT of 11 July 2014. The map coverage is between 125°E-150°E and 25°N-50°N. Every 2 TECU is displayed as a color level. VTEC contours are clearly exhibited in Fig. 4a. VTEC decreased gradually with the increase of latitude. From 12:00 UT to 15:00 UT, a series of northwest-southeast (NW-SE) wave stripes were formed. These structures had wavelengths of several hundred kilometers, propagated to southwest at ~ 150 m/s and recovered gradually after 16:00 UT, which is believed to be a typical nighttime medium-scale traveling ionospheric disturbances (MSTID) (Tsugawa et al. 2007). The example in Fig. 4 shows that PSGM can obtain VTEC measurements for the whole Japan region, and PSGM can at least reflect the structure and evolution of ionospheric disturbances that greater than 2.0 TECU level.

Cross-validation
Though there are some comparisons between PSGM and PBGM in "VTEC result" section, they cannot be used as the exact evidence to explain the accuracy of the method as they all derived from the same observation data. To further illustrate the accuracy of PSGM, all GEONET stations are randomly divided into two groups for cross-validation. The stations of the two groups are completely different. The PSGM-VTEC results of these two groups are named GROUPA-VTEC and GROUPB-VTEC, respectively. Figure 5 shows the comparison of fitting results on 4 example days, the average difference at 135°E 35°N between GROUPA-VTEC and GROUPB-VTEC is only 0.36 TECU. If all the grids with both sets of VTEC results are counted in, the average VTEC difference between the two groups is 0.46 TECU and only 1.92% grid VTEC differences are higher than 2.0 TECU. Cross-validation further proves the accuracy of PSGM algorithm. It should be noted that since the number of stations of each group is reduced by half, the number of available IPP neighbor pairs is also greatly reduced, so the fitting accuracy of total stations should be higher than that of any group.

RMSE result
Parameter RMSE in Formula (8) can reflect how well the data is fitted. The RMSE level is the possible difference between the two measurements in the same grid. Figure 6 shows the RMSE result in 2014. On the whole, PSGM-RMSE is the largest in summer and smallest in winter. The average PSGM-RMSE of 364 days is 0.40 TECU. The maximum PSGM-RMSE is 0.73 TECU on 13 May and the minimum RMSE is 0.26 TECU on 21 February. This means that PSGM is appropriate to be applied to the study of ionospheric disturbances at 1.0 TECU level theoretically.
As a contrast, the average PBGM RMSE of 364 days is 0.64 TECU. The maximum PBGM RMSE is 1.12 TECU on 16 March and the minimum RMSE is 0.39 TECU on 18 January. PSGM-RMSE is about 37% less than PBGM RMSE. On the other hand, the biggest season of PBGM RMSE has changed to spring, which is caused by the larger TEC spatial gradient in spring. The mean value of the rate of TEC index (ROTI) (Pi et al. 1997) of all observation in a day is introduced as the measurement of ionospheric disturbing level. Figure 6b examined the relationship between PSGM-RMSE and mean ROTI. It can be seen that PSGM-RMSE increases with the increase of mean ROTI. The correlation between ROTI and PSGM-RMSE is a clue that the accuracy of PSGM can be affected by ionospheric disturbances.
PSGM can also be applied to a short period rather than an entire day. Figure 7 shows the PSGM-RMSE obtained  Li et al. Earth, Planets and Space (2020) 72:14 by fitting the observation data of 4 h in the daytime and the nighttime, respectively. In summer, the daytime PSGM-RMSE is smallest, but the nighttime PSGM-RMSE is largest. And in winter, PSGM-RMSE increased in daytime and decreased at night. It is interesting that the relationship between the daytime RMSE and nighttime RMSE is opposite in winter and summer. On the whole, though the daytime VTEC is much larger, the average RMSE at night (0.33 TECU) is slightly larger than that during the day (0.29 TECU), which is also consistent with the objective reality that more ionospheric disturbances happened at night.

Shell height parameter optimization
At the beginning of "VTEC result" section, we assume that the height of ionospheric shell is fixed at 400 km to simplify the calculation. But actually, the height of the ionospheric shell always changes with time and place (Wang et al. 2016;Lu et al. 2017). For PSGM, the change of shell height parameter is easier to map IPPs to different grids under different satellite elevation than PBGM. In turn, this feature can also be used to optimize the ionospheric shell height parameter. When the shell height parameter changes, not only the grid mapping of the IPP will change, but also the zenith angle χ, so PSGM-RMSE usually does not change smoothly with the shell height parameter. In addition, due to the narrow station distribution of GEONET, when the shell height parameter is too large, the number of neighbor IPP pairs composed of different satellites will be greatly reduced, and the PSGM-RMSE will always tend to decrease as shell height parameter increases. Even though, at most of the time, we can still find a shell height parameter between 250 and 600 km which can make PSGM-RMSE minimum. Figure 8 shows the relationship between RMSE and shell height parameter of the 6 example days 01:00-05:00 UT in 2014. We performed calculation every 10 km of 250-600 km according to the shell height parameter and keep all other conditions unchanged. The shell height parameter of minimum PSGM-RMSE in these 6 days is 340 km, 390 km, 390 km, 450 km, 440 km and 360 km. The reason why we choose the observation data of 4 h at noon rather than the entire day is that the ionosphere height changes the least and the TEC value is the highest during this period. Figure 9 shows the variation of the optimum height (OH) parameters estimated by the minimum RMSE method during 01:00-05:00 in 2014. The overall trend is the lowest in winter and the highest in summer, which is consistent with the result estimated by Zhao and Zhou (2018) using single station DCB method. The OH in summer is often over 500 km; while in winter it is often as low as 300 km but never less than 250 km. The parabola in the figure is the trend of the shell height parameter variation fitted by polynomial, which is consistent with the trend of hmF2 in IRI. PSGM-RMSE under the OH parameter is about 4.4% less than that of 400 km. It shows that there is a potential to further improve the fitting accuracy if the variable shell height parameter is used instead of the constant shell height parameter.

Summary
The accuracy of VTEC estimation using dual-frequency GNSS system is a very meaningful research topic. In this paper, the RMSE of VTEC fitting residuals is analyzed by using the GEONET observation and the PSGM. The PSGM is based on two assumptions. The first assumption that the integer ambiguities and hardware delay biases of dual-frequency GPS phase observation are invariant within the same phase connected arc can help to get rid of the dependence on low-precision pseudorange measurements. The second assumption that the VTEC is identical in the same 0.1° × 0.1° grid can minimize the prerequisite for spatio-temporal smoothing.
The trends of PSGM, PBGM, UNIBE and IRI are consistent in both diurnal and seasonal variations. PSGM can reduce the impact of large TEC spatial gradient. The PSGM-VTEC can be used to detect some ionosphere phenomena such as VTEC enhancement and MSTID.
The cross-validation error is about 0.46 TECU in average. The PSGM-RMSE in summer is larger than that in winter and increases with the increase of mean ROTI. It is found that under the simple isotropic fixed height ionospheric shell model, the PSGM-RMSE is as low as 0.40 TECU in average in the year of 2014. This progress will refresh the understanding of the upper limit of the fitting accuracy for thin shell model. The trend of PSGM-RMSE is different in daytime and nighttime. In daytime, PSGM-RMSE is the largest in winter and the smallest in summer, and it is the opposite at night.
PSGM excludes most error sources and minimizes the impact of spatial variation, so that the difference between the assumed shell height and the actual shell height can be reflected in the change of RMSE. By comparing the PSGM-RMSE corresponding to different shell height parameters, the OH parameter that can minimize the RMSE can be found. In 2014, the OH in the daytime is the lowest in winter and the highest in summer.