A Gaussian Sum Filter Approach for DGNSS Integrity Monitoring

Youngsun Yun; Ho Yun; Doyoon Kim; Changdon Kee

doi:10.1017/S0373463308004906

A Gaussian Sum Filter Approach for DGNSS Integrity Monitoring

Published online by Cambridge University Press: 02 October 2008

Youngsun Yun ,

Ho Yun ,

Doyoon Kim and

Changdon Kee

Show author details

Youngsun Yun*: Affiliation:
(Seoul National University, Republic of Korea)
Ho Yun: Affiliation:
(Seoul National University, Republic of Korea)
Doyoon Kim: Affiliation:
(Seoul National University, Republic of Korea)
Changdon Kee: Affiliation:
(Seoul National University, Republic of Korea)
*: (E-mail: zoro@snu.ac.kr)

Article contents

Abstract
INTRODUCTION
DGNSS PSEUDORANGE ERROR MODEL
SIGMA OVERBOUNDING
GAUSSIAN SUM FILTER
INTEGRITY-MONITORING ALGORITHMS
SIMULATION
CONCLUSIONS
References

Rights & Permissions

Abstract

With conventional snapshot RAIM algorithms, it is difficult to detect small errors and simultaneous multiple faults. Assuming that we know the system dynamics, filtering algorithms, such as the Kalman filter, can provide better integrity-monitoring performance than the snapshot algorithms because the filter reduces the noise level of measurements. However, because the Kalman filter presumes that measurement noise and disturbance follow the Gaussian distribution, its performance might degrade if the assumption is wrong. To address this problem, we propose a fault detection and exclusion algorithm using Gaussian sum filters. Because GNSS measurement noise does not follow the Gaussian distribution perfectly, the Gaussian sum filter can estimate the posterior distribution more accurately; therefore it has better integrity-monitoring performance. This paper describes the detailed algorithms and shows simulation results to evaluate the integrity-monitoring performance of the algorithms. The proposed algorithms detect about 30% smaller faults and generate 35% lower protection levels than the conventional methods. The results show that the proposed algorithms can provide better accuracy and availability performance.

Keywords

Integrity monitoring RAIM Gaussian sum filter Non-Gaussian error distribution

Type: Research Article
Information: The Journal of Navigation , Volume 61 , Issue 4 , October 2008 , pp. 687 - 703

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

1. INTRODUCTION

For safety-critical applications of global navigation satellite systems (GNSSs), such as aircraft and missile navigation systems, it is important to be able to detect and exclude faults that could adversely affect accuracy and integrity so that the navigation system can operate continuously without any degradation in performance. For high-accuracy systems, the fault detection and exclusion (FDE) function needs to be able to detect and exclude smaller and multiple biases. Because it is difficult to detect small errors and simultaneous multiple faults using the conventional snapshot receiver autonomous integrity-monitoring (RAIM) algorithms, various filtering methods have been studied to reduce the measurement noise level or to integrate a GNSS with other sensors, so that the navigation system can estimate its position more accurately and reliably. Because filters reduce the noise level of measurements using previous information, they can provide a better integrity-monitoring performance than snapshot algorithms. Moreover, a filtering method can calculate measurement residuals individually, so that it can detect simultaneous multiple faults and exclude faulty measurements more easily than conventional RAIM algorithm. However, because most filters, such as Kalman filters, presume that the measurement error and disturbance follow a Gaussian distribution, their performance can degrade if this assumption is wrong. Because the GNSS measurement error does not follow a Gaussian distribution perfectly, the Kalman filter approach has to use an inaccurate error model that may cause performance degradation.

To address this problem, we propose a new integrity-monitoring algorithm using Gaussian sum filters. It is known that Gaussian sum filters and particle filters have an almost identical GNSS integrity-monitoring performance (Yun et al., Reference Yun, Kim and Kee2006). However, particle filters need high computational power and are therefore difficult to implement in real time systems. A Gaussian sum filter is proposed as an alternative to the particle filter method. Gaussian sum filters have been researched over the last few decades as an alternative for solving nonlinear/non-Gaussian problems (Daniel and Sorenson, Reference Alspach and Sorenson1972; Fruhwirth, Reference Fruhwirth1997). The proposed algorithm estimates the distribution of measurement residuals individually from the posterior density and detects exceedingly large residuals to satisfy a false alarm rate. This algorithm can detect faults based on accurate estimation of the individual posterior distribution and can detect and exclude faults almost simultaneously so that the system can exclude a faulty measurement easily. With a non-Gaussian measurement error, a Gaussian sum filter can estimate the distribution of the state more accurately than a Kalman filter and therefore has better integrity-monitoring performance. In addition, if the system is nonlinear, the performance will also be better. Gaussian sum filters for nonlinear systems have undergone much research (Hanebeck et al., Reference Hanebeck, Briechle and Rauh2003), but the advantage of these filters with a non-Gaussian error distribution has not been investigated. Therefore, our work focused on the effect of a non-Gaussian error distribution of the DGNSS measurement on the integrity-monitoring performance, which was assessed by simulation. The proposed algorithm was implemented in a WADGPS system with three conventional integrity-monitoring algorithms and its performance was evaluated by comparison with other algorithms. We implemented four integrity-monitoring algorithms in this work: the snapshot χ² method (SSX2), the Kalman filter χ² method (KFX2), the Kalman filter individual residual (KFIR) method, and the Gaussian sum filter individual residual (GSFIR) method. Our proposed algorithm for the integrity monitor function based on the Gaussian sum filter is described in sections 4 and 5 and the other methods will be treated briefly.

2. DGNSS PSEUDORANGE ERROR MODEL

This section and the next describe the true error model and the sigma overbounding method, which are based on the work of Shively and Braff (Reference Shively and Braff2000).

In most GNSS applications, including the Space-Based Augmentation System (SBAS) and the Ground-Based Augmentation System (GBAS), the pseudorange measurement error is assumed to follow a Gaussian distribution. However, this is not true. In general, the core part of the error distribution can be characterized well using a Gaussian distribution, but the tail part of the distribution is heavier than that of a Gaussian distribution. The heavier tail is due to ground-reflected multipaths or to systematic receiver/antenna errors (Sayim et al., Reference Sayim, Pervan, Pullen and Enge2002). This type of distribution can be represented by a Gaussian core-Laplacian tail (GCLT) probability density function (PDF). In this paper, we take the GCLT model as the receiver's true error distribution to simulate a heavy-tailed model:

(1)

${\rm p}\lpar n\rpar \equals \openup3\left\{\! {\matrix{ {{1 \over {\sqrt {2\pi \sigma } }}\exp \left( \displaystyle{{{n^{\setnum{2}} } \over {2\sigma ^{\setnum{2}} }}} \right)\comma } \tab {\left\vert n \right\vert \les n_{tr} } \cr {p_{tr} \times \exp \left( { \minus \displaystyle{{p_{tr} } \over {c_{tr} }}\left\vert {n \minus n_{tr} } \right\vert} \right)\comma } \tab {\left\vert n \right\vert \gt n_{tr} } \cr} } \right.$

where n _tr=2·58σ is the transition point from the Gaussian to the Laplacian distribution. The other constants are defined as follows:

(2)

$p_{tr} \equals {\rm p}_{N\lpar \setnum{0}\comma \sigma \rpar } \lpar n_{tr} \rpar \times PUF\comma \quad c_{tr} \equals \int\nolimits_{n_{{tr}} }^{\infty } {{\rm p}_{{N\rm \lpar \setnum{0}\comma }\sigma {\rm \rpar }} } \lpar n\rpar dn \times TUF$

PUF and TUF provide confidence in consideration of the volume of data. The narrower curve shown in Figure 1 represents the PDF.

Figure 1. The Gaussian–Laplacian model and the Gaussian overbound model.

3. SIGMA OVERBOUNDING

Because the measurement error does not follow a Gaussian distribution, sigma overbounding methods are used to simplify implementation and to prevent integrity risks (Shively and Braff, Reference Shively and Braff2000; Rife et al., Reference Rife, Pullen and Pervan2004). Overbounding is a procedure used to find the inflation factors, which are the ratio of the overbounded and the observed sigmas.

(3)

$\sigma _{over} \equals INF \cdot \sigma$

The inflation factor is determined using equation (4):

(4)

$INF \equals {M \over {5.95}}\comma \int\nolimits_{M}^{\infty } {{\rm p}_{true} \lpar x\rpar dx \equals } \int\nolimits_{\setnum{5}.\setnum{95}}^{\infty } {{\rm p}_{N\lpar \setnum{0}\comma \setnum{1}\rpar } \lpar x\rpar dx \equals 3.12 \times 10^{ \minus \setnum{9}} }$

where p_true (x) is the observed true PDF, which was taken to represent the GCLT distribution in our work. The denominator value of 5·95 is a multiplier used for SBAS users to calculate a defined fault-free protection level (RTCA, 2001). This method provides users with an overbound broadcast sigma that they can use to estimate the fault-free protection levels to bound position errors with the required probability. Figure 2 shows a pictorial representation of this method, which shows that the one-sided tail probabilities of the Gaussian and the GCLT distributions have the same observed standard deviation. Figure 1 shows the PDF values of the true error model (GCLT) and the overbounded Gaussian (OG) model. The OG model had a broader distribution, so that the tail was conservative enough to bound the heavy-tailed GCLT model. This method can satisfy integrity performance, however its residual distribution is so conservative that it calculates a protection level (PL) that is too large. The system compares the value of the PL to a predefined alert limit (AL) for an expected operation, for example, precision approach, non-precision approach and enroute, to decide whether or not the navigation system will use the solution. If PL is larger than AL, the system will not use the solution for reasons of safety, therefore the availability of the system can be degraded. To address this problem we propose a new integrity-monitoring algorithm using the Gaussian sum filter.

Figure 2. Set inflation factor.

4. GAUSSIAN SUM FILTER

We can approximately express an arbitrary PDF with a weighted sum of Gaussian distributions, which is called a Gaussian mixture model (GMM). When we assume that the measurement error follows the model, we can use several parallel Kalman filters that deal with each Gaussian component to estimate the distribution of states more accurately. This method is called a Gaussian sum filter (Fruhwirth, Reference Fruhwirth1997).

First, let us assume that a j-th measurement error distribution can be expressed by a sum of N _comp Gaussian distributions:

(5)

${\rm p}_{GM} \lpar n^{\hskip1 j} \rpar \equals \sum\limits_{i \equals \setnum{1}}^{N_{{{\rm comp}}} } {\gamma ^{ \lt\hskip-1 i\hskip-1 \gt } {\rm p}_{\rm N} \lpar n^{\hskip1j} \semi {\rm \ }m^{\hskip1j\comma \lt\hskip-1 i\hskip-1 \gt } \comma v^{\hskip1j\comma \lt\hskip-1 i\hskip-1 \gt } \rpar}$

where $\sum\limits_{i \equals \setnum{1}}^{N_{{comp}} } {\gamma ^{ \lt\hskip-1 i\hskip-1 \gt } \equals 1 \mathop{\rm p}\nolimits} _{\rm N} \lpar x\semi {\rm \ }m\comma C\rpar$ : Gaussian probability density of x, m is mean, C is variance; $\left(\ {} \right)^{ \lt\hskip-1 i\hskip-1 \gt}$ : i-th component of Gaussian mixture model.

If we assume that all N _M measurements are independent, a distribution of the measurement vector is obtained by the following equation:

(6)

${\rm p}\lpar {\bf n}\rpar \equals \prod\limits_{j \equals \setnum{1}}^{N_{M} } {{\rm p}_{GM} \lpar n^{\hskip1 j} \rpar \equals \prod\limits_{j \equals \setnum{1}}^{N_{M} } {\left( {\sum\limits_{i \equals \setnum{1}}^{N_{{{\rm comp}}} } {\gamma ^{ \lt\hskip-1 i\hskip-1 \gt } {\rm p}_{\rm N} \lpar n^{\hskip1j} \semi {\rm \ }m^{\hskip1j\comma \lt\hskip-1 i\hskip-1 \gt } \comma v^{\hskip1j\comma \lt\hskip-1 i\hskip-1 \gt } \rpar } } \right)}}$

and it can be expressed by a sum of multivariate normal distributions:

(7)

${\rm p}\lpar {\bf n}_{k} \rpar \equals \sum\limits_{i \equals \setnum{1}}^{N_{{mode}} } {\varpi _{k}^{ \lt\hskip-1 i\hskip-1 \gt } {\rm p_{N}} \lpar {\bf n}_{k} \semi {\bf M}_{k}^{{\rm \lt }i{\rm \gt }} \comma {\bf R}_{k}^{{\rm \lt }i{\rm \gt }} \rpar } \comma$

where N _mode is the number of mode, $N_{\rm mode} \equals \lpar N_{{\rm comp}} \rpar ^{N_{M} } \comma \quad \sum\limits_{i \equals \setnum{1}}^{N_{{{\rm mode}}} } {\varpi _{k}^{ \lt\hskip-1 i\hskip-1 \gt } } \equals 1$ .

With this measurement error distribution, we implement the Kalman filters for every Gaussian component. First, after time propagation using equation (8), the GSF has the prior distribution as equation (9):

(8)

$\tilde{\bf x}_{k}^{ \lt j \gt } \equals {\bf F}_{k} {\rm \hat{\bf x}}_{k \minus \setnum{1}}^{ \lt j \gt } \comma \quad \tilde{\bf P}_{k}^{ \lt j \gt } \equals {\bf F}_{k} \hat{\bf P}_{k \minus \setnum{1}}^{ \lt j \gt } {\bf F}_{k}^{T} \plus {\bf Q}_{k \minus {\rm \setnum{1}}} \comma$

(9)

${\rm p\lpar }{\bf x}_{k} \left\vert {{\bf z}_{{\rm \setnum{1}\colon }k \minus \setnum{1}} } \right.{\rm \rpar \equals }\sum\limits_{j{\rm \equals \setnum{1}}}^{N_{{{\rm post}}} } {\pi _{k}^{ \lt j \gt } {\rm p}_{\rm N} \lpar {\bf x}_{k} \semi {\rm \ }\tilde{\bf x}_{k}^{ \lt j \gt } \comma \tilde{\bf P}_{k}^{ \lt j \gt } \rpar } \comma$

where:

${\rm \tilde{x}}_{\rm k}$ : predicted state vector at time k
${\rm \hat{x}}_{\rm k}$ : measurement updated state vector at time k
${\rm z}_{{\rm \setnum{1}\colon k}}$ : measurement set up to time k
F_k: state transition model which is applied to the previous state x_k−1.
P: estimation error covariance matrix
Q: process noise covariance matrix

$\hskip11\sum\limits_{j \equals \setnum{1}}^{N_{{{\rm post}}} } {\pi _{k}^{ \lt\hskip1 j\hskip1 \gt } \equals 1}\hfill$

From equation (7), the measurement distribution at the current epoch is:

(10)

${\rm p\lpar }{\bf z}_{k} \left\vert {{\bf x}_{k} } \right.{\rm \rpar \equals }\sum\limits_{i \equals \setnum{1}}^{N_{{{\rm mode}}} } {\varpi _{k}^{ \lt\hskip1 i\hskip1 \gt } {\rm p}_{\rm N} \lpar z_{k} \semi {\rm \ }{\bf H}_{k} {\bf x}_{k} \comma {\bf R}_{k}^{ \lt\hskip-1 i\hskip-1 \gt } \rpar } \comma$

where H is the observation matrix, R is the observation noise covariance matrix.

GSF updates the weights of the Gaussian components, the states and the covariance matrix by using equations (11–13):

(11)

$\omega _{K}^{ \lt i\comma j \gt } \propto \varpi _{k}^{ \lt\hskip-1 i\hskip-1 \gt } \pi _{k}^{ \lt\hskip-1 i\hskip-1 \gt } {\rm p}\lpar {\bf z}_{k} \semi {\rm \ }{{\bf H}_{k} {\tilde{\bf x}}_{_{k} }^{ \lt j \gt } \comma {\bf R}_{_{k} }^{ \lt\hskip-1 i\hskip-1 \gt } \plus {\bf H}_{k} {\tilde{\bf P}}_{_{k} }^{ \lt j \gt } {\bf H}_{_{k} }^{T} \rpar \comma$

where $\sum\limits_{i \equals \setnum{1}}^{N_{{\bmod {\rm e}}} } {\sum\limits_{j \equals \setnum{1}}^{N_{{{\rm post}}} } {\omega _{k}^{ \lt i\comma j \gt } \equals 1}}$ ,

(12)

${\hat{\bf x}}_{k}^{ \lt i\comma j \gt } \equals {\tilde{\bf x}}_{k}^{ \lt j \gt } \plus {\hat{\bf P}}_{k}^{ \lt i\comma j \gt } {\rm H}_{k}^{T} \lpar {\bf R}_{k}^{ \lt\hskip-1 i\hskip-1 \gt } \rpar ^{ \minus \setnum{1}} \lpar {\bf z}_{k} \minus {\bf H}_{k} {\tilde{\bf x}}_{k}^{ \lt j \gt } \rpar \comma$

(13)

${\hat{\bf p}}_{\rm k}^{{\rm \lt }i\comma j{\rm \gt }} \equals \lsqb \lpar {\tilde{\bf p}}_{k}^{ \lt j \gt } \rpar ^{ \minus \setnum{1}} \plus {\bf H}_{k}^{T} \lpar {\bf R}_{k}^{ \lt\hskip-1 i\hskip-1 \gt } \rpar ^{ \minus \setnum{1}} {\bf H}_{k} \rsqb ^{ \minus \setnum{1}}.$

Now we have the estimation of the current states as follows.

(14)

${\rm p\lpar }{\bf x}_{k} {\rm \vert }{\bf z}_{\setnum{1}\colon k} {\rm \rpar } \equals \sum\limits_{i{\rm \equals }\setnum{1}}^{N_{{{\rm mode}}} } {\sum\limits_{j \equals \setnum{1}}^{N_{{{\rm post}}} } {\omega _{k}^{ \lt\hskip-1 i\comma\hskip-1 j\hskip-1 \gt } {\rm p}_{\rm N} {\rm \lpar }{\bf x}_{k} {\rm \semi }\,{\hat{\bf x}}_{k}^{ \lt\hskip-1 i\comma\hskip-1 j\hskip-1 \gt } {\rm \comma }{\hat{\bf p}}_{k}^{ \lt\hskip-1 i\comma\hskip-1 j \hskip-1\gt } {\rm \rpar}}}$

5. INTEGRITY-MONITORING ALGORITHMS

GNSS-based navigation systems must have an integrity-monitoring system that contains two main functions: 1) detection and exclusion of satellite faults and 2) estimation of the uncertainty of the position solutions. First, the system calculates decision variables and compares these to thresholds that have been set to satisfy the integrity requirements for the desired operation. If any decision variable exceeds the threshold, the system concludes that the corresponding measurement has a fault, and excludes the detected fault. Second, because the system cannot know the true position estimation errors, it estimates the uncertainty of the solution, which is known as the protection level.

5.1. Fault Detection and Exclusion

The proposed algorithm is a residual method that uses the posterior density of a state, which is estimated from the filters and the measurement error distribution. A residual is defined as the difference between an expected and an actual measurement, which will be described in detail in the following subsections. Residual methods using the snapshot algorithm (SSX2) and the Kalman filter (KFX2) calculate the decision variables using the squared sum of the residuals, which is known as the weighted squared sum of the error (WSSE) that follows a χ² distribution whose degrees of freedom are determined by the number of measurements, as in equation (19).

(15)

$WSSE \equals {{{\bf r}^{T} {\bf r}} \over {{\bf C}_{r} }}\sim \left\{ {\matrix{ {\chi ^{\setnum{2}} \lpar 0\comma N_{M} \minus 4\rpar \comma } \tab {{\rm SSX2}} \cr \hskip-14{\chi ^{\setnum{2}} \lpar 0\comma N_{M} \rpar \comma } \tab {{\rm KFX2}} \cr} } \right.$

This type of algorithm can be found in literature (Parkinson and Axelrad, Reference Parkinson and Axelrad1988; McGraw et al., Reference McGraw, Murphy, Brenner, Pullen and Van Dierendonck2000).

The filtering approaches can examine the individual residuals, but not the combined one, because they have a previously estimated position. In the KFIR, the algorithm compares a measurement residual to the multiple of the residual standard deviation. Because our proposed method, GSFIR, assumes a non-Gaussian error, it does not use the residual standard deviation. Instead, it estimates the residual distribution, and then decides whether the residual is located in a normal region or not. These individual residual algorithms have to estimate the residual distribution and determine the individual detection threshold for each residual.

The expected measurement was calculated from the propagation of the previous epoch's state vector as:

(16)

${\tilde{\bf x}}_{k} \equals E\lsqb {\rm f}_{k} \lpar {\hat{\bf x}}_{k \minus \setnum{1}} \comma {\bf v}_{k \minus \setnum{1}} \rpar \rsqb \equals {\rm f}_{k} \lpar {\hat{\bf x}}_{k \minus \setnum{1}} \comma 0\rpar.$

The expected measurement is a projection of the propagated state from the position domain to the measurement domain. Under normal conditions, the expected and the actual measurement distributions lie close to each other but when there is a measurement fault they are separated by a distance. In our work the distance between the two measurements is referred to as the residual.

(17)

${\bf r}_{k} \equals {\rm h}_{k} \lpar {\tilde{\bf x}}_{k} \comma 0\rpar \minus {\bf z}_{k}$

The Gaussian sum filter has the residual components:

(18)

${\bf r}_{k}^{ \lt i\comma j \gt } \equals {\bf H}_{k} {\tilde{\bf x}}_{k}^{ \lt\hskip-1 i\hskip-1 \gt } \minus {\bf z}_{k}^{ \lt j \gt } \comma$

the associated weights:

(19)

$\omega _{k}^{ \lt i\comma j \gt } \equals \omega _{k}^{ \lt\hskip-1 i\hskip-1 \gt } \varpi _{k}^{ \lt j \gt } \comma$

and the covariances:

(20)

${\bf C}_{r\comma k}^{ \lt i\comma j \gt } \equals {\bf H}_{k} {\tilde{\bf P}}_{k}^{ \lt\hskip-1 i\hskip-1 \gt } {\bf H}_{k}^{T} \plus {\bf R}_{k}^{ \lt j \gt }.$

Therefore, the residual distribution from the Gaussian sum filter is expressed as follows:

(21)

${\rm p}_{{\bf R}_{k} } \lpar {\bf r}_{k} \rpar \equals \sum\limits_{l \equals \setnum{1}}^{N_{{post}} } {\omega _{k}^{ \lt l \gt } {\rm p}_{\rm N} \lpar {\bf r}_{k} \semi \,{\bf r}_{k}^{ \lt l \gt } \comma {\bf C}_{r\comma k}^{ \lt l \gt } \rpar}$

A sample of a residual distribution under nominal conditions is shown in Figure 3. The expectation is near zero, and the large residuals have low probabilities. When the residual is in the region between the two thresholds, known as the normal region, the system declares that there is no fault in the measurement. Otherwise, it determines the measurement as having a fault, and excludes it for continuous operation.

Figure 3. PDF of a residual and thresholds.

In real implementation, the expectation of the PDF is located at the calculated residual and not at zero. Therefore, we calculate a one-sided tail probability below or above zero as a decision variable using equation (22) and compare it with the individual false alarm rate that will be defined in the next subsection for each measurement.

(22)

$pr_{k}^{j} \equals \left\{ {\matrix{ {\int\nolimits_{ \minus \infty }^{\setnum{0}} {{\rm p}_{R_{k}^{\hskip1j} } \lpar r_{\hskip-1k}^{{\hskip1j}}\hskip-2 \prime \rpar dr_{\hskip-1k}^{{\hskip1j}}\hskip-2 \prime \comma } } \tab {r_{\hskip-1k}^{\hskip1j} \ges 0} \cr {\int\nolimits_{\setnum{0}}^{\infty } {{\rm p}_{R_{k}^{\hskip1j} } \lpar r_{\hskip-1k}^{{\hskip1j}}\hskip-2 \prime \rpar dr_{\hskip-1k}^{{\hskip1j}}\hskip-2 \prime \comma } } \tab {r_{\hskip-1k}^{\hskip1j} \lt 0} \cr} } \right.$

Because the KFIR method assumes a Gaussian error, the PDFs are all Gaussian and the residual distribution can be expressed simply by a standard deviation, so that the thresholds are determined as the multiples of the standard deviations. The threshold is set to satisfy a specified false alarm rate when a χ² method based on the Gaussian distribution, SSX2 or KFX2, is employed.

(23)

$\int\nolimits_{T_{{\chi ^{\setnum{2}} }} }^{\infty } {{\rm p}_{\chi ^{\setnum{2}} } \lpar x\rpar dx \equals P_{FA} }$

However, the individual residual test requires a threshold for each measurement and therefore an individual false alarm rate, p _fa, needs to be determined first. We used the multivariate normal distribution and the residual covariance matrix estimated from the Kalman filter, and assuming that the individual false alarm rates for all residuals are the same, we have the following equation:

(24)

$P_{FA} \equals 2\left( {1 \minus {1 \over {\sqrt {\vert {\rm C}_{r} \vert \lpar 2\pi \rpar ^{N_{M} } } }}\int\nolimits_{ \minus \infty }^{K\sigma _{\setnum{1}} } {\int\nolimits_{ \minus \infty }^{K\sigma _{\setnum{2}} } { \ldots \int\nolimits_{ \minus \infty }^{K\sigma _{{N_{M} }} } {\exp \left(\!\! \minus {1 \over 2}{\bf r}^{T} {\rm C}_{r} {\bf r}\right) d{\bf r}} } } } \right)$

where $\sigma _{i} \equals {\rm C}_{r} \lpar i\comma i\rpar$ . Because K is the only unknown in this equation, we can calculate it by iteration. With the determined K, the individual false alarm rate is calculated as follows:

(25)

$P_{fa\comma i} \equals 2\int\nolimits_{k\sigma _{i} }^{\infty }\!\! {P_{R_{i} } \lpar r_{i} \rpar dr_{i} \equals } 2\int\nolimits_{k}^{\infty } \!\!{P_{N\lpar \setnum{0}\comma \setnum{1}\rpar } \lpar n\rpar dn \equals P_{fa} }$

Because in Section 2 we assumed that the measurement error followed a symmetrical distribution, the proposed algorithms declare a fault present if a one-sided tail probability of a residual $pr_{k}^{j}$ exceeded the value of P _fa/2.

5.2. Protection Level

SBAS users must calculate both fault-free and faulted protection levels (XPL ₀ and XPL ₁, respectively). The faulted protection levels indicate the position error protection assuming that a fault measurement exists. SBAS assumes that a fault is induced from the reference receiver errors. However, we assumed that the fault includes all possible errors contained in the measurements that are used for position estimation. The formulas used to calculate the protection levels based on a snapshot method assuming a Gaussian error were derived as follows. (The Kalman filter and the proposed filtering methods can obtain the protection levels using a similar procedure, which will be discussed later).

In GNSSs, the measurement equation can be expressed using the linear function:

(26)

${\bf z} \equals {\bf Hx} \plus {\bf n}.$

The state vector is estimated using the weighted least squares method and the estimation error can be obtained from equation (27):

(27)

$\delta \hat{\bf x} \equals {\bf H}_{W}^{ \plus } \delta {\bf z}\comma \quad {\bf H}_{W}^{ \plus } \equals \lpar {\bf H}^{\rm T} {\bf WH}\rpar ^{ \minus \setnum{1}} {\bf WH}^{\rm T} \comma$

where ${\bf W} \equals \lpar diag\vert \sigma _{n_{\setnum{1}}}^{\setnum{2}} \ldots \sigma _{n_{j}}^{\setnum{2}} \ldots \sigma _{n_{{N_{M} }} }^{\setnum{2}} \vert \rpar ^{ \minus \setnum{1}}.$

When a fault T occurs in measurement j, the protection levels are calculated using the following equations:

(28)

$\eqalign{{\rm HPL}_{\rm \setnum{0}} \tab \equals k_{ffmd\comma h} \sigma _{h} \equals HPL_{\setnum{0}\comma n} {\rm \comma }\,{\rm HPL}_{\rm \setnum{1}}^{\rm j} \equals k_{md} \sigma _{h} \plus \vert T_{j} \sqrt {{\bf H}_{w}^{ \plus } \lpar 1\comma j\rpar ^{\setnum{2}} \plus {\bf H}_{w}^{ \plus } \lpar 2\comma j\rpar ^{\setnum{2}} } \vert\cr\tab \equals HPL_{\setnum{1}\comma n} \plus HPL_{b}^{j} \comma$

(29)

${\rm VPL}_{\rm \setnum{0}} \equals k_{ffmd\comma v} \sigma _{v} \equals VPL_{\setnum{0}\comma n} {\rm \comma }\,{\rm VPL}_{\rm \setnum{1}}^{\rm j} \equals k_{md} \sigma _{v} \plus \vert T_{j} \cdot {\bf H}_{w}^{ \plus } \lpar 3\comma j\rpar \vert \equals VPL_{\setnum{1}\comma n} \plus VPL_{b}^{j} \comma$

where $\sigma _{h}^{\setnum{2}} \equals \lpar {\bf H}^{{\rm T}} {\bf WH}\rpar _{\setnum{1}\comma \setnum{1}}^{ \minus \setnum{1}} \plus \lpar {\bf H}^{\rm T} {\bf WH}\rpar _{\setnum{2}\comma \setnum{2}}^{ \minus \setnum{1}}, \sigma _{v}^{\setnum{2}} \equals \lpar {\bf H}^{\rm T} {\bf WH}\rpar _{\setnum{3}\comma \setnum{3}}^{ \minus \setnum{1}}, k_{ffmd\comma h} \equals 6.0, k_{ffmd \comma v} \equals 5.33$ and k _md=3.29

The standard deviation multipliers, k _ffmd,h, k _ffmd,ν, and k _md, are defined in WAAS MOPS (RTCA, 2001). The value of the fault T is equal to the threshold for each FDE algorithm because it is the maximum fault that the system cannot detect. The filtering methods can estimate the protection levels with the above procedure, however, the values of XPL _n and XPL _b^j are determined differently for each algorithm. The Kalman filter method can obtain the standard deviations, σ_h and σ_v, from the estimated covariance matrix and T _j is calculated from the residual standard deviation and the individual false alarm rate. The Gaussian sum filter methods can estimate them from the estimated posterior density and the residual density, respectively as follows:

(30)

$\int\nolimits_{ \minus \infty }^{k_{{md}} } {{\rm p}_{Gaussian} \lpar n\rpar dn \equals \int\nolimits_{ \minus \infty }^{VPL_{n} } {{\rm p}_{{\rm \hat{x}}_{{\rm k}} \lpar \setnum{3}\rpar } \lpar x_{\setnum{3}} \rpar dx_{\setnum{3}} } }$

(31)

$P_{fa} \equals 2\int\nolimits_{T_{j} }^{\infty } {{\rm p}_{R_{k}^{j} } \lpar r_{k}^{j} \rpar dr_{k}^{j} }$

6. SIMULATION

In this section, the proposed algorithm is verified by performing a real time test. For the real time test, the proposed algorithm was implemented in the user system of Korean WADGPS Test Bed (KWTB) (Kim et al., Reference Kim, Yun, Park, Jeon, Sohn and Kee2006).

6.1. Simulation Setup

A Trimble 4000ssi receiver with a quartz crystal oscillator (Trimble, 1996) was used for the user system. The state variable was as follows:

(32)

${\bf x} \equals \lsqb e\ n\ u\ B\ \rmDelta B\rsqb ^{T}$

The propagation matrix is defined as in equation (33).

(33)

${\bf F} \equals \left[ {\matrix{{ {\bf I}{\rm \lpar 3} \times {\rm 3\rpar }} \tab {\bf 0} \cr {\bf 0} \tab {{\bf F}_{\rm C} } \cr} } \right]\comma \quad {\bf F}_{\rm C} \left[ {\matrix{ 1 \tab {\rmDelta t} \cr 0 \tab 1 \cr} } \right]$

Δt is the sampling time (1 s).

In the process noise covariance matrix Q, the elements associated with user position are set to be a very small value, close to 0, because the user does not move in this test.

(34)

${\bf Q} \equals {\rm E\lsqb vv}^{\rm T} \rsqb \equals \left[ {\matrix{ {0.001 \times {\bf I}_{\setnum{3} \times \setnum{3}} } \tab {\bf 0} \cr {\bf 0} \tab {{\bf Q}_{C} } \cr} } \right]$

The elements associated with clock error are represented as in equation (35) (Brown and Hwang, Reference Brown and Hwang1997).

(35)

${\bf Q}_{C} \equals \left[ {\matrix{ {S_{f} \rmDelta t \plus {{S_{g} \rmDelta t^{\setnum{3}} } \over 3}} \tab {S_{g} \rmDelta t} \cr {S_{g} \rmDelta t} \tab {{{S_{g} \rmDelta t^{\setnum{2}} } \over 2}} \cr} } \right]$

Because the element describing the clock drift was added in the state vector, the measurement matrix H was also edited as follows:

(37)

${\bf H} \equals \left( {\matrix{ {{\bf e}_{u}^{\setnum{1}^{T} } } \hfill \tab { \minus 1} \hfill \tab 0 \hfill \cr {{\bf e}_{u}^{\setnum{2}^{T} } } \hfill \tab { \minus 1} \hfill \tab 0 \hfill \cr \hfill\vdots\hfill \hfill \tab \hfill\vdots \hfill \tab 0 \hfill \cr {{\bf e}_{u}^{N_{M}^{\ \ \ T} } } \hfill \tab { \minus 1} \hfill \tab 0 \hfill \cr} } \right).$

6.2. Analysis of Pseudorange Error and Error Modelling

In order to find the pseudorange error distribution of the WAD correction data, we saved the pseudorange error after applying the WAD correction. The error level is set to be 3·09σ, because 3·09σ includes 99·9% of the random variable in Gaussian distribution. According to SBAS requirement, all of these correction errors must be smaller than this error level.

The histogram in Figure 4 (Top left) shows the probability density of the correction error collected in 3 wide-area reference stations (WRS). As expected, the correction error distribution is similar to the Gaussian distribution, but it has a heavier tail. The presented Gaussian distribution has a standard deviation that has been calculated from correction data that is UDRE and GIVE.

Figure 4. KWTB correction error. Top left – Probability density function; Top right – QQ-plot; Bottom – Sigma overbound.

In order to perform integrity monitoring with the filtering algorithm, an error model needs to be set up. Figure 4 (Top right) is a Quantile–Quantile (QQ) plot of the collected correction error distribution. This figure shows that the distribution becomes more Gaussian-like as it gets closer to a straight line, and a larger slope means a larger standard deviation. The collected data distribution is Gaussian-like between −2 and 2. However, outside of that interval it is different from Gaussian distribution. Thus the error model has to be expressed as a Gaussian mixture model. It is hard to determine the precise tail distribution because of the limited number of data. Therefore, the tail distribution of the error model is set to follow the tail distribution of the GCLT model, so that the error model follows data distribution in the core, and the GCLT model in the tail. In Figure 4 (Top right), the Gaussian mixture model follows the collected data in the interval of −1·7~1·7, and the GCLT model outside the interval, and bounds collected data in whole interval.

The simulation sets PUF=0·52 and TUF=1 in equation (2) to follow the convoluted three-receiver error model. The standard deviation of overbounded Gaussian distribution was calculated as in equations (3) and (4). The calculated inflation factor was 2·2, and the chosen coefficients of the Gaussian mixture model are in equation (38). As this model is applied, Figure 4 (Bottom) shows the GCLT model, the Gaussian mixture model and the overbounded Gaussian model using the inflation factor.

(38)

$\matrix{ {\gamma ^{ \lt\hskip-1 \setnum{1}\hskip-1 \gt } \equals 0{\cdot}81\comma } \tab {m^{ \lt\hskip-1 \setnum{1}\hskip-1 \gt } \equals 0{\cdot}0\comma } \tab {v^{ \lt\hskip-1 \setnum{1}\hskip-1 \gt } \equals 0{\cdot}4^{\setnum{2}} } \cr {\gamma ^{ \lt\hskip-1 \setnum{2}\hskip-1 \gt } \equals 0{\cdot}19\comma } \tab {m^{ \lt\hskip-1 \setnum{2}\hskip-1 \gt } \equals 0{\cdot}0\comma } \tab {v^{ \lt\hskip-1 \setnum{2}\hskip-1 \gt } \equals 1{\cdot}8^{\setnum{2}} } \cr \hskip10{\gamma ^{ \lt\hskip-1 \setnum{3}\hskip-1 \gt } \equals 0{\cdot}0008\comma } \tab {m^{ \lt\hskip-1 \setnum{3}\hskip-1 \gt } \equals 0{\cdot}0\comma } \tab {v^{ \lt\hskip-1 \setnum{3}\hskip-1 \gt } \equals 2{\cdot}8^{\setnum{2}} } \cr}$

6.3. Real-time Test

A simulation was performed to assess our proposed integrity monitoring algorithms based on the Gaussian sum filter. To investigate the advantage of the non-Gaussian assumption, three conventional methods using the Gaussian assumption were also implemented in parallel. These methods are shown in Table 1.

Table 1. Four methods implemented in simulation.

To evaluate the integrity-monitoring performance we injected bias into the pseudorange, as in Figure 5. The first four biases were injected as step shape in order to evaluate the integrity performance to a sudden system fault. After that, decreasing bias was injected to find what bias size each algorithm could detect. Four increasing biases with different slopes were injected, so that we could find out how it copes with diverging error. This fault scenario takes about 1000 s to complete.

Figure 5. Fault scenario.

6.4. Measurement Domain Results

In this section we assess the suitability of each method for each bias type. The decision variable is calculated and compared with the threshold to determine whether there is a fault or not. If it is declared that a fault exists, the number of fault measurements needs to be verified. If the method finds that there are two or more fault measurements, it has failed to detect and exclude the fault, because we injected bias into only one measurement. Lastly, position error and protection level are verified, so that we can estimate the user performance.

Fault scenario 1 – step bias

SSX2 and KFX2 succeed in detecting the first two biases, at 20 and 15 m, but cannot detect a bias of 5 and 10 m at all. In the case of the 15 m bias, SSX2 could detect the fault but a false alarm occurred because it determined that two measurements were faulty. KFIR detected and excluded the 20 and 15 m biases successfully, but detected the 10 m bias only in 2 epochs. GSFIR detected and excluded 20, 15, and 10 m biases well. No method could detect the 5 m bias.

Fault scenario 2 – decreasing bias

As shown in Figure 6, SSX2 detected faults in the first few epochs, and KFX2 detected faults until 15 s. After excluding a faulty measurement, if the number of measurements that remained normal is less than four, the position cannot be calculated and if less than five, the decision variable cannot be calculated. The number of visible satellites during this test is six or seven and we injected bias into one measurement, therefore the parts for which the decision variable equals to zero are those in which a normal measurement is declared as faulty. KFIR succeeded in detecting fault until 17 s, failed to detect in one epoch and detected again until 22 s. GSFIR detected until 26 s, failed to detect in a few epochs and detected again until 38 s. Through this result the minimum detectable bias of each method can be calculated, as shown in Table 2. As expected GSFIR detected the smallest bias.

Figure 6. Decision variables and threshold: Top left – KFX2; Top right – SSX2; Bottom – KFIR and GSFIR.

Table 2. Minimum detectable bias for each method.

Fault scenario 3 – increasing bias

Each method begins to detect the ramp failure at the time indicated in Table 3. Also, in this fault scenario GSFIR can detect the smallest bias and can detect fault most quickly.

Table 3. The time required to detect ramp failure.

6.5. Position Domain Results

Figure 7 shows the horizontal position error and protection level after exclusion. The GSFIR method can detect and exclude the highest number of faults therefore it has the smallest error level. The GSFIR method provides the smallest protection level, as expected. Table 4 shows the numerical values of this result. The horizontal and vertical protection levels of GSFIR are, respectively, 30% and 35% smaller than those of KFIR.

Figure 7. Horizontal position error and protection level after exclusion.

Table 4. Position error and protection level after exclusion.

Figure 8 shows the fault detection threshold for the residuals of KFIR and GSFIR that affect directly the performance of fault detection and the decision of protection level. As shown in Table 5, the threshold of the GSFIR method is about 3 m smaller than that of KFIR. Therefore, GSFIR can cope with a smaller bias and also can provide a smaller protection level.

Figure 8. Fault detection threshold for residual.

Table 5. Average of fault detection threshold for residual.

Under the assumption that the system is used in aircraft navigation, the availability was calculated as applied in APV-II (Approaches with vertical guidance II) and the CAT-I (Category I) precision approach. As seen in Table 6, for APV-II, all filtering methods have a very high availability; in this case there is little difference between KFIR and GSFIR. However, when applied to the stricter requirements of CAT-I, GSFIR was superior to other methods. The availability of KFIR is 50% higher than that of KFX2, and that of GSFIR is 35% higher than that of KFIR. It can be said that although the difference in error distribution is very small, it is very critical at strict requirement.

Table 6. Availability of each method.

6.6. Result Analysis

The reason that the GSFIR method outperformed other methods can be explained as follows.

The SSX2 data showed the worst performance due to the high noise level, which is straightforward to understand. The KFX2 and KFIR methods use the same filtering procedures, but have different decision variables: WSSE and individual residuals, respectively. For a fixed false alarm rate, the χ² method has a larger threshold, because it needs to consider the noise of the other fault-free measurements contained within the WSSE calculations. In contrast, because the individual methods only consider a single residual at a time, they have smaller thresholds. This difference is responsible for the better detection performance and availability of the KFIR method.

The KFIR and GSFIR methods use the same FDE algorithm, but they use the overbounded Gaussian and non-Gaussian error model, respectively. The impact of the different error models can be described as follows. In equations (28) and (29), the important factors are the detection threshold, T, and the position error noise, VPL _n. When T is small, the system can detect smaller biases and smaller values of T and VPL _n lead to lower protection levels that result in higher availability. Figure 9 shows the one-sided tail probability of the residuals of the KFIR method (dotted) and the GSFIR method (solid). The fault detection threshold is the residual at which the tail probability is equal to P_fa/2, as defined previously. In Figure 9, the thresholds for the KFIR and GSFIR methods differ by D1, due to the overbounding effect, and the difference of the residual values at which the tail probabilities were equal to the missed detection rate caused the discrepancy in the noise part. A projection of the difference, D2, to the position domain is shown as the protection level difference.

(39)

$T\vert _{KFIR} \equals T\vert _{GSFIR} \plus D1$

(40)

$VPL_{n} \vert _{KFIR} \equals VPL_{n} \vert _{GSFIR} \plus \sum\limits_{i \equals \colon }^{N_{M} } {{\bf H}_{W}^{ \plus } \lpar 3\comma i\rpar \cdot D2\lpar i\rpar}$

Equations (39) and (40) show that the GSFIR method has a smaller minimum detectable bias and both equations indicate that the GSFIR method has smaller protection levels.

Figure 9. The effect of the overbounded sigma.

7. CONCLUSIONS

In this work, we have developed an integrity-monitoring algorithm that tests measurement residuals individually using a Gaussian sum filter. The simulation results showed that our proposed algorithm detected smaller measurement biases and generated protection levels smaller by about 35%, which implies that systems with the algorithm can provide better integrity-monitoring performance than other methods that assume an overbounded Gaussian measurement. The data show that the conservative sigma overbounding methods degrade the integrity-monitoring performance. Inversely, using a true error model or more accurate overbounding methods, systems can be more available without loss of integrity or continuity.

Although our proposed algorithm requires a high computational capacity, the real time process was possible by using a conventional PC. Therefore it can be used in many sophisticated systems such as airliners, military aircraft, GNSS reference stations, etc. Moreover it can be used for development and performance assessment of sigma overbounding methods in conventional systems.

In the future, more computationally efficient algorithms need to be developed for our method to be widely used. In addition, because our algorithm tests the residuals individually, it can be used to detect simultaneous multiple faults, which is another interesting topic for future research.

ACKNOWLEDGMENTS

This research was supported in part by the Institute of Advanced Machinery and Design, Institute of Advanced Aerospace Technology at Seoul National University, and the BK-21.

References

REFERENCES

Brown, R. G., Hwang, P. Y. C. (1997). Introduction to Random Signals and Applied Kalman Filtering. John Wiley & Sons.Google Scholar

Alspach, D. L. and Sorenson, H. W. (1972). Nonlinear Bayesian estimation using Gaussian sum approximations. IEE Transactions on Automatic Control, Vol. AC-17, No. 4, 439–447.Google Scholar

Fruhwirth, R. (1997). Track fitting with non-Gaussian noise. Computer Physics Communications, Vol. 100, 1–16. Elsevier Science.Google Scholar

Hanebeck, U. D., Briechle, K., and Rauh, A. (2003). Progressive Bayes: A New Framework for Nonlinear State Estimation. Proceedings of SPIE, vol. 5099. AeroSense Symposium, 256–267.CrossRef Google Scholar

Kim, D., Yun, Y., Park, B., Jeon, S., Sohn, Y. and Kee, C. (2006). Development and preliminary test results of Korean WADGPS test bed using NDGPS infrastructure in Korea. ION GPS/GNSS 2006, Fort Worth, TX.Google Scholar

McGraw, G. A., Murphy, T., Brenner, M., Pullen, S. and Van Dierendonck, A. J. (2000). Development of the LAAS accuracy models. Proceedings of the ION GPS 2000.Google Scholar

Parkinson, B. W. and Axelrad, P. (1988). Autonomous GPS integrity monitoring using the pseudorange residual. The Journal of Navigation, 35, 255–274.Google Scholar

Rife, J., Pullen, S. and Pervan, B. (2004). Core overbounding and its implications for LAAS integrity. Proceedings of the ION GNSS 2004, Sept.Google Scholar

RTCA Special Committee 159 Working Group 2 (2001). Minimum operational performance standards for global positioning system/wide area augmentation system airborne equipment. RTCA Document Number DO-229C.Google Scholar

Sayim, I., Pervan, B., Pullen, S. and Enge, P. (2002). Experimental and theoretical results on the LAAS sigma overbound. Proceedings of the ION GPS 2002.Google Scholar

Shively, C. A. and Braff, R. (2000). An overbound concept for pseudorange error from the LAAS ground facility. MP00W0000138, The MITRE Corporation.Google Scholar

Trimble Navigation Limited. (1996). Series 4000 Receiver Reference.Google Scholar

Yun, Y., Kim, D. and Kee, C. (2006). A fault detection and exclusion algorithm using particle filters for non-Gaussian GNSS measurement noise, 2006. International Symposium on GPS/GNSS, Jeju, Korea.Google Scholar