M7-T6 — Navigation Integrity: Fault Detection and Protection Levels

Learning objectives

After this lesson, you should be able to:

• Separate accuracy from integrity and state the three quantities an integrity requirement fixes.
• Build the residual test from the least-squares projector and set its threshold from a false-alert probability.
• Compute the smallest detectable bias of each measurement and the position error it causes.
• Identify a measurement that cannot be checked and say what it implies.

The intuition

UAV M3-T13: Innovation monitoring tests whether a filter's innovations are as large as predicted, and UAV M3-T16: Sensor fault detection and isolation covers detection and isolation in general. This lesson specializes the idea to a position fix and goes on to what the others do not compute: the bound that a fault which was not detected can still reach.

Accuracy is not integrity

A fix with 1 m of noise per range is accurate to about a metre. If one beacon has a fault, a biased clock, a reflected signal or a wrong surveyed position, the fix is pulled away and nothing in its statistics shows it. A system used for landing needs more than a typical accuracy: it needs a protection level, a bound on the position error that holds with a stated probability, even allowing for a fault that was missed, and an alert when the bound exceeds what the operation allows. Three numbers define the requirement: the probability of a false alert, the probability of a missed detection, and the alert limit.

Redundancy is the detector

With six pseudoranges and four unknowns there are two more equations than needed. A least-squares fit leaves a residual in each, and if all measurements are consistent the residuals are small and random; the sum of their squares, divided by the noise variance, follows a chi-square distribution with 2 degrees of freedom. A fault in one measurement raises the sum. The test compares it with a threshold chosen so that a fault-free fix sets off a false alert only once in a thousand. With seven beacons the degrees of freedom are three and the threshold on the normalized sum of squares is 16.27.

The smallest fault that can be seen

A fault is detected with probability $1 - p_{md}$ only if it is large enough. For each beacon the smallest detectable bias depends on how much of its error appears in the residual, the diagonal element $S_{ii}$ of the residual projector: a beacon whose error is mostly absorbed by the position solution leaves little in the residual. For Hopper's six ground beacons the smallest detectable bias is 9.4 to 10.5 m, and the position error that bias causes is 2.8 to 4.2 m horizontally. The protection level is that worst case plus the noise term: 6.2 m horizontally.

The beacon that cannot be checked

The balloon beacon is the only one above the vehicle, and the only one that fixes the height. If its range is wrong the position solution adapts the height to match it exactly and the residual is zero: $S_{ii} = 0$. Its fault cannot be detected, whatever the threshold. A bias of 377 m would be needed for even the formula to give a detectable residual, and it would cause a vertical error of 410 m. The beacon that improved the geometry most in the last lesson is the one that the redundancy cannot check, and an integrity-critical system needs a second one.

Integrity needs redundancy. The residual test finds faults that the other measurements contradict, the smallest detectable fault sets the protection level, and a measurement with no redundancy behind it cannot be protected at all.

Definitions

Integrity

The ability of a system to warn the user when it should not be used, with a stated probability of missing a fault.

False-alert probability $P_{fa}$, missed-detection probability $P_{md}$

The two error rates of the fault test.

Residual projector $\mathbf{S}$

$\mathbf{S} = \mathbf{I} - \mathbf{H}(\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}\mathbf{H}^{\mathsf{T}}$: it maps measurement errors to least-squares residuals, with rank n − 4.

Slope

For each measurement, the ratio between the position error caused by a bias in it and the part of that bias that appears in the test statistic.

Protection level

A bound on the position error that holds with the required probability even with an undetected fault.

Derivation: the test statistic and the minimum detectable bias

From the residual of the least-squares fix to the bias at which the test fails to detect.

1. $$\mathbf{r} = \mathbf{z} - \mathbf{H}\hat{\mathbf{x}} = \bigl(\mathbf{I} - \mathbf{H}(\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}\mathbf{H}^{\mathsf{T}}\bigr)\mathbf{z} = \mathbf{S}\mathbf{z}$$

   The least-squares residual is the part of the measurement errors that the fit cannot absorb. S is symmetric and idempotent, $\mathbf{S}^2 = \mathbf{S}$, a projector.

2. $$\mathbf{r}^{\mathsf{T}}\mathbf{r} = \mathbf{z}^{\mathsf{T}}\mathbf{S}\mathbf{z},\quad \mathrm{rank}(\mathbf{S}) = n - 4$$

   The sum of squared residuals. A projector of rank $n-4$ leaves $n-4$ independent directions of the Gaussian noise, each with variance $\sigma^2$.

3. $$T = \mathbf{r}^{\mathsf{T}}\mathbf{r}/\sigma^2 \sim \chi^2_{n-4}$$

   A sum of $n-4$ squares of standard normal variables. The threshold is the quantile that leaves $P_{fa}$ in the tail: 16.27 for 3 degrees of freedom and $P_{fa} = 10^{-3}$.

The minimum detectable bias

1. $$\mathbf{z} = \boldsymbol\nu + b\,\mathbf{e}_i \;\Rightarrow\; \mathbf{r} = \mathbf{S}\boldsymbol\nu + b\,\mathbf{S}\mathbf{e}_i$$

   A bias b on measurement i adds a deterministic part to the residual.

2. $$T \sim \chi^2_{n-4}(\lambda),\quad \lambda = \dfrac{b^2\,\mathbf{e}_i^{\mathsf{T}}\mathbf{S}\mathbf{e}_i}{\sigma^2} = \dfrac{b^2 S_{ii}}{\sigma^2}$$

   The statistic is noncentral chi-square, with noncentrality equal to the squared size of the biased residual over the noise variance.

3. $$\Pr\bigl(\chi^2_{n-4}(\lambda) > T_{\mathrm{th}}\bigr) = 1 - P_{md} \;\Rightarrow\; \lambda = 48.1$$

   Solve for the noncentrality that gives a detection probability of 0.999 with 3 degrees of freedom and a threshold of 16.27.

4. $$b_{\min,i} = \sigma\sqrt{\lambda/S_{ii}},\quad \delta\mathbf{p}_i = \mathbf{a}_i\,b_{\min,i}$$

   The smallest bias that is detected with the required probability, and the position error it causes through the row $\mathbf{a}_i$ of the solution matrix. If $S_{ii} = 0$ no bias is detected.

The residual contains only what the other measurements contradict. The detector sees a bias b at the fraction S_ii of its size, so the smallest detectable bias is sigma times the square root of lambda over S_ii, and a measurement with S_ii = 0 can never be checked.

A worked example

Hopper's range with seven beacons: the six ground beacons and the balloon of M7-T1: The Position Fix: Pseudoranges, Least Squares and Dilution of Precision, Hopper at (40, 20, 60) m, 1 m of range noise.

1. Degrees of freedom and threshold. Seven pseudoranges and four unknowns leave 3 degrees of freedom. For $P_{fa} = 10^{-3}$ the threshold on $T$ is 16.27. For $P_{md} = 10^{-3}$ the noncentrality is $\lambda = 48.1$.
2. Weights. The residual weights $S_{ii}$ of the ground beacons are 0.43 to 0.55. The balloon beacon has $S_{ii} = 0.000$.
3. Smallest detectable bias. For the ground beacons 9.4 to 10.5 m (beacon 5 is best, beacon 0 worst). The position error a bias of that size causes is 2.8 to 4.2 m horizontally, the worst being beacon 0 at 4.2 m.
4. Protection level. Bias term 4.16 m plus a noise term of 3.29 times the largest horizontal standard deviation (2.06 m) gives a horizontal protection level of 6.2 m, against a 1σ noise error of 0.83 m.
5. The balloon. Its formula gives a minimum detectable bias of 377 m and a vertical error of 410 m for it: effectively undetectable. A fault of any size on that beacon goes through to the height.
6. Monte Carlo. With a bias on beacon 0 equal to half the minimum (5.3 m) the test detects in 36.8% of 3000 trials; at the minimum (10.5 m) 99.9%, as designed; at twice the minimum 100%. Without a fault the false-alert rate is 0.10% (target 0.1%).

import math
import random
import gnckit as g

beacons = [(-300.0, -200.0, 2.0), (350.0, -250.0, 5.0), (400.0, 180.0, 0.0), (-250.0, 300.0, 8.0), (50.0, 420.0, 3.0), (60.0, -430.0, 1.0), (-30.0, 10.0, 150.0)]
truth = (40.0, 20.0, 60.0)
sigma = 1.0
use = [0, 1, 2, 3, 4, 5, 6]          # six ground beacons and a tethered balloon

def geometry(idx):
    H = []
    for i in idx:
        b = beacons[i]; r = math.dist(truth, b)
        H.append([(truth[0] - b[0]) / r, (truth[1] - b[1]) / r, (truth[2] - b[2]) / r, 1.0])
    return H
H = geometry(use)
n = len(H)
Ht = g.mat_T(H)
Q = g.mat_inv(g.mat_mul(Ht, H))
A = g.mat_mul(Q, Ht)                                  # position and clock from pseudoranges: x = A z
Sproj = g.mat_sub(g.mat_eye(n), g.mat_mul(H, A))      # residual projector S = I - H A, rank n - 4
dof = n - 4
print(f'{n} pseudoranges, 4 unknowns (3 position, 1 clock): {dof} degrees of freedom in the residual')

# detection: the sum of squared residuals over sigma^2 is chi-square with dof degrees of freedom; for dof = 3 a closed form of the tail does not exist, so use the series
def chi2_sf(x, k):
    # regularized upper incomplete gamma Q(k/2, x/2), by the series of the lower function
    a, xx = k / 2.0, x / 2.0
    term = 1.0 / a; s = term
    for j in range(1, 400):
        term *= xx / (a + j); s += term
    lower = s * math.exp(-xx + a * math.log(xx) - math.lgamma(a))
    return 1.0 - lower
def chi2_isf(p, k):
    lo, hi = 0.0, 100.0
    for _ in range(100):
        mid = (lo + hi) / 2
        if chi2_sf(mid, k) > p: lo = mid
        else: hi = mid
    return (lo + hi) / 2
def ncx2_sf(x, k, lam):
    # noncentral chi-square survival function: Poisson mixture of central ones
    tot = 0.0
    for j in range(0, 200):
        w = math.exp(-lam / 2 + j * math.log(lam / 2) - math.lgamma(j + 1)) if lam > 0 else (1.0 if j == 0 else 0.0)
        tot += w * chi2_sf(x, k + 2 * j)
    return tot
pfa, pmd = 1e-3, 1e-3
thr = chi2_isf(pfa, dof)
print(f'threshold for a false-alert probability {pfa}: SSE / sigma^2 > {thr:.3f}')
lo, hi = 0.0, 200.0
for _ in range(80):
    mid = (lo + hi) / 2
    if ncx2_sf(thr, dof, mid) < 1 - pmd: lo = mid           # detection probability 1 - pmd needs the survival function above 1 - pmd
    else: hi = mid
lam = (lo + hi) / 2
print(f'noncentrality needed to detect with probability {1 - pmd}: lambda = {lam:.2f}')

print('beacon   residual weight S_ii   smallest detectable bias (m)   horizontal slope   position error from that bias (m)')
rows = []
for i in range(n):
    sii = Sproj[i][i]
    bmin = sigma * math.sqrt(lam / sii)
    slope = math.sqrt(A[0][i] ** 2 + A[1][i] ** 2)
    rows.append((i, sii, bmin, slope, slope * bmin, abs(A[2][i]) * bmin))
    print(f'  {i}       {sii:8.3f}             {bmin:10.2f}                {slope:8.3f}           {slope * bmin:8.2f} horizontal, {abs(A[2][i]) * bmin:8.2f} vertical')
checkable = [r for r in rows if r[1] > 1e-6]               # a beacon with S_ii = 0 leaves no residual: its fault cannot be seen at all
for r in rows:
    if r[1] <= 1e-6:
        print(f'beacon {r[0]} has S_ii = 0: it is the only beacon that fixes the height, a fault there leaves no residual and cannot be detected')
worst = max(checkable, key=lambda r: r[4])
k_noise = 3.29                                              # two-sided 1e-3 quantile of the normal distribution
noise_h = k_noise * sigma * math.sqrt(max(Q[0][0], Q[1][1]))
print(f'horizontal protection level over the detectable beacons: bias term {worst[4]:.2f} m (beacon {worst[0]}) + noise term {noise_h:.2f} m = {worst[4] + noise_h:.2f} m')
print(f'for comparison the 1-sigma horizontal noise error is {math.sqrt(sigma ** 2 * (Q[0][0] + Q[1][1])):.2f} m, and the 1-sigma vertical {math.sqrt(sigma ** 2 * Q[2][2]):.2f} m')
# Monte Carlo of detection
rng = random.Random(2)
def detect_rate(i, bias, trials=3000):
    hits = 0
    for _ in range(trials):
        z = [rng.gauss(0, sigma) for _ in range(n)]
        z[i] += bias
        r = g.mat_vec(Sproj, z)
        if sum(c * c for c in r) / sigma ** 2 > thr: hits += 1
    return hits / trials
print('Monte Carlo of the detection probability with a bias on the worst beacon:')
for bias in (0.5 * worst[2], worst[2], 2 * worst[2]):
    print(f'  bias {bias:7.2f} m: detected in {detect_rate(worst[0], bias) * 100:5.1f}% of 3000 trials')
print(f'  no fault: false alerts in {detect_rate(0, 0.0) * 100:.2f}% of 3000 trials (target {pfa * 100:.1f}%)')

7 pseudoranges, 4 unknowns (3 position, 1 clock): 3 degrees of freedom in the residual
threshold for a false-alert probability 0.001: SSE / sigma^2 > 16.266
noncentrality needed to detect with probability 0.999: lambda = 48.10
beacon   residual weight S_ii   smallest detectable bias (m)   horizontal slope   position error from that bias (m)
  0          0.433                  10.53                   0.395               4.16 horizontal,     4.76 vertical
  1          0.530                   9.53                   0.321               3.06 horizontal,     0.00 vertical
  2          0.467                  10.15                   0.378               3.83 horizontal,     0.63 vertical
  3          0.496                   9.85                   0.337               3.32 horizontal,     3.39 vertical
  4          0.526                   9.56                   0.305               2.92 horizontal,     1.38 vertical
  5          0.548                   9.37                   0.296               2.78 horizontal,     1.95 vertical
  6          0.000                 377.11                   0.006               2.33 horizontal,   409.52 vertical
horizontal protection level over the detectable beacons: bias term 4.16 m (beacon 0) + noise term 2.06 m = 6.22 m
for comparison the 1-sigma horizontal noise error is 0.83 m, and the 1-sigma vertical 1.25 m
Monte Carlo of the detection probability with a bias on the worst beacon:
  bias    5.27 m: detected in  36.8% of 3000 trials
  bias   10.53 m: detected in  99.9% of 3000 trials
  bias   21.07 m: detected in 100.0% of 3000 trials
  no fault: false alerts in 0.10% of 3000 trials (target 0.1%)

In robotics

Integrity is the difference between a navigation sensor and a navigation system certified to land aircraft.

• Receiver autonomous integrity monitoring. Aviation GNSS receivers run exactly this test on their pseudoranges and compute protection levels against an alert limit, which for approach operations is tens of metres horizontally and a few tens of metres vertically. If the protection level exceeds the limit, the receiver warns the pilot.
• Fault isolation. Detection says something is wrong; to continue, the faulty measurement must be identified and removed. That needs one more redundant measurement (here 8 beacons), and each exclusion costs a degree of freedom.
• Filters and solution separation. A filter has no simple redundancy of the form above. The usual method is to run several filters, each with one measurement excluded, and to compare their solutions with the full one: the separation is the detector and its bound the protection level.
• Design with integrity in mind. The geometry of the last lesson now has a second requirement: not only a low DOP but every important measurement with a nonzero residual weight. A balloon beacon needs a twin.
• Quality of the noise model. The threshold assumes the assumed $\sigma$ is correct. If the real noise is twice as large the false-alert rate is far above 0.1%, and the model of the noise is as important as the test.

The lab M7-L3: Ride Through the Outage asks for the outage and the faulty beacon. The innovation tests that this lesson builds on are in UAV M3-T13: Innovation monitoring.

Common misconceptions

Key takeaways

1. Integrity is a bound on the error with a stated probability, including undetected faults; accuracy is the typical error.
2. The residual of a least-squares fix is $\mathbf{r} = \mathbf{S}\mathbf{z}$ with $\mathbf{S}$ a projector of rank n − 4; $\mathbf{r}^{\mathsf{T}}\mathbf{r}/\sigma^2$ is chi-square with n − 4 degrees of freedom.
3. The threshold is set by the false-alert probability: 16.27 for 3 degrees of freedom and 10⁻³.
4. The smallest detectable bias is $\sigma\sqrt{\lambda/S_{ii}}$: 9.4 to 10.5 m for Hopper's ground beacons, giving 2.8 to 4.2 m of horizontal error.
5. The horizontal protection level is 6.2 m, against a 0.83 m 1σ noise error.
6. The balloon beacon has residual weight 0: its faults are invisible and cause up to 410 m of vertical error.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• P. D. Groves, Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, 2nd edition, Artech House, 2013, chapter 17 on integrity monitoring. Paid.
• R. G. Brown, "A baseline GPS RAIM scheme and a note on the equivalence of three RAIM methods", Navigation 39(3), 1992. Paid.
• RTCA DO-229, Minimum Operational Performance Standards for GPS/WAAS Airborne Equipment, the definition of protection levels used in aviation. Paid (RTCA); summaries are free on the FAA website.
• ESA Navipedia, RAIM , a free summary of receiver autonomous integrity monitoring. Free to read; the licence is stated on the site.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Before
  M7-T5: The Aiding Sensors of Aircraft and Spacecraft

  Give the measurement equation and dominant error of a barometric and a radar altimeter, an air-data system, a magnetometer, a sun sensor, a star tracker and terrain-relative navigation.

• Lab
  M7-L3: Ride Through the Outage

  The fix is lost for 60 s and then returns with one faulty beacon. Complete the measurement update of a loosely and a tightly coupled filter, compare their drift, and reject the fault.

Other courses on this site treat neighbouring topics from their own angle:

• UAV M2-T1: GNSS fundamentals
• UAV M2-T3: GNSS receivers and positioning
• UAV M2-T4: GNSS errors and multipath
• UAV M2-T6: GNSS jamming and spoofing
• UAV M3-T6: Error-State Kalman Filter
• UAV M3-T5: Extended Kalman Filter
• UAV M3-T11: Multi-sensor fusion architectures
• UAV M3-T13: Innovation monitoring
• UAV M3-T16: Sensor fault detection and isolation
• UAV M6-T5: PX4 state estimation and EKF2
• UAV M2-T11: Visual-Inertial Odometry
• SLAM M8-T5: GPS-Denied Localization
• Linear Algebra M5-T6: Weighted Least Squares
• Linear Algebra M6-T5: Singular Value Decomposition
• Probability for Robotics (modules 7, 8 and 12: Kalman, extended and error-state filters)

This course leaves the following to the pages that already teach them:

• Constellations, signals, error sources, multipath, jamming and spoofing
• The Kalman, extended and error-state filter derivations
• Visual, LiDAR and radar odometry, VIO and SLAM as aiding sources
• Innovation tests (NIS, NEES) as such