Introduction

In Satellite-Based Augmentation Systems (SBASs), such as the US Wide-Area Augmentation System (WAAS), the ionospheric delay corrections are broadcast in terms of ionospheric grid delay (IGD) and grid ionospheric vertical error (GIVE) on GPS L1 frequency (1575.42 MHz) at each ionospheric grid point (IGP) (Chao et al. 1995; Shah and Desai 2018). The user receives the information to correct their own ionospheric delays so as to improve positioning accuracy (Andrei et al. 2009). Since the initial operation of WAAS in 2003, IGD estimation algorithms, including the planar method and the kriging method based on the thin-shell model, have been successively applied to improve the availability and service performance of WAAS. It is also proved that the kriging method can better represent the changes of IGD than the planar method (Sparks et al. 2010; Walter et al. 2018; Sparks et al. 2011a, b).

A further study on WAAS IGD estimation algorithms has been carried out. For example, Blanch et al. (2004a) took advantage of three-dimensional ideas to expand the conventional single-layer thin-shell kriging method to a multi-layer thin-shell kriging method to calculate IGDs. The results show that the new kriging algorithm can reduce the ionospheric delay error by 30–50%, while the planar method can only reduce the error by 15–30%. Furthermore, Blanch et al. (2004b) applied this multi-layer thin-shell kriging algorithm to estimate the ionospheric delay in low-latitude regions so that the impact of ionospheric disturbances on slant ionospheric errors can be reduced by nearly 30% in relation to the planar fit. Compared to the thin-shell model, where a mapping function has to be used for the conversion of total electron content (TEC) from the slant to the vertical direction, the ionospheric tomography overcomes the limitations of 2D ionospheric delay modeling and can realize 3D or even 4D ionospheric electron density inversion, and thus resulting in accurate TEC estimation (Howe et al. 1998; Bust and Mitchell 2008); especially, ionospheric tomography is widely used to research and monitor large-scale ionospheric structures at both ionospherically quiet times and storm times (Bust et al. 2007).

As one of the advanced tomographic algorithms, the Multi-Instrument Data Analysis System (MIDAS), proposed by Mitchell and Spencer (2003), can incorporate various data to carry out the inversion of 3D ionospheric electron density and then calculate the ionospheric delay. For example, Global Positioning System (GPS) data during the major storm of July 2000 over Europe were selected by Meggs et al. (2004) to reconstruct the vertical TEC (VTEC) from MIDAS without the use of the conversion of slant TEC into VTEC by a mapping function in the thin shell model. They also found that the inversion offered a more accurate estimate of ionospheric delay and provided improvements over the thin shell in TEC mapping at middle latitudes. Dos Santos Prol et al. (2018) developed a new tomographic method using both ionosonde and radio-occultation (RO) measurements to generate a new ionospheric background in the region of Brazil. Compared with global ionospheric maps in 6 GNSS stations that were not used in tomography, the new method has shown an improvement of 59% in TEC and 31% in the single-frequency precise point positioning (PPP).

For the first time by virtue of electron density inversions, the tomographic method (MIDAS) alone and further combined with the kriging approach (referred to as “kriging-combined tomographic method” hereafter) are proposed innovatively to estimate WAAS IGDs in the range of 70° W-125° W and 20° N-50° N over the US during ionospheric disturbed periods from September 7 to 9, 2017. The planar method and the kriging method are also used to make comparison. Then, those IGDs are validated by independent GPS dual-frequency data at six selected WAAS reference stations. Furthermore, the feasibility and advantage of tomography in estimating IGP delays at the strongest storm time between 0:00–3:00UT and 12:00–15:00UT on September 8 are also examined.

Ionospheric disturbance

It is found that geomagnetic storms with perturbed geomagnetic indices greatly affect ionospheric disturbances, and they are highly correlated. For example, the higher the Kp and the lower the Dst, the more intense the ionospheric disturbance. Therefore, in this study, geomagnetic indices Kp and Dst which are obtained from the Geomagnetic Data Service Center are used to analyze the characteristics of ionospheric disturbances. Variations in Kp and Dst during September 7–9, 2017, are illustrated in Fig. 1.

Fig. 1
figure 1

Variations in Kp and Dst during September 7–9, 2017

Full size image

Figure 1 shows that Kp values throughout the day on September 9, 2017, were below 3, and Dst values were in the range of -30nT-0nT, indicating that the geomagnetic activity was low. However, on September 7 and 8, the ionospheric activity changed significantly. For example, Dst began to gradually decrease starting at 21:00 UT on September 7 and Kp directly increased to more than 7. Until 01:00UT on September 8, Dst reached the minimum of -124 nT, and Kp increased to the maximum of 8, indicating that the ionospheric disturbance reached its strongest point. From then on, Dst gradually increased and Kp decreased, but the ionospheric conditions were still active as shown in Fig. 1; Dst was below -50nT and Kp kept above 4 on September 8. Generally, as Kp values increase from 0 to 9 the intensity of geomagnetic activity becomes stronger. On the contrary, as Dst values decrease lower, especially below zero, the level of geomagnetic activity becomes stronger.

Ionospheric grid delay estimation methods

In this section, four ionospheric grid delay estimation methods, i.e., the planar method, the kriging method, the tomographic method, and the kriging-combined tomographic method, will be described in detail.

Ionospheric delay

The slant total electron content (STEC) can be represented as the linear integral of the electron density along the satellite-to-receiver transmission path, which is given by:

$${\text{STEC}} = \int\limits_{{s^{^{\prime}} }} {N_{e} } {\text{d}}s$$
(1)

where \(N_{e}\) is the ionospheric electron density and \(s\) is the satellite to a receiver transmission path. The vertical ionospheric delay \(d_{{{\text{ion}}}}\) is given by:

$$d_{{{\text{ion}}}} { = }\frac{40.3}{{f^{2} }}{\text{VTEC}}$$
(2)

where \(f\) denotes the carrier frequency of the GPS signal; VTEC is the vertical total electron content at the ionospheric pierce point (IPP). IPP refers to the intersection of the satellite-to-ground receiver transmission signal and the ionospheric thin shell as described below.

The ionospheric thin shell model assumes that free electrons in the earth’s atmosphere are concentrated on a thin spherical shell with a thickness of 0 and a fixed height from the earth’s surface (Mannucci et al. 1993; Wan et al. 2016). The altitude of the thin shell is assumed to be 350 km in this study.

STEC can be typically derived from the following equations:

$${\text{STEC = VTEC}} \times {\text{MF}} + Q_{{{\text{GPS}}}}$$
(3)
$${\text{MF}} = \cos \left( {\arcsin \left( {\frac{{R_{E} }}{{R_{E} + H}}\sin \alpha } \right)} \right)$$
(4)
$$Q_{{{\text{GPS}}}} = \frac{1}{40.308}\frac{{\left( {f_{1}^{2} f_{2}^{2} } \right)}}{{\left( {f_{1}^{2} - f_{2}^{2} } \right)}}\left( {c*\vartriangle \tau } \right)$$
(5)

where \({\text{MF}}\) denotes the obliquity factor converting STEC to VTEC and \(Q_{{{\text{GPS}}}}\) is the error related to differential code biases (DCBs) of GPS satellites and receivers in equation (3);\(R_{{\text{E}}}\) stands for the radius of the earth (6378 km), \(H\) is the fixed height of the thin shell (350 km), and \(\alpha\) represents the elevation angle in (4); \(f_{1}\) and \(f_{2}\) are the GPS operating frequencies (Hz), \(c\) is the velocity of light (m/s), and \(\vartriangle \tau\) are DCB data obtained from the International GNSS Service (IGS) website ftp://igs.ign.fr/pub/igs/products/mgex/dcb/ in equation (5).

Planar method

The delay of each ionospheric grid point to be estimated can be fitted by the IPP near the IGP, and a local plane fitting equation is needed to determine the ionospheric delay. The plane fitting is given as follows:

$$I\left( x \right) = \mathop a\nolimits_{0} + \mathop a\nolimits_{1} \mathop x\nolimits_{{\text{E}}} + \mathop a\nolimits_{2} \mathop x\nolimits_{{\text{N}}} + r\left( x \right)$$
(6)

where \(x = \left( {\mathop x\nolimits_{{\text{E}}} ,\mathop x\nolimits_{{\text{N}}} } \right)\), \(\mathop x\nolimits_{{\text{E}}}\) and \(x_{N}\) are the longitude and latitude of IPP, \(I(x)\) is the vertical ionospheric delay at the IPP located in \(x = (x_{E} ,x_{N} )\), \(\mathop a\nolimits_{0}\), \(\mathop a\nolimits_{1}\) and \(\mathop a\nolimits_{2}\) are the coefficients of the fitted equation, and \(r(x)\) is the residual error estimation of the ionospheric delay at IPPs in the fitting domain.

Kriging method

It is a common spatial estimation method in geostatistics, which is suitable for the least mean square estimation of spatial data. Compared with the planar method, the kriging method based on the thin shell model can predict unknown points by using observations of random fields and can describe the random distribution characteristics of the ionosphere. The following are the specific calculation steps of the kriging method.

Step1: According to the principle of spatial statistics, points closer together have more similar attributes than points farther away. The spatial similarity is measured by creating a semivariogram. The IGP ionospheric delay is estimated by the IPP delay in the fitting domain (Sparks et al. 2011a, b; Ren and Yin 2020). Hence, it is essential to determine the fitting radius and calculate the distance between all IPPs and IGPs in the fitting area. These pairs of IPPs and IGPs must be grouped according to the step length in order to obtain a vertical ionospheric delay estimation variance to measure the difference. Finally, an admissible and widespread model for the semivariogram (Webster and Oliver 2008) is considered, i.e., the exponential model, which is given in Blanch (2004c):

$$\gamma \left( h \right) = \mathop C\nolimits_{0} + \mathop C\nolimits_{1} \left( {1 - \mathop e\nolimits^{{ - {\raise0.7ex\hbox{$h$} \!\mathord{\left/ {\vphantom {h a}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$a$}}}} } \right)$$
(7)

where \(h\) denotes the step length, \(\gamma \left( h \right)\) represents the semivariogram model, \(C_{0}\) is the nugget, that is, the magnitude of the intercept, and \(C_{0} = 0.04\;{\text{m}}^{2}\), \(C_{0} + C_{1}\) is the sill value that the semivariogram model attains at the range (the value on the y-axis), \(C_{1} = 2m^{2}\), and \(a\) is the variation distance. The above constants are given in Blanch (2004c).

Step2: The ionospheric delay of IGP is estimated by the ionospheric delay of IPPs in domain fitting given by:

$$\mathop z\nolimits^{*} \left( {\mathop x\nolimits_{0} } \right) = \sum\limits_{i = 1}^{n} {\mathop \lambda \nolimits_{i} } Z\left( {\mathop x\nolimits_{i} } \right)$$
(8)

where \(\mathop z\nolimits^{*} \left( {\mathop x\nolimits_{0} } \right)\) denotes the estimated value at position \(x_{0}\), \(Z(x_{i} )\) is the measured value at position \(x_{i}\), \(n\) is the number of measurements used to estimate the process, \(\lambda_{i}\) represents the weight of \(Z(x_{i} )\), and \(\sum\nolimits_{i = 1}^{n} {\lambda_{i} } = 1\) is to ensure unbiased estimation during calculation. \(\lambda_{i}\) can be determined according to Blanch (2004c).

Ionospheric tomography

MIDAS is chosen as the ionospheric tomographic method to estimate IGDs both at quiet times and storm times. First, the ionosphere over the inversion area is divided into voxels along latitude, longitude, and altitude, and the integral expression (1) is converted into a matrix form to derive the TEC. Second, the electron density in each voxel is inverted according to:

$${\mathbf{Hx = Z}}$$
(9)

where \({\mathbf{H}}\) is the vertical intercepts of the transmitted rays through the ionospheric grid, \({\mathbf{x}}\) denotes the ionospheric electron density, and \({\mathbf{Z}}\) represents the slant TEC along the propagation path.

The number of transmitted rays distributed within the ionospheric grid is greatly affected by the number of ground observation stations. In fact, the uneven distribution and sparsity of ground stations cause some voxels in the discretized ionosphere to have no signal observation information, resulting in a serious rank deficit of matrix \({\mathbf{H}}\) and the non-unique solution of (9). In MIDAS, this problem is usually solved by constructing the mapping matrix \({\mathbf{M}}\). Generally, equation (9) can be transformed into the following:

$${\mathbf{HMy = Z}}$$
(10)

where \({\mathbf{y}}\) is the solution of the transformation basis that can be derived from singular value decomposition (SVD).

MIDAS applies a variety of functions to construct the mapping matrix \({\mathbf{M}}\), for example, a set of empirical orthogonal functions (EOFs) is computed to map the radial profile, and a voxel based grid is bounded in latitude, longitude, and altitude (Spencer and Mitchell 2007). The international reference ionosphere (IRI) model 2012 (Bilitza et al. 2014) or Chapman functions (Rishbeth and Garriott 1969) are generally used to create EOFs. Two EOFs are generated from the IRI 2012 model in this work. Then, the solution of the electron density can be obtained from:

$${\mathbf{x = My}}$$
(11)

The estimation of the ionospheric delay of IGP in WAAS by MIDAS can follow two steps as below:

Step 1: Set the inversion area in the range of 10° N-62° N, and 50° W-130° W, which is usually larger than the study WAAS area, to keep good quality of tomographic results. The selection of MIDAS grids is highly dependent on the distribution of ground GPS data. In order to obtain accurate tomographic results, different MIDAS grids are chosen with respect to different coverage of ground data used for the inversion (Yin et al. 2017). In this study, we attempt to use a similar number and distribution of ground stations to those of WAAS to do the inversion so as to evaluate the performance of tomographic IGDs and WAAS kriging IGDs under similar circumstances. Therefore, a sparse configuration of receivers is selected for reconstruction over the study region, and the step in latitude, longitude, and altitude is set as 4°, 4°, and 40 km, respectively, and the time resolution is 30 min.

Step 2: Estimate the IGP ionospheric delay based on the distance difference between tomographic grid points and IGPs. The tomographic grid points are selected as close to the IGPs as possible. Therefore, in the study area, only the case where the difference between the latitude and longitude is not larger than 2° is considered. Suppose the difference between the longitude and latitude of the to-be-estimated point in the study area and its nearby tomographic grid point is within 1°. Then, the delay value of the tomographic grid point closest to the IGP is directly selected as the IGP delay. If the difference of the two points is within [1°, 2°], the mean value of the ionospheric delay corresponding to all tomographic grid points within the range of 2° difference is selected as the IGP delay.

Kriging-combined tomography

Since MIDAS grids are chosen differently from WAAS grids, as shown in Fig. 2, it is essential to add the kriging step to the tomography for the estimation of IGDs. In order to take advantage of both features of ionospheric tomography and the kriging approach, a method of tomography combined with kriging is proposed innovatively for the WAAS application. First, tomographic delays over the inverted area are obtained with MIDAS, then the inversion grid points are treated as the IPPs, and finally, the IGDs are estimated by the kriging method.

Fig. 2
figure 2

Distribution of WAAS stations (solid blue circles), six reference stations (red triangles), and 84 IGPs (small black points) in the study area

Full size image

Results and discussions

WAAS IGDs only in the range of 70° W-125° W and 20° N-50° N in the USA from September 7 to 9, 2017, are investigated in the study. The WAAS grid is divided as 5° × 5° in latitude and longitude, so there are a total of 84 IGPs in the study area. Figure 2 shows the distribution of those IGPs (small black points) and 20 WAAS stations (solid blue circles) used to estimate IGDs by the planar and kriging methods, as well as six reference stations (red triangles) selected for comparison.

TEC is first calculated by GPS dual-frequency carrier phase observations (Ren and Yin 2020). Then, the planar method, the kriging method, the tomographic method, and the kriging-combined tomographic method are applied to calculate IGP delays. Finally, the estimated results of the aforementioned four methods are compared with ‘true’ IGDs, which are calculated from dual-frequency carrier phase data at six evenly distributed WAAS reference stations.

By calculating the IPP vertical ionospheric delay at each reference station from September 7 to 9, the distances between all IPPs and nearby IGPs are calculated. The ionospheric delay at the IPPs closest to the IGP to be estimated is selected to evaluate mean absolute percentage error (MAPE) and root mean squared error (RMSE). These IGD estimates by four methods are compared with ‘true’ IGDs at six sites for 24 h of three days. The calculation of RMSE and MAPE is given by:

$$RMSE = \sqrt {{\raise0.7ex\hbox{${\sum\limits_{{i,j = 1}}^{{24}} {\left( {{\text{ }}x_{i} - {\text{ }}x_{j} } \right)} ^{2} }$} \!\mathord{\left/ {\vphantom {{\sum\limits_{{i,j = 1}}^{{24}} {\left( {{\text{ }}x_{i} - {\text{ }}x_{j} } \right)} ^{2} } {24}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${24}$}}}$$
(11)
$${\text{MAPE = }}\frac{1}{{24}}\sum\limits_{{i,j = 1}}^{{24}} {\frac{{\left| {{\text{ }}x_{i} - {\text{ }}x_{j} } \right|}}{{{\text{ }}x_{j} }}} * {\text{100}}\%$$
(12)

where \(x_{i}\) are the estimated IGDs of four estimation methods and \(x_{j}\) are the’true’ IGDs.

The IGD comparison results at six stations on three days are shown in Tables 1, 2, 3, 4, 5, 6. The six stations are coded with ZDC1, ZLA1, ZSE1, ZDC1, YWG1, and ZDV1, located in the east, west, west, south, and north of the study area, respectively, as illustrated in Fig. 2 with red triangles.

Table 1 IGD comparisons at YWG1 on September 7–9, 2017
Full size table
Table 2 IGD comparisons at ZDC1 on September 7–9, 2017
Full size table
Table 3 IGD comparisons at ZDV1 on September 7–9, 2017
Full size table
Table 4 IGD comparisons at ZHU1 on September 7–9, 2017
Full size table
Table 5 IGD comparisons at ZLA1 on September 7–9, 2017
Full size table
Table 6 IGD comparisons at ZSE1 on September 7–9, 2017
Full size table

From the results listed in the six tables, it can be seen that the overall IGD errors of the two tomographic methods are smaller than those of the planar method and the kriging method, and errors of the planar method in both RMSE and MAPE are the worst of four methods at six stations on three days. RMSEs of the planar method, the kriging method, the tomographic method, and the kriging-combined tomographic method at most stations decrease sequentially. The two tomographic methods do not show too much improvement compared with the kriging method. However, regarding MAPEs, two tomographic methods present a significant advantage over the planar and the kriging methods, especially on September 8 when the ionosphere became disturbed. MAPEs at all six stations are decreased by 2–4 times with two tomographic methods compared to the two methods based on the thin shell model, indicating the tomographic technique is better than the thin shell model at storm times.

Based on the analysis in the aforementioned, the ionosphere was relatively quiet on September 7 and 9 but became greatly disturbed on September 8. Hence, the IGD errors during the periods of strong disturbance will be evaluated together with TEC maps in the following section.

In order to further explore the impact of ionospheric disturbances on the estimation results of different methods, this section compares and analyzes estimation results at 6 stations during the most disturbed periods, that is, 0:00–3:00UT and 12:00–15:00UT on September 8. Since TEC variations between 12:00 and 15:00UT were not obvious, we focus the investigation on the tempo-spatial variations in the ionosphere over 0:00–3:00UT.

Figure 3 shows the distribution of GPS stations used for tomography and IPP trajectories at the height of 350 km between ground stations and GPS satellites over the study area from 1:00 to 2:00 UT on September 8, 2017. To keep consistent with WAAS stations, 32 stations, including 20 WAAS stations as shown in Fig. 2 and additional 12 CORS stations, are chosen for the reconstruction so that the coverage of IPPs is not too poor to make the inversion as needed.

Fig. 3
figure 3

Distribution of 32 GPS stations for the inversion area (solid red circles) and IPP trajectories at the height of 350 km between ground stations and GPS satellites over the study area from 1:00 to 2:00 UT on September 8, 2017. These IPP tracks are plotted to show the coverage of ground GPS data over the certain period

Full size image

Meanwhile, tomographic TEC maps between 0:00 and 3:00UT on three days are illustrated in Fig. 4. The left panels show TEC maps from 0:00 to 3:00UT on September 7, the central panels show TEC maps from 0:00 to 3:00UT on September 8, and the right panels show TEC maps from 0:00 to 3:00UT on September 9. When compared with these TEC maps, it is obvious that TEC increased a lot and large TEC gradients from northeast to southwest over time occurred on September 8 rather than on September 7 and 9.

Fig. 4
figure 4

Hourly TEC maps from 0:00 to 3:00 UT (top to bottom) on September 7–9 (left to right panels). The color bar shows the TEC range in TECU

Full size image

Next, when considered in combination with the locations of six reference stations in Fig. 2, TEC at ZSE1 and ZDV1 changed a lot between 0:00 and 2:00UT with a drop from about 23–25 TECU directly to about 4–6 TECU due to the large TEC gradient. Meanwhile, TEC at ZDC1 and YWG1 presented a large drop from 0:00 to 1:00UT. On the contrary, the impact of ionospheric activities at ZLA1 and ZHU1 was relatively small as TEC at both stations reduced smoothly over 3 h.

IGD estimation errors of four methods at each of the six stations between 0:00—3:00 UT and 12:00—15:00 UT on September 8 are shown in Fig. 5 regarding RMSE and MAPE. As the previous results of three days, errors of the planar method are the largest among those of the four methods, while errors of the kriging method are smaller than planar errors, confirming its advantage over the planar fitting in the operational WAAS. Furthermore, the tomographic method proposed in this study shows a great improvement over the planar and kriging methods since RMSEs and MAPEs of two tomographic methods at all six stations are much lower than those of the planar and kriging methods. For example, compared with the kriging method, MAPEs of the kriging-combined tomography method are reduced by 77% at YWG1, 72% at ZHU1, 66% at ZDV1, 55% at ZLA1, 44% at ZDC1, and 27% at ZSE1.

Fig. 5
figure 5

Estimation errors of four methods at the six stations between 0:00–3:00 UT and 12:00–15:00 UT on September 8, 2017. The top panel shows the variation in RMSE, and the bottom panel shows the variation in MAPE

Full size image

Obviously, the differences of MAPEs, as well as RMSEs between the kriging method and the kriging-combined tomography method, are small at the reference station ZSE1. Besides, MAPEs of all four methods at this station are greatest in all six stations as shown in Fig. 5, indicating that the accuracy of the four estimation methods is relatively poor. This may be related to the fact that the reference station ZSE1 is located in the western boundary of the study area, where a relatively sharp TEC gradient was experienced during the ionospheric disturbance. On the contrary, the station ZLA1 is also located along the west coast, but IGD errors of the tomographic methods there are not very large due to smooth TEC variations over time.

The kriging-combined tomographic method has a greater improvement at the remaining four stations than the kriging method alone. The RMSE and MAPE values of station ZHU1 are the smallest because the location of ZHU1 is in the southern part of the study area, where TEC changes relatively smoothly, as shown in Fig. 4. While ZDV1 is located in the middle of the study area, the RMSE value there is larger than that of ZHU1, resulting from large TEC gradients over the station. Since TEC also dropped quickly in one hour at YWG1 and ZDC1, the IGD estimation errors are a bit large.

For the comparison between the two tomographic methods, it can be found that the RMSEs of the kriging-combined tomography methods are lower than those of the tomographic method by 0.23 m, 0.37 m, 0.19 m, and 0.08 m at the reference stations YWG1, ZDC1, ZDV1, and ZHU1, respectively. In the mean of time, errors of the tomographic method are much less than those of the kriging method and the planar method, indicating that the kriging-combined tomography method can achieve the greatest improvement among these methods in reducing the IGD errors. However, at ZLA1 and ZSE1 the errors of the kriging-combined tomography are slightly worse than those of the tomographic method.

Figure 6 shows the mean error of the four estimation methods at all six stations between 0:00–3:00UT and 12:00–15:00UT on September 8. It can be clearly seen from the figure that RMSEs of the planar method, the kriging method, the tomographic method, and the kriging-combined tomographic method decrease in order. Compared with the kriging method, RMSE and MAPE of the kriging-combined tomographic method are reduced from 2.25 m to 1.53 m and 3.54% to 1.93%, respectively. However, in terms of MAPEs, the kriging-combined tomographic method seems not to work as well as the tomographic method by an increase of 0.2%, inferring that the tomographic method alone with less computation load may replace the kriging-combined tomographic method during the strong ionospheric disturbance.

Fig. 6
figure 6

Mean error of estimation results of different methods at all six stations between 0:00–3:00UT and 12:00–15:00 UT on September 8, 2017

Full size image

Conclusion

In order to estimate SBAS ionospheric grid delays, the planar method and the kriging method based on the thin-shell model have been applied separately in operational systems, such as WAAS. However, these methods may deteriorate against ionospheric tomography during disturbed ionospheric conditions. In this study, for the first time by virtue of electron density inversions, a tomographic method, as well as a kriging-combined tomographic method, is proposed innovatively to apply for estimating IGDs over part of the WAAS region during ionospheric disturbances on September 7–9, 2017. The planar method, kriging, tomographic, and kriging-combined tomographic methods are implemented based on WAAS GPS data to estimate the ionospheric grid point delays for these three days. The results of four methods are then compared with those of independent dual-frequency observations at six reference stations.

In general, the estimated errors of the planar method, the kriging method, the tomographic method, and the kriging-combined tomographic method are sequentially reduced. It is quite clear that the RMSEs and MAPEs of the planar method are the largest at six stations on all three days, indicating that the planar method is the worst among the four methods. Furthermore, two tomographic methods are performed better than those based on the thin-shell model, especially when the ionosphere was strongly disturbed on September 8. As far as MAPEs are concerned, the improvement in the two tomographic methods is 2–4 times that of the methods based on the thin-shell model, showing the significant advantage of tomography during ionospheric storms. At the same time, it is found that the improvement in the IGD estimate with the kriging-combined tomography is relatively small compared with the tomographic method alone.

In addition, the performance of four estimation methods over strong ionospheric disturbed times is evaluated. Two periods from 0:00 to 3:00 UT and from 12:00 to 15:00 UT on September 8 were chosen to conduct the study. The results show that the performance of the two tomographic methods is still better than the methods based on the thin shell model in estimating the ionospheric delay. Besides, the kriging-combined tomographic method still presents the best improvement in reducing IGD errors. For example, the RMSEs of the kriging-combined tomography methods are lower than those of the tomographic method by 0.23 m, 0.37 m, 0.19 m, and 0.08 m at the reference stations YWG1, ZDC1, ZDV1, and ZHU1, respectively. In comparison with the kriging method, the RMSE and MAPE of the kriging-combined tomographic method are reduced from 2.25 m to 1.53 m, and 3.54% to 1.93%, respectively.

However, in terms of MAPEs, the kriging-combined tomographic method seems not to work as well as the tomographic method by an increase of 0.2%, inferring that the tomographic method alone with less computation load might replace the kriging-combined tomographic method under strong ionospheric disturbance. When analyzing together with the change of TEC gradients between 0:00 and 3:00UT on September 8, it is noted that the kriging-combined tomography performs better than tomography alone where there is small TEC gradients, e.g., at ZHU1. Whereas the TEC gradient changes sharply, e.g., at ZSE1, the kriging-combined tomography tends to be slightly worse than the tomographic method alone. Finally, it should be clarified that the above performance of tomographic and kriging-combined tomographic methods is limited to a sparse configuration of receivers that is similar to that of WAAS, and if smaller grids could be used in MIDAS taking advantage of highly populated ground stations, better performance might be achieved and additional interpolation might not be needed.