A Difference Test Method for Early Detection of Slowly Growing Errors in GNSS Positioning

Shaojun Feng; Washington Ochieng

doi:10.1017/S037346330700433X

A Difference Test Method for Early Detection of Slowly Growing Errors in GNSS Positioning

Published online by Cambridge University Press: 09 August 2007

Shaojun Feng and

Washington Ochieng

Show author details

Shaojun Feng: Affiliation:
(Imperial College London)
Washington Ochieng*: Affiliation:
(Imperial College London)
*: (Email: w.ochieng@imperial.ac.uk)

Article contents

Abstract
INTRODUCTION
OVERBOUNDING OF A DISTRIBUTION
DIFFERENCE TEST FOR DETECTING SLOWLY GROWING ERRORS
RESULTS
DISCUSSIONS AND CONCLUSIONS
References

Rights & Permissions

Abstract

Satellite based navigation system integrity monitoring is essential both for liability and safety critical services. The threats to system integrity are the potential failure modes that could occur at system, operational and user sensor levels. Research on identification, characterisation and modelling of such failure modes, suggests that failures that grow slowly over time (e.g. ramp error less than 2 m/s), also referred to as slowly growing errors (SGE) are the most difficult to detect early. Conventional snapshot algorithms detect SGEs when the test statistic crosses the threshold usually after relatively long periods of time. A recent concept based on the average of residuals over time has been found to have a significant weakness. This paper proposes a “difference test” algorithm capable of detecting SGEs early. In this method, the test statistic is the difference between the norm of current residuals and the norm of the residuals at a previous epoch. The distribution of the test statistic is over-bounded by a normal distribution whose parameters are derived from two Chi-distributions. Results show that the new algorithm results in a significant reduction in detection time compared to conventional methods.

Keywords

GNSS System Integrity Slowly growing errors

Information

Type: Research Article
Information: The Journal of Navigation , Volume 60 , Issue 3 , September 2007 , pp. 427 - 442

DOI: https://doi.org/10.1017/S037346330700433X [Opens in a new window]
Copyright: Copyright © The Royal Institute of Navigation 2007

1. INTRODUCTION

Integrity is a measure of trust, which can be placed in the correctness of information supplied by the total system. Integrity includes the ability of a system to provide timely and valid warnings to the user. This is vital for liability and safety critical applications. Integrity risk (loss of integrity) of a navigation system is defined as the probability that a user will experience a position error (PE) larger than the alert limit (AL) without an alert being raised within the specified time-to-alert (TTA) at any instant in time and at any location in the coverage area. Loss of integrity can happen in three ways (Ochieng, Reference Ochieng2003; Tiemeyer, Reference Tiemeyer2002): an unsafe condition is not detected, an unsafe condition is detected but the isolation function removes the wrong satellite, and an unsafe condition is detected but annunciation takes longer than the TTA.

The integrity of a global navigation satellite system (GNSS) can be monitored at the independent network level (e.g. the European geostationary navigation overlay service, EGNOS), system level (e.g. Galileo), and user level. The main approach at the user sensor level is receiver autonomous integrity monitoring (RAIM). An extension of RAIM which accommodates the requirement for additional sensors to improve the performance of a GNSS such as GPS, is referred to as user level autonomous integrity monitoring (UAIM) (Feng et al Reference Feng, Ochieng, Walsh and Ioannides2006a).

Integrity monitoring at the user level can be divided into two parts: determination of the protection level (PL), and failure detection and exclusion (FDE). The PL is the upper bound of the positioning error, which is compared against the alert limit (AL) to determine if an alert should be raised. The FDE algorithm is used to detect unsafe conditions (including abnormal error in measurements). Once detected and if the exclusion function is available, then the algorithm excludes the failed measurement(s) and uses the remaining ones for positioning. If the remaining measurements cannot meet the minimum requirements (either geometry or number of visible satellites), an alert is raised because the conditions necessary for the execution of the detection function are inadequate (this is usually referred to as RAIM unavailability).

In general, failure detection algorithms are based on a number of assumptions, the most important of which is that residual errors in the measurements are normally distributed. Failure detection consists of three main steps: the construction of a test statistic; the characterisation of the test statistic and the determination of a threshold to reflect the user requirement in terms, for example, of probability of false alert; and decision making. One of the key features in the design of any integrity algorithm is its sensitivity to various types of failure modes. For positioning using GNSS, a significant error such as clock jump can be easily detected by a snapshot algorithm, while the most difficult types of failure to detect are those that grow slowly over time (e.g. clock drift at less than 2 m/s) (Lee et al Reference Lee and O'Laughlin1999). Such failures are of particular concern in filtering based UAIM because the Kalman filter tends to adapt to them. This results in the positioning solution being contaminated with the consequence of misleading information both in terms of accuracy and integrity. For this reason, early detection of this type of failure is vital.

Previous research on the detection of slowly growing errors (SGEs) has explored the suitability of the estimation of the mean of the residuals over time from the innovation of the Kalman Filter (Diesel et al Reference Diesel and Luu1995), culminating in the so-called Extrapolation Method (Lee et al Reference Lee and O'Laughlin1999). The process of estimating the mean involves storing measurements and residuals in buffers for a period of 30 minutes or more and estimating the mean of residuals over a selected time interval. The test statistic used is then the sum-squared of the estimated mean residuals which has a chi-square distribution if there is no failure and a non-central chi-square distribution otherwise (Diesel et al Reference Diesel and Luu1995). By taking a relatively long averaging interval, the covariance of the averaged residual should be small in the absence of failure and large otherwise. This property, in theory, makes it possible to detect SGEs. However, there is an issue with the degrees of freedom of the distribution if the number of satellites changes during the period of time used in the averaging process. In this case, the threshold determined from the distribution with either the current or an average number of visible satellites has the potential to result in a higher number of missed detections, false alerts or both. Any decision based on the threshold for example, does not reflect accurately the continuity risk (from which the probability of false alert is derived). Put simply, the averaging technique suffers from the fact that in some instances, the distributions used to determine the threshold and the test statistic may not match. It is important that the distribution used to determine the threshold should at least overbound the distribution of the test statistic. Given this significant weakness, the averaging concept is not discussed further in this paper.

The other approach is known as the ratio test which is based on the F distribution and is normally used to test whether there is a significant difference between two independent variants. The ratio test may be applicable to the comparison of two variants, one of which has to be fault-free. If the failure (ramp error) occurs before the test, the ratio test cannot give a correction decision. Therefore, the ratio test is not discussed further in this paper.

Given the background above, this paper proposes a “difference test” algorithm for the detection of slowly growing errors. A new test statistic is derived from the difference between the norm of the residuals at the current epoch and a previous epoch. The distribution that overbounds the test statistic is defined and justified by simulations. Section 2 of the paper gives background statistical information on overbounding. Section 3 develops the “difference test” method including the methods proposed for the determination of the test statistic, the distribution and the threshold. Section 4 looks at the application of the algorithm to early detection. Section 5 presents and discusses the results. The paper is concluded in Section 6.

2. OVERBOUNDING OF A DISTRIBUTION

The distribution of a test statistic could be non-regular which cannot be fully described by a typical distribution (e.g. normal, student etc.). Overbounding is to ensure that the probability of a known distribution is higher than the actual distribution at least for a part of a distribution such as the tail. There are normally two types of overbounding: the Probability Density Function (PDF) and Cumulative Distribution Function (CDF) (DeCleene, Reference DeCleene2000; Ober et al Reference Ober, Farnworth and Breeuwer2001). The other techniques of overbounding are largely based on these e.g. excess-mass functions (Rife et al Reference Rife, Walter and Blanch2004).

In general, PDF overbounding can be expressed as:

(1)

$p_{T} \lpar x\rpar \gt p_{A} \lpar x\rpar \quad \forall x \in \lpar a\comma b\rpar$

where p _T(x) is the PDF of a typical distribution, p _A(x) is the PDF of the actual distribution, a and b define the interval. If a=−∞, then p _T(x) overbounds p _A(x) on the left tail. If b=∞, then p _T(x) overbounds p _A(x) on the right tail.

Similarly, in general, CDF overbounding can be expressed as

(2)

$P_{T} \lpar x\rpar \gt P_{A} \lpar x\rpar$

(3)

$1 \minus P_{T} \lpar x\rpar \gt 1 \minus P_{A} \lpar x\rpar$

where P _T(x) is the CDF of a typical distribution, P _A(x) is the CDF of the actual distribution. Expression (2) stands for P _T(x) overbounding P _A(x) on the left tail while expression (3) stands for P _T(x) overbounding P _A(x) on the right tail.

The CDF overbounding and PDF overbounding are correlated. If p _A(x) is overbounded by p _T(x) on the left tail (a=−∞), then P _A(x) is overbounded by P _T(x) on the left tail (x∈(−∞,b)). If p _A(x) is overbounded by p _T(x) on the right tail (b=∞), then P _A(x) is overbounded by P _T(x) on the right tail (x∈(a,∞)).

In GNSS integrity monitoring, the normal distribution is often used to overbound non-Gaussian distributed errors (e.g. multipath errors). Therefore, appropriate normal distribution in terms of overbounding is a key research issue. Suppose we have two normal distributions N ₁∈N(μ₁,σ₁) and N ₂∈N(μ₂,σ₂). Three typical scenarios are shown in Figures 1a–c where one can bound the other at some intervals. The intersection points A and B are the critical points for one distribution to overbound the other. Figure 1a demonstrates that N ₂ overbounds N ₁ in PDF at (−∞,a) and (b,∞) under the conditions of μ₁=μ₂ and σ₁<σ₂. This can be referred to as standard deviation overbounding. Figure 1b demonstrates that N ₂ overbounds N ₁ in PDF at (a,∞) under the conditions of μ₁<μ₂ and σ₁=σ₂. This can be referred to as mean overbounding. Figure 1c demonstrates that N ₂ overbounds N ₁ in PDF at (−∞,a) and (b,∞) under the conditions of μ₁≠μ₂ and σ₁<σ₂.

Figure 1a.

Overbounding of same mean and different standard deviations.

Figure 1b.

Overbounding of same standard deviation and different means.

Figure 1c.

Overbounding of different means and different standard deviations.

The critical values a and b can be resolved from the following equation (Abramowita et al, Reference Abramowita and Stegun1970)

(4)

${1 \over {\sigma _{\setnum{1}} \sqrt {2\pi } }}e^{{{ \minus \lpar x \minus \mu _{\setnum{1}} \rpar ^{\setnum{2}} } \over {\setnum{2}\sigma _{\setnum{1}} ^{\setnum{2}} }}} \equals {1 \over {\sigma _{\setnum{2}} \sqrt {2\pi } }}e^{{{ \minus \lpar x \minus \mu _{\setnum{2}} \rpar ^{\setnum{2}} } \over {\setnum{2}\sigma _{\setnum{2}} ^{\setnum{2}}}}}$

Let m=σ₁²−σ₂², p=2σ₂²μ₁−2σ₁²μ₂, $q \equals \sigma _{\setnum{1}}^{\setnum{2}}\hskip2\mu _{\setnum{2}}^{\setnum{2}} \minus \sigma _{\setnum{2}}^{\setnum{2}}\hskip2 \mu _{\setnum{1}}^{\setnum{2}} \minus 2\sigma _{\setnum{1}}^{\setnum{2}} \hskip2\sigma _{\setnum{2}}^{\setnum{2}} \,{\rm In}\left( {{{\sigma _{\setnum{1}} } \over {\sigma _{\setnum{2}} }}} \right)$ , d=p ²−4mq.

If σ₂>σ₁, and d⩾0, the intersection points (i.e. A,B in Figure 1) can be determined to be:

(5)

$\left\{\openup2 {\matrix{ {a \equals {\displaystyle{{ \minus p \minus \sqrt d } \over {2m}}}} \cr {b \equals {\displaystyle{{ \minus p \plus \sqrt d } \over {2m}}}} \cr} } \right.$

N(μ₂,σ₂) is used to overbound N(μ₁,σ₁) in the tails, with the overbounding K factor being:

(6)

$\eqalign{\tab K\left\{\hskip-5 {\matrix{ { \lt K_{\setnum{1}} } \tab {{\rm \lpar left\ tail\rpar }} \cr { \gt K_{\setnum{2}} } \tab {{\rm \lpar right\ tail\rpar }} \cr} } \right. \cr \tab {\rm \ \ } \cr}$

where $K_{\setnum{1}} \equals {\textstyle{{x_{\setnum{1}} \minus \mu {}_{\setnum{2}}} \over {\sigma _{\setnum{2}} }}}$ , and $K_{\setnum{2}} \equals {\textstyle{{x_{\setnum{2}} \minus \mu {}_{\setnum{2}}} \over {\sigma _{\setnum{2}}}}}$ which reflect the overbounding probability.

In the case of σ₁=σ₂ (i.e. m=0), the only intersection point is $x \equals \minus {q \over {p}}$ and the corresponding K factor is:

(7)

$K \equals {\textstyle{{x \minus \mu {}_{\setnum{2}}} \over {\sigma _{\setnum{2}} }}}.$

The analysis above suggests that a normal distribution can be overbounded on the left tail by another normal distribution with a larger sigma value and a left offset in the mean. The same is true for overbounding of the right tail by another normal distribution with a larger sigma value and a right offset in the mean.

3. DIFFERENCE TEST FOR DETECTING SLOWLY GROWING ERRORS

Failure detection in RAIM is based on statistical consistency checks using redundant measurements. There are two different RAIM schemes for use with measurements; snapshot and filtering (Brown, 1996). In the snapshot scheme only the current redundant measurements are used to check measurement consistency. In snapshot methods the test statistic is derived from the relationship between measurements and position given by:

(8)

$y \equals Gx \plus \varepsilon$

where

y is the difference between the actual measured parameter (e.g. pseudorange) and the predicted parameter based on a nominal user position and clock bias. It is a n×1 vector. n is the number of measurements.
G is the design matrix (i.e. the observation matrix in a local horizontal coordinates system, e.g. ENU: East-North-Up).
x is the vector of the unknown parameters (i.e. three components of true position plus a user clock bias).
ε=[ε₁ ε₂ … ε_n]^T is the measurement error vector caused by receiver noise, wave propagation, ephemeris etc. The elements of the error vector are assumed to be independent zero-mean Gaussian random variables with the same variance σ_m.

The square root of the sum of squared errors (SSE) is taken as a test statistic, which can be calculated by the least squares or parity methods. The SSE can be obtained from the residual vector (v) by:

(9a)

$SSE \equals v^{T} v$

where v=[I−G(G ^TG)⁻¹G ^T]ε, I is an identity matrix (Feng et al Reference Feng, Ochieng, Walsh and Ioannides2006b). In a weighted mode where each measurement has a different error (σ), the SSE is expressed by:

(9b)

$SSE \equals v^{T} Wv$

Where W is the inverse of the measurement covariance matrix. Based on the assumption that the measurements errors are normally distributed with zero-mean, the SSE follows the Chi-square distribution with (n−4) degrees of freedom (Parkinson et al, Reference Parkinson and Axelrad1988). Therefore, the range residual parameter (norm of the residual), $\sqrt {SSE}$ , follows a Chi-distribution with (n−4) degrees of freedom.

In the case of filtering methods, the test statistic is derived from the “innovation” of the Kalman filter as follows:

(10)

$r \equals z \minus Hx_{ \minus}$

where r is the innovation vector, z is the current measurement vector, H is the observation matrix, and x_ is the state of the previous epoch. The test statistic is:

(11)

$T \equals \sqrt {r^{\it T} Wr} \equals \sqrt {SSE}$

where W=(HP ⁽⁻⁾H ^T+R)⁻¹, P ⁽⁻⁾ is the predicted covariance, and R is the measurement noise matrix in the Kalman filter.

Based on the same assumption as snapshot methods for the measurements errors, the SSE follows a Chi-square distribution with degrees-of-freedom of (n) (Diesel et al, Reference Diesel and Luu1995; Lee et al, Reference Lee and O'Laughlin1999). Therefore, the range residual parameter (norm of the residual), $\sqrt {SSE}$ , follows a Chi-distribution with (n) degrees of freedom. One important feature is that the statistical distributions of SSE and $\sqrt {SSE}$ are completely independent of the satellite geometry for any (n) based on the assumption above under fault-free conditions.

Neither the snapshot nor filtering methods are designed for detecting slowly growing error. The concept proposed based on the mean of the residual to detect SGEs has been shown in the introduction section to have a fundamental weakness related to the change of the number of satellites. This section proposes a new approach based on the difference of the norm of the residuals between two epochs.

3.1. Test statistic

Normally GPS exhibits long term stability with normal residual measurement noise. However, in the presence of SGEs, the GPS residuals exhibit a rate of growth. Therefore, a test statistic based on the difference between the conventional test statistics at different epochs has the potential to enable the detection of SGEs (as illustrated in Figure 4). This paper proposes a “difference test” algorithm to check whether there is a significant difference between the current residuals and residuals at a previous epoch. This time period depends on the characteristics of the different failure mode. The test statistic is thus the difference between the two sets of residuals. In the basic mode, the test statistic can be expressed as:

(12a)

$T_{\rmDelta t} \equals \sqrt {{\textstyle{{SSE_{t} } \over {\sigma _{m} }}}} \minus \sqrt {{\textstyle{{SSE_{t \minus \rmDelta t} } \over {\sigma _{m}}}}}$

In the weighted mode, the test statistic can be expressed as:

(12b)

$T_{\rmDelta t} \equals \sqrt {SSE_{t} } \minus \sqrt {SSE_{t \minus \rmDelta t}}$

where t is current time and Δt is the time interval selected. The corresponding degrees of freedom are denoted as dof ₂ and dof ₁ respectively. σ_m is the standard deviation of the measurement error which is a constant value determined with the prior-knowledge of the measurement characteristics. Any change in the σ_m value based on previous information could be contaminated by the failures we are trying to detect. The test statistic is based on a moving window and always captures the difference between the two edges of the window. As described above, the norm of the residual follows a Chi-distribution. Therefore, the test statistic (T _Δt) is in fact the difference of two Chi-distributed variants. In order to determine the threshold, the distribution of the test statistic must be known or determined.

3.2. Distribution and threshold determination

The mean (μ), standard deviation (σ), skewness (γ₁) and kurtosis (γ₂) of a Chi-distribution can be expressed as:

(13)

$\mu \equals {{\sqrt 2 \rmGamma \lpar {\textstyle{1 \over 2}}\lpar n \plus 1\rpar \rpar } \over {\rmGamma \lpar {\textstyle{1 \over 2}}n\rpar}}$

(14)

$\sigma ^{\setnum{2}} \equals {{2\left( {\rmGamma \lpar {\textstyle{1 \over 2}}n\rpar \rmGamma \lpar {\textstyle{1 \over 2}}n \plus 1\rpar \minus \rmGamma ^{\setnum{2}} \lpar {\textstyle{1 \over 2}}\lpar n \plus 1\rpar \rpar } \right)} \over {\rmGamma ^{\setnum{2}} \lpar {\textstyle{1 \over 2}}n\rpar }} \equals n \minus \mu ^{\setnum{2}}$

(15)

$\gamma _{\setnum{1}} \equals {\textstyle{\mu \over {\sigma ^{\setnum{3}} }}}\lpar 1 \minus 2\sigma ^{\setnum{2}} \rpar$

(16)

$\gamma _{\setnum{2}} \equals {\textstyle{2 \over {\sigma ^{\setnum{2}} }}}\lpar 1 \minus \mu \sigma \gamma _{\setnum{1}} \minus \sigma ^{\setnum{2}} \rpar \minus 3$

where $\rmGamma \lpar n\rpar \equals \int_{\setnum{0}}^{\infty } {t^{n \minus \setnum{1}} } e^{ \minus t} dt$ is a Gamma function. Therefore, the theoretical mean of the difference of two Chi-distributed variants can be expressed as:

(17)

$\mu _{d} \equals \mu _{\setnum{2}} \minus \mu _{\setnum{1}} \equals {{\sqrt 2 \rmGamma \lpar {\textstyle{1 \over 2}}\lpar dof_{\setnum{2}} \plus 1\rpar \rpar } \over {\rmGamma \lpar {\textstyle{1 \over 2}}dof{}_{\setnum{2}}\rpar }} \minus {{\sqrt 2 \rmGamma \lpar {\textstyle{1 \over 2}}\lpar dof_{\setnum{1}} \plus 1\rpar \rpar } \over {\rmGamma \lpar {\textstyle{1 \over 2}}dof_{\setnum{1}} \rpar }}$

where, dof ₁ and dof ₂ are the degrees-of-freedom of two Chi-distributions respectively.

The variance of the difference of two Chi-distributed variants can be expressed through error propagation as:

(18)

$\sigma _{d}^{\setnum{2}} \equals \sigma _{\setnum{1}}^{\setnum{2}} \plus \sigma _{\setnum{2}}^{\setnum{2}}$

where σ₁² and σ₂² are the value of expression (14) with (n) equal to dof ₁ and dof ₂ respectively.

The variance of the difference of two Chi-distributions converges from 0·3634 (dof=1) to 0·5 as a result of increasing degrees of freedom. Therefore, the variance of the test statistic (T _Δt) converges to 1·0 as a result of increasing both dof ₁ and dof ₂. The minimum value is 0·7268 if both dof ₁ and dof ₂ take the value of 1.

Skewness is a measure of the asymmetry of the probability distribution. Roughly speaking, a distribution has a positive skew (right-skewed) if the right (higher value) tail is longer that the left tail and a negative skew (left-skewed) if the left (lower value) tail is longer than the right tail. The value of the skewness is zero for normal distribution where the probability distribution is symmetrical. The skewness of the distribution of the test statistic (T _Δt) is a function of the two degrees of freedom (dof ₂, dof ₁ and dof ₂−dof ₁).

Kurtosis is the degree of peakedness of a distribution. For the normal distribution, the kurtosis is 0. It is difficult to figure out the analytical solution in terms of kurtosis for the distribution of the test statistic (T _Δt). However, it can be approximated by an empirical expression using simulation based studies.

3.2.1. Simulation based analysis

This section relies on the assumption of normally distributed GPS residual errors and a conventional Chi-distributed test statistic to attempt to figure out the distribution for the new test statistic based on the difference between two variants of Chi-distribution. Hence the simulation presented is purely mathematical and scalable, and does not need to be based on GPS data.

A number of Monte Carlo simulations were carried out to develop the distribution. The simulations were based on 100 000 samples each of Chi-distribution for the difference (dof ₂−dof ₁) of two sets of degrees-of-freedom ranging from 1 to 20. Figure 2a shows the residual of the mean of the simulated test statistic compared to the theoretical value (μ_d, expression 17). The values are very small with a minimum of −0·0092 and a maximum of 0·0086. The colour scheme (also in Figures 2b, 2c and 2e) emphasises the variation in three dimensions between the two sets of degrees-of-freedom and the parameter under investigation. Figure 2b shows the standard deviation (the square root of variance) of the simulated test statistic. It demonstrates that the standard deviation increases with both dof ₂ and dof ₁ with a minimum of 0·8506 and a maximum of 0·9982 both within the upper bound of 1. This increase is expected due to properties of measurement noise addressed in expressions (14) and (18).

Figure 2a.

The residual of mean.

Figure 2b.

The standard deviation.

Figure 2c.

The skewness.

Figure 2d.

The relation of skewness and dof ₂−dof ₁.

Figure 2e.

The kurtosis.

Figure 2c shows the skewness of the simulated test statistic samples. Skewness is strongly related to the value of dof ₂−dof ₁ with a minimum of −0·1974 and a maximum of 0·2047. The larger the absolute value of (dof ₂−dof ₁) the larger the skewness, and vice versa, while the sign determines whether the distribution is left or right skewed. The absolute value of the skewness tends to be large when either dof ₁ or dof ₂ has a small value (i.e. less than 5). The maximum value of the skewness is 0·0748 and the minimum value is −0·0671 when dof ₂⩾5 and dof ₁⩾5 as shown in Figure 2d. In general, the skewness is close to that of the normal distribution.

Figure 2e shows the kurtosis of the simulated test statistic samples. There is no significant offset to the normal distribution. The figure suggests that the kurtosis has a characteristic that is inverse to the standard deviation (see Figure 2b).

3.2.2. Determination of the distribution for detection

Based on the analysis of the mean, the standard deviation, the skewness, and the kurtosis, the distribution of the test statistic (T _Δt) can be approximated by a normal distribution N(μ_d,σ_d) if the value of skewness is zero. If the skewness is considered, N(μ_d,σ_d) overbounds the right tail if the skewness is negative, and overbounds the left tail if the skewness is positive as shown in Figure 3a. Therefore, the test statistic can be overbounded by a normal distribution N(μ_d,1) if the value of skewness is close to zero where (dof ₂−dof ₁) is small in value. This reflects actual GNSS positioning in most environments where the number of visible satellites is around seven with a small possibility of significant deviation in the short term.

Figure 3a.

Distribution case 1.

Figure 3b.

Distribution case 2.

To have a general and conservative expression of the normal distribution that overbounds the test statistic, the skewness is used for the determination of the mean value while the standard deviation is taken to be 1 since all the values of σ_d are less than one. Hence, the proposed mean value is determined as:

(19)

$\mu _{D} \equals \mu _{d} \minus \gamma _{\setnum{1}d}$

Where μ_d is the theoretical mean, $\gamma _{\setnum{1}d} \equals \lpar \lpar \gamma _{\setnum{12}} \rpar ^{{\textstyle{\setnum{1} \over \setnum{3}}}} \minus \lpar \gamma _{\setnum{11}} \rpar ^{{\textstyle{\setnum{1} \over \setnum{3}}}} \rpar \sigma _{d} \sol 2$ is the conservative offset factor with γ₁₂ and γ₁₁ being the skewness corresponding to dof ₂ and dof ₁ respectively. The offset applied to the mean value ensures that the absolute value of μ_D is always larger than μ_d while keeping the sign the same as that of μ_d. Therefore, the test statistic is overbounded by both the mean and standard deviation. A decision threshold can then be determined from the normal distribution N(μ_D,1) by taking account of the required navigation performance (RNP), specifically, the integrity and continuity risk from which the probability of missed detection (P _MD) and the probability of false alert (P _FA) can be derived. (Note that RNP is a concept endorsed by the International Civil Aviation Organisation (ICAO), and is a statement of the navigation performance necessary for operation within a defined airspace) (Ochieng et al, Reference Ochieng2003).

The relationship between threshold T and P _FA, is (Feng et al, Reference Feng, Ochieng, Walsh and Ioannides2006b):

(20)

$P_{FA} \equals \int\limits_{T}^{\infty } {f\left( x \right){\rm \ }} dx$

where f(x) is the probability density function of normal distribution in the fault-free case. The threshold can be calculated for a known P _FA.

3.3. Early detection and identification of slowly growing error

One important factor in the test statistic constructed in section 3.1 is the time interval Δt. The choice of time interval depends on the rate of error growth and the length of the data buffer designed. Obviously, the slower the error growth rate, the longer the time interval required to detect the error.

To detect and identify SGEs, a multiple window scheme is developed below where a number of test statistics (e.g. three in Figure 4, T _Δt1, T _Δt2, T _Δti) are used for different time intervals (Δt ₁, Δt ₂, Δt _i). For each test statistic, a threshold can be derived from the distribution defined in section 3.2 (i.e. (TH _Δt1,TH _Δt2,TH _Δti)). Figure 4 demonstrates the three test statistics case.

Figure 4.

Formations of multiple test statistics.

The scenarios for failure detection are as follows:

● If T _Δt3>TH _Δt3, while T _Δt2<TH _Δt2 and T _Δt1<TH _Δt1, a very slow ramp error is detected;
● If T _Δt3>TH _Δt3, and T _Δt2>TH _Δt2 while T _Δt1<TH _Δt1, a slow ramp error is detected;
● If all the test statistics are larger than the corresponding thresholds, a relatively fast growing error is detected;
● If all the test statistics are less than the corresponding thresholds, no growing error is detected.

This difference test can be used to detect a failure significantly earlier compared to the conventional methods. Results demonstrate that it is able to detect a SGE about 20 seconds earlier than the conventional method as shown in Figure 6. This is crucial for safety critical applications such as aviation where the time-to-alert ranges from 5 minutes for the En-route phase of flight to 6 seconds for Category I precision approach with even more stringent requirements expected for Category II and Category III (ICAO, 2006).

Further reduction in the time taken to detect a failure can be achieved by applying the difference test to a number of sequential epochs. This is based on the hypothesis that if there is no failure, the test statistics are independent for sequential epochs. Therefore, for a given threshold, the joint probability of false alert is:

(21)

$P_{FA} \equals P_{FA\setnum{1}} \times P_{FA\setnum{2}} \times \ldots \times P_{FAn}$

where P _FA is the joint probability of false alert which reflects the user defined RNP, P _FAi (i=1 … n) are the probabilities of false alert for each sequential epoch. If the values of P _FAi are taken to be the same, the following equation holds:

(22)

$P_{FA\setnum{1}} \equals P_{FA\setnum{2}} \equals \ldots \equals P_{FAn} \equals \lpar P_{FA} \rpar ^{{\textstyle{\setnum{1} \over n}}}$

Figure 5 illustrates early detection where n=3. The single epoch detection requires a large threshold (TH ₂) to reflect the RNP (e.g.1×10⁻⁶/sample false alert), while detection based on three sequential epochs requires a small threshold (TH ₁ for 1×10⁻²/sample false alert) to reflect the same RNP parameters. If a measurement is faulty, the above hypothesis is breached. A failure is claimed if all the three sequential epochs detect a failure instead of waiting until the test statistic is large enough to exceed the single epoch threshold TH ₂.

Figure 5.

Early detection scheme.

The maximum number of sequential epochs is based on the K factor using the algorithm in section 2 and the parameters determined in this section above. The corresponding probability of false alert is the maximum value to be taken as P _FAi for each sequential epoch. Therefore

(23)

$n_{\max } \equals floor\lpar \log \lpar P_{FA} \rpar \sol \log \lpar P_{FAi} \rpar \rpar$

where floor(x) means rounding to the nearest lower integer. The n _max is therefore used to determine the actual P _FAi from expression (22). Another benefit is that the scheme can give the confidence level at the first early detection point (e.g. 97·85% at point 1) and subsequent detection points (e.g. 99·95% at point 2).

4. RESULTS

Trial data (pseudorange measurements) contaminated by a simulated slowly growing error are used to verify the proposed methods. The trial was carried out near Aberporth, United Kingdom on 13 August 2005. The positioning algorithm uses a Kalman filter.

To demonstrate the performance of the algorithms, different time intervals for the “difference test” and different ramp errors are applied to PRN 21 at the 1500th second from epoch 0. The required probability of false alert is 3·333×10⁻⁷/sample. Figure 6 shows the test statistic (difference test) at a time interval of 360 seconds for a ramp error of 0·2 m/s. The conventional test statistic is also shown in the figure. The thresholds for the “difference test” are the same if the differences in the number of visible satellites are the same, and vice versa. The comparison of the conventional method and the difference test shows that the latter is able to detect the failure significantly earlier than the former.

Figure 6.

Comparison of the conventional method and the difference test.

Figure 7 shows the test statistic (difference test) at time intervals of 120, 240, and 360 seconds for a ramp error of 0·2 m/s. It compares the conventional method with the difference test implementing the new early detection scheme. The algorithm with the shortest time interval (120 s) is not sensitive to this error (test 1 in the figure). On the other hand the algorithm with the longer time intervals (i.e. 240 s and 360 s) is sensitive to this error (test 2 and test 3 in the figure), and can detect the failure much earlier (about 180 s) than the conventional method.

Figure 7.

Comparison of the conventional method and the difference test with the early detection scheme.

Figure 8 shows the test statistic (difference test) at time intervals of 300, 600, and 900 seconds for a ramp error of 0·2 m/s. All the three tests are able to detect the failure much earlier than the conventional methods. In comparison with the results shown in Figure 7, it is clear that the algorithm proposed is sensitive to the time interval used.

Figure 8.

Comparison at different time intervals.

Figure 9 shows the test statistic (difference test) at time intervals of 300, 600, and 900 seconds for a ramp error of 0·05 m/s. Only the algorithm with the longest time interval (test 3) is able to detect the error earlier than the conventional method. The algorithm with a medium time interval is at the critical point where the test statistic is around the threshold (test 2). The algorithm with the shortest time interval does not detect the error. In comparison with the results shown in Figure 8, it is clear that the algorithm is sensitive to the rate of the ramp as well as emphasising the need for a good understanding of failure modes of integrity.

Figure 9.

Comparison at different rates of error growing.

Table 1 presents the key quantities used in the comparison of the conventional method and the difference test approach. This table also shows the detectability in relation to time interval and failure mode (in terms of rate of ramp).

Table 1.

Comparison of results.

5. DISCUSSIONS AND CONCLUSIONS

The results in section 4 demonstrate the efficiency and sensitivity of the new algorithms in the early detection of the most difficult type of failure, the slowly growing errors. The length of the time interval used is the key factor that affects the sensitivity of the algorithm in terms of the minimum ramp error and the initialization time (i.e. the time it takes to start to execute a RAIM function). The multiple time interval algorithm can deal with the sensitivity issues completely and is also able to identify the ramp at least at a qualitative level (fast or slow).

The main contribution of this paper is the new “difference test” RAIM algorithm and the development of the distribution of the new test statistic, used to calculate the threshold. A slowly growing error detection algorithm is proposed based on the “difference test” with multiple time intervals. The “difference test” can detect a potential failure significantly earlier than conventional methods. It can be enhanced by the proposed early detection scheme to enable a further gain in the reduction in the time taken to detect a failure.

In terms of overall integrity monitoring, the first task is the characterisation of various failure modes, followed by implementation of appropriate algorithms. In this, a minimum configuration might involve the use of the difference test algorithm with carefully selected “differencing” intervals together with conventional methods in order to detect failures of different characteristics. Care should be exercised in the cases of variable rate or intermittent SGEs, with an example being multipath, which has a cyclic behaviour. In such cases the performance of the proposed algorithm will depend not only on the characteristics of the failure modes and the length of the time interval, but also the number of intervals.

ACKNOWLEDGEMENTS

This paper is the result of research carried out within SPACE project. SPACE (Seamless Positioning in All Conditions and Environments) is an EPSRC-funded collaborative research project. It is being undertaken by a consortium of eleven UK university and industrial groups. The university groups are the CAA Institute of Satellite Navigation at the University of Leeds, the Department of Civil and Environmental Engineering at Imperial College London, the Department of Geomatic Engineering at University College London (UCL), and the Institute of Engineering Surveying and Space Geodesy at University of Nottingham, and the industrial partners are the Civil Aviation Authority, EADS Astrium, Leica Geosystems, Nottingham Scientific Ltd, Ordnance Survey, QinetiQ and Thales Research and Technology.

References

REFERENCES

Abramowita, M., Stegun, I. A. (1970). Handbook of mathematical functions with formulas, graphs and mathematical tables. Dover Publications, Inc., New York.Google Scholar

DeCleene, B. (2000), Defining Pseudorange Integrity-Overbounding, Proceedings of ION GPS 2000, Salt Lake City, UT, 19–22, September, 1916–1924Google Scholar

Diesel, J., and Luu, S., (1995) GPS/IRS AIME: Calculation of thresholds and protection radius using Chi-square methods, Proceedings of ION GPS 1995, Palm Springs, CA, 12–15, September, 1959–1964Google Scholar

Feng, S., and Ochieng, W. Y., (2006a), User Level Autonomous Integrity Monitoring for Seamless Positioning in All Conditions and Environments, Proceedings of the European Navigation Conference, Manchester, 7–10, May, 2006Google Scholar

Feng, S., Ochieng, W. Y., Walsh, D. and Ioannides, R., (2006b), A measurement domain receiver autonomous integrity monitoring algorithm, GPS Solutions Journal, 10 (2), 85–96CrossRef Google Scholar

ICAO (2006), Annex 10 Aeronautical Telecommunications, Volume I Radio Navigation Aids (Sixth Edition), 23/11/2006, 3–70.Google Scholar

Lee, Y. C., O'Laughlin, D. G., (1999). Performance analysis of a tightly coupled GPS/inertial system for two integrity monitoring methods, Proceedings of ION GPS 1999, Nashville, Tennessee, 14–17, September, 1187–1197Google Scholar

Ober, P. B., Farnworth, R., Breeuwer, , et al. , (2001), SBAS Integrity Verification, Proceedings of ION GPS 2001, Salt Lake City, UT, 11–14Google Scholar

Ochieng, W. Y. et al. . (2003). GPS integrity and potential impact on aviation safety. The Journal of Navigation, 56, 51–65CrossRef Google Scholar

Parkinson, B. W., and Axelrad, P., (1988), Autonomous GPS Integrity Monitoring Using the Pseudorange Residual, NAVIGATION: Journal of the Institute of Navigation, 35(2), pp255–274Google Scholar

Rife, J., Walter, T., Blanch, J., (2004), Overbounding SBAS and GBAS Error Distributions with Excess-Mass functions, Proceedings of the 2004 International Symposium on GNSS/GPS, Sydney, AustraliaGoogle Scholar

Tiemeyer, B. (2002). Performance evaluation of satellite navigation and safety case development. http://137.193.200.177/ediss/tiemeyer-bernd/inhalt.pdf Google Scholar