## 1. Introduction

The distribution of electrons in the ionosphere is of interest to scientists and also to engineers working on applications such as earth–space communication systems, which must transmit through the ionosphere, and skywave systems, which make use of ionospheric refraction. Electron content is commonly examined using total electron content (TEC) mapping. This mapping ﬁnds use in other applications such as studying the evolution of magnetic storms, which have, in the past, had profound effects on satellite communication systems and on other critical ground-based systems. Information on the electron content of the ionosphere can be collected using the global positioning system (GPS) and by examining the phase and amplitude changes which occur in paths between transmitting satellites and ground-based receivers. These data can then be processed in order to create maps of the ionospheric TEC.

As the number of paths between GPS ground stations and satellites is relatively low, producing TEC maps is an exercise in reconstruction from sparse data. Recent research has mainly focused on methods, such as tomography, that provide time-dependent volumetric reconstructions [1], [2]. However, when the data points are too sparsely distributed, these techniques are underconstrained and do not produce meaningful results. In ionospheric studies, problems relating to sparsity are particularly prevalent in historic data sets. For example, in 1992, there were only 25 receiver sites operated by the International GPS Service (IGS) in the U.S. [3], by 1996, there were over 75, and now, there are over 500. Therefore, while the issues due to undersampling have largely disappeared for TEC imaging systems utilizing modern GPS data, they still remain for older data and regularly arise in other geoscience applications [4], [5]. Consequently, interpolation methods still have an important role to play in ionospheric studies. The most commonly used interpolation technique for TEC-mapping studies are kriging and cubic B-spline. In addition, there are many other interpolation methods for geophysical data that have received little recent attention from the ionospheric imaging community.

## 2. Basic concepts on interpolation techniques

Interpolation methods can be divided into two categories, local and global, depending upon the locality of the points which are used to derive a given output point. Local techniques make use of a definition of locality to compute output values; only data which fall within a given point’s local neighborhood are used to calculate output values. Global techniques use a weighted sum of all data to compute output values, and for large numbers of input points, an approximation is generally used. When a new datum is added to a globally interpolated field, the whole field must be recalculated, whereas for a locally interpolated field, only those positions within the neighborhood of the added datum need to be recalculated. These two points tend to favor the use of local techniques.

### 2.1. Cubic B-spline

This part introduces a procedure for multi-dimensional local ionospheric model. The model consists of a given reference part and an unknown correction part expanded in term of B-spline functions. This approach is used to compute regional models of Vertical Total Electron Content (VTEC) based on the International Reference Ionosphere (IRI 2012) and GPS observations from terrestrial Global Navigation Satellite System (GNSS) reference stations. The approach can be used for local and global modelling of different ionospheric parameters. In this paper, the focus lies on the three-dimensional local modelling of the VTEC depending on horizontal position ϕ,λ and on time t. The unknown VTEC is separated into a reference part VTEC_{ref} taken from IRI 2012 and a correction part ΔVTEC which is modelled by a series expansion in tensor products of three systems of 1-D normalized endpoint-interpolating B-splines with and unknown coefficients d_{k1},_{k2},_{k3} [9]. The basic observation equation reads

with

Each 1-D basis function system consists of *K=*2^{J}+2 single B-spline functions *J=*3. The number *K* of the functions depends on the level *J* of the spline.

### 2.2. Ordinary kriging

The kriging technique is a linear interpolator that belongs to the best linear unbiased estimator family estimators. Thus, the main purpose of the kriging technique is to estimate a certain unknown variable

In order to apply the ordinary kriging technique, it is necessary to assume that the random function

After differencing Eq. (5) with respect to

Therefore, in order to get the weights

### 2.3. Semivariogram and covariance function

Previously to use the kriging technique it is necessary to determine the semivariogram (or alternatively two-times it, variogram). This is a function that describes the spatial correlation among the data used in the interpolation, which knowledge is important since it is used as the main input of the kriging algorithm (7). The semivariogram function is computed by means of doing the squared difference between pairs of observations at a semivariogram function is computed by means of doing the squared difference between pairs of observations at a ﬁxed distance

Once the experimental semivariogram is computed, the next step is to adjust this experimental semivariogram to a theoretical one

For the definition of the *covariance function* the stationarity of the first two moments (mean and covariance) of the random function is essential.

Some properties of the covariance function are:

But unlike the semivariogram it can also take negative values.

which is bounded by

In general, the reverse is not true, because the semivariogram is not necessarily bounded. Thus, the hypothesis of second-order stationarity is less general than the intrinsic hypothesis (for the monovariate case) and unbounded semivariogram models do not have a covariance function counterpart.

### 2.4. Different semivariogram models

The use of a semivariogram in a Kriging procedure requires continuous semivariogram values for every distance *Kriging system* where the values of an experimental semivariogram can lead to negative Kriging variances.

There are several reasons to favor the semivariogram instead of the covariance function. The semivariogram is a more general tool than the covariance. Another reason is more of practical interest: The semivariogram, unlike the covariance function, does not depend on the existence of a mean value. In practice, the mean is not known in most cases and has to be estimated out of the data, which also adds a bias. Therefore, the semivariogram is often preferred to the covariance function.

Using *h* to represent lag distance, *a* to represent (practical) range, and *c* to represent sill, the three most frequently used models (fig. 1) are:

Spherical:

Exponential:

Gaussian:

These three models are shown below:

Most of semivariograms are defined through several parameters:

**Sill:** The semivariance value at which the semivariogram levels off. Also used to refer to the “amplitude” of a certain component of the semivariogram. For the plot above, “sill” could refer to the overall sill (1.0) or to the difference (0.8) between the overall sill and the nugget (0.2). Meaning depends on context.

**Range:** The lag distance at which the semivariogram (or semivariogram component) reaches the sill value. Presumably, autocorrelation is essentially zero beyond the range.

**Nugget:** In theory the semivariogram value at the origin (0 lag) should be zero. If it is significantly different from zero for lags very close to zero, then this semivariogram value is referred to as the nugget. The nugget represents variability at distances smaller than the typical sample spacing, including measurement error.

## 3. Results

### 3.1. VTEC maps

The September 28th Coronal Mass Ejection (CME) impacted Earth’s magnetic field at 22:20 UT, September 30, 2012 sparking strong Geomagnetic storms at high latitudes. The Bz component of the interplanetary magnetic field (IMF) sharply deviated to -35 nT during the impact. Geomagnetic K-index reached Kp=7 levels on October 1st at 03:00 UT (fig.2) and NOAA/SWPC issued G3 (Strong) Geomagnetic Storm Level alert which is slightly higher than initially predicted.

In our research we used the GPS data of 18 EPN stations (tab.1) near to the EGNOS Ranging and Integrity Monitoring Stations (RIMS) network at the mid-latitude. After analyzing the geographic location of IPPs (fig.3) for all the observational epochs, a region located between 30° - 60° latitude and -40° - 45° longitude was selected to produce the local ionosphere maps each 15 minutes on a 2.5°x 2.5° grid using two interpolation methods (Kriging and cubic B-spline).

The TEC values obtained at IPPs were interpolated using kriging and the cubic B-spline model in order to create high-resolution regional maps of the ionosphere. The results, produced using the above methods, are compared and analyzed in the following section. Temporal evaluation of TEC distribution over study area on the first disturbed day (30 September) is presented in Fig. 4 via the series of TEC maps. When producing the maps we used 15 min averages of TEC data. This approach provides detailed analysis of the ionospheric response to the storm. A significant feature in latitudinal variations of the ionosphere was the presence of the trough. In Fig. 4, one can see that the trough first occurred east and after that was heading due west.

The ionosphere was modeled for the period of 18 hours (5:00 to 23:00 UT) during three days. The ionosphere maps obtained using Kriging and cubic B-spline methods present good agreement. The maximum TEC was observed around 12:00 UT during the first disturbed day (30 September). This might be explained by active geomagnetic conditions (see Fig. 2).

### 3.2. Semivariogram modelling for a single location

The main purpose of the GPS is to determine the position and velocity of a fixed or mobile object, placed over or near the earth surface, using the signals of the 32 satellites on earth orbit. GPS is a complex and expensive constellations of 32 satellites distributed in 6 orbital planes, at 20,200 km altitude, with an orbit inclination of 55 degrees and an approximately 12 hour period. The first step in interpolation using kriging is the formation of a semivariogram [6], [7].

As it has been mentioned before, one of the steps to solve the kriging equations is to compute an experimental semivariogram, in order to adjust a theoretical one. Thus, in principle, one semivariogram should be computed for each realization of the kriging equations (7). For any set of data we can compute the semivariance for each pair of points. These values can then be plotted against the lag distance as a scatter diagram, called the “semivariogram cloud” by Chauvet (1982) [10].

Figure 5 contains all of the information on the spatial relations in the data to lag. In principle, we could fit a model to it to represent the regional semivariogram, but in practice it is almost impossible to judge from it if there is any spatial correlation present, what form it might have, and how we could it. A more sensible approach is to average the semivariances for each of a few lags and examine the results [11]. Nevertheless, the semivariogram cloud shows the spread of values at the different lags, and it might enable us to detect outliers or anomalies. The tighter this distribution is, the stronger is the spatial continuity in the data.

In present work the semivariogram calculations were provided for different magnetic/ionospheric conditions separately. Correlation distance determined from the semivariogram for disturbed conditions is about 10 degrees. Survey data in two dimensions are often unevenly distributed. Each pair of observations is separated by a potentially unique lag in both distance and direction. To obtain averages containing directional information we must group the separations by direction as well as by distance. We choose a lag interval, the multiples of which will form a regular progression of nominal lag distances as in the one-dimensional case. We then choose a range in distance, usually equal to the lag interval. We also choose a set of direction and a range in direction. When all comparisons have been made the experimental semivariogram will consist of the set of averages for the nominal lags in both distance and direction. We can extend further by computing the average experimental semivariogram over all directions.

The semivariogram is sensitive to outliers and to extreme values in general. If the extreme is near the margin of the region then it will contribute to fewer comparisons than if it is near the centre. The end point on a regular transect, for example, contributes to the average just once for each lag, whereas points near the middle contribute many times. If data are unevenly scattered then the relative contributions of extreme values are even less predictable. The result is that the experimental semivariogram is not inflated equally over its range, and this can add to its erratic appearance.

## 4. Conclusions

This work demonstrates the concept and practical examples of mapping of regional ionosphere, based on GPS observations from permanent stations near to the EGNOS Ranging and Integrity Monitoring Stations (RIMS) network. Interpolation/prediction techniques, such as kriging (KR) and the cubic B-spline, which are suitable for handling multi-scale phenomena and unevenly distributed data, were used to create total electron content (TEC) maps. Their computational efficiency (especially the B-spline) and the ability to handle undersampled data (especially kriging) are particularly attractive. The data sets have been collect into strong geomagnetic storm at September 2012. TEC maps have a spatial resolution of 2.5° and 2.5° in latitude and longitude, respectively, and a 15-minutes temporal resolution. The time series of the TEC maps can be used to derive average monthly maps describing major ionospheric trends as a function of time, season, and spatial location.