Optimal Geostatistical Methods for Interpolation of the Ionosphere: A Case Study on the St Patrick’s Day Storm of 2015

Geostatistical Analyst is a set of advanced tools for analysing spatial data and generating surface models using statistical and deterministic methods available in ESRI ArcMap software. It enables interpolation models to be created on the basis of data measured at chosen points. The software also provides tools that enable analyses of the data variability, setting data limits and checking global trends, as well as creating forecast maps, estimating standard error and probability, making various surface visualisations, and analysing spatial autocorrelation and correlation between multiple data sets. The data can be interpolated using deterministic methods providing surface continuity, and also by stochastic techniques like kriging, based on a statistical model considering data autocorrelation and providing expected interpolation errors. These properties of Geostatistical Analyst make it a valuable tool for modelling and analysing the Earth’s ionosphere. Our research aims to test its applicability for studying the ionosphere, and ionospheric disturbances in particular. As raw source data, we use Global Navigation Satellite Systems (GNSS)-derived ionospheric total electron content. This paper compares ionosphere models (maps) developed using various interpolation methods available in Geostatistical Analyst. The comparison is based on several indicators that can provide the statistical characteristics of an interpolation error. In this contribution, we use our own method, the parametric assessment of the quality of estimation (MPQE). Here, we present analyses and a discussion of the modelling results for various states of the ionosphere: On the disturbed day of the St Patrick’s Day geomagnetic storm of 2015, one quiet day before the storm and one day after its occurrence, reflecting the ionosphere recovery phase. Finally, the optimal interpolation method is selected and presented.


Introduction
In recent decades, Global Navigation Satellite Systems (GNSS) have found several applications in a broad range of geosciences. Along with the ongoing development into GNSS technology, the offered accuracy has increased, opening up new research possibilities. GNSS signals are primarily applied to provide the user's position [1,2]. This technology is also increasingly used to monitor deformations in the Earth's crust [3]. However, precise GNSS positioning and GNSS-based geodetic and geodynamic studies require accurate corrections of ionospheric delay [4,5]. In this respect, another field of studies based on GNSS data is atmosphere remote sensing [6,7]. Therefore, satellite data are also often used in modelling the ionosphere. There are a number of various modelling methods that differ in terms of accuracy and reliability [8]. A comprehensive review of the most popular GNSS-derived ionosphere models is provided in [9,10].
A predominant group of GNSS-derived global ionospheric maps (GIMs) are elaborated by using mathematical methods [11], such as spherical harmonic expansion (SHE). This approach is often used for two reasons. Firstly, because of using noisy phase-smoothed pseudorange data. Secondly, it is often used due to the heterogeneous resolution of the observations. This approach reflects long-term effects of spatial and temporal changes of total electron content (TEC) well. However, it is characterized by a loss of information about the local structures in the ionosphere. The obtained ionosphere maps have low spatial and temporal resolution, and only allow for a general overview of the ionosphere response to a geomagnetic storm.
In order to overcome the irregular data problem, several authors have proven that the use of the kriging technique is also possible. This approach takes into account the spatial correlation among the data to be interpolated. Odijk [4] and Orus et al. [12] conducted research on the application of kriging technique to global ionosphere mapping. In Orus et al., ordinary kriging was applied to improve the ionospheric maps from the Technical University of Catalonia (UPC), computed with GPS data. The authors show that new UPC kriging GIMs (UQRG) are characterised by lower Root Mean Square (RMS) error in the observed slant TEC (sTEC) than the original UPC GIMs and IGS GIMs. However, this model remains at a resolution-namely, 2.5 • × 5.0 • × 15 min-that does not provide fine details of the ionosphere.
Kriging is increasingly used in regional ionosphere modelling. Stanislawska et al. [13] applied this technique in regional vertical TEC (vTEC) estimation over the European area. The analyses were carried out for a few quiet and disturbed days in September 1999. Other researchers also successfully used this interpolation technique in regional ionosphere studies. Deviren and Arikan [14] presented the mapping algorithm based on universal kriging with linear trend for midlatitude regions and ordinary kriging for other regions. According to Shukla et al. [15], for the Indian region, it is clearly more suitable to use ordinary kriging in place of planar fit to estimate the ionospheric delay. They also found out that ordinary kriging performed better than the bilinear interpolation technique.
Due to the dynamic nature of changes in the disturbed ionosphere, the accuracy and the resolution of the developed models is still often insufficient to analyze the storm time effects in detail [16] as well as to support precise positioning applications, especially those requiring fast ambiguity resolution [10]. Therefore, the development of high-accuracy models with higher spatial and temporal resolution is still a popular research topic.
In this contribution, we aim at providing examples of evaluation and validation of various TEC interpolation techniques-namely statistical and deterministic methods. Therefore, Geostatistical Analyst was used to provide a set of ionospheric grids that were subsequently evaluated in terms of their accuracy and reliability. Recently, a new approach based on processing precise carrier phase GNSS data was proposed by Krypiak-Gregorczyk et al. [17]. This approach results in precise ionospheric TEC estimations in an irregular cloud of ionosphere piercing points (IPP). This cloud of points has to be provided to users in the form of a regular grid. However, there is a large number of available interpolation methods, and it is difficult to choose the most suitable one that would be appropriate to a range of geographical regions and ionospheric conditions. Therefore, in this initial contribution, we aim to test different geostatistical methods offered by Geostatistical Analyst for modelling the ionosphere over the European region (e.g., ionospheric mid-latitudes). All the geostatistical methods available in the Geostatistic Wizard of ESRI ArcMap10.4, divided into deterministic and stochastic methods, were used in the work. Geostatistics is a group of algorithms, based on generalised least-squares regression, enabling solutions of both deterministic and stochastic models [18]. It has a solid mathematical theoretical foundation developed by Matheron [19,20]. An advantage of the deterministic method is the lack of any requirement as to the normality of the data distribution of the analysed feature, because it assumes a non-random character of the studied phenomenon and describes the model as one fixed function defined in space. In contrast, stochastic methods assume a random nature of the value of an environmental variable, which is a stationary Gaussian stochastic process [21,22]; as a result, krigings allow geostatistical simulations [23].
The applicability of geostatistical methods has been well documented in Earth Sciences [23][24][25]. However, it was only in the 1980s that these methods started to be applied on a more regular basis in disciplines other than geology, where they originate. Geostatistics is a branch of spatial data statistics that not only includes the attributes in the analysis, but also their location in space and/or time [26]. Its specific nature results from the fact that it analyses and models continuous variables (attributes), i.e., variables with a determined (though usually unknown) value at each point of the analysed field, or, in the case of qualitative variables, a determined state [27]. This method is founded on probability theory, specifically the random functions theory [28,29].
Since the early 1990s, geostatistics has mostly been used to describe spatial structure and to estimate the values of attributes in non-measured fields. Geostatistics takes into consideration the spatial distribution of attributes. It is used to analyse how attribute values depend on the distance between the points at which particular values have been measured. That is why it is easier to statistically describe and interpret the spatial distribution of an analysed attribute. We can distinguish geostatistical methods that make it possible to analyse not only quantitative, but also qualitative data. Detailed descriptions of the geostatistical methods applied in the research are available in the literature of the subject [23,24,30]. These methods are based on research conducted in mathematics and life sciences [20,22,31], hence they have solid theoretical bases. Moreover, geostatistical methods can be used to evaluate precisely the uncertainty of the estimates of the attribute value in areas where the attributes were not actually measured. Geostatistical Analyst is a set of advanced tools for analysing spatial data and generating surface models through the use of statistical and deterministic methods. It enables the interpolation of models on the basis of data measured at chosen points. Geostatistical Analyst enables analyses of data variability, setting data limits and checking global trends, as well as creating forecast maps, estimating standard error and probability, making various surface visualisations, including contour lines (isolines), and analysing spatial autocorrelation and correlation between multiple data sets. The data can be interpolated using deterministic methods providing surface continuity, and also by determining an estimation of the most probable value, while the kriging techniques, which are based on a statistical model considering data autocorrelation, evaluate the expected error. These properties of Geostatistical Analyst make it a valuable tool for modelling and analysing the Earth's ionosphere. Hence, our research aims to test its applicability for studying the ionosphere, and ionospheric disturbances in particular. The need to analyse the applied interpolation methods and the effects of selecting the wrong interpolation methods for spatial data visualisations has been pointed out by Ogryzek [32]. The article uses the method of the parametric assessment of the quality of estimation (MPQE) proposed by Ogryzek [33]. Here, the validation is made on the basis of comparing methods by using root mean square prediction error (RMSE) in particular [30,34].

Methodology
For providing ionospheric TEC data for further analysis, dual-frequency multi-GNSS data (GPS + GLONASS) from the Polish active geodetic network (ASG-EUPOS) and European EUREF Permanent Network (EPN) stations were used [35]. The presented regional ionosphere model is computed using exclusively precise, absolute (non-differenced) carrier phase GNSS measurements, several orders of magnitude more precise than pseudorange ones. The geometry-free (LGF) carrier phase linear combination is used to eliminate geometry-related observational errors-receiver and GNSS satellite clock biases, tropospheric delays, etc. Unfortunately, this combination, along with information about ionospheric delays, also includes carrier phase ambiguities and carrier phase hardware delays that have to be estimated.
The University of Warmia and Mazury in Olsztyn (UWM) model used in this study is based on the single layer model (SLM) approach [36] (Figure 1). In the first processing step, carrier phase biases (ambiguities + hardware delays) are estimated.
For each continuous satellite arc [17]. In the second step, the obtained biases are used, together with dual-frequency, multi-GNSS carrier phase observations, to calculate the ionospheric TEC at IPP Single layer model geometry [36]. z-satellite's zenith distance at the receiver's location; z'-satellite's zenith distance at the ionospheric pierce point; R-the mean Earth radius; Hthe height of the single layer.
In the modelling process, the semivariogram of the approximated function has been studied using nugget, spherical, Gaussian, power, exponential and linear base models. We defined not only the range, but also the direction to specify the sector type ( Figure 2). The choice of the appropriate base model was carried out using the MPQE approach proposed by Ogryzek [33], which is based on: ME (mean prediction error), RMSE (root mean square prediction error), ASE (average standard error), MSE (mean standardised prediction error), and RMSSE (root mean square standardised prediction error). In the present study, for the geostatistical methods the abovementioned errors received different weightings: RMSE = 60%, ME = 10%, ASE = 10%, MSE=10%, RMSSE = 10%. Accurately matching the model to the source data may not generate the lowest errors, which is why it is important to validate the results prior to choosing the optimal model. In practice, choosing the best method is problematic, due to a large amount of quality statistics of the estimation and the lack of clear selection criteria [33]. The MPQE method uses an optimisation algorithm based on estimation parameters from the validated cross validation (CV) and subsets validation (SV) model. Single layer model geometry [36]. z-satellite's zenith distance at the receiver's location; z'-satellite's zenith distance at the ionospheric pierce point; R-the mean Earth radius; H-the height of the single layer.
In the modelling process, the semivariogram of the approximated function has been studied using nugget, spherical, Gaussian, power, exponential and linear base models. We defined not only the range, but also the direction to specify the sector type ( Figure 2). The choice of the appropriate base model was carried out using the MPQE approach proposed by Ogryzek [33], which is based on: ME (mean prediction error), RMSE (root mean square prediction error), ASE (average standard error), MSE (mean standardised prediction error), and RMSSE (root mean square standardised prediction error). In the present study, for the geostatistical methods the abovementioned errors received different weightings: RMSE = 60%, ME = 10%, ASE = 10%, MSE=10%, RMSSE = 10%. Single layer model geometry [36]. z-satellite's zenith distance at the receiver's location; z'-satellite's zenith distance at the ionospheric pierce point; R-the mean Earth radius; Hthe height of the single layer.
In the modelling process, the semivariogram of the approximated function has been studied using nugget, spherical, Gaussian, power, exponential and linear base models. We defined not only the range, but also the direction to specify the sector type ( Figure 2). The choice of the appropriate base model was carried out using the MPQE approach proposed by Ogryzek [33], which is based on: ME (mean prediction error), RMSE (root mean square prediction error), ASE (average standard error), MSE (mean standardised prediction error), and RMSSE (root mean square standardised prediction error). In the present study, for the geostatistical methods the abovementioned errors received different weightings: RMSE = 60%, ME = 10%, ASE = 10%, MSE=10%, RMSSE = 10%. Accurately matching the model to the source data may not generate the lowest errors, which is why it is important to validate the results prior to choosing the optimal model. In practice, choosing the best method is problematic, due to a large amount of quality statistics of the estimation and the lack of clear selection criteria [33]. The MPQE method uses an optimisation algorithm based on estimation parameters from the validated cross validation (CV) and subsets validation (SV) model. Accurately matching the model to the source data may not generate the lowest errors, which is why it is important to validate the results prior to choosing the optimal model. In practice, choosing the best method is problematic, due to a large amount of quality statistics of the estimation and the lack of clear selection criteria [33]. The MPQE method uses an optimisation algorithm based on estimation parameters Sensors 2020, 20, 2840 5 of 12 from the validated cross validation (CV) and subsets validation (SV) model. The indicators (parameters) that have been analysed are the characteristics of the statistical error of the interpolation: RMSE and ME for each tested interpolation model. For each sample (analysis day), the MPQE values are based on the data of 45% RMSE CV and 5% ME CV (parameters from the validated cross validation) and 45% of the RMSE SV and 5% ME SV data (parameters from the validated subsets validation (SV) model).
To choose the optimal method, we use a weighting scheme: 45% RMSE (CV) + 5% ME (CV) and 45% for RMSE (SV) 5% ME (SV), i.e., for each model: where IDW is inverse distance weighting, GPI is global polynomial interpolation and OK is ordinary kriging. The method with the lowest weighted prediction error is the optimal method.

Experiment
In order to test the applicability of the Geostatistical Analyst tools to the ionosphere interpolation, a test period of three days, characterised by very different geomagnetic and ionospheric activity, was selected. This period includes three days, from 16  Examples of IPP locations for their measurements collected at 11.10 UT on 17 March 2015 are presented in Figure 3. The indicators (parameters) that have been analysed are the characteristics of the statistical error of the interpolation: RMSE and ME for each tested interpolation model. For each sample (analysis day), the MPQE values are based on the data of 45% RMSE CV and 5% ME CV (parameters from the validated cross validation) and 45% of the RMSE SV and 5% ME SV data (parameters from the validated subsets validation (SV) model).
To choose the optimal method, we use a weighting scheme: 45% RMSE (CV) + 5% ME (CV) and 45% for RMSE (SV) 5% ME (SV), i.e., for each model: where IDW is inverse distance weighting, GPI is global polynomial interpolation and OK is ordinary kriging. The method with the lowest weighted prediction error is the optimal method.

Experiment
In order to test the applicability of the Geostatistical Analyst tools to the ionosphere interpolation, a test period of three days, characterised by very different geomagnetic and ionospheric activity, was selected. This period includes three days, from 16  Examples of IPP locations for their measurements collected at 11.10 UT on 17 March 2015 are presented in Figure 3.

Validation
Here, we aim to evaluate and validate various TEC interpolation techniques-statistical and deterministic methods. The procedure of validation of the map: Stage 1: Preliminary data analysis; Stage 2: Mapping by different interpolation methods; Stage 3: Execution of validation; Stage 4: Comparison of estimation assessment parameters; Stage 5: Selection of the optimal geostatistical method.
The results obtained from the analyses (even comparing raw RMSE and ME from cross-validation) showed that the polynomial interpolation (LPI) method is the best among those tested, as the RMSE and ME values are a single study, although not always. It often happened that the OK method had better results than the LPI method. However, on average, LPI is the method that gives lower error values. Additionally, geostatic methods (krigings) require the implementation of a complex process, the analysis of errors of many models and the choice of the one that is closest to the criterion of estimation quality assessment. The choice of an appropriate model type in kriging methods is not unambiguous; the final decision is made by an analyst with extensive knowledge in this area. For example, for data during a storm (Table 1): For each daily data set, we produced 72 maps at 20 min intervals. The validation was carried out to confirm, in a manner consistent with the assumptions, whether the applied interpolation procedures led to the expected results (Table 2). It is clearly visible from Table 1 that, during the St Patrick's Day storm, the MPQE increased by half. The lowest values were obtained one day after the storm. The worst results were obtained for global polynomial interpolation (GPI) and for simple kriging (SK). A day before the storm, the results are worse than the LPI results by 30-40%, on the day of the storm the difference is 70-80%, and one day after the storm it is more than 100-120%. The lowest MPQE was obtained on each day using the local polynomial interpolation (LPI) method. The MPQE value increased on the day of the storm, and then abruptly fell one day after the storm, during the recovery phase. The ordinary kriging (OK) method lies within +/-0.05% statistical significance of MPQE, but only a day before and during the storm. On the day after the storm, the OK method is worse than the LPI method by 21%. This means that only the local polynomial interpolation method should be used after the storm, but in both calm and stormy periods either of the two methods (LPI or OK) can be used to model the ionosphere. The average results for the analysis performed during the three tested days reveal that the best method was LPI, although each method has a statistically significant difference in results. Generally speaking, the results obtained for RMSE should not be averaged over several days, and modelling should not be made without recognising the RMSE errors. The differences in deviations from errors are at the level from 6% to 70%. The averaged MPQE results are close to errors obtained on a calm day before the storm, which makes it possible to conclude that LPI is the best method on calm days, even though the research undertaken on a calm day fails to show statistically significant difference of MPQE for the OK method. Figure 4 shows examples of the state of the ionosphere over Europe on 16 March 2015 (a), 17 March 2015 (b) and 18 March 2015 (c) at 11.20 UTC derived with the LPI method. A total of 72 maps were analysed for each interpolation method. According to Krivoruchko [18], RMSE is an indicator (parameter) recommended for comparing the results of spatial analysis using various interpolation methods. Its value should be as small as possible. For each sample (day of analysis), the MPQE (RMSE 100% data and RMSE 90% data) was the lowest if the LPI method was applied. In addition, errors generated during storms are the largest, while in calm periods, the errors drop on average by 65%. An important part of the analysis of the results is the analysis of partial results, where at times the OK method generated the smallest errors, or where the differences between errors generated in the OK and LPI methods were below the statistical significance. The similarities and differences in error generating in the ordinary kriging and polynomial interpolation methods, and the possibilities of obtaining similar errors as a result, have been described by Matheron [31]. Zhou et al. [25] suggested ordinary kriging as an alternative to estimation through the polynomial interpolation method.  .20 UTC derived with the LPI method. A total of 72 maps were analysed for each interpolation method. According to Krivoruchko [18], RMSE is an indicator (parameter) recommended for comparing the results of spatial analysis using various interpolation methods. Its value should be as small as possible. For each sample (day of analysis), the MPQE (RMSE 100% data and RMSE 90% data) was the lowest if the LPI method was applied. In addition, errors generated during storms are the largest, while in calm periods, the errors drop on average by 65%. An important part of the analysis of the results is the analysis of partial results, where at times the OK method generated the smallest errors, or where the differences between errors generated in the OK and LPI methods were below the statistical significance. The similarities and differences in error generating in the ordinary kriging and polynomial interpolation methods, and the possibilities of obtaining similar errors as a result, have been described by Matheron [31]. Zhou et al. [25] suggested ordinary kriging as an alternative to estimation through the polynomial interpolation method.  An important part of the analysis of the results is the analysis of partial results, where at times the OK method generated the smallest errors, or where the differences between errors generated in the OK and LPI methods were below the statistical significance. The similarities and differences in error generating in the ordinary kriging and polynomial interpolation methods, and the possibilities of obtaining similar errors as a result, have been described by Matheron [31]. Zhou et al. [25] suggested ordinary kriging as an alternative to estimation through the polynomial interpolation method. Own study). Figure 5 shows a comparison for a stormy day at the same time for three different methods. The differences are clearly visible, which proves the need for model development and validation to select the most suitable one for a particular dataset.

Conclusions
Global maps of ionospheric total electron content are produced by interpolating GNSS-derived TEC measurements. These maps are produced to test data acquisition, monitoring facilities, and mapping techniques. The TEC mapping can provide accurate ionospheric calibrations to navigation systems. These maps are also used to monitor ionospheric weather, and to nowcast ionospheric storms, which often occur in response to activities in solar wind and the Earth's magnetosphere and thermosphere. The objective of the presented analysis was to evaluate the accuracy of the interpolation methods, available in the analysis software ArcGIS 10.2.1, for the ionospheric TEC modelling for the area of Europe. Then, a number of interpolation methods were tested, both deterministic-inverse distance weighting (IDW), global polynomial interpolation (GPI), radial basis function (RBF), local polynomial interpolation (LPI), and geostatistical-ordinary kriging (OK), simple kriging (SK), universal kriging (UK) and disjunctive kriging (DK). The result is 72 maps per day representing the state of the ionosphere for the area of Europe. The interpolation precision indicates the average size of the resulting errors. Based on the case studies conducted, among the deterministic methods, the most accurate method is characterised by LPI, while among the geostatistical methods it is OK. This may suggest that the ionospheric TEC provided in an irregular cloud of IPPs is a process that is determined by the ionospheric physical properties, and to a lesser extent by random factors. However, when analysing quiet days only, the results show that geostatistical models may perform better. It means that one may choose different interpolation

Conclusions
Global maps of ionospheric total electron content are produced by interpolating GNSS-derived TEC measurements. These maps are produced to test data acquisition, monitoring facilities, and mapping techniques. The TEC mapping can provide accurate ionospheric calibrations to navigation systems. These maps are also used to monitor ionospheric weather, and to nowcast ionospheric storms, which often occur in response to activities in solar wind and the Earth's magnetosphere and thermosphere. The objective of the presented analysis was to evaluate the accuracy of the interpolation methods, available in the analysis software ArcGIS 10.2.1, for the ionospheric TEC modelling for the area of Europe. Then, a number of interpolation methods were tested, both deterministic-inverse distance weighting (IDW), global polynomial interpolation (GPI), radial basis function (RBF), local polynomial interpolation (LPI), and geostatistical-ordinary kriging (OK), simple kriging (SK), universal kriging (UK) and disjunctive kriging (DK). The result is 72 maps per day representing the state of the ionosphere for the area of Europe. The interpolation precision indicates the average size of the resulting errors. Based on the case studies conducted, among the deterministic methods, the most accurate method is characterised by LPI, while among the geostatistical methods it is OK. This may suggest that the ionospheric TEC provided in an irregular cloud of IPPs is a process that is determined by the ionospheric physical properties, and to a lesser extent by random factors. However, when analysing quiet days only, the results show that geostatistical models may perform better. It means that one may choose different interpolation methods when modelling the quiet and disturbed ionosphere.
The results obtained from the analyses (even comparing raw RMSE and ME from cross-and subset validation) showed that the LPI method, among the subjects, generates the lowest errors. It often happened that the OK method had better results than the LPI method. However, on average, LPI is the method that gives lower error values. Additionally, geostatic methods (kriging) require the implementation of a complex process, the analysis of errors of many models, and the selection of the one that is closest to the criterion of estimation quality assessment. The choice of an appropriate model type in kriging methods is not unambiguous; the final decision is made by an analyst with extensive knowledge in this area. For example, during the modelling process they have to perform many iterations for one series of data to compare results obtained by different approximation methods and (several) functions. However, in some cases, especially after the end of the storm, lower RMSE and ME parameters could be obtained by cracking methods. Nevertheless, taking into account that geostatistical methods require additional knowledge from the data analyst and should rather be used on smaller plots and require more time and, in addition, that the results do not show significant statistical differences, the LPI can be considered as the optimal method, where selecting only a number of sectors is a substitute for geostatistical modelling. Since this is an initial study, more focus should be made on the validation methods in the future. In the further analysis, interpolation methods should be looked at using a control set (validation layer) due to the large number of observations (every 20 min), or by using the cross-validation method of the leave-one-out type, consisting in the sequential elimination of data from each point and, performing interpolation at this point on the basis of the remaining data.

Conflicts of Interest:
The authors declare no conflict of interest.