InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Earth and Planetary Sciences » Oceanography and Atmospheric Sciences » "Mitigation of Ionospheric Threats to GNSS: an Appraisal of the Scientific and Technological Outputs of the TRANSMIT Project", book edited by Riccardo Notarpietro, Fabio Dovis, Giorgiana De Franceschi and Marcio Aquino, ISBN 978-953-51-1642-4, Published: July 17, 2014 under CC BY 3.0 license. © The Author(s).

Chapter 11

Regional Ionosphere Mapping with Kriging and B-spline Methods

By Oksana Grynyshyna-Poliuga, Iwona Stanislawska and Anna Swiatek
DOI: 10.5772/58776

Article top


Spherical, Exponential and Gaussian semivariograms models
Figure 1. Spherical, Exponential and Gaussian semivariograms models
The values of Kp index during the experiment
Figure 2. The values of Kp index during the experiment
Example location of IPPs and their TEC values at 12:00-13:00 UT, September 30, 2012, when seven satellites were simultaneously observed.
Figure 3. Example location of IPPs and their TEC values at 12:00-13:00 UT, September 30, 2012, when seven satellites were simultaneously observed.
Ionospheric TEC maps created with Kriging (a) and cubic B-spline (b) interpolations methods
Figure 4. Ionospheric TEC maps created with Kriging (a) and cubic B-spline (b) interpolations methods
An example of semivariograms for disturbed conditions with logarithmic presentation. Crosses are averaged data, lines – fitted data, red − full data set, and blue crosses – data created artificially
Figure 5. An example of semivariograms for disturbed conditions with logarithmic presentation. Crosses are averaged data, lines – fitted data, red − full data set, and blue crosses – data created artificially

Regional Ionosphere Mapping with Kriging and B-spline Methods

Oksana Grynyshyna-Poliuga1, Iwona Stanislawska1 and Anna Swiatek1

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 finds 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 VTECref 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 dk1,k2,k3 [9]. The basic observation equation reads




Each 1-D basis function system consists of K=2J+2 single B-spline functions ΦkJ(x) equally distributed on the unit interval for 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 (Z*) as a linear combination of the known values (Zi):

ωi being the weights computed by the kriging equations (7), that are applied to each value Z(xi)=Zi.

In order to apply the ordinary kriging technique, it is necessary to assume that the random function Zi belongs to the stationary random functions family, which means that the mean values and the standard deviation of Zi have to be independent of the location. Moreover, the unbiased condition over the weights (iωi=1) is imposed. Then, the variance is minimized with the help of the Lagrange multipliers in order to impose the unbiased condition (8) for details:


E[(Z*Z)2] being the Z* variance that can be expressed as a function of the semivariogram:


After differencing Eq. (5) with respect to λ and ωi, and equating to 0, the ordinary kriging equations are obtained in compact form:


Therefore, in order to get the weights (ωi), the following equation, that is expressed in matrix notation, has to be solved:



Ω being the vector that contains the weights ωi and the Lagrange multiplier λ. Γ is the matrix that contains the semivariogram estimations for the known values and locations, and Γ0 is the vector that contains the semivariogram estimations for the unknown values but with known locations.

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 fixed distance dl±Δdl/2, as it is show in (1), see [8] for details:


γ being the experimental semivariogram, m(dl) the number of pairs of observations at a distance dl, Zi and Zj are the observation values that correspond to points at xi and xj at a distance |xixj|=dl.

Once the experimental semivariogram is computed, the next step is to adjust this experimental semivariogram to a theoretical one γ(dl), which must verify several mathematical conditions [8] in order to be applied in the kriging equations (7). It has to be noted that such theoretical semivariograms are classified in well-known semivariogram families, which try to take into account the most number of semivariogram types for families examples.

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:

  • The covariance function is bounded and its absolute value does not exceed the variance

  • Similar to the semivariogram, it is an even function

But unlike the semivariogram it can also take negative values.

  • The covariance function divided by the variance is called the correlation function


which is bounded by

  • Furthermore, the semivariogram function can be deduced from a covariance function 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.

  • A covariance is a positive and negative definite functions.

2.4. Different semivariogram models

The use of a semivariogram in a Kriging procedure requires continuous semivariogram values for every distance |h|. Of course, this cannot be provided by the experimental semivariogram since only discrete measurements can be realized in practice. Fitting the experimental semivariogram by an appropriate semivariogram function helps to overcome this problem. Using a theoretical semivariogram also guarantees that the semivariance of any linear combination of sample values is positive. This is important for setting up a 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: g(h)={c(1.5(ha)0.5(ha)3)c

Exponential: g(h)=c(1exp(3ha))

Gaussian: g(h)=c(1exp(3h2a2))

These three models are shown below:


Figure 1.

Spherical, Exponential and Gaussian semivariograms models

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.


Figure 2.

The values of Kp index during the experiment

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).


Figure 3.

Example location of IPPs and their TEC values at 12:00-13:00 UT, September 30, 2012, when seven satellites were simultaneously observed.

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).

Station Latitude (°N) Longitude (°E)
Maartsbo, Sweden60.617.3
Kirkkonummi, Finland60.224.4
Suldrup, Denmark56.89.7
Daresbury, UK53.3-2.6
Hailsham, UK50.90.3
Potsdam, Germany52.413.1
CBKA, Poland52.221.1
Saint-Mande, France48.82.4
Zimmerwald, Switzer.46.97.5
Vigo, Spain42.0-8.8
Toulose, France43.61.5
Santa Cruz, Portugal39.4-31.1
Cascais, Portugal38.7-9.4
Malaga, Spain36.7-4.4
Palma de Mallorca, Spain39.52.6
Moto, Italy36.915.0
Athens, Greece38.023.9
Ankara, Turkey39.932.8

Table 1.

The coordinates of the selected stations


Figure 4.

Ionospheric TEC maps created with Kriging (a) and cubic B-spline (b) interpolations methods

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].


Figure 5.

An example of semivariograms for disturbed conditions with logarithmic presentation. Crosses are averaged data, lines – fitted data, red − full data set, and blue crosses – data created artificially

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.


This research work is undertaken in the scope of the TRANSMIT ITN (, funded by the Research Executive Agency within the 7th Framework Program of the European Commission, People Program, Initial Training Network, Marie Curie Actions – GA No. 264476.


1 - C. N. Mitchell and P. S. J. Spencer, “A three-dimensional time-dependent algorithm for ionospheric imaging using GPS,” Ann. Geophys., vol. 46, no. 4, pp. 687–696, 2003.
2 - J. M. Pallares, G. Ruffini, and L. Ruffini, “Ionospheric tomography using GNSS reflections,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 2, pp. 321–326, Feb. 2005.
3 - E. Brockmann and W. Gurtner, “Combination of GPS solutions for densification of the European network: Concepts and results derived from 5 European associated analysis centers of the IGS,” in Proc. EUREF Workshop, Ankara, Turkey, May 1996.
4 - M. Liao, T. Wang, L. Lu, W. Zhouzhou, and D. Li, “Reconstruction of DEMS from ERS-1/2 TANDEM data in mountainous area facilitated by SRTM data,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 7, pp. 2325–2335, Jul. 2007.
5 - M. Gianinetto and P. Villa, “Rapid response flood assessment using minimum noise fraction and composed spline interpolation,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 10, pp. 3204–3211, Oct. 2007.
6 - N. A. C. Cressie, Statistics for Spatial Data. Hoboken, NJ: Wiley, 1991.
7 - H. Omre, “The variogram and its estimation,” in Geostatistics for Natural Resources Characterization. Amsterdam, The Netherlands: Reidel, 1984, pt. 1, pp. 107–125.
8 - Cressie, N.A.C., 1993. Statistics for Spatial Data, revised ed. Wiley, New York, EUA.
9 - Schmidt M, Dettmering D, Mößmer M, Wang Y, Zhang J (2011) Comparison of spherical harmonic and B-spline model for VTEC. Radio Sci 46:RSOD11. doi:10.1029/2010RS004609
10 - P. Chauvet, (1982). The Variogram Cloud. In Proceedings of the 17th APCOM Symposium. Colorado School of Mines, Golden, CO, pages 757-764.
11 - R. Webster, M.A. Oliver: Geostatistics for Environmental Scientists. Wiley, Chichester, 2001. Hardbound, 271 pp., ISBN 0471965537.