Escher A.C., Ecole Nationale de l'Aviation Civile, France

Macabiau C., Ecole Nationale de l'Aviation Civile, France

Roturier B., Service Technique de la Navigation Aerienne, France

Martin N., Thales Avionics, France

**GPS/IRS HYBRIDIZATION: DEFINITION OF EXCLUSION RADIUS USING SOLUTION
SEPARATION METHOD**

##### ABSTRACT

When hybridizing GPS with IRS, the integrity of GPS signals has to be
checked so that slowly growing errors on GPS measurements don't affect
inertial calibration.

Several solutions exist to perform fault detection and exclusion of GPS
signal.

Some of these depend only on GNSS information but the capacity of RAIM
(Receiver Autonomous Integrity Monitoring) is limited. This algorithm
greatly depends on the constellation geometry (at least 5 SV are needed
for detection) and doesn't have sufficient availability for stringent
phases of flight. AAIM (Aircraft Autonomous Integrity Monitoring)
algorithms are taking into account additional onboard information such
as IRS. The Solution Separation method can be implemented for RAIM or
AAIM.

The aim of this paper is to describe the technique inspired from
Solution Separation used to compute isolation radius in addition to the
classical detection protection level, and to present some simulation
results we obtain by implementing this algorithm in a tightly coupled
Kalman.

After reviewing civil aviation requirements and defining the tightly
coupled Kalman filter that was used, the detection and isolation
satellite failure method is described. This method is inspired from the
Solution Separation method which is generalised to N subsets of N-1
Kalman filters used to define the isolation radius as presented in
\[2\].

Finally simulated results of the algorithm for NPA are shown.

##### 

##### INTRODUCTION

Tight integration of GPS and inertial information allows enhancing GNSS
satellite FDE (fault detection and exclusion) and continuity performance
(calibrated IRS may be used for coasting when GNSS information is
unavailable).

Solution Separation method has first been introduced by Brenner in 1995
\[1\]. Assuming that only one failure may occur at a time, fault
detection is performed by estimating the separation between the solution
given by a main filter and solution from different sub-filters that each
exclude measurement from one satellite. This separation is an estimate
of the impact of the satellite failure on the IRS position error
correction using GPS measurements.

The FDE method presented in this paper is an extension of Brenner's
solution separation as it has been defined in \[2\]. By introducing one
more rank of sub-sub-filters (each one is excluding measurement from 2
satellites) the algorithm succeeds in isolating the faulty satellite and
proposes a computation of exclusion level.

The aim of this paper is to present some simulation results obtained by
implementing the FDE algorithm with false alarm rate and missed
detection probability corresponding to NPA phase of flight.

A point is first made on civil aviation assumptions concerning tightly
integrated GPS/inertial systems. The FDE process is described and the
method for computing protection level is defined. Finally simulated test
results are presented.

##### 

##### CIVIL AVIATION ASSUMPTIONS

**[FDE requirements]{.underline}**

Appendix R to RTCA/DO 229C \[3\] clarifies FDE requirements to GPS/IRS
application for en-route to non-precision approach. (Table 1)

  ----------------------------------------------- -----------------------
  Missed detection probability                    0.001

  False detection rate (SA off)                   1/3 10^-6^/hr

  Probability (p~MI~) of exceeding HPL            10^-7^/hr

  Rare normal performance rate                    \<10^-7^/sample

  Failed exclusion probability                    0.001

  Continuity                                      1/3 10^-6^/hr
  ----------------------------------------------- -----------------------

Table 1. FDE requirements.

The rare normal performance rate is the probability that the horizontal
error exceeds the horizontal alert limit in a fault free case ie in the
case of noise only. In case of SA off its value is 10^-7^/sample.
Demonstration in \[3\] shows that the contribution of the rare normal
performance limit $\text{HPL}_{0}$ which depends on this probability is
essential in integrated systems. In the system configuration presented
here $\text{HPL}_{0}$ is the rare normal performance limit of the
primary Kalman filter.

**[FDE performance parameters]{.underline}**

Nominal fault detection performance involves three parameters (test
metric, decision threshold and protection level$\text{PL}$).

The test metric is an observed quantity that is compared to a decision
threshold. The decision threshold is chosen on the basis of statistical
characteristics of the test metric so that false alert doesn't occur
more than a false alert rate$P_{\text{FA}}$. This rate is defined for
each phase of flight and depends on the continuity requirement
(continuity loss-of-function per hour or phase duration) imposed by ICAO
and the GNSS measurements time correlation.

The protection level ($\frac{\text{HPL}}{\text{VPL}}$) is the position
error that the FDE algorithm guarantees will not be exceeded without
being detected by the fault detection function. It depends on the values
of missed detection probability $P_{\text{MD}}$ and$P_{\text{FA}}$.

As it will be described in the next paragraph, exclusion of the faulty
satellite is performed by examining each subset of $N - 1$ satellites
(where $N$ is the numerous of tracked satellites) and searches the
subset without a fault detection condition. Therefore fault exclusion
performance involves some similar parameters: a test statistics,
decision threshold and exclusion level. The exclusion level is a radius
in the horizontal/vertical $(\frac{\text{HEL}}{\text{VEL}})$ plane that
guarantees that any errors on or beyond this radius due to a satellite
failure will be eliminated because of the exclusion function. It is
given as the largest $\text{PL}$ of the subsets of $N - 1$ satellites
computed in function of failed exclusion probability $P_{\text{FE}}$
and$P_{\text{FA}}$.

##### SOLUTION SEPARATION METHOD

**[Fault Detection]{.underline}**

Fault Detection algorithm consists in maintaining a primary Kalman
filter which incorporates the N measurements of the whole satellite
constellation, and as many Kalman sub-filters as the number of
satellites (i.e. N sub-filters). The primary filter $F_{\text{00}}$
provides IRS correction. Sub-filters which each of these incorporate the
measurements from N-1 satellites are dedicated to detection only. They
are noted$F_{0n}\ ,n = 1\ldots N$. Figure 1.

Figure 1. Detection filters hierarchy.

Solution separation between $F_{\text{00}}$ and each $F_{0n}$ is
estimated: it is the difference in the horizontal/vertical plane between
the state vector estimated by $F_{\text{00}}$ and the state vector
estimated by$F_{0n}$.

Fault detection is performed by monitoring the separation between the
main-solution and each of sub-solutions and comparing it to a computed
detection threshold that depends on the separation statistics
and$P_{\text{FA}}$. When a failure occurs at least one solution
separation will exceed this threshold. This method also guarantees a
protection radius against any failure (even slow ramp) in straight
relationship with the threshold and the value of the
required$P_{\text{MD}}$.

**[Fault Exclusion]{.underline}**

It assumes that there is only one satellite failure at a time.

In order to perform exclusion for each sub-filter N-1 sub-sub-filters
are also maintained. Each of these is excluding the measurement excluded
by its "parent" sub-filter, and the measurement from a different
satellite. These sub-sub-filters are noted
$F_{\text{nm}}\ ,\ n = 1\ldots N,\ m = 1\ldots N,\ n \neq m\text{.}$Figure
2.

![](media/image1.wmf){width="2.4925196850393703in"
height="0.35in"}![](media/image2.wmf){width="1.353543307086614in"
height="0.3637795275590551in"}![](media/image3.wmf){width="0.4507874015748031in"
height="0.35in"}

Figure 2. Exclusion filters hierarchy.

When there is a detection, the separation between each sub-sub-filter
and its parent sub-filter is computed and compared to an exclusion
threshold that depends on the expected separation statistics
and$P_{\text{FA}}$. If for sub-filter$F_{0n}$, there exists one solution
separation such that one separation between it and sub-sub-filter
$F_{\text{nm}}$ exceeds the threshold, it can't be the faulty satellite.
But if there is only one sub-filter $F_{0n}$ for which all separations
between its solution and its sub-sub-filters $F_{\text{nm}}$ are under
the threshold, satellite $n$ is the faulty satellite.

##### GPS/IRS KALMAN HYBRIDIZING FILTERS

All of the Kalman filters are running in the same way.

For any of them, the implemented system is composed of three units: an
inertial unit (Inertial measurement Unit + Inertial Reference System), a
GNSS receiver (GPS measurements) and an integration process (Kalman
filters) that also performed FDE function. Figure 3.

Figure 3. GPS/IRS architecture

All measurements and simulations have been made with Matlab.

**[Inertial unit]{.underline}**

It provides to the integration process Matlab-generated inertial
position, velocity and attitude angles and gyrometric and accelerometric
measurements.

![](media/image4.wmf){width="0.5200787401574803in"
height="0.35in"}![](media/image5.wmf){width="0.49921259842519683in"
height="0.4437007874015748in"}

Figure 4. Inertial Unit

**Inertial Measurement Unit: IMU**

This unit generates realistic gyro and accelero measurements at
$\text{100}\ \text{Hz}$rate from the data of an aircraft trajectory and
attitude angles evolution. These measurements are respectively:

$a_{\text{m/i}}^{m}$ non-gravitational acceleration of the mobile
relative to the inertial (absolute) frame in the mobile frame

$\omega_{\text{m/I}}^{m}$angular rate of the mobile relative to the
inertial frame in the mobile frame

Sensor noises and biases are also modelled in order to simulate
different inertial sensors quality.

**Inertial Reference System: IRS**

IMU measurements are processed by Strapdown inertial navigation to
provide inertial solution at$\text{100}\ \text{Hz}$. The navigation
function gives the mobile positions and velocities relative to the earth
frame in the navigation frame. Strapdown inertial navigation scheme is
described in figure 5

![](media/image6.wmf){width="0.4925196850393701in"
height="0.35in"}![](media/image7.wmf){width="0.7952755905511811in"
height="0.4582677165354331in"}![](media/image8.wmf){width="0.2692913385826772in"
height="0.2074803149606299in"}![](media/image8.wmf){width="0.2692913385826772in"
height="0.2074803149606299in"}![](media/image9.wmf){width="0.6035433070866142in"
height="0.35in"}![](media/image10.wmf){width="0.5200787401574803in"
height="0.35in"}![](media/image10.wmf){width="0.5200787401574803in"
height="0.35in"}![](media/image11.wmf){width="0.41732283464566927in"
height="0.3515748031496063in"}![](media/image11.wmf){width="0.41732283464566927in"
height="0.3515748031496063in"}

Figure 5. IRS platform computations

$B_{n}^{m}$is the coordinate transformation matrix that rotates the
mobile body frame $(m)$ acceleration measurements to the local north,
east and down navigation frame $(n)$. Its components are computed from
the integration of gyro measurements.

The attitude determination function determines $B_{n}^{m}$ and the
attitude of the mobile body frame to the navigation frame that is to say
Euler angles, where $\phi,\ \theta,\ \psi$ stand for roll, pitch, yaw
respectively. They are the rotation angles that allow passing from the
$m -$frame to the $n -$frame.

**Attitude determination method**

The attitude algorithm used is the quaternion method ($\int_{}^{}q$)
which offers better numerical and stable characteristics than Euler or
Direct Cosine methods.\[4\]

The single rotation that allows transformation from the $m -$frame to
the $n -$frame is represented by the quaternion$q$. It evolves in
accordance with a simple differential equation that is solved using
Edwards' method. See \[4\]

**[GPS Receiver]{.underline}**

Pseudo range measurements (PRs) are generated at $1\text{Hz}$ . Random
noise is added on each satellite measurement and is composed of a
correlated noise due to inospheric delay ($b_{\text{iono}}$) and
white-Gaussian noise due to noise process ($b_{\text{PR}}$).

Figure 6. GPS measurements generator

**[GPS/IRS Kalman filtering system]{.underline}**

As shown in figure 3, the output position of the whole system is the
corrected inertial position. The role of Kalman filters is to estimate
inertial errors using GPS measurements in order to correct inertial
outputs as done in figure 7. This may be done in an open loop manner or
in a closed loop manner.

The dynamical evolution of the system is given by inertial error model
equations. The measurement vector consists of the difference between two
PRs to each satellite (GPS PR and predicted PR computed with inertial
data).

![](media/image12.wmf){width="0.5389763779527559in"
height="0.2374015748031496in"}![](media/image13.wmf){width="1.338976377952756in"
height="0.15551181102362205in"}![](media/image14.wmf){width="0.9098425196850394in"
height="0.14291338582677166in"}![](media/image15.wmf){width="0.547244094488189in"
height="0.24645669291338582in"}![](media/image16.wmf){width="4.5275590551181105e-2in"
height="8.543307086614173e-2in"}![](media/image17.wmf){width="1.198031496062992in"
height="0.2637795275590551in"}Figure 7. Open loop GPS/IRS hybridization

**Filter error state model**

Each component $\delta x$of the Kalman filter state vector stands for
the difference between the true value $x$ and the measured
$\overset{\sim}{x}$(or computed$\hat{x}$) value.

The state vector consists of a 17 error state variables

$$\Delta X = \lbrack\rho,\text{δv},\text{δp},\text{δω},\text{δf},\text{δb}\rbrack^{T}$$

where

  ----------------------------------------------- --------------------------------------------------
  $$\rho$$                                        ***Attitude error vector***

  $$\text{δϕ}\ \left( \varepsilon_{N} \right)$$   roll error

  $$\text{δθ}\ \left( \varepsilon_{E} \right)$$   pitch error

  $$\text{δψ}\ \left( \varepsilon_{D} \right)$$   yaw error

  $$\text{δv}$$                                   ***Inertial velocity error vector***

  $$\text{δv}_{N}$$                               north velocity error

  $${\delta v}_{E}$$                              east velocity error

  $${\delta v}_{D}$$                              vertical velocity error

  $$\text{δr}$$                                   ***Inertial position error***

  $$\delta L$$                                    latitude error

  $$\delta G$$                                    longitude error

  $$\delta h$$                                    altitude error

  $$\text{δω}$$                                   ***Gyro drift*** $m -$***frame***

  $$\text{δω}_{x}$$                               x-axis gyro drift $(m)$

  $$\text{δω}_{y}$$                               y-axis gyro drift $(m)$

  $$\text{δω}_{z}$$                               z-axis gyro drift $(m)$

  $$\text{δf}$$                                   ***Accelero bias*** $m -$***frame***

  $${\delta f}_{x}$$                              x-axis accelero bias $(m)$

  $${\delta f}_{y}$$                              y-axis accelero bias $(m)$

  $${\delta f}_{z}$$                              z-axis accelero bias $(m)$

  $$\text{δb}$$                                   ***Receiver clock error vector***

  $$b$$                                           receiver clock bias

  $$\overset{˙}{b}$$                              receiver clock bias drift
  ----------------------------------------------- --------------------------------------------------

Table 1. Kalman filter state vector.

Many different inertial error models are available in literature. They
are actually equivalent \[5\]. In our system the inertial navigation
error model applied is

$\begin{array}{r}
\overset{˙}{\rho} = \text{δω} + \omega_{\text{In}}^{n} \times \rho \\
\delta\overset{˙}{v} = \text{δf} + \left( \omega_{\text{Ie}}^{n} + \omega_{\text{In}}^{n} \right) \times \text{δv} + \rho \times a_{\text{mI}}^{n} \\
\delta\overset{˙}{r} = \text{δv} + \omega_{\text{en}}^{n} \times \text{δr}
\end{array}$ (1)

Using inertial data as nominal trajectory these non-linear equation are
linearized and lead to the matrix presentation of the dynamical
evolution equation of the linearized Kalman filter \[6\]:

$\overset{˙}{X} = F \cdot X + v$ (2)

where $F$ is the state transition matrix

$F = \begin{bmatrix}
F_{\text{ρρ}} & F_{\text{ρv}} & F_{\text{ρr}} & F_{\text{ρω}} & 0 & 0 & 0 \\
F_{\text{vρ}} & F_{\text{vv}} & F_{\text{vr}} & 0 & F_{\text{vf}} & 0 & 0 \\
F_{\text{rρ}} & F_{\text{rv}} & F_{\text{rr}} & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & F_{\text{ωω}} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & F_{\text{ff}} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}$ (3)

and the noise state vector$w$is a zero-mean Gaussian white noise vector,
whose components are all independent.

Let $\Phi$ be the discret form of$F$.

Let be $Q$ the state noise covariance matrix:

$Q = E\left\lbrack v \cdot v^{T} \right\rbrack$ (4)

**Filter error measurement model**

Each component of the measurement vector $z$ at filter's input is the
difference between two PRs to each satellite. One is the measured PR
input from the GPS receiver; the other is the predicted PR computed on
the basis of the satellite positions obtained from the GPS receiver and
the user location as calculated by the IRS.

The coefficients of the measurement matrix $H$are direction cosine
computed from the GPS navigation equation linearization.

Let $R$ be the measurement noise covariance matrix: all measurement
noises are independent.

The size of the $z -$vector and the$R$ and$H -$matrices depend on the
number of tracked satellite.

**Kalman filters implementation**

Let $\text{ij}$be the subscript for the Kalman filter$F_{\text{ij}}$.
The maximum number of tracked satellites is$N = \text{10}$.

There is one primary filter, N sub-filters, and N\*(N-1) sub-sub-filters
that are running in parallel.

The estimated state vector for each are

primary filter: $\Delta{\hat{x}}_{\text{00}}$

sub-filter: $\Delta{\hat{x}}_{\text{0n}}\ ,\ n = 1,N$

sub-sub-filter:$\Delta{\hat{x}}_{\text{nm}}\ ,\ n = 1,n;\ m = 1,N;\ m \neq n$

The estimation error covariance matrix for each is

primary filter: $P_{\text{00}}$

sub-filter: $P_{\text{0n}}\ ,\ n = 1,N$

sub-sub-filter:$P_{\text{nm}}\ ,\ n = 1,n;\ m = 1,N;\ m \neq n$

where

$$P_{\text{ij}} = E\left\lbrack \text{δx}_{\text{ij}} \cdot \text{δx}_{\text{ij}^{T}} \right\rbrack,\ \text{δx}_{\text{ij}} = \Delta{\hat{x}}_{\text{ij}} - \text{Δx}$$

$\Delta x$is the true error between estimated IRS and true position:
$\Delta x = {\hat{X}}_{\text{IRS}} - X$

The Kalman Gain matrix for each is

primary filter: $K_{\text{00}}$

sub-filter: $K_{\text{0n}}\ ,\ n = 1,N$

sub-sub-filter:$K_{\text{nm}}\ ,\ n = 1,n;\ m = 1,N;\ m \neq n$

Because GPS and IRS data are not available at the same rate, the error
model is updated at $\text{100}\text{Hz}$and Kalman filters measurements
are updated at$1\text{Hz}$.

So the Kalman corrections are available every second.

Next section details implementation of FDE algorithm.

**[FAULT DETECTION AND HORIZONTAL PROTECTION LEVEL]{.underline}**

Fault detection is performed by estimating the separation between the
primary filter estimate and each of the N sub-filters estimates in the
horizontal plane.

At each estimation time the discriminator for the n^th^ sub-filter is
based on Solution Separation vector between the primary filter and the
n^th^ sub-filter:

$\begin{array}{r}
\text{dx}_{\text{0n}^{+}}(k) = \Delta{\hat{x}}_{\text{00}^{+}}(k) - \Delta{\hat{x}}_{\text{0n}^{+}}(k) \\
\  = \ \text{δx}_{\text{00}^{+}}(k) - \text{δx}_{\text{0n}^{+}}(k),\ n = 1\ldots N
\end{array}$ (5)

whose statistics are described by the covariance matrix

$\text{dP}_{\text{0n}^{+}}(k) = E\left\lbrack \left( \text{dx}_{\text{0n}^{+}}(k) \right) \cdot \left( \text{dx}_{\text{0n}^{+}}(k) \right)^{T} \right\rbrack$
(6)

Since $\text{dx}_{\text{0n}}$and $\text{δx}_{\text{00}}$ are statically
dependent, let us form the 34x1 error vector

$\text{δx}^{+}(k) = \left\lbrack \text{δx}_{\text{00}^{+}}(k),\ \text{Δx}_{\text{0n}^{+}}(k) \right\rbrack^{T}$
(7)

whose statistics are described by

$\begin{array}{r}
\text{dP}_{\text{0n}}^{\text{dual}^{+}}(k) = E\left\lbrack \left( \text{δx}^{+}(k) \right) \cdot \left( \text{δx}^{+}(k) \right)^{T} \right\rbrack \\
\  = \begin{bmatrix}
P_{\text{00}^{+}}(k) & \text{cross}\ \text{corr}\text{.} \\
\text{cross}\ \text{corr}\text{.} & \text{dP}_{\text{0n}^{+}}(k)
\end{bmatrix}
\end{array}$ (8)

Given the fact that

$\begin{array}{r}
\text{dx}_{\text{0n}^{+}}(k) = \Delta{\hat{x}}_{\text{00}^{+}}(k) - \Delta{\hat{x}}_{\text{0n}^{+}}(k) \\
\  = \delta{\hat{x}}_{\text{00}^{+}}(k) - \Delta{\hat{x}}_{\text{0n}^{+}}(k) + \text{Δx}
\end{array}$ (9)

and

$\text{δx}^{-}(k) = \underset{\Pi}{\overset{\begin{bmatrix}
\Phi & 0 \\
0 & \Phi
\end{bmatrix}}{\underbrace{}}} \cdot \text{δx}^{+}(k - 1) + \begin{bmatrix}
v \\
0
\end{bmatrix}$ (10)

one can show that

$\text{dP}_{\text{0n}}^{\text{dual}^{+}}(k) = \Sigma \cdot \text{dP}_{\text{0n}}^{\text{dual}^{-}}(k) \cdot \Sigma^{T} + \Gamma \cdot R \cdot \Gamma^{T}$
(11)

with

$\Sigma = \begin{bmatrix}
I - K_{\text{00}} \cdot H_{\text{00}} & 0 \\
\left( K'_{\text{0n}} - K_{\text{00}} \right) \cdot H_{\text{00}} & I - K'_{\text{0n}} \cdot H_{\text{00}}
\end{bmatrix}$ (12)

$\Gamma = \begin{bmatrix}
K_{\text{00}} \\
K_{\text{00}} - K'_{\text{0n}}
\end{bmatrix}$ (13)

$$\text{δP}_{\text{0n}}^{\text{dual}^{-}}(k) = \Pi \cdot \text{δP}_{\text{0n}}^{\text{dual}^{+}}(k - 1) \cdot \Pi^{T} + \begin{bmatrix}
Q & 0 \\
0 & 0
\end{bmatrix}$$

> \(14\)

$K'_{\text{0n}}$ is constructed from $K_{\text{0n}}$so that it operates
on the full set of measurements by extending the 17x(N-1) matrix
$K_{\text{0n}}$ sub-solution matrix with a n^th^ zeroed column.

**Detection threshold** $D_{0n}$ **vs. test statistics** $d_{0n}$

The statistic of the horizontal separation
vector$X_{H} = \text{dx}_{\text{0n}^{+}}\left\lbrack \text{7:8} \right\rbrack$,
is described by
$L_{\text{0n}} = \text{dP}_{\text{0n}^{+}}\left\lbrack \text{7:8,7:8} \right\rbrack$.
Because separations on North and East axes are correlated,
$L_{\text{0n}}$is projected in an orthogonal plane so that,

$L_{\text{0n}} = P_{} \cdot \Delta \cdot P_{{}^{T}}$ (15)

and $X_{} = P_{{}^{T}} \cdot X_{H}$ is a Gaussian vector whose
covariance matrix is the diagonal matrix $\Delta$ (eigenvalues matrix
of$L_{\text{0n}}$).

One of the two eigenvalues is dominating. Let $\lambda^{\text{dP}}$ be
the approximate variance ($\Delta(1,1)\ \text{or}\ \Delta(2,2)$) of the
error in the horizontal plane.

For each sub-filter the decision threshold is set so that $d_{0n}$will
exceed $D_{0n}$ with the probability$P_{\text{FA}}$. Since all
sub-filters have the same chance of false alarm

$D_{0n} = \sqrt{\lambda^{\text{dP}}} \cdot Q^{- 1}\left( \frac{P_{\text{FA}}}{2N} \right)$with
$Q^{- 1}(u) = \text{erfc}(u)$

and $d_{0n} = \left| X_{}(1) \right|\ \text{or}\ \left| X_{}(2) \right|$
(16)

depending on
$\lambda^{\text{dP}} = \Delta(1,1)\ \text{or}\ \lambda^{\text{dP}} = \Delta(2,2)$

The test for detecting one failure satellite is

$H_{0}$ no detection: $d_{0n} \leq D_{0n}$

$H_{1}$ detection: $d_{0n} > D_{0n}$ for at least one $n,\ n = 1,N$

**Horizontal protection level HPL**

By definition the $\text{HPL}$is the error bound that contains the
primary filter error with a probability of $1 - p_{\text{MD}}$ when
$d_{0n} = D_{0n}$.

Since for each sub-filter

$\text{δx}_{\text{00}^{+}}(k) = \text{dx}_{\text{0n}^{+}}(k)\  + \text{δx}_{\text{0n}^{+}}(k),\ n = 1\ldots N$
(17)

at detection (for the n^th^ subset)

$\text{HPL}_{n} = D_{0n} + a_{0n}$ (18)

where $a_{0n}$ is an upper bound of the n^th^ sub-filter horizontal
error ${\Delta x}_{\text{0n}^{+}}\left\lbrack \text{7:8} \right\rbrack$
whose distribution is described by
$L = P_{\text{0n}^{+}}\left\lbrack \text{7:8,7:8} \right\rbrack$.

Let $\lambda^{P_{0n}}$be the maximum eigenvalue of $L$

$a_{0n} = \sqrt{\lambda^{P_{0n}}} \cdot Q^{- 1}\left( 1 - p_{\text{MD}} \right)$
(19)

Besides we have to consider the rare normal performance that is to say
the contribution of the primary filter $\text{HPL}_{0}$

$\text{HPL}_{0} = \sqrt{\lambda^{P_{\text{00}}}} \cdot Q^{- 1}\left( \frac{P_{\text{ff}}}{2} \right)$
(20)

where $P_{\text{ff}}$ is the rare normal performance rate and
$\lambda^{P_{\text{00}}}$ the maximum eigenvalue of
$P_{00^{+}}\left\lbrack \text{7:8,7:8} \right\rbrack$.

The horizontal protection level is

$\text{HPL} = \text{max}\left\{ \text{HPL}\text{,max}\left( \text{HPL}_{n} \right) \right\}\ ,\ n = 1,N$
(21)

**[FAULT EXCLUSION AND HORIZONTAL EXCLUSION LEVEL]{.underline}**

Fault exclusion is accomplished in the same manner than detection but
one layer down in the filters hierarchy.

The discriminator of each of the sub-sub-filter is the separation
between each sub-sub-filter estimate solution and its parent's
sub-filter estimate solution in the horizontal plane. For the n^th^
sub-filter they are

$\begin{array}{r}
\text{dx}_{\text{nm}^{+}}(k) = \Delta{\hat{x}}_{\text{0n}^{+}}(k) - \Delta{\hat{x}}_{\text{nm}^{+}}(k) \\
\  = \ \text{δx}_{\text{0n}^{+}}(k) - \text{δx}_{\text{nm}^{+}}(k),\ m = 1\ldots N,m \neq n
\end{array}$ (22)

whose statistics are described by the covariance matrix

$$\begin{array}{r}
\text{dP}_{\text{nm}^{+}}(k) = E\left\lbrack \left( \text{dx}_{\text{nm}^{+}}(k) \right) \cdot \left( \text{dx}_{\text{nm}^{+}}(k) \right)^{T} \right\rbrack \\
E\left\lbrack \left( {\delta x}_{\text{0n}^{+}}(k) - {\delta x}_{\text{nm}^{+}}(k) \right)\left( {\delta x}_{\text{0n}^{+}}(k) - {\delta x}_{\text{nm}^{+}}(k) \right)^{T} \right\rbrack \\
P_{\text{0n}^{+}}(k) - P_{\text{nm}}^{\text{cross}^{+}}(k) - \left( P_{\text{nm}}^{\text{cross}^{+}}(k) \right)^{T} + P_{\text{nm}^{+}}(k)
\end{array}$$

\(23\)

**Detection threshold** $D_{\text{nm}}$ **vs test statistics**
$d_{\text{nm}}$

For each sub-sub-filter the decision threshold is set so that (assuming
that each filter $F_{\text{nm}}$has the same chance for a false alarm)

$D_{\text{nm}} = \sqrt{\lambda^{\text{dP}_{\text{nm}}}} \cdot Q^{- 1}\left( \frac{P_{\text{FA}}}{2(N - 1)} \right)$
(24)

where $\lambda^{\text{dP}_{\text{nm}}}$ the maximum eigenvalue of
$L_{\text{nm}} = P_{\text{nm}^{+}}\left\lbrack \text{7:8,7:8} \right\rbrack$

The test statistic associated to $D_{\text{nm}}$ is

$d_{\text{nm}} = \left| X_{m_{}}(1) \right|\ \text{or}\ \left| X_{m_{}}(2) \right|$
(25)

where
$X_{m\bot = P_{m\bot{}^{T} \cdot \text{dx}_{\text{nm}^{+}}\left\lbrack \text{7:8} \right\rbrack}}$

with$L_{\text{nm}} = P_{m\bot \cdot \Delta \cdot P_{m\bot{}^{T}}}$

Satellite $r$is excluded as the failed satellite if and only if

$d_{\text{rm}} < D_{\text{rm}}$ for all $m \neq r$

and

$d_{\text{nm}} \geq D_{\text{nm}}$ for at least one $m \neq n$for all
$n \neq r$

**Horizontal exclusion level HEL**

By definition $\text{HEL}$ is the error bound that contains the primary
filter error with a probability of $1 - p_{\text{FE}}$ when the failure
is excluded.

"For a failed exclusion to occur, the solution of one of the sub-filters
$\Delta{\hat{x}}_{\text{0n}}$which does contain the failed satellite
must be separated from one of its sub-sub-filter's solution
$\Delta{\hat{x}}_{\text{nm}}$ by less than the threshold $D_{\text{nm}}$
plus the sub-sub-filter position error" \[2\]

Let $a_{\text{nm}}$ be the sub-sub-filter horizontal estimation error
bound and $\lambda^{P_{\text{nm}}}$ be the maximum eigenvalue of
$P_{\text{nm}^{+}}\left\lbrack \text{7:8,7:8} \right\rbrack$

$a_{\text{nm}} = \sqrt{\lambda^{P_{\text{nm}}}} \cdot Q^{- 1}\left( 1 - p_{\text{FE}} \right)$
(26)

For the nth sub-filter $\text{HEL}_{n}$is given by

$\text{HEL}_{n} = \text{max}\left\{ D_{\text{nm}} + a_{\text{nm}} \right\}\ ,\ m = 1,N$
(27)

And the horizontal exclusion level is

$\text{HEL} = \text{max}\left\{ \text{HEL}_{n} \right\}\ ,\ n = 1,N$
(28)

**[RESULTS OF SIMULATIONS]{.underline}**

The whole system has been implemented in Matlab language. This section
presents the firsts results obtained in testing the FDE algorithm for
detection and isolation of failures causing slowly increasing ramp
error.

**Conditions for the simulation**

GPS constellation

-   Mask angle 5° (10 SVs, no constellation change)

GPS receiver

-   Total receiver noise 12.5m (1σ)

-   Correlated noise: 1srt order Markov process
    ($\tau = \text{15}\text{min,}\sigma = 4\text{.}8m$ )

Inertial error sources

  --------------------------------------- ------------ ------------------
  *Error source (one sigma)*              *Good*       *Low cost*

  Gyro constant bias: °/hr                0.01         1

  Gyro white noise:                       3.10^-5^     0.001
  °/s/$\sqrt{\text{Hz}}$                               

  Accelero constant bias: μg              50           500

  Accel. white noise:                     6.4          64
  μg/$\sqrt{\text{Hz}}$                                
  --------------------------------------- ------------ ------------------

Table 2. Good / low cost IRS error sources

GPS/IRS hybridization is performed in an open-loop manner for good IRS,
and closed and open loop manner for low cost IRS.

FDE algorithm

-   No. of independent tests/hr=4

-   $$P_{\text{FA}} = \frac{1}{3}\text{.}\frac{\text{10}^{- 6}}{\text{hr}} = \frac{1}{\text{12}}\text{.}\frac{\text{10}^{- 6}}{\text{test}}$$

-   $$P_{\text{ff}} = \frac{1}{3}\text{.}\frac{\text{10}^{- 8}}{\text{test}}$$

-   ramp error size (m/s) (time 400s)

> $$0\text{.}1,\ 0\text{.}2,\ 0\text{.}5,\ 1,\ 2,\ 5$$
>
> on the worst case satellite (SV 4).

Trajectory: duration 600s

-   The simulated trajectory is an approach to Toulouse-Blagnac Airport
    from AGN VOR. It includes a turn at time 180s and a descent at time
    280 s.

![](media/image18.png){width="3.2194444444444446in"
height="1.6229166666666666in"}

![](media/image19.png){width="3.2194444444444446in"
height="1.6229166666666666in"}

Figure 8.Simulated trajectory

-   During the turn, heading angle is varying from 098° to 145°. Roll
    angle is varying with a rate of 1 (180°/60s). Pitch angle is also
    varying since the aircraft is supposed to keep its altitude.

-   During the descent only pitch angle is varying according trajectory
    slope.

**Results**

Figure 9 shows Horizontal Protection Level for RAIM and AAIM (HPL), and
Horizontal Exclusion Level for AAIM (HEL) using Solution Separation
algorithm.

AAIM HPL and HEL convergence occurs after 250s of simulation. HPL and
HEL are much smaller (under 30m) than RAIM HPL (70m). AAIM HEL is
greater than HPL by approximately 5m. AAIM HPL and HEL obtained using
closed loop are the same than those obtained using open loop.

Figure 10 shows how detection is achieved on channel 4 in presence of a
failure when GPS measurement noise is only due to white-Gaussian noise.
Before the failure occurs at 400s the discriminator is zero-mean and is
lower than the decision threshold. The failure detection occurs when the
discriminator exceeds the decision threshold. The IRS quality seems to
have an influence: HPL and HPE obtained hybridizing GPS with a good IRS
are in the order of 10m under those obtained hybridizing GPS with a low
cost IRS. But no extensive Monte-Carlo simulation has been performed.

![](media/image20.png){width="3.2194444444444446in"
height="1.8854166666666667in"}

Figure 9. Horizontal protection levels.

![](media/image21.png){width="3.2194444444444446in"
height="1.9597222222222221in"}

Figure 10. Discriminator deviation due satellite failure causing ramp

The following tables show examples of detection and exclusion
performance with only white measurement noise.

In these examples, AAIM algorithm is detecting satellite failure earlier
than RAIM algorithm and the delay between detection and isolation
doesn't exceed 20 s (table 3). Besides, for all the simulations the
faulty satellite was detected and isolated before the horizontal error
grows over around 13m (table 4, 5).

In the case of open-loop, the low cost IRS drifts tend to make the whole
system derivate. In the case of closed loop bad estimates of inertial
errors due to ramp failure tend to miss-calibrate the inertial platform.
Thus position error drift is not only due to satellite failure but also
inertial drift. In these cases detection occurs earlier, these results
might yield us to think that the false alarm rate could increase.

  -------------------- ------------ ------------ ------------ ------------
  ***Failure           RAIM         AAIM                      
  introduced at                                               
  400s***                                                     

                                    **open                    **closed**
                                    loop**                    

                                    **good**     **Low cost   
                                                 IRS**        

  *Ramp: 0.1(m/s)*                                            

  *detection*          \> 600s      566 s        594s         \> 600s

  *SV4 isolation*      no simu.     586 s        \> 600s      \> 600s

  *Ramp: 0.5 (m/s)*                                           

  *Detection*          565 s        451 s        439s         443s

  *SV4 isolation*      no simu.     458 s        445s         446s

  *Ramp: 1 (m/s)*                                             

  *Detection*          500 s        435 s        433s         428s

  *SV4 isolation*      no simu.     439 s        435s         432s

  *Ramp: 5 (m/s)*                                             

  *Detection*          415          411 s        412s         411s

  *SV4 isolation*      no simu.     413 s        413s         412s
  -------------------- ------------ ------------ ------------ ------------

Table 3. RAIM and AAIM times of detection and isolation.

  -------------------- ------------ ------------ ------------ ------------
  *Ramp size*          RAIM         AAIM                      

                                    **Open                    **closed**
                                    loop**                    

                                    **good**     **Low cost   
                                                 IRS**        

  *0.1(m/s)*           \>7.1 m      5.2m         9.5m         \> 6m

  *0.5 (m/s)*          34.1m        5.3m         7m           5.3m

  *1 (m/s)*            32.5m        5.7m         12m          10.5 m

  *5 (m/s)*            8.5m         5.8m         10.6m        9 m
  -------------------- ------------ ------------ ------------ ------------

Table 4. Error size at detection time in the horizontal plane (good IRS)

  -------------------- ---------------- ---------------- ----------------
  *Ramp size*          AAIM                              

                       **Open loop**                     **closed**

                       **good IRS**     **Low cost IRS** 

  *0.1(m/s)*           6.1m             \> 9.5m          \> 6m

  *0.5 (m/s)*          6.5m             9.5m             5.8m

  *1 (m/s)*            7m               13m              11.2 m

  *5 (m/s)*            7.2m             12.5m            10.5 m
  -------------------- ---------------- ---------------- ----------------

Table 5. Error size at AAIM isolation time in the horizontal plane (good
IRS)

One main advantage of this algorithm is that it is able to detect any
type of failures. But it may become a drawback if the process is bad
modelled in the Kalman filter. Simulations were made assuming the
measurement noise was not only white noise. Many false detections (only
due to measurement noise) and false isolations were observed. Indeed in
that case correlated noise was interpreted by the FDE algorithm as a
failure.

Figure 11 and 12 show AAIM HPL with only white measurement noise
generated and mixed measurement noise. In the simulation context of
figure 12, results show a false detection alert at time $\text{404}s$
whereas the ramp error on measurement is not large enough to be detected
in the case illustrated in figure 10. In fact, in figure 12, correlated
noises make the GPS measurement look like being affected by ramp failure
which is certainly the reason of this alert. Figure 13 shows the
miss-modelling impact on the discriminator computed for sub-filter
$F_{\text{04}}$ that is no more a zero-mean variable due to correlated
noise addition on range measurement. On figure 10 the discriminator
deviation is only due to satellite failure causing a ramp.

![](media/image22.png){width="3.2194444444444446in" height="1.89375in"}

Figure 11. White noise measurement noise observation

![](media/image23.png){width="3.2194444444444446in" height="1.89375in"}

Figure 12. White + correlated noise measurement observation

![](media/image24.png){width="3.2194444444444446in"
height="1.9993055555555554in"}

Figure 13. Discriminator deviation due to correlated noise

Solution to take into account the correlated noise on GPS measurement is
to estimate the evolution of this process. The Kalman filters state
vector is extended by adding as many states as the number of tracked
satellites. The i^th^ supplemental state stands for GPS range bias error
for i^th^ satellite.

**[CONCLUSION]{.underline}**

Solution Separation method is a very appealing Fault detection and
Exclusion algorithm because of its rigorousness. This paper presents an
open-loop GPS/IMU system being investigated by the ENAC in order to
start appreciating the Solution Separation algorithm ability of fault
isolation. Simulations were performed assuming NPA requirements and
satellite failure causing a ramp from 0.1m/s to 5m/s. AAIM HPL and HEL
are the order of 20-30 m and are lower than RAIM HPL (closed to 70m).
Isolation of the faulty satellite occurs less than 20 seconds after
fault detection arise. Low cost IRS is affected by inertial drift that
pollutes the hybridized solution causing false alarm in open-loop
hybridization. Finally the whole FDE algorithm performance seems to
greatly depend on Kalman filter miss-modelization. Thus implementation
improvements have to be taken account such as the correlation
measurement noise in the Kalman filter model and will be the issue of
further works.

##### REFERENCES

1\. - M. BRENNER, Integrated GPS/inertial detection availability,
Journal of The Institute of Navigation, Sept. 1995.

2 - K. VANDERWERF, FDE Using Multiple Integrated GPS/Inertial Kalman
Filters in the Presence of Temporally and Spatially Correlated
Ionospheric Errors, ION GPS 2001, 11-14 Sept. 2001, Salt Lake City

3\. - Requirements and Test Procedures for Tightly Integrated
GPS/Inertial Systems Appendix R to DO-229C,.

4.- J. C. RADIX, Systèmes Inertiels à Composants Liés "Strapdown",
Cépadues Educs, 1991.

5.- R. DA, Investigation of a Low-Cost and High-Accuracy GPS/IMU System,
Proceedings of ION GPS, Sept. 1996.

6.- FARRELL J. A., BARTH M., The Global Positionning System & Inertial
navigation, McGraw-Hill, 1998.