The Reliability Evaluation of GNSS Observations in the Presence of Ionospheric Perturbations

1. Introduction

Global Navigation Satellite Systems (GNSS) observations are intrinsically uncertain and inaccurate. In fact, the influence of some phenomenon on the accuracy of the GNSS observations can be relatively easily reduced or removed. However, some other random and deterministic phenomena occurring in the GNSS signal propagation path, like ionospheric perturbations, are very difficult (or perhaps impossible) to predict, detect and model.

The ionospheric perturbations are described as fast and random variability of plasma density in the ionosphere. All of those irregularities can produce diffraction and refractions effects causing signal fading. Such power drops can affect the operations of GNSS receivers denying the signal acquisition or worsening the signal tracking, leading in some cases to a loss of lock. Therefore, they are especially harmful to real-time kinematic applications, with an autonomously working single-receiver. In such circumstances one of the most important requirements standing behind the autonomously working GNSS receiver is to ensure a high level of trust of correctness of the observations used by the positioning algorithm. In such applications verification and confirmation of a high-level reliability for the GNSS observations is a very critical issue. Due to that reason the quality monitoring of the observations affected by the ionospheric perturbations can play a crucial role to enhance reliability of positioning solution.

In presence of GNSS signal perturbations one of the most important requirements is effective integrity monitoring of GNSS observations. The user-level integrity monitoring scheme, the so-called Receiver Autonomous Integrity Monitoring (RAIM) should work independently of any external tools. The RAIM is a powerful technique to check consistency of positioning solution, can play strategic role in reliable positioning in the presence of any irregular perturbation of the GNSS observations. The main task of RAIM is to provide to the user up-to-date and valid warnings information when the system’s performance exceeds a user specified tolerance.

In the paper the weighted least-squares-residuals (WLSR) method for reliability control of the GNSS observations is applied. In this method, the process of the Fault Detection and Exclusion (FDE) is performed by statistical tests. The algorithm assumes that the GNSS user estimate a single epoch (instantaneous) navigation solution by performing the Weighted Least-Squares (WLS) estimation. It has to be emphasized that the approach requires at least one redundant measurement. Selected GNSS data-sets from one continuously operated GNSS station located at high latitude, where ionospheric disturbances occur more frequently, have been used for the analysis and for evaluation of applicability of the proposed algorithm.

2. Receiver Autonomous Integrity Monitoring RAIM

In order to verify the integrity of the positioning solution, in the second half of 1980’s, a concept of RAIM was formalized. Since that time a number of definitions of the receiver autonomous integrity were proposed. To start a discussion about RAIM one of the essential definitions should be quoted

Concept Paper 1 (WP/43), AWOP (All Weather Operations Panel) Working Group Meeting, Kobe, Japan, February/March 1994, "Required navigation performance (RNP) - Considerations for the Approach, Landing and Departure Phases of the Flight".

: "Integrity is that quality which relates to the trust which can be placed in the correctness of the information supplied by the total system. Integrity risk is theprobability of an undetected failure of the specified accuracy. Integrity includes the ability of a system to provide timely warnings to the user when the system should not be used for the intended operation."

The fundamental part of integrity monitoring and reliability assurance is application of a selected Fault Detection and Exclusion (FDE) algorithm. Here the least-squares-residuals method (LSR) [10,11] for reliability monitoring has been taken into consideration. The approach was improved by including individual weighting of the code pseudorange measurements [11]. The LSR is one of the most frequently used RAIM method classified as snapshotschemein which only current observation epochs containing redundant number of measurements are processed. In opposite to this, in the sequentialschememeasurements from the previous epochs are also taken into account (eg. Kalman filter).

2.1. Weighted Least-Squares (WLS) estimation

The main observables used in the GNSS positioning solution are distances between satellite $\mathit{s}$and receivers $\mathit{r}$, derived from signal TOA/TOF (Time of Arrival / Time of Flight). In these investigations only the code-phase observations have been used. The extended nonlinear equation for GPS code measurements$\mathbf{ }{\mathit{r}}_{\mathit{r}}^{\mathit{s}}$, expressed in meters, can be written as:

$$P_r^s={\rho }_r^s+c\left({\delta t}_r-\mathrm{ }{\delta t}^s\right)+I_r^s+Z_r^s+c\left({HD}_r-{HD}^s\right)+{ϵ}_r^s$$E1

where ${\rho }_r^s$is the geometrical distance between satellite $s$and receiver $r$. The other parameters are: ${\delta t}_r$and ${\delta t}^s$receiver and satellite clock corrections; ${HD}_r$and ${HD}^s$-code signal delays in the hardware of the receiver and satellite; $I_r^s$-ionospheric refraction; $Z_r^s$-tropospheric refraction; ${ϵ}_r^s$-observation noise. The nonlinear observations equations must be linearized around of approximate initial values of the coordinates $x_0$, and the receiver clock correction, using the Taylor series expansion, to solve for the parameters using the LS adjustment. Vector of the estimated parameters contains corrections$\mathrm{ }\mathrm{\Delta }x$, $\mathrm{\Delta }y$, $\mathrm{\Delta }z$to the approximated position $x_0$and the receiver’s clock correction ${\delta t}_r$. The linearized observation model in the matrix notation for the code-phase observation can be written as:

$$∆\rho =A∆x+\mathrm{ }∆\epsilon$$E2

Where $∆\rho$is the misclosure vector, defined as the difference between "observed-calculated" (o-c) code measurements; $A$is the geometry or so-called, designed matrix; $\mathrm{\Delta }x$are four unknowns parameters (three corrections to the coordinates and correction to the receiver clock); $∆\epsilon$is the vector containing measurement noise.

After the estimation process the user’s coordinates are obtained by correcting approximated position $x_0$using estimated incremental vector $∆\widehat{x}$:

$$\widehat{x}=x_0+∆\widehat{x}$$E3

Since $\mathrm{\Delta }\rho$has some unknown and random errors, the equation should be treated as a stochastic model. Then $\mathrm{\Delta }x$can be estimated using the weighted least-squares approach. Required weight matrix $\mathrm{\Sigma }$is assumed to be known, and estimation of x can be obtained from:

$$∆\widehat{x}={\left(A^T{\mathrm{\Sigma }}^{-1}A\right)}^{-1}{\mathrm{\Sigma }}^{-1}{A}^T\mathrm{\Delta }\rho$$E4

where $A^T{\mathrm{\Sigma }}^{-1}A$is a non-singular (invertible matrix) variance-covariance matrix of estimated parameters.

Currently, in the stochastic approach many different observation weighting models have been used. One of the simplest and most commonly used is based on the satellite elevation angle. However, in the case of ionospheric perturbations, which are not directly dependent on the elevation angle (see Figure 3), applicability of this weighting approach can be useless. Due to that in this analysis one of the signal quality parameters which described all of the imperfections in the signal has been used. Signal Quality is usually represented as signal-to-noise ratio (SNR) or as carrier-to-noise ratio (CNR). Both of those parameters are essential to assess the performance of GPS receiver and they are directly related to the precision of code-phase and carrier-phase pseudorange observations [9]. SNR is obtained at the correlator output and is described as a ratio of the signal power $S_{corr}$to the noise power $N_{corr}$of the modulated signal. CNR is obtained at the receiving antenna and is described as a ratio of the signal power $C_{ant}$to the noise power $N_{ant}$of the modulated signal. Due to the fact that signal and noise power are amplified (between antenna and correlator output) by approximately the same factor we can assume that ratio of those parameters is also almost the same (equation 5).

$$CNR=\frac{C_{ant}}{N_{ant}}\mathrm{ }\approx \mathrm{ }\frac{S_{corr}}{N_{corr}}=SNR$$E5

The matrix $\mathrm{\Sigma }$, being the weight matrix, describes the noise characteristics related to the measurements:

$$\mathrm{\Sigma }=\left|\begin{array}{ccc}{\sigma }_1^2& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\sigma }_i^2\end{array}\right|$$E6

The diagonal components of the weight matrix (6) are the variances of the code measurement ${\sigma }_i^2.$This variance model (the so called $sigma-ϵ$) for weighting of the GPS observations has been proposed by [3] and is defined as:

$${\sigma }_i^2=a+b\mathrm{ }\bullet {10}^{\frac{-C/N_0}{10}}$$E7

Where $C/N_0$is the carrier-to-noise power density ratio expressed in $dBHz$unit. For the code-phase pseudorange observations the $SNR$can be used instead of $C/N_0$as well. The SNR describes the ratio of the signal power and noise power in a given bandwidth, expressed in $dB$unit. The parameters $a$and $b$have to be chosen according to the local environment. In the paper [6,7] the authors proposed the following values of the parameters:

• for heavily degraded signal condition $a=0.01\mathrm{ }\left[\frac{m^2}{s^2}\right]$and $b=25\mathrm{ }\left[\frac{m^2}{s^2}Hz\right]$;

• for lightly degraded signal condition $a=10\mathrm{ }\left[\frac{m^2}{s^2}\right]$and $b=150\mathrm{ }\left[\frac{m^2}{s^2}Hz\right]$.

RAIM FDE techniques have been developed for reliability monitoring based on statistical tests with the aim to detect and exclude faulty measurement. This process is used for checking consistency of the measurements. It is carried out by means of a statistical hypothesis test of the residuals of a least squares estimation of the GNSS position. Method presented in this paper, the so called weighted least-square-residual approach uses weighted sum of square of the errors(WSSE) as the test value. The WSSE can be explained as a quantity used in describing discrepancy between the data and an estimation model. It has to be emphasized that the approach requires at least one redundant observation. Faulty observations on the input to the navigation solution will be then excluded [14]. In the WLS method the estimated code-pseudorange residuals are:

$$\widehat{v}=A\widehat{x}-∆\rho =-R∆\rho$$E8

where, the so-called redundancy matrix R is:

$$R=\mathrm{ }{C}_{\widehat{v}}{\mathrm{\Sigma }}^{-1}$$E9

The trace of the R represents degree of freedom for the established model. $C_{\widehat{v}}$matrix in the equation (9) is the variance-covariance matrix of the estimated residuals:

$$C_{\widehat{v}}=\mathrm{\Sigma }-A{\left(A^T{\mathrm{\Sigma }}^{-1}A\right)}^{-1}{A}^T$$E10

To obtain normally distributed $N\left(\mathrm{0,1}\right)$observations, the residuals $\widehat{v}$, must be standardized:

$Z_i=\left|\sqrt{\frac{{\widehat{v}}_i}{{\left(C_{\widehat{v}}\right)}_{ii}}}\right|,i=1\dots \text{.}n$E11

The FDE is based on statistical tests for outlier detection using null and alternative hypotheses. The presented scheme of FDE includes a global test for detection of presence of the unacceptable error, followed by a local test for exclusion of the faulty measurement.

Decision	Situation
	$H_0$is true; no blunder presents	$H_0$is false; blunder presents
$H_0$is reject; blunder detect	Type I error: probability of false alarm$\alpha$; (significance level)	Correct decision: probability of rejecting a true blunder$\left(1-\beta \right)$(power of test)
$H_0$is accept; no blunder detect	Correct decision: probability of making a correct decision$\left(1-\alpha \right)$(confidence level)	Type II error: probability of missed detection$\beta$

Table 1.

Statistical hypothesis testing with the errors classification

Two types of erroneous decision can occur [5] (Table I). The type Iof erroneous decision occurs if the null hypothesis rejects the correct observation. Error of type IIappears if the faulty observation is accepted. In the Table 1, the so-called ”significance level”$\alpha$denotes probability of committing of the type Ierrors and the $\left(1-\alpha \right)$is probability of making the correct decision and is also called ”confidence level”. The probability of committing type IIerrors is denoted by$\beta$where $\left(1-\beta \right)$is the probability of rejecting a true blunder observation, also called “power of the test”.

The threshold values for the statistical tests

The threshold values for the global and the local tests must be predefined based on following parameters:

• $\alpha$-is the false alarm probability of the global test

• ${\alpha }_0$-is the false alarm probability of the local test,

• ${\beta =\beta }_0$-is the missed detection probability, should be this some for the local and the global test.

The above parameters $\alpha$, $\beta$and ${\alpha }_0$are related by the following formula [1].

$$\lambda ={\left({\lambda }_0\right)}^2={\left(n_{1-\frac{{\alpha }_0}{2}}+n_{1-\beta }\right)}^2$$E12

$${\chi }_{\beta ,n-p,\lambda }^2={\chi }_{1-\alpha ,n-p}^2$$E13

Where $\lambda$is the non-centrality parameters of the non-centrally chi-square distribution ${\chi }^2$. Thus, only two of them can be independently chosen. Furthermore, definition of values of the selected parameters is essential for consideration the following consequences:

• The large value ${\alpha }_0$implies a smaller critical value of the local test $n_{1-\frac{{\alpha }_0}{2}}$, causing exclusion of a higher number of correct observations.

• The large value of ${\beta }_0$causes higher probability of missed detection, it means that more erroneous observation will be accepted as correct one.

Global Test-for the detection of faulty observation

Global test is used for evaluation if the set of GNSS observations include an erroneous observation. When the measurement errors are zero-mean normally distributed $N\left(0,\mathrm{\Sigma }\right)$, the testing follows the central chi-square distribution ${\chi }^2$. In such a case the test threshold value is defined by the inverse chi-square cumulative distribution function (CDF). This value is tested against the test parameters ${\widehat{v}}^T{\mathrm{\Sigma }}^{-1}\widehat{v}$, the so-called weighted sum of square error(WSSE). In the set of observations with distribution $N\left(0,\mathrm{\Sigma }\right)$the hypothesis is tested as:

$$\begin{array}{c}{H}_0:\mathrm{ }{\widehat{v}}^T{\mathrm{\Sigma }}^{-1}\widehat{v}\le {\chi }_{1-\alpha ,n-p}^2\\ H_a:\mathrm{ }{\widehat{v}}^T{\mathrm{\Sigma }}^{-1}\widehat{v}>{\chi }_{1-\alpha ,n-p}^2\end{array}$$E14

The expression $(n-p)$is degree of freedom, where n and p are number of observations and estimated parameters respectively. In a case when the result exceeds the threshold value ${\chi }_{1-\alpha ,n-p}^2$, the null hypothesis $H_0$is rejected and consequently alternative hypothesis is performed (see Fig.1). If $H_0$is rejected and $H_a$is accepted, some inconsistency of the observations exists. In such a case the so called local test, assuming that only one blunder is present in the data set, should be applied.

Local Test-for the identification of faulty observation

A failed global test indicates that there is at least one erroneous observation. In such situation, the FDE algorithm starts execution of a local test to identify and reject that observation. The standardized residuals from the equation 10 are used to conduct the local test. Those residuals are compared with the ${\alpha }_0$-quantile of the standard normal distribution $N\left(0,1\right)$. The standardized residuals are normally distributed [15] with zero expectation if the $H_{0,i}$is correct and with non-zero expectation value otherwise. The local test is based on the following hypothesis:

$$\begin{array}{c}{H}_{0,i}:\mathrm{ }\left|z_i\right|\le n_{1-\frac{{\alpha }_0}{2}}\\ H_{a,i}:\mathrm{ }\left|z_i\right|>n_{1-\frac{{\alpha }_0}{2}}\end{array}$$E15

where the probability ${\alpha }_0$was divided equally to the obtained right-tailed and left-tailed test. The null hypothesis $H_{0,i}$denotes that i-th observation is not a blunder. Thus, when the $H_{0,i}$is rejected then the alternative hypothesis $H_{a,i}$is used to recognize if the threshold value is exceeded. The excluding process of the k-th erroneous observation is based on the following formulation [15]:

$$H_{a,k}:\mathrm{ }{Z}_k\le Z_i\mathrm{ }\forall \mathrm{ }i,\mathrm{ }\wedge \mathrm{ }{Z}_k>n_{1-\frac{{\alpha }_0}{2}}\mathrm{ }$$E16

After exclusion of an erroneous observation and repetition of navigation solution, the statistical test should be repeated (Figure 1). As it was mentioned before, the process can be repeated until no more errors are detected or until the condition of redundant observation has stops being met.

Figure 1.

Fault Detection and Exclusion Scheme

3. Reliability tests with data collected in the presence of ionospheric perturbations

The investigations have been conducted using long time series of GNSS observations gathered at high latitudes 67.5 N LCKI (Kiruna/Sweden), see in the Fig. 2, where ionospheric perturbations occur more frequently and are stronger. The analysis were performed using high rate (1Hz) GPS data from the two selected days with low (2012 July, 05-Day of the Year-DOY: 187) and high (2012 July, 15-Day of the Year-DOY: 197) ionospheric perturbations. Some selected results of our investigations are presented in the figures below.

Figure 2.

Localization of the GNSS receiver (LCKI, Kiruna/Sweden), visibility of satellites and sTEC (CODE/GIM from the Center for Orbit Determination in Europe, Astronomical Institute of the University of Bern). Slant TEC (in TECU unit) has been interpolated at the ionospheric pierce point (IPP) for each observed satellites. DOY:197

3.1. Influence of ionospheric perturbations onto quality of GNSS positioning

The ionospheric irregularities are correlated with the solar radiation and the Earth’s geomagnetic field. Both of these phenomenon cause variability of electrons density in space and time. Since the ionospheric perturbations are directly associated with geomagnetic storms, for identification of those phenomena the planetary Kp-index has been used. This index represents irregular disturbances of the geomagnetic field caused by solar particle radiation [2]. The Figure 3 shows planetary Kp-index values during low-and high ionospheric perturbations.

Figure 3.

Planetary Kp-index: DOY187 (left) and DOY 197, (right) data LCKI – (Kiruna/Sweden), source: http://wdc.kugi.kyoto-u.ac.jp/kp/index.html)

Some other useful parameter used for monitoring the changes of the density of ions is continuous analysis of the time derivative of TEC [13] (ROT, rate of change of Total Electron Content). This parameter can be used for describing of direct correlation between ionospheric perturbations and GPS observations. The formula for ROT is written as:

$$ROT=\frac{{TEC}_t^s-{TEC}_{t-1}^s}{t_t-t_{t-1}}$$E17

where $s$is the visible satellite and $t$is the observation epoch. The TEC is derived from the dual frequency carrier phase measurements, and with usage of the ROT the "problem" of the carrier phase ambiguity fixing can be avoid.

Figure 4.

The relation between ROT and elevation angles presented for satellite PRN13 (DOY197), LCKI (Kiruna/Sweden)

Figure 5.

The relation between SNR (L2 frequency) and elevation angles presented for satellite PRN13 (DOY197), LCKI (Kiruna/Sweden)

In order to emphasize influences of ionospheric perturbation onto precision and accuracy of GNSS position displayed (see Figure 6), the residuals of GNSS coordinates are presented for the both selected days. This analysis shows some influences of the ionospheric disturbances on accuracy and precision of the GNSS positioning. The results obtained from the data recorded during strong ionospheric perturbations show higher degradation of precision and accuracy then the results obtained during low ionospheric perturbations. Some of the results are unacceptable by the applications where the high precision and reliability is required.

The effects of the ionospheric perturbations are also clearly seen in the Figure 7 where the presented carrier-phase residuals during high ionospheric perturbations are much bigger then carrier-phase residuals during the low ionospheric perturbations, Figure 8). The residuals "observed-calculated" (o-c) can be used to analyse observational noise level and to identify outliers in observations. Here the (o-c) values are used as an indicator describing the noise level of the carrier phase observations recorded under strong atmospheric disturbances.

Figure 6.

The confidence ellipse for the residuals of North and East coordinates, DOY:187 (upper figure) and DOY:197 (lower figure) (confidence level: 95%)-LCKI (Kiruna/Sweden)

Figure 7.

Carrier-phase residuals (o-c) per satellite DOY:187, LCKI (Kiruna/Sweden)

Figure 8.

Carrier-phase residuals (o-c) per satellite DOY:197, LCKI (Kiruna/Sweden)

An initial investigation was conducted in order to evaluate applicability of FDE method in order to improve performance of absolute positioning approaches in the presence of ionospheric perturbations. In this test the certainty levels for the reliability monitoring were predefined with the following values: the false alarm probability was set to ${\alpha }_0=10\mathrm{\%}$and the probability of missed detection was set to $\beta =20\mathrm{\%}$.

In the Figure 9 solution of Single Point Positioning method are presented. This test has been performed using a set of simulated data. A simulated data approach has been used in order verify correctness of the implemented. The red points in the left graph (Figure 9) have been identified as unreliable solution. For those points inconsistency of GPS observations the have been detected. The right graph (Figure 9) shows solutions of the some set of GSP observations using the FDE. From this figure we can see that applicability of FDE gives ability to detect and to remove the erroneous observation and consequently to improve final solution.

Figure 9.

The Single Point Positioning solution without FDE (left) and with FDE (right)

4. Summary and conclusions

The paper considered integrity monitoring at the user-level for GNSS positioning applications. In order to improve performance of absolute point positioning algorithms, applicability of the fault detection and exclusion methods have been investigated. The suggested approach can be also used to support mitigation of ionospheric threats in GNSS real-time positioning solutions. Reliability testing of the fault detection and exclusion method has been performed with the weighted least-square residuals approach. Due to the strong influence of ionospheric perturbations onto the signal noise, the stochastic model has been defined using the weight matrix containing SNR (signal-to-noise ratio) values.

In this analysis only code-phase pseudo-ranges have been processed using a real-time single-epoch approach. It has been shown that the snapshot scheme allows unequivocal identification of blunder observations in a real-time “single epoch” standard point positioning approach and it can support mitigation of ionospheric perturbation influences. Due to strong influence of the geometry of the satellite constellation onto positioning quality, analysis of the DOP (dilution of precision) parameters will be taken into consideration as well. For the purpose of future analysis a stochastic model related to information about ionospheric perturbations will be developed and applied.

Acknowledgments

The research is funded by the FP7 People Programme through the Marie Curie Initial Training Network TRANSMIT-Training Research and Applications Network to Support the Mitigation of Ionospheric Threats.

In addition, acknowledgements are given to Space Weather Application Center Ionosphere (SWACI) operated by German Aerospace Center (DLR) in Neustrelitz for providing the SWACI data sets.