Applicability Analysis of VTEC Derived from the Sophisticated Klobuchar Model in China

Although the Klobuchar model is widely used in single-frequency GPS receivers, it cannot effectively correct the ionospheric delay. The Klobuchar model sets the night ionospheric delay as a constant, i.e., it cannot reflect temporal changes at night. The observation data of seventeen International Global Navigation Satellite System Service (IGS) stations within and around China from 2011 provided by the IGS center are used in this study to calculate the Total Electron Content (TEC) values using the Klobuchar model and the dual-frequency model. The Holt–Winters exponential smoothing model is used to forecast the error of the 7th day between the Klobuchar model and the dual-frequency model by using the error of the former six days. The forecast results are used to develop the sophisticated Klobuchar model when no epochs are missing, considering that certain reasons may result in some of the observation data being missing and weaken the relationship between each epoch in practical applications. We study the applicability of the sophisticated Klobuchar model when observation data are missing. This study deletes observation data of some epochs randomly and then calculates TEC values using the Klobuchar model. A cubic spline curve is used to restore the missing TEC values calculated in the Klobuchar mode. Finally, we develop the sophisticated Klobuchar model when N epochs are missing in China. The sophisticated Klobuchar model is compared with the dual-frequency model. The experimental results reveal the following: (1) the sophisticated Klobuchar model can correct the ionospheric delay more significantly than the Klobuchar model; (2) the sophisticated Klobuchar model can reflect the ionosphere temporal evolution, particularly at night, with the correct results increasing with increasing latitude; and (3) the sophisticated Klobuchar model can achieve remarkable correction results when N epochs are missing, with the correct results being nearly as good as that of the dual-frequency model when no epochs are missing.


Introduction
Satellite signals passing through the ionosphere may be influenced by the ionosphere.The ionospheric error of a signal may vary from meters to tens of meters.This error is one of the main positioning errors in the Global Navigation Satellite System (GNSS) [1].The Total Electron Content (TEC) is the key element that causes ionospheric delay and is an important parameter for satellite navigation [2].Because the ionosphere is a dispersive medium, when electromagnetic wave signals pass through the ionosphere, the delay is related to the frequency of the signal and electron density distribution.The ionospheric delay error can be eliminated via the ionosphere-free combination of two or more frequencies [3,4].However, this requires receivers that can accept a dual-frequency signal and even multi-frequency signals from a satellite.Regarding Global Positioning System (GPS) commercial users, their receivers are mainly single-frequency receivers.For a single-frequency GPS receiver, different types of models or algorithms can be used to correct the ionospheric delay.These models are applied to the receivers and can be classified as follows: (1) Broadcast Ionospheric Models (BIMs), which include the Klobuchar model (BKM) in GPS [5], the NeQuick model in Galileo [6], the BeiDou satellite navigation System (BDS) ionospheric model [7], and the ionospheric grids that are broadcast by the Satellite-Based Augmentation System (SBAS) [8]; (2) Empirical Ionospheric Models (EIMs), such as the International Reference Ionosphere (IRI) model [9] and the Parameterized Ionospheric Model (PIM) [10]; and (3) Global Ionospheric Maps (GIM) provided by each International Global Navigation Satellite System Service (IGS) analysis center [11].These models use observational data to establish a model to calculate the TEC all over the world.
The BIM can calculate ionospheric delay data easily for a broadcast ephemeris, satellite elevation, and satellite azimuth only.Because BIMs can compute the ionospheric delay easily and conveniently, they are widely used in single-frequency receivers.In particular, BKM is frequently used, although its correcting accuracy is not high [12][13][14].The EIM can calculate the ionospheric delay using parameters of the previous few days that can influence the ionosphere, such as the activity of the sun and geomagnetic activity [15].However, these parameters are not easily obtained by the users.Moreover, the EIM cannot provide the real-time ionospheric delay.This limits the usage of the EIM in single-frequency receivers.The GIMs are provided by the IGS center.The final product may be delayed by almost 14 days.Even the rapid product is delayed by two days, and the mean accuracy is 2-8 TECu worldwide [16,17].Alternatively, some IGS analysis centers can provide a two-day ahead prediction product [18].However, this product cannot satisfy the requirement of high-accuracy positioning if the area is lacking an IGS station [19].We conclude that BIMs are models that are easy to use, effective, and widely used in practice [20].
BKM was proposed for the GPS system by J. A. Klobuchar in 1987.It is used to calculate the ionospheric delay and can reflect a delay of 50%-60% in the mid-latitude.It is best applied in mid-latitude.The model uses the eight parameters in the broadcast ephemeris to calculate the ionospheric delay.Moreover, the model treats the delay values as the positive part of a cosine function during the daytime.The peak of the ionosphere is set at 14:00 Local Time (LT).The delay at night is set as a constant 5 ns.The model's accuracy is limited due to the following two reasons.First, a table of predefined values is used to update the eight coefficients in the GPS system.These values are predefined for a particular condition in the solar cycle.They are not set to update every day.Second, the diurnal variation of the ionosphere is calculated using the eight parameters in the broadcast ephemeris, and setting the delay values at night as a constant makes the model fail to model the temporal changes, particularly at night [21].Scholars have attempted to find a way to solve these two primary problems.
For the first approach, the Center for Orbit Determination in Europe (CODE) [22] fits the eight parameters again and releases them to the internet.Single-frequency users can use this approach to address the baseline easily and achieve a higher accuracy of ionospheric delay.Yuan et al. [23] refined the eight parameters using data from hundreds of IGS stations and the Crust Movement Observation Network of China (CMONOC).This approach improved the accuracy by nearly 15% in the China region when compared with BKM.Moreover, the refined eight parameters method of BKM is also used in BDS regional ionospheric correction [24].The Indian Regional Navigation Satellite System (IRNSS) [25] broadcasts to the single-frequency users a set of coefficients updating measurements over India.Although the method to refine the eight parameters using satellite observation data can achieve significant correction results in real applications, it only considers the diurnal variation during the daytime.However, the variation at night may still result in a 20%-30% positioning error and is not suitable over China.
To correct the delay at night, Filjar et al. [26] established a Klobuchar-like model to correct the delay at night.They used a linear fit to simulate the delay at night.Moreover, the changing of amplitude takes seasonal change into account.This approach can better correct the ionospheric delay in the south Adriatic.Moreover, Han et al. [27] established a Klobuchar-like model with fourteen parameters.Two of the parameters are used to calculate the vertical ionospheric delay at night.This approach provides a certain degree of improvement at night.These methods to improve BKM can validly correct the delay at night.However, their overall accuracy must be further improved.
With the improvement of technology, single-frequency users are eager to improve their positioning accuracy.The methods mentioned above can better correct the ionospheric delay overall, but they are limited in their use in the receivers because of their low correction accuracy and complex methods.To address these issues, this study treats the bias between BKM and the dual-frequency observation model (DM) as a time series.The time-series analysis method is used to predict the bias.The predicted bias is used to improve BKM.The observational data of the IGS stations within and around China are used to test the applicability of the sophisticated BKM.Moreover, this study will discuss the applicability of the model in China in the case of some missing epochs.

Dual-Frequency Observation Model
The GPS-TEC can be obtained via pseudo-range and carrier phase measurements.It is an important way to estimate the electron density in navigation systems.By computing the differential of the code and carrier phase measurements, dual-frequency GPS receivers can provide integral information regarding the ionosphere.Thus, the TEC derived from the dual-frequency receivers is proposed as an input to an assimilative model of the ionosphere.This is why the TEC data derived from the dual-frequency receivers have been utilized for this study and treated as true values.At the same time, we can learn the following from Ref. [28].The observation data from the pseudo-range are unambiguous and can yield absolute TEC values.However, they are affected by large noise, which limits their widespread use in practice.On the other hand, the observation data from the carrier phase are ambiguous and yield relative TEC values, although they contain less noise.This limits their widespread use in real engineering practice, as well as most scientific applications.To solve this problem, the pseudo-range noise is reduced by smoothing pseudo-range data using carrier phase measurements.This technique is called carrier phase smoothing [29].
Hatch initially proposed a widely used algorithm to smooth the code pseudo-range using carrier phase measurements [30].An epoch-dependent smoothing weight factor (SWF) is further introduced to improve this algorithm [31].A modification is further made to this algorithm by reducing the SWF by a constant value with the increase of the number of GPS data epochs used for smoothing [32].Compared to the TEC data not being smoothed, the smoothed TEC data typically show an improvement of a few TECu in precision.The algorithm of this model is described as follows: In Equation (1), P(t i ) and L(t i ) are the code pseudo-range and carrier phase measurements at epoch t i , respectively; ω i is the smoothing weight factor, with an initial value of 1 and varying with [ 0 ∼ 1]; P(t i ) sm is the resultant smoothed code measurement at epoch t i .In this investigation, the weight factor ω i decreases with the extension of smooth time.The decrease rate can be confirmed using experimental results.The decrease rate is usually 0.01 when i > 100 and ω i = 0.01 and never changes again.
To evaluate the accuracy of the smoothed code pseudo-range measurement, the expanded analytical form of Equation (1) must be developed.Starting from the initial condition of P(t 1 ) sm = ω 1 P(t 1 ) at the first epoch, the following expanded expression can be derived for the smoothed code measurement at epoch t i : Once the smoothed code pseudo-range measurements at the L1 and L2 band are obtained, ionospheric total electron content can be derived as follows [33]: In Equation (3), TEC(t i ) is the TEC at epoch t i ; P(t i ) 1,sm and P(t i ) 2,sm are the smoothed code pseudo-range measurements for L1 and L2 frequencies, respectively; γ = ( f 1 / f 2 ) 2 = (1575.41/1227.6) 2 ; differential code biases (DCBs) are inter-frequency electronic range delays produced in the instrument of the receivers (DSB R ) and the satellite (DCB S ) and can be removed via the method presented in Ref. [34].
Moreover, the TEC derived from the dual-frequency is slant TEC (STEC) along the signal path between the satellite and the receiver on the ground.To free the STEC from dependence on the ray path geometry from the satellite to the receiver through the ionosphere, it is necessary to convert STEC to vertical TEC (VTEC).At the same time, VTEC is a more important parameter for characterizing the TEC over a given receiver position and is used as a good indicator of the overall ionization of the Earth's ionosphere [35].Hence, Equation ( 4) is used to convert STEC to VTEC by assuming that the ionosphere is equivalent to a thin shell encircling the earth, with its center the same as that of the Earth [36].In terms of the zenith angle χ at the Ionospheric Piecing Point (IPP) and the zenith angle χ at the receiver position on the ground, the relationship between STEC and VTEC can be given as: where Here, R e is the Earth's radius in kilometers; h m is the height of the maximum electron density at the F2 peak, which usually ranges from 250 km to 350 km at mid-latitudes and from 350 km to 500 km at the equator [37,38].In this study, the height of the maximum electron density, h m = 350 km, was chosen because most of China is at mid-latitudes and, at this height, the ionosphere is assumed to be spatially uniform.

Klobuchar Model
BKM is formed using the statistical properties of the variation of the ionosphere.The algorithm of BKM [5] is given in Equation ( 6): In Equation ( 6), I(t) denotes the vertical ionospheric delay value, in units of seconds; A 1 denotes the constant value at night, with a value of 5 ns; A 2 denotes the amplitude of the cosine function during the day and can be calculated from α n in the broadcast ephemeris (Equation ( 7)); A 3 denotes the initial phase and is always the time of the extremum value (here, the time is 14:00 LT (50,400 s)); A 4 denotes the cycle of the cosine function and can be calculated from β n in the broadcast ephemeris (Equation ( 8)); and t denotes the local time of the intersection point between the ionosphere and the line between the satellite and the receiver: In Equations ( 7) and ( 8), φ M denotes the geomagnetic latitude of IPP, in units of semicircles.

Holt-Winters Exponential Smoothing Model to Sophisticated Klobuchar Model
To improve the usage of BKM over China, we assume the differences of a certain day between BKM and DM influenced by the previous few days; the recent data influences the day's data more than the forward data.For this reason, the Holt-Winters exponential smoothing model was chosen to sophisticate BKM.
The Holt-Winters exponential smoothing model [39] was developed by Holt in 1957.The model has been widely used in many other areas [40][41][42].The model assumes that all known data can influence the data that are to be predicted.However, the recent data influence the prediction data more than the forward data.The influence varies in a geometric series as time changes.The non-seasonal model expression is as follows: In Equation ( 9), S t is the stable composition at time t; X t is the observation at time t; b t is the trend composition at time t; α, β ∈ [0, 1] are smoothing parameters; ∆ t i+m denotes prediction values in the latter m period; and m is the prediction period.Moreover, we suggest that readers read Refs.[40][41][42] to learn more about the Holt-Winters model.
After performing many repeated experiments, this study finally uses the former six days of observation data provided by the IGS stations (ftp://cddis.gsfc.nasa.gov) to calculate the VTEC values using BKM and DM.First, the differences between BKM and DM at each epoch of the previous six days are calculated.Second, the Holt-Winters exponential smoothing model is used to forecast the differences at each epoch of the seventh day.Last, the prediction data of each epoch are used to develop the sophisticated Klobuchar model (SKM).The algorithm used is given by Equation (10): Here, Equation (11) shows how ∆ t i is derived from the Holt-Winters model.∆ t denotes the differences between BKM and DM at each epoch of the previous six days; m = 1; and B t is the trend composition at time t.Both ∆ t and B t can be derived from Equation (6).
Moreover, the correct values ∆ t are divided into ∆ t 1 and ∆ t 2 , which vary as time changes.Here, we divided one day into two time periods, as in BKM.The first time period is |t − denotes the predicted biases of the first; and ∆ t 2 denotes the predicted biases of the second.At the same time, the total number of ∆ t 1 and ∆ t 2 in a day is 2880, combining them as a time series.Meanwhile, in practice, ∆ t 1 and ∆ t 2 can be updated once a day.The other algorithms are the same as those of the Klobuchar model.

Results and Discussion
To establish the sophisticated Klobuchar model over China, this study chooses ten IGS stations within China and seven IGS stations around China (Figure 1).The observation data are from the day of the year (DOY) 095 to 101 in 2011.The sophisticated method of each station is mentioned in Section 2.3.Moreover, this study divides one day into six time periods.Each period contains 4 h (480 epochs).We computed statistical parameters over time periods of 4 h at each station.For each time period, we compute the average bias (Equation ( 12)) between BKM or SKM and the DM at each station, the average corrective rate (Equation ( 13)) between BKM or SKM and the DM at each station, and the root mean square error (RMSE, Equation ( 14)) at each station to test the corrected result.The unit of bias is TECu, the unit of the corrective rate is percentage (%), and the unit of RMSE is TECu.At the same time, a continuous curvature gridding algorithm [43] mapping technique is used to obtain the VTEC over all of China from the isolated point.It is a method that is widely used in earth science and can retain the basic information of each isolated point.

Result of Modeling in China
The chosen data are used to establish SKM (mentioned in Section 2.3) in the China region.To understand the improvement of SKM, statistical analysis is performed on the average bias between BKM or SKM and DM for DOY 101, 2011; the result is shown in Figure 2. Hence, Figure 2a-c are the distributions of the average values of DM, BKM and SKM in the China region for DOY 101, 2011, respectively.Figure 2d shows that the bias between BKM and DM in China is between −4 TECu and 7 TECu.Among these biases, the largest bias is approximately 6.76 TECu, which appears at the XIAN

Result of Modeling in China
The chosen data are used to establish SKM (mentioned in Section 2.3) in the China region.To understand the improvement of SKM, statistical analysis is performed on the average bias between BKM or SKM and DM for DOY 101, 2011; the result is shown in Figure 2. Hence, Figure 2a-c are the distributions of the average values of DM, BKM and SKM in the China region for DOY 101, 2011, respectively.Figure 2d shows that the bias between BKM and DM in China is between −4 TECu and 7 TECu.Among these biases, the largest bias is approximately 6.76 TECu, which appears at the XIAN station.The lowest bias is approximately −3.38 TECu, which appears at the KUNM station.Figure 2d also indicates that BKM has a better correction accuracy in the mid-latitude, but it still needs to be improved; otherwise, the accuracy cannot meet the requirement of high accuracy positioning.Figure 2e shows that the bias of SKM in China is between 0 TECu and 1 TECu.Among these biases, the largest bias is approximately 0.94 TECu, which appears at the GUAO station.The lowest bias is approximately 0.59 TECu, which appears at the KUNM station.Figure 2e indicates that the biases in China improved significantly compared to Figure 2d.Thus, the approach is effective in improving BKM.Moreover, this approach can significantly improve the correction result in the ionosphere.

Result of Modeling in China, Considering Missing Epochs
In Section 3.1, we did not change the sample interval of observation data, with the value remaining at 30 s.Thus, a certain relationship exists between different epochs.However, in practice, some situations may cause the receivers to be unable to accept the signal from the satellite.This weakens the relationship between different epochs and reduces the accuracy of the correcting result, although SKM can achieve a rather good, correct result.This study will consider the applicability of the model for the case of N epochs missing in China.
To determine the accuracy that SKM can achieve while considering epochs missing in the observation data, we consider the following cases: random deletion of 480 epochs (480), 960 epochs (960), 1440 epochs (1440), 1920 epochs (1920), and 2400 epochs (2400).The elimination of random epochs is carried out for each day.For example, random deletion of 480 epochs (480) means that these random 480 epochs are missing for all six previous days.Moreover, cubic spline interpolation is a special case of spline interpolation that is used very often to avoid the problem of Runge's phenomenon.It gives an interpolating polynomial that is smoother and has a smaller error than some other interpolating polynomials.Therefore, a cubic spline curve [44] is used to restore the missing TEC values calculated by BKM.Next, the TEC values in China are determined using the method

Result of Modeling in China, Considering Missing Epochs
In Section 3.1, we did not change the sample interval of observation data, with the value remaining at 30 s.Thus, a certain relationship exists between different epochs.However, in practice, some situations may cause the receivers to be unable to accept the signal from the satellite.This weakens the relationship between different epochs and reduces the accuracy of the correcting result, although SKM can achieve a rather good, correct result.This study will consider the applicability of the model for the case of N epochs missing in China.
To determine the accuracy that SKM can achieve while considering epochs missing in the observation data, we consider the following cases: random deletion of 480 epochs (480), 960 epochs (960), 1440 epochs (1440), 1920 epochs (1920), and 2400 epochs (2400).The elimination of random epochs is carried out for each day.For example, random deletion of 480 epochs (480) means that these random 480 epochs are missing for all six previous days.Moreover, cubic spline interpolation is a special case of spline interpolation that is used very often to avoid the problem of Runge's phenomenon.It gives an interpolating polynomial that is smoother and has a smaller error than some other interpolating polynomials.Therefore, a cubic spline curve [44] is used to restore the missing TEC values calculated by BKM.Next, the TEC values in China are determined using the method described in Section 2.3.
The statistical mean bias is determined for DOY 101, 2011 between BKM or SKM and DM when considering N epochs missing for each of the previous six days in China.The respective results are shown in Figures 3 and 4.This may be caused by the Equatorial Ionization Anomaly [45].
We can learn from Figures 3 and 4 that the correction result of SKM for China is significantly improved when compared to that of BKM.However, this approach cannot reflect the ionospheric correction result of each time period.Thus, this study counts the correction results from the IRKJ station (LT = UT + 7), the CHAN station (LT = UT + 8), the LHAZ station (LT = UT + 6) and the PIMO  Figure 4 suggests that the biases of SKM change little with the increase in the number of missing epochs, regardless of how much of the epoch is missing for each of the previous six days.The biases change in the range of −5 TECu to 1 TECu.In this range, the largest is approximately 0.50 TECu, which was observed at the GUAO station.The lowest is approximately −4.33 TECu, which was observed at the KUNM station.The biases of other stations range from approximately −1 TECu to 0.50 TECu.The correction results of SKM are significantly improvement compared to that of BKM, regardless of how much of each epoch is missing when compared with the results of Figure 3.Moreover, the absolute biases in LHAZ, KUNM and TWTF are larger than those of any other stations.This may be caused by the Equatorial Ionization Anomaly [45].
We can learn from Figures 3 and 4 that the correction result of SKM for China is significantly improved when compared to that of BKM.However, this approach cannot reflect the ionospheric correction result of each time period.Thus, this study counts the correction results from the IRKJ station (LT = UT + 7), the CHAN station (LT = UT + 8), the LHAZ station (LT = UT + 6) and the PIMO (LT = UT + 8) station when N epochs are missing for each of the previous six days.Figure 5 shows the statistical results.(LT = UT + 8) station when N epochs are missing for each of the previous six days.Figure 5 shows the statistical results.Figure 5 indicates that SKM fits the data of DM well in both the low-latitude and the mid-latitude when no epoch is missing.This shows that SKM can better reflect the temporal variation of the ionosphere, particularly at night.Moreover, the biases between SKM and DM are approximately ±3 TECu in many situations.SKM has a relatively high accuracy with DM.The biases increase as the latitude decreases.In the low-latitude, the biases are larger than that in the mid-latitude.When N epochs are missing in the observation data, the ionospheric correcting accuracy of SKM can also fit the data of DM well in both the low-latitude and the mid-latitude.Moreover, SKM can better reflect the temporal variation, particularly at night.This indicates that the correction accuracy of SKM with N epochs missing is not lower than that of SKM with no epochs missing and can better reflect the temporal variation of the ionosphere, despite the missing epoch.Otherwise, it has little change during the maximum; this behavior may be caused by the exquisite different temporal variations of the ionosphere near the maximum.This influences the forecasting results of the bias and makes the fitted values of SKM approximations near the maximum.
To further show the correction results of SKM with no epochs missing and with N epochs missing, this study counts the average correction results of BKM, SKM with no epoch missing and SKM with N epochs missing for each time period for the IRKJ station, the CHAN station, the LHAZ station and the PIMO station.The statistical results are shown in Table 1.
Table 1 indicates that the correction results of SKM are significantly improved compared to those of BKM when no epochs are missing.Moreover, the ionospheric correction results in the mid-latitude are better than those in the low-latitude.Moreover, the ionospheric correction results increase with the increase of the latitude.In the low-latitude area (PIMO), the corrective rate of SKM with no epochs missing improved by nearly 30% over BKM, while the RMSE improved by nearly 3 TECu during the day.At night, the corrective rate improved by nearly 30%, while the RMSE improved by nearly 5 TECu.In the mid-latitude area (IRKJ, CHAN, LHAZ), the corrective rate of SKM with no epochs missing improved by nearly 20% compared to BKM, while the RMSE improved by nearly 6 TECu Figure 5 indicates that SKM fits the data of DM well in both the low-latitude and the mid-latitude when no epoch is missing.This shows that SKM can better reflect the temporal variation of the ionosphere, particularly at night.Moreover, the biases between SKM and DM are approximately ±3 TECu in many situations.SKM has a relatively high accuracy with DM.The biases increase as the latitude decreases.In the low-latitude, the biases are larger than that in the mid-latitude.When N epochs are missing in the observation data, the ionospheric correcting accuracy of SKM can also fit the data of DM well in both the low-latitude and the mid-latitude.Moreover, SKM can better reflect the temporal variation, particularly at night.This indicates that the correction accuracy of SKM with N epochs missing is not lower than that of SKM with no epochs missing and can better reflect the temporal variation of the ionosphere, despite the missing epoch.Otherwise, it has little change during the maximum; this behavior may be caused by the exquisite different temporal variations of the ionosphere near the maximum.This influences the forecasting results of the bias and makes the fitted values of SKM approximations near the maximum.
To further show the correction results of SKM with no epochs missing and with N epochs missing, this study counts the average correction results of BKM, SKM with no epoch missing and SKM with N epochs missing for each time period for the IRKJ station, the CHAN station, the LHAZ station and the PIMO station.The statistical results are shown in Table 1.  1 indicates that the correction results of SKM are significantly improved compared to those of BKM when no epochs are missing.Moreover, the ionospheric correction results in the mid-latitude are better than those in the low-latitude.Moreover, the ionospheric correction results increase with the increase of the latitude.In the low-latitude area (PIMO), the corrective rate of SKM with no epochs missing improved by nearly 30% over BKM, while the RMSE improved by nearly 3 TECu during the day.At night, the corrective rate improved by nearly 30%, while the RMSE improved by nearly 5 TECu.In the mid-latitude area (IRKJ, CHAN, LHAZ), the corrective rate of SKM with no epochs missing improved by nearly 20% compared to BKM, while the RMSE improved by nearly 6 TECu during the day.At night, the corrective rate improved by nearly 30%, while the RMSE improved by nearly 8 TECu.With some epochs missing, the ionospheric correction results of SKM are nearly identical to that of SKM with no epochs missing in both the low-latitude and the mid-latitude, regardless of how many epochs are missing.Among them, the corrective rate changes between −0.5% and 0.5%, with few of the periods changing by more than −1% to 1%.The RMSE changes between −0.1 TECu and 0.1 TECu, with a few of the periods changing by more than −0.2 TECu to 0.2 TECu.
The corrective rate of the low-latitude is better than that of the mid-latitude; however, the low-latitude RMSE is worse than that of the mid-latitude during the maximum values' time period.This phenomenon is shown in Table 1.This may occur because of the direct influence of the sun.With this influence, the VTEC over stations in the low-latitude areas is larger than that in the mid-latitude areas.This makes the biases between BKM and SKM larger than that in the mid-latitude.Thus, the corrective rate in the low-latitude is higher, whereas the RMSE value is smaller.Otherwise, a negative correction rate value occurs at night, as indicated in Table 1.This may be caused by the constant ionosphere delay value at night.It can also indicate that the constant ionospheric delay value at night cannot accurately reflect the temporal variation at night.We can also observe from Table 1 that the ionospheric correction results of SKM with N epochs missing are nearly identical to those of SKM with no epochs missing; however, a slight difference remains.This indicates that although the cubic spline curve method has its advantage of being time invariant, the invariance results still exhibit little difference from the observation data.When an epoch is missing, the difference may result in the correction results of SKM with an epoch missing and SKM with no epochs missing to exhibit a bias.The bias may be random.
Recently, many GPS/GNSS users have used dual-frequency receivers or even multi-frequency receivers, and the Continuously Operating Reference Stations (CORS) has been widely applied in many cities.This situation makes it possible and easy to obtain dual-frequency observation data.With these data, it is convenient to calculate the TEC values.Moreover, the corrective bias can be updated once a day.Thus, when users use single-frequency receivers to perform static measurements with a short baseline, observation data are missing because of the sudden change of the surroundings or solar activity.The cubic spline curve method can be used to repair the observation data.With this approach, the missing TEC values can be restored.SKM can be used as a post-correction treatment to address the ionospheric delay value to improve the baseline calculating precision.Moreover, when the single-frequency receivers log in to CORS, the single receivers can obtain the correction value in real time.This approach can reduce the workload of receivers and improve the work efficiency.

Conclusions
This study used the observation data from DOY 095 to 101 in 2011 provided by IGS stations within and around the China region to obtain an improved model.In this approach, BKM and DM are used to calculate the TEC values of each day.The TEC value of the previous six days is regarded as input data.The Holt-Winters exponential smoothing model is used to forecast the biases between BKM and DM of the seventh day.The forecast results are used to improve BKM.Moreover, some epochs may be missing in practice because of particular situations, such as sudden change of the surroundings or the activity of the sun.Thus, some epochs are randomly deleted.After we compared SKM with no epochs missing to SKM with an epoch missing with DM via experiments and analyzed the results, we drew the following conclusions: (1) SKM with no epoch missing and SKM with an epoch missing can better fit the trend of DM.SKM can better reflect the temporal changes of the ionosphere, especially at night.(2) The correction results of the ionospheric delay increase with increasing latitude.All correction results of the ionospheric delay are improved compared to those of BKM during the day.However, at night, the correction results are significantly improved compared to those of BKM.This finding indicates the validity of this method to improve the BKM.(3) The cubic spline curve method is used to repair the missing epoch data when the observation data have an epoch missing in some situations.Then, SKM is improved.It can also yield a remarkable result.The correction results of ionospheric delay are nearly identical to those of SKM with no epochs missing, with only a slight difference remaining.This difference may be random.
The analysis presented here shows the results of China over one day.Our conclusions may be limited by the time period and the area.We plan to extend the same analysis to a larger time period and a larger area to validate extensively the SKM algorithm performances.Moreover, we will analyze the impact caused by the update of the eight GPS broadcast coefficients, particularly for the case when there is large variation of these coefficients.
), i = 1 denotes BKM; i = 2 denotes SKM; K i j denotes the VTEC value of BKM or SKM at the j-th epoch; D j denotes the VTEC value of DM at the j-th epoch; K i j − D j denotes the absolute error between the two models at the j-th epoch; n denotes the initial epoch; and m denotes the final epoch.ISPRS Int.J. Geo-Inf.2017, 6, 75 6 of 14 ), i = 1 denotes BKM; i = 2 denotes SKM; denotes the VTEC value of BKM or SKM at the j-th epoch; denotes the VTEC value of DM at the j-th epoch; | − | denotes the absolute error between the two models at the j-th epoch; n denotes the initial epoch; and m denotes the final epoch.

Figure 1 .
Figure 1.The distribution of IGS stations within and around China used to establish the sophisticated model.

Figure 1 .
Figure 1.The distribution of IGS stations within and around China used to establish the sophisticated model.

14 Figure 2 .
Figure 2. (a) distribution of average VTEC value of DM in the China region for DOY 101, 2011; (b) distribution of average VTEC value of BKM in the China region for DOY 101, 2011; (c) distribution of average VTEC value of SKM in the China region for DOY 101, 2011; (d) distribution of average bias between BKM and DM in the China region for DOY 101, 2011; and (e) distribution of average bias between SKM and DM in the China region.The color denotes the degree of the bias for DOY 101, 2011.

Figure 2 .
Figure 2. (a) distribution of average VTEC value of DM in the China region for DOY 101, 2011; (b) distribution of average VTEC value of BKM in the China region for DOY 101, 2011; (c) distribution of average VTEC value of SKM in the China region for DOY 101, 2011; (d) distribution of average bias between BKM and DM in the China region for DOY 101, 2011; and (e) distribution of average bias between SKM and DM in the China region.The color denotes the degree of the bias for DOY 101, 2011.

Figure 3 .
Figure 3. Distribution of mean biases for DOY 101, 2011 between BKM and DM in China when (a) 480 epochs; (b) 960 epochs; (c) 1440 epochs; (d) 1920 epochs; and (e) 2400 epochs are missing for each of the previous six days.The color denotes the degree of bias.

Figure 3
Figure 3 indicates that the distributions of biases between BKM and DM with different epochs missing in China are nearly identical.The biases are between −4 TECu and 7 TECu.Among these biases, the largest bias is approximately 6.80 TECu, which appears at the XIAN station.The lowest bias is approximately −3.40 TECu, which appears at the KUNM station.

Figure 3 .
Figure 3. Distribution of mean biases for DOY 101, 2011 between BKM and DM in China when (a) 480 epochs; (b) 960 epochs; (c) 1440 epochs; (d) 1920 epochs; and (e) 2400 epochs are missing for each of the previous six days.The color denotes the degree of bias.

Figure 4
Figure 4 suggests that the biases of SKM change little with the increase in the number of missing epochs, regardless of how much of the epoch is missing for each of the previous six days.The biases change in the range of −5 TECu to 1 TECu.In this range, the largest is approximately 0.50 TECu, which was observed at the GUAO station.The lowest is approximately −4.33 TECu, which was observed at the KUNM station.The biases of other stations range from approximately −1 TECu to 0.50 TECu.The correction results of SKM are significantly improvement compared to that of BKM, regardless of how much of each epoch is missing when compared with the results of Figure 3.Moreover, the absolute biases in LHAZ, KUNM and TWTF are larger than those of any other stations.This may be caused by the Equatorial Ionization Anomaly[45].We can learn from Figures3 and 4that the correction result of SKM for China is significantly improved when compared to that of BKM.However, this approach cannot reflect the ionospheric correction result of each time period.Thus, this study counts the correction results from the IRKJ station (LT = UT + 7), the CHAN station (LT = UT + 8), the LHAZ station (LT = UT + 6) and the PIMO

Figure 3
Figure 3 indicates that the distributions of biases between BKM and DM with different epochs missing in China are nearly identical.The biases are between −4 TECu and 7 TECu.Among these biases, the largest bias is approximately 6.80 TECu, which appears at the XIAN station.The lowest bias is approximately −3.40 TECu, which appears at the KUNM station.Figure 4 suggests that the biases of SKM change little with the increase in the number of missing epochs, regardless of how much of the epoch is missing for each of the previous six days.The biases change in the range of −5 TECu to 1 TECu.In this range, the largest is approximately 0.50 TECu, which was observed at the GUAO station.The lowest is approximately −4.33 TECu, which was observed at the KUNM station.The biases of other stations range from approximately −1 TECu to 0.50 TECu.The correction results of SKM are significantly improvement compared to that of BKM, regardless of how much of each epoch is missing when compared with the results of Figure 3.Moreover, the absolute ISPRS Int.J. Geo-Inf.2017, 6, 75 10 of 14

Figure 5 .
Figure 5. Daily variation of vertical TEC at some stations.Top left: CHAN station.Top right: IRKJ station.Bottom left: PIMO station.Bottom right: LHAZ station.Blue lines represent the VTEC values of BKM; black lines represent the VTEC values of SM; red lines represent the VTEC values of SKM with no epochs missing; green lines represent the VTEC values of SKM with 480 epochs missing; purple lines represent the VTEC values of SKM with 960 epochs missing; bottle green lines represent the VTEC values of SKM with 1400 epochs missing; brown lines represent the VTEC values of SKM with 1920 epochs missing; pink lines represent the VTEC values of SKM with 2400 epochs missing.

Figure 5 .
Figure 5. Daily variation of vertical TEC at some stations.Top left: CHAN station.Top right: IRKJ station.Bottom left: PIMO station.Bottom right: LHAZ station.Blue lines represent the VTEC values of BKM; black lines represent the VTEC values of SM; red lines represent the VTEC values of SKM with no epochs missing; green lines represent the VTEC values of SKM with 480 epochs missing; purple lines represent the VTEC values of SKM with 960 epochs missing; bottle green lines represent the VTEC values of SKM with 1400 epochs missing; brown lines represent the VTEC values of SKM with 1920 epochs missing; pink lines represent the VTEC values of SKM with 2400 epochs missing.

Table 1 .
The statistics of the corrective rate and RMSE for each time period with BKM, SKM and SKM with N epochs missing for each of the previous six days.