Using Ordinary Kriging for the Creation of Scintillation Maps

This paper introduces Ordinary Kriging as a method for the creation of maps of scintillation indices S4 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 er‐ ror estimate of the map values, could ultimately lead to their direct implementation in the positioning’s stochastic model.


Introduction
This paper introduces Ordinary Kriging as a method for the creation of maps of scintillation indices S 4 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 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].
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),  Mitigation of Ionospheric Threats to GNSS: an Appraisal of the Scientific and Technological Outputs of the TRANSMIT Project where a positioning computation could be affected by multiple satellite-receiver links under scintillation.

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.

Scintillation indices
According to [8] a simplified model for one particular signal as seen by the receiver can be expressed as (1): 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(t) = 2PδPs(t)cos (ωt + ϕ + δϕ) + n ( t ) (2) The following indices, S 4 (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 4 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].

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]. 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 σ ϕ T 2 (on GPS L1) as (4).
Mitigation of Ionospheric Threats to GNSS: an Appraisal of the Scientific and Technological Outputs of the TRANSMIT Project The Conker model for phase scintillation (on L1) is given as (5), (6) and (7).
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.

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 4 and σ ϕ60 (standard deviation over 60 seconds) maps, which avoid some of these complications.

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

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.

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

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).
The Ordinary Kriging system estimation variance is given by (12).
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.

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 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 gener-ated 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.

Overcoming some of the limitations
The limitations discussed in paragraph 5.2 are: 1. The thin shell model

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

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.

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] Using Ordinary Kriging for the Creation of Scintillation Maps http://dx.doi.org/10.5772/58781 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.
• 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.  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.

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