An End-to-end Approach for Near Real Time Ionosphere Monitoring Over Mid-latitudes from GPS Data using Krigning Interpolation and IGS Products
Elhadi Takka, Microwaves Laboratory, Polytechnic Military School, Algeria; Aichouch Belhadj Aissa, USTHB University, Algeria; Hocine Kimouche, Microwaves Laboratory, Polytechnic Military School, Algeria
Ionospheric variability and disturbances affect the performance of Global Navigation Satellite Systems (GNSS). Different applications request increasingly for real-time information about the state of the ionosphere, predictions and warnings. This study presents an end-to-end approach for ionosphere monitoring in near-real time using hourly GPS data over mid-latitudes to routinely map Total Electron Content (TEC). The generated TEC maps provide valuable information about the ionosphere behavior and its evolution.
2.1-Receiver DCBs estimation
The developed software package starts from receiver DCBs (Differential Code Bias) estimation of non-IGS stations using rapid GIM (Global Ionosphere Model) from IGS (International GNSS Service) and nighttime GPS observations. We assume that day-to-day DCBs variations are negligible and we choose Rapid IGS GIM with one-day latency as a priori model. Basing on this assumption, satellite DCBs and IGS stations DCBs are extracted from GIM and we consider them as valid for the day of interest. Only measurements observed during nighttime (2000-0300LT) are used considering daytime ionosphere has larger variability that may generate larger uncertainty in the DCB estimation. Then, receiver DCBs of non-IGS stations are estimated every measurement epoch over nighttime by introducing for each satellite-receiver pair measurement the corresponding vTEC (vertical TEC) extracted from GIM using interpolation operated at both spatial and temporal domains to obtain the vTEC value at the given location and time. The median value is chosen to determine the daily receiver DCBs.
2.2-Refinement of TEC samples
After carrier-phase smoothing of code measurements on the L_1 (1,575.42 MHz) and L_2 (1,227.60 MHz) frequencies, and generating of vTEC measurements, these latter may still be contaminated with persistent errors may be caused by a poor receiver DCBs estimation or the related error to ionospheric thin-layer assumption and Modified single-layer model (MSLM) mapping function especially when there are large geomagnetic gradients. These contaminated measurements can affect drastically the quality of TEC maps. It is clear that low quality inputs leads to maps of limited quality. Thus, to maximize the accuracy in estimating the spatial distribution of TEC, we compare these measurements with their corresponding values from CODE'S (Centre of Orbit Determination in Europe) predicted ionosphere maps using interpolation operated at both spatial and temporal domains to obtain the vTEC value at the given location and time. All bad samples are rejected according to a pre-determined threshold. We note, with much less samples than originally planned it is possible to produce quality maps as it is showed in literature that the meridional correlations of ionospheric day-to-day variability shows a strong correlation with lengths about 7 degrees in mid-latitudes and the zonal correlation lengths are approximately 20 degrees at mid-latitudes. This explains, the no need as many samples to produce quality maps.
TEC mapping is based on an automatic spatial interpolation algorithm using ordinary kriging and a planar de-trending by decomposing data samples into a deterministic trend component and an autocorrelated random component. kriging takes into account the spatial correlation among the sample data using experimental semivariogram. For the computation of the latter an appropriate lag distance is chosen automatically such as that the number of pairs is greater than a fixed threshold. The spatial variation is assumed omnidirectional and the resulted semivariogram is isotropic. The experimental semivariogram needs to be well adjusted by a theoretical model. The developed software package apply an automatic fitting by testing several parametric models namely spherical, exponential, Gaussian, Matern family and M. Stein’s parameterization model. During fitting process, the variogram parameters namely nugget, sill and the effective range are determined using Weighted Least Squares (WLS) estimator. For Matern models different degrees of smoothness of the process varying from 0.1 to 10 are tested to select their optimum value. The algorithm is iterative with stopping criterion based on citeria: 1) the maximum number of iterations has been reached and 2) the distance between the experimental variogram and the model is small enough. The algorithm iterates over the candidate semivariogram models and picks the model that shows the smallest Sum of Squared Error (SSE) between semivariogram model and sample semivariogram. Once the appropriate model is determined, we need to choose whether we use local or global kriging. Local kriging uses a set of IPP (Ionosphere Pierce Points) in the neighborhood of IGP (Ionosphere Grid Points) while global kriging uses all IPP. If the correlation range is less than the linear dimension of the data sample extent local kriging is chosen otherwise global kriging is used. For the local kriging the maximum radius to use to select a set of IPP is set dynamically to the correlation range value. This kriging algorithm can be repeated for any desired temporal interval, and vTEC maps are obtained automatically with high spatial resolution.
The performance of kriging model is assessed by the 10–fold cross-validation where the original sample data is split randomly into 10 equal subsample. Each single subsample is retained as the validation data for testing the model, and the remaining 9 subsamples are used as training data. The cross-validation process is then repeated 10 times and the residual error can be inferred from observed and interpolated values. Accuracy is assessed by calculating the Root Mean Square Error (RMSE) over all sample data.
The data processing software package is implemented partly in Matlab and partly in R. Matlab is used for GPS data preprocessing, receiver DCBs estimation and to extract good measurements of vTEC over IPPs using CODE'S predicted ionosphere maps. These vTEC data were interpolated using R extension packages such as automap and gstat. The automatic mapping system results in two maps: a map of the kriging prediction and a map of the kriging standard error.
The developed algorithm is applied over the west-side of the Mediterranean Sea [20°W-15°E, 25°N-50°N] during January-March 2014 using hourly data from 6 RGPA (Réseau GPS Permanent Algérien) stations and 27 IGS stations to generate hourly-TEC fixed grids of 1° x1° .
3.1-validation of DCBs estimation
To investigate the performance of the proposed method, DCBs of IGS receiver stations estimated during days 1-59 of the year 2014 by the proposed method were compared with their corresponding DCBs released by IGS in their respective final ionospheric products. It has been shown that the daily difference of all stations is comprised between -0.4ns and 0.4ns. The RMS of all differences is lower than 0.32 ns. This confirms the reliability of the daily receiver DCBs used method.
3.2-Performance of data refinement
For each vTEC over IPP the corresponding value from CODE’s predicted ionosphere maps is extracted using interpolation operated at both spatial and temporal domains to obtain the vTEC value at the given location and time and the difference is compared to a threshold. A threshold of ±1TECu (TEC unit) reject generally more than 60% of samples. During the day 7 of the year 2014, the mean derived RMSE from a cross validation of TEC maps genertaed from refined vTEC samples using a threshold of ±1TECu and raw vTEC samples is 0.08±0.05TECu and 0.65±0.34TECu respectively. This comparison shows that data refinement leads for suitable results.
3.3- validation of kriging technique
During one day, the structure of the experimental semivariogram can be very different during the hours of the day and they can change significantly. The spatial correlation of TEC decreases with distance. Fitting of experimental variogram over 59 days, shows that Matern and M. Stein’s parameterization models were selected respectively 45% and 43% of realizations where Exponential and Gaussian models were used respectively 8% and 4% of realizations. It can be observed that Spherical model has never been selected. This was expected as from Matern Family all models can be generated. We notice the nugget effect, varies during the day especially during sunrise and noon.
Under quite conditions, the RMSE obtained from a cross-validation procedure of derived maps were comprised between 0.04 and 0.19TECu and the mean RMSE of derived maps is 0.10 TECU.
The second test uses CODE GIM post-processed product to get an accuracy estimate of the different TEC maps. The performance can be obtained by computing the bias, standard deviation and RMS during days 1-59 of the year 2014. The comparison shows a good agreement with mean differences of 0.81± 0.59 TECU.
The developed software is a Matlab and R package used to estimate receiver DCBs from dual-frequency GPS observations and rapid IGS GIM products. Also it is used to process in near-real time hourly GPS data to produce hourly vTEC grids with spatial resolution of 1° x 1°. To improve quality TEC maps, the measured vTEC over IPPs are compared to their corresponding values from CODE'S predicted ionosphere maps following a predefined threshold to select a quality inputs for interpolation TEC mapping is based on an automatic spatial interpolation algorithm using ordinary kriging and a planar de-trending. During the automatic fitting of the experimental semivariogram several parametric models namely spherical, exponential, Gaussian, Matern family and M. Stein’s parameterization model are tested to select the appropriate model and to estimate its related parameters . A test is applied to a mid-latitude regional GPS Network containing 6 RGPA stations and 27 IGS stations distributed non-uniformly over Western Europe and North Africa to assess performance of the software package and kriging interpolation method.