Using Ordinary Kriging for the Creation of Scintillation Maps

1. Introduction

This paper introduces Ordinary Kriging as a method for the creation of maps of scintillation indices S₄ and σ_ϕ which could be used to provide the user information on the accuracy of specific receiver-satellite links. Sufficiently accurate scintillation maps, also providing an error estimate of the map values, could ultimately lead to their direct implementation in the positioning’s stochastic model.

Small scale (in the order of a few hundred meters, time varying irregularities in electron density in the ionosphere can cause scattering and diffraction of the radio signals, passing through these irregularities [10].

Figure 1.

Scintillation Principle [2]

Figure-1 provides a graph illustrating scintillation by looking at the incident wave plane (plane of equal phase) which is arriving at a layer with rapidly changing irregularities; the emerging wave front is no longer a plane due to scattering (signal hitting the structures and going off into multiple directions) and diffraction (part of a wave hitting a structure and going forward as a circular wave-front causing an interference pattern). The result for the observer is a signal which is rapidly changing in phase (phase scintillation) and intensity (amplitude scintillation) [2].

Figure 2.

Scintillation Occurrence [3]

The morphology and climatology of scintillations have been studied extensively over the last half century. A brief summary of this is given in figure-2 [3], which shows a geographical overview of scintillation occurrences during periods of maximum and minimum solar activity.

From figure-2, we can see that the main areas of scintillation occurrence are in the equatorial zone (after sunset) and at high latitudes near the auroral oval (at night-time). During solar maximum conditions the scintillation effects are much more pronounced than during solar minimum conditions [3].

The effects of scintillations on the GPS receiver tracking ability, range from a complete loss of lock on the signal (satellite no longer available for positioning purposes), to degradation in the accuracy of the measured satellite – receiver range, which if not de-weighted in the stochastic positioning model, could degrade the accuracy of the computed position [1].

The creation of scintillation maps [11] using Ordinary Kriging is thus particularly valuable for areas with the highest frequency of scintillation occurrence (equatorial and high-latitude), where a positioning computation could be affected by multiple satellite-receiver links under scintillation.

2. Metrics of scintillation

The magnitude of scintillation occurrences (phase and amplitude) can be quantified through scintillation indices. Furthermore some parameters describing the effect of scintillation on the performance of receiver tracking loops, e.g. the Phase Locked Loop (PLL) can be obtained from the receivers’ high rate phase and intensity data, if available to the user. These indices and receiver tracking performance parameters are the possible (combined) inputs for the scintillation maps and are normally only provided by scintillation monitoring receivers, which are specifically designed to maintain lock on the satellite even under high scintillation levels.

The following paragraphs provide an overview of the indices and receiver tracking parameters. The last paragraph concludes with a discussion on the preliminary choice for the map content.

2.1. Scintillation indices

According to [8] a simplified model for one particular signal as seen by the receiver can be expressed as (1):

$$r\left(t\right)=\sqrt{2P}s\left(t\right)\cos \left(\omega t+\mathrm{ }\varphi \right)+n\mathrm{ }\left(t\right)$$E1

In which P is the received signal power (watts), ω is the carrier frequency (rad/s), s (t) is the normalised transmitted signal and n (t) is noise. Scintillation causes a perturbation on both the amplitude (δP) and phase (δϕ) and thus the signal under scintillation can be modelled as (2).

$$r\left(t\right)=\sqrt{2P\delta P}s\left(t\right)\cos \left(\omega t+\mathrm{ }\varphi +\delta \varphi \right)+n\mathrm{ }\left(t\right)$$E2

The following indices, S₄ (amplitude scintillation), σ_ϕ (phase scintillation) and SI (signal intensity), are frequently used to describe the severity of scintillation occurrence. A short description [8] is given here:

1. Signal Intensity (SI): Intensity of the received signal discarding noise. SI can be calculated from the difference between Narrow Band Power (NBP) and Wide Band Power (WBP) of the receiver measured over the same interval (within one bit of the navigation message-20 milliseconds). For the practical calculation of SI, see [12].

2. S₄ index: describing the severity of the amplitude fades caused by scintillation. The S4 index is defined as the standard deviation of the detrended power fluctuation δP, which has a Nakagami-m probability density function. For the practical calculation of S4, Total S4, see [12].

3. σ_ϕ index: describing the phase variations caused by scintillation. The σ_ϕ index is defined as the standard deviation of the de-trended phase. The de-trending is performed to remove slow variations due to e.g. satellite movement. The abbreviation σ_ϕ60 means that the standard deviation is taken over a sixty second period. For the practical calculation of σ_ϕ, see [12].

2.2. PLL and DLL tracking performance parameters

The influence of scintillation (amplitude and phase) can be observed in the receiver through various parameters. Amplitude scintillation has an effect on the Carrier-to-Noise density ratio (C/No) [13], which when calculated over a short timespan decreases with the amplitude fades observed in amplitude scintillation, see figure-3 [9]. Phase scintillation has an influence on the spectral parameters of the phase PSD (spectral slope) and T (spectral strength). Both are calculated from the de-trended phase spectrum [12].

Figure 3.

Amplitude Scintillation Effect on C/No [9]

The actual influence of amplitude and phase scintillation on the receiver tracking performance can be expressed by the tracking jitter variance (or its standard deviation) for the PLL and DLL. The tracking jitter variance may be directly estimated at receiver level or derived (in post processing) from the high rate data when combining this high rate data with the scintillation sensitive Conker model [4, 5]. As an example (3) shows the standard thermal noise formula for the DLL jitter variance without scintillation, whereas under amplitude scintillation the Conker model formulates the DLL jitter variance ${\mathit{\sigma }}_{{\mathit{\varphi }}_{\mathit{T}}}^2$ (on GPS L1) as (4).

$${\sigma }_{{\varphi }_T}^2=\frac{B_{n}d\left(1+\frac{1}{\eta {\left(c/n_0\right)}_{L1-C/A}}\right)}{{2\left(c/n_0\right)}_{L1-C/A}}$$E3

$${\sigma }_{{\varphi }_T}^2=\frac{B_{n}d\left(1+\frac{1}{\eta {\left(c/n_0\right)}_{L1-C/A}\left(1-{2S}_4^2\right)}\right)}{{2\left(c/n_0\right)}_{L1-C/A}\left(1-S_4^2\right)}$$E4

The Conker model for phase scintillation (on L1) is given as (5), (6) and (7).

$${\sigma }_{{\varphi }_{\epsilon }}^2={{\sigma }_{{\varphi }_S}^2+\sigma }_{{\varphi }_T}^2+{\sigma }_{{\varphi }_{,\mathrm{ }osc}}^2$$E5

$${\sigma }_{{\varphi }_T}^2=\frac{B_n\left(1+\frac{1}{2\eta {\left(c/n_0\right)}_{L1-C/A}\left(1-{2S}_4^2\right)}\right)}{{\left(c/n_0\right)}_{L1-C/A}\left(1-S_4^2\right)}$$E6

$${\mathrm{\sigma }}_{{\mathrm{\varphi }}_{\mathrm{S}}}^2=\frac{\mathrm{\pi }\mathrm{T}}{\mathrm{k}{\mathrm{f}}_{\mathrm{n}}^{\mathrm{p}-1}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\left[2\mathrm{k}+1-\mathrm{p}\right]\mathrm{\pi }}{2\mathrm{k}}}$$E7

Since scintillation monitor receivers (e.g. from Septentrio and NovAtel) are outputting their high rate data and/or scintillation indices files containing the spectral parameters (Septentrio), this data combined with the Conker model can be used to calculate the tracking jitter variance for the different GNSS signals (e.g. L1, L2, etc.). As can be seen in (4-7) the Conker model is dependent on the design of the receiver tracking implementation.

2.3. The preliminary choice for the map inputs

The preferred ‘type’ for the interpolated maps would be DLL and PLL tracking jitter maps, generated for all GNSS signal types. These could potentially (if the maps are accurate enough) be directly used as an input for the stochastic model of a receiver within the area of the map. As presented in the previous paragraph, the tracking jitter (variance) is a function of some obtainable receiver parameters (within the network). Unfortunately some of these parameters (e.g. p, T and C\No) are heavily receiver dependent. It will also require the map generation process to include a mapping function to verticalize the slant jitter values and vice versa. To demonstrate the ability of our Ordinary Kriging concept this paper presents S₄ and σ_ϕ60 (standard deviation over 60 seconds) maps, which avoid some of these complications.

3. Short introduction to Ordinary Kriging

The following paragraphs provide a short introduction to Ordinary Kriging which is used for the map interpolation [14].

3.1. Ordinary Kriging (OK)

In Ordinary Kriging a point is estimated as a weighted average from the points around it. It is an unbiased estimator because a constraint is set to the weights to sum up to one.

The Ordinary Kriging algorithm is given in (8). The weights (which sum up to one) are calculated from the underlying geospatial structure which can be found by means of an empirical variogram, described in the next section.

$$Z_{OK}^{\mathrm{*}}\mathrm{ }\left({\mathrm{x}}_0\right)=\sum _{\alpha =1}^n{w}_{\alpha }Z\mathrm{ }\left(x_{\alpha }\right)$$E8

3.2. The Variogram and covariance function

The geospatial relationship (dissimilarity between points over distance / direction) within the area of interest can be described with an empirical Variogram. The inter-distances between the sample-points are calculated, see figure-4. The values making up the variogram are then calculated, through ordering the samples with inter-distances falling within a distance bin (h1, h2, etc.); see figure-5, calculating the value as an average of the square root of the difference in the samples γ (h) (9).

$$\gamma \left(h\right)=\frac{1}{2n}\sum _{\alpha =1}^n{\left(Z\left(X_{\alpha }+h\right)-Z\left(X_{\alpha }\right)\right)}^2$$E9

Figure 4.

Sample points on the ionospheric shell and their inter-distances

The empirical variogram is subsequently replaced with a matching theoretical variogram (e.g. an exponential, Gaussian or spherical variogram model) to filter outliers and to make it practical to calculate the variogram values for different distances, necessary to be used in the Ordinary Kriging system.

Figure 5.

Empirical Variogram Calculation

The variogram can be described with the following characteristics; the nugget, the range and a possible sill, see figure-6. When the variogram has a sill (bounded variogram), it can be replaced with a covariance function through the relationship defined by (10).

Figure 6.

Variogram Characteristics

$$\gamma \left(h\right)=b-C\mathrm{ }\left(h\right)$$E10

When we also order the samples making up the empirical variogram by direction we are able to find a possible directional dependency (checking for anisotropy).

3.3. The Ordinary Kriging system

To calculate the weights in (8), an Ordinary Kriging system can be setup, which minimizes the estimation variance and constrains the weights to sum up to one (11).

$$\left[\begin{array}{cccc}\gamma \left(X_1-X_1\right)& \cdots & \gamma \left(X_1-X_n\right)& 1\\ ⋮& \ddots & ⋮& ⋮\\ \gamma \left(X_n-X_1\right)& \cdots & \gamma \left(X_n-X_n\right)& 1\\ 1& \cdots & 1& 0\end{array}\right]\left[\begin{array}{c}{w}_1^{OK}\\ ⋮\\ w_n^{OK}\\ {\mu }_{OK}\end{array}\right]\mathrm{ }=\left[\begin{array}{c}\gamma \left(X_1-X_0\right)\\ ⋮\\ \gamma \left(X_n-X_0\right)\\ 1\end{array}\right]$$E11

The Ordinary Kriging system estimation variance is given by (12).

$${\sigma }_{OK}^2={\mu }_{OK}-\gamma \left(X_0-X_0\right)+\sum _{\alpha =1}^n{w}_{\alpha }^{OK}\gamma \left(X_{\alpha }-X_0\right)$$E12

Intrinsic Multivariate Correlation is the direct correlation between two or multiple variables, independent from their spatial correlation function. Additional variables having an intrinsic correlation structure with the variable of interest can be directly used as a secondary input, for the creation of the map. Additional variables which have a spatial correlation with the primary variable can also be incorporated through co-Kriging.

4. Ordinary Kriging applied to scintillation maps

4.1. The basic concept

The scintillation data is obtained from a network of stations. All sampled scintillation indices, for all stations and all satellites, with an elevation above the elevation mask are projected on the Scintillation Thin Shell (a model where all scintillations take place at the same height above the Earth’s surface, with the height often taken as 350km). A comprehensive dataset (one year’s worth of network data including all levels of ionospheric disturbances) is used to retrieve the underlying geospatial structure through the calculation of variogram (s). The overall variogram is calculated by combining the values calculated at each epoch or by combining an epoch with another epoch with a fixed time offset, to only take one variable (distance) into account for its calculation. Ordinary Kriging is applied to calculate values for all grid points in the map. Only sample points within the ‘range’ of the Variogram (neighbourhood) are taken into account to calculate the interpolation weights and perform the interpolation. Maps are created epoch by epoch for 24 hours during a day of interest, and can be shown as a film (evolving map during the day).

4.2. The limitations

The limitations of the basic concept described above are threefold; the ionospheric single shell assumption, the limited number of ionospheric data samples in an epoch, and a spatial correlation which is not constant within the area of interest (non-stationarity within the map).

The scintillation thin shell assumption (all scintillations take place at the same height above the Earth’s surface) projects the scintillation occurrences at their respective IPPs (Ionospheric Pierce Point; the projection onto the thin shell model). The latitudinal and longitudinal position of these IPPs depends on the azimuth and elevation of the satellites, the used ellipsoid and the height of the thin shell model. However in reality these scintillation events are taking place mostly throughout the F-Region of the ionosphere (250 – 400km above the Earth’s surface). Depending on the elevation (worse with lower elevation) the projection on the thin shell model will cause an error in the position of the scintillation occurrence. For the estimation of the geospatial correlation (creation of the variograms) this will not lead to significant errors as long as the variograms are based on a large enough dataset and the scintillation thin shell model uses a representative height for the scintillation occurrences (mean). However for the generated maps even a representative shell height will still leave a remaining inaccuracy on the position of the scintillation occurrences.

The use of one epoch of data samples leaves us a very limited sample-set to work with. If we assume 10 stations with on average 8 usable satellites (satellites above the elevation mask) we still only have 80 samples over a significantly large map area. Ordinary Kriging used as an interpolation method will statistically provide the best interpolation, but this does not mean that the interpolated values are correct. The more densely sampled the map is, the more accurate it can present the actual situation. Furthermore scintillations do not present a smooth phenomenon; they arise from small scale structures (a few hundred meters) and Ordinary Kriging interpolation therefore offers a more smoothed climatological (probabilistic) model of what is actually happening.

Our scintillation maps are generated every epoch, and the local time of the sampled data changes during the day, i.e. the Earth is turning underneath a sun fixed system. Since scintillation occurrences depend on local time and latitude, the spatial correlation will change during the day and will also vary for different latitude zones. Since scintillation activity is also dependent on solar activity, the geospatial correlation will also be different for quiet and active days.

4.3. Overcoming some of the limitations

The limitations discussed in paragraph 5.2 are:

1. The thin shell model

2. Limited amount of sample points and smoothing

3. The non-stationarity within the map

The thin shell model has to be applied to allow to make a two dimensional map of the scintillation occurrences. Furthermore from our measurements we are also not able to determine the height where the scintillations occur. However if we manage to choose a representative height (average height of scintillation occurrences) we can minimize the effects on the geo-statistical interpolation model and we can also minimize the average error introduced on the position of the scintillation occurrences.

The limited amount of sample points from just one epoch can be improved in two ways: the inclusion of previous epochs and through inclusion of secondary data with intrinsic correlation with the primary data. The inclusion of previous epochs can be achieved by projecting these epochs on where they would be after the Earth has rotated over the time difference with the current epoch. The de-correlation over time has then to be taken into account as well and is done by introducing a measurement error which increases with time. This method also allows a form of prediction.

The use of secondary parameters could be achieved from the inclusion of other non-scintillation monitoring stations, e.g. Rate of Change of TEC; which are correlated with scintillations.

The non-stationarity within the area of the map can be overcome with a method introduced in [6]. The map has to be split in different local time and latitudinal regions. For each of these regions a variogram is established for different disturbance levels. For each grid-point of the map a distance weighted ‘localised variogram’ is established from the surrounding regional variograms, see figure-7 and 8. The Kriging for the grid point proceeds as normal, including sample points which are within the range of this ‘local variogram’.

Figure 7.

Division into zones

Figure 8.

Distance weighted combined variogram

4.4. Validation of the results

Validation of the results is performed by leaving one network station (located approximately in the middle of the map) out of the map computation (sample points not used for the Ordinary Kriging interpolation). The measured scintillation data for this station is compared with the interpolated values from the map at the location of the stations IPPs.

For each scintillation index, the map will provide two values: the interpolated value and the error estimate for this interpolated value. The Kriging interpolation method is considered to be satisfactory when the following is observed:

• Overall the map value and the actual value are in agreement with each other.

• When the difference between the map value and the actual value is small, then the interpolation error estimate is small as well.

• When the difference between the map value and the actual value gets larger, then the interpolation error estimate gets larger as well.

5. Results and conclusions

5.1. The dataset

The variograms to establish the spatial and temporal correlation are generated from the complete 2012 archive of the CHAIN network [7], see figure-9.

Figure 9.

The CHAIN Network [7]

On average this network consists of 8 to 10 active stations (depending on maintenance etc.). The dataset for the creation of the S4 and σ_ϕ60 maps is a small subset (one day, 9th of March 2012) of the 2012 dataset. The calculation is performed only with these stations (sparse network) and thus no secondary data has been used in the Kriging interpolation. An elevation mask of 20 degrees has been used for the generation of the maps.

5.2. Example Results

Figure-10 till 13 show some examples of the obtained geospatial correlation structure. The overall structure used consists of:

• Local hours: 0-6, 3-9, 6-12 etc.

• Latitude: 45-55, 55-65, etc.

• Kp: 0-2, 3-5, 6+

The figures clearly demonstrate the variability of the correlation structure in time (during the day), latitude and disturbance level.

Figure-14 and 15 show an example σ_ϕ60 map (one epoch) of the obtained geospatial correlation structure and its corresponding interpolation error map (Ordinary Kriging error). The maps will normally be calculated for each epoch during a complete day and can then be presented as a video file, showing the evolution of patches with higher scintillation. When they are used as an input for the stochastic model of a receiver’s position computation they can be stored as a grid file similar to an IONEX (Ionosphere map exchange format) file.

Figure-16 and 17 show a comparison between the obtained map values and the actual values from station TALC which is left out from the map calculations. The map values show compliance with the trend in the actual σ_ϕ60 values and often follow the shape of the peaks, but do not manage to follow the extreme values (smoothing). In areas which are lacking enough (representative) observations they do not manage to follow the shape at all.

5.3. Conclusions and future work

The Ordinary Kriging algorithm seems a promising method to produce scintillation maps. When comparing the map with the directly obtained values it manages to follow the general trend and often follows the same pattern of rise and decline in the observed scintillation. As the Kriging interpolation algorithm has a smoothing effect, it does not manage to follow the extreme values. However one could take this smoothing effect into account when using the maps.

When the map values do not follow the observed pattern, it is due to a lack of (representative) data in the area of interest. The method to improve on this (including secondary data) has been discussed and will be part of the future work on this topic.

Figure 10.

Local Time: 0-6, Lat.: 45-55, Kp: 0-2, Nugget 0.1e-4, Range 725km, Sill 0.9e-4, Power 2.0

Figure 11.

Local Time: 3-9, Lat.: 45-55, Kp: 0-2, Nugget 0.1e-4, Range 725km, Sill 0.9e-4, Power 2.0

Figure 12.

Local Time: 0-6, Lat.: 45-55, Kp: 3-5, Nugget 0.2e-4, Range 600km, Sill 1.6e-4, Power 2.0

Figure 13.

Local Time: 0-6, Lat.: 45-55, Kp: 3-5, Nugget 0.15e-4, Range 600km, Sill 1.2e-4, Power 2.0

Figure 14.

Local Time: 0-6, Lat.: 55-65, Kp: 0-2, Nugget 0.1e-4, Range 400km, Sill 0.65e-4, Power 1.0

Figure 15.

Local Time: 0-6, Lat.: 55-65, Kp: 3-5, Nugget 0.2e-4, Range 400km, Sill 1.2e-4, Power 1.0

Figure 16.

Phi60 map 9th of March 2012

Figure 17.

Phi60 Kriging Error map 9th of March 2012

Figure 18.

PRN 13 Kriging – Observed at station TALC

Figure 19.

PRN 17 Kriging – Observed at station TALC

Acknowledgments

The research described in this paper focuses on a part of the research undertaken for the EU Marie Curie Actions project called “TRANSMIT” (Training Research and Applications Network to Support the Mitigation of Ionospheric Threats). The topic of the research is the development of tools to mitigate ionospheric scintillation effects in both real time and in post mission GNSS positioning computations. We would like to thank the European Commission for their support and making projects like “TRANSMIT” possible.

We would also like to thank the University of New Brunswick allowing us to use a comprehensive scintillation dataset from their CHAIN network (Canadian High Artic Ionospheric Network) [7].