GB2380819A

GB2380819A - navigation system integrity monitoring and fault detection process

Info

Publication number: GB2380819A
Application number: GB0124144A
Authority: GB
Inventors: John Demetrios Loizou
Original assignee: VEGA GROUP PLC
Current assignee: VEGA GROUP PLC
Priority date: 2001-10-08
Filing date: 2001-10-08
Publication date: 2003-04-16
Also published as: GB0124144D0

Abstract

An integrity monitoring and fault detection process for navigation systems (eg GPS) which includes a weighted Total Least Squares (TLS) analysis of an over-determined positioning problem, used to generate a test statistic indicating the self consistency of the linear model used in the TLS analysis. In addition, for a given over-determined positioning problem a fault detection threshold is determined for the test statistic, as a function of a specified probability of false alarm. Furthermore, appropriate weightings may be determined to be used in Total Least Squares solutions of over-determined positioning problems, such that emphasis can be given to the required dimensions to be protected, such as Vertical Protection Limits. Also, for a given over-determined positioning problem, a limit on position errors such as a Vertical Protection Limit can be determined, as a function of specified false alarm and missed detection probabilities.

Description

NAVIGATION SYSTEM INTEGRITY MONITORING AND FAULT DETECTION PROCESS This invention relates to a process for the provision of navigation system integrity by detection of outliers in a set of range measurements.

Navigation system integrity refers to the ability of the system to provide timely warnings to users when the system should not be used for navigation. In systems such as the Global Positioning System (GPS), satellites broadcast their positions (as a function of time) along with ranging signals which are processed by the receiver to estimate the user's position. Where more than the minimum number of measurements necessary to provide a position solution are available to the user the algebraic problem to be solved is said to be"over-determined", and outliers can be detected by evaluating the self-consistency of the measurements. An outlier is a measurement whose value stands out from other measurements and predictions to the extent that it arouses suspicions that it was created by a different mechanism from the main set of measurements, such as a fault.

In general, in an over-determined system there is no unique solution which provides an estimate of the user position exactly matching all of the measured ranges. A linear model is generally employed and the least squares solution is obtained to find the best estimate of the user's position. For each measured range there will be some difference between the actual measurement and the distance between the estimated user position and the measurement predicted by the linear model. This difference is known as the"residual"for that measurement. By comparing each residual, or some scalar test statistic derived from the residuals, with some threshold value, outliers can be detected. Such techniques are known generically as Receiver Autonomous Integrity Monitoring (RAIM) methods. The principles of a number of RAIM methods are described in Chapter 5 of"Global Positioning System: Theory and Applications Volume II", edited by B. W. Parkinson and J. J. Spilker Jr (ISBN 1-56347-107-8), hereby incorporated by reference.

A fundamental element of most RAIM methods is that the magnitude of the detected faulty measurement should be proportional to the position error as well as to the value of the test statistic. Thus, in the event of a fault on one of the measurements, there is a relationship between the magnitude of the position error and the value of the test statistic.

When specifying system integrity requirements, two key parameters are the Probability of Missed Detection (P MD) and the Probability of False Alarm (P FA). The magnitude of the test statistic threshold used to declare an alarm from a set of measurements is typically a function of P FA; false alarms are undesirable in an integrity monitoring system, and as

the false alarm specification becomes more stringent (i. e. as P FA reduces) the test statistic threshold increases. The magnitude of the position error protected by the integrity monitoring system has two components which are functions of P FA and P MD respectively. The first component is the magnitude of the position error proportional to the required test statistic threshold. Added to this is the second component, which is a function of the allowable probability of a missed detection (given a fault on a ranging measurement) and the probability density function of position errors.

The sum of these components provides a position protection limit for the integrity system.

Where three-dimensional position integrity is required (for example, for aviation applications), a Horizontal Alarm Limit (HAL) and Vertical Alarm Limit (VAL) are typically specified as operational requirements. The integrity system therefore has two distinct functions: to generate alarms in the event of detected outliers, and to generate alarms in the event that integrity is not available to meet the specified protection limits.

In navigation systems such as GPS, there may be errors in both the linear model used to solve the positioning problem (for example, due to errors in the broadcast position of the satellites) as well as in the measured ranges. For RAIM methods where the test statistic is derived solely from the self-consistency of the vector of range residuals the available position protection limit is often large relative to the actual position errors caused by erroneous range measurements.

An object of the present invention is to provide an integrity monitoring process for navigation systems.

A further object of the present invention is to provide a fault detection process for navigation systems.

Another object of the present invention is to provide a process for the determination of appropriate weightings to be used in Total Least Squares solutions of over-determined positioning problems, such that emphasis can be given to required dimensions to be protected, such as Vertical Protection Limits (VPL) or Horizontal Protection Limits (HPL).

Yet another object of the present invention is to provide a process wherein for a given over-determined positioning problem a limit on position errors such as a Vertical Protection Limit or Horizontal Protection Limit can be determined.

According to the present invention a process, including a weighted Total Least Squares (TLS) analysis of an over-determined positioning problem, is used to generate a test statistic indicating the self consistency of the linear model used in the TLS analysis, formed from the ratio of the norms of mismatches of the observation model and measurement data.

In addition, for a given over-determined positioning problem a fault detection threshold is determined for the test statistic, as a function of a specified Probability of False Alarm, by determination and application of a residual vector comprising a critical bias on one element only, with all other elements set to zero.

Furthermore, appropriate weightings may be determined to be used in Total Least Squares solutions of over-determined positioning problems, which reduces the dispersion of the slopes on a chart of notional position error against test statistic, such that emphasis can be given to required dimensions to be protected, such as Vertical Protection Limits or Horizontal Protection Limits.

Also, for a given over-determined positioning problem, a limit on position errors such as a Vertical Protection Limit or Horizontal Protection Limit can be determined, as a function of a specified Probability of False Alarm (P FA) and Probability of Missed Detection (P MD), by determination and application of a residual vector that causes the largest position error for a given test statistic threshold.

All or part of this process may be embodied as re-programmable or embedded software within navigation equipment, or in any other manner as appropriate to the specific application. Although the integrity monitoring and fault detection process would typically be undertaken in real-time to meet specified integrity requirements, parts of the process such as the evaluation of appropriate weighting matrices as discussed in the subsequent example may be performed off-line, with the results of this analysis used in the real-time integrity monitoring and fault detection process.

A specific embodiment of the invention will now be described by way of example.

The steps of the navigation system integrity monitoring and fault detection process, as applied to an implementation for aircraft positioning to protect vertical position such as during precision or non-precision approach, using measurements from navigation satellites such as the GPS constellation, are as follows:

1. Produce the observation matrix G. G is an n x 4 matrix, where n is the number of ranging measurements available (in the case of a satellite navigation system, n IS the number of satellites in view). As is well known, the first three columns of the observation matrix consists of a set of unit vectors from the presumed user position (about which the observation matrix is linearised) to the assumed position of each satellite, and the fourth column contains "1", signifying the existence of a clock offset.

This is a common step in satellite navigation.

2. For each observation i (where i equals 1 to n), estimate the standard deviation of the error, s (i), and the variance, s (i) 2. In a satellite navigation system, s (i) corresponds to the assumed User Equivalent Range Error (UERE). Typically, UERE will be specified as a function of elevation angle, therefore this step will usually involve calculating the angle of elevation of each satellite in view, and estimating s (i) as a function of this angle, by calculation, from a look-up table or by some other means. This is a common step in satellite navigation.

3. Produce the weighting matrix W, which is the inverse of the covariance matrix. For simplification, and in the absence of any information to the contrary, it may be assumed that the error sources on each observation are uncorrelated with each other. Therefore, in a satellite navigation system, the diagonal elements are the inverses of the variances corresponding to each satellite, with all off-diagonal terms set to zero. This is a common step in satellite navigation.

4. Produce the 4 x n matrix K, the weighted pseudo-inverse (also known as the Moore- Penrose inverse) of G. This is a common step in satellite navigation.

5. Produce an n x n matrix P that is the linear transformation that takes the range measurements into the resulting residual vector: P=GK This is a common step in satellite navigation RAIM methods.

6. For each observation i (where i equals 1 to n), using standard Least Squares Residuals (LSR) RAIM methods, evaluate the slope of positioning error against test statistic. For this satellite navigation aircraft positioning example, where errors in the vertical dimension are to be protected, this slope for each satellite i is represented by the term Vslope (i) and is evaluated by:

Where: IK (3i) 1 is the absolute value of the third row, column i, of the K matrix.

P (ii) is the value of row i, column i, of the P matrix.

If horizontal position rather than vertical position is to be protected, the parameter to be evaluated is Hslope (i), in which case IK (3i) 1 above should be replaced by [K (1 i) 2 + K (2i) 2]1/2 7. Find the observation with the greatest value of this slope. In this example, this is henceforth referred to as Sat (crit), and the value associated with it is Vslope (max).

Find also the standard error associated with this satellite, referred to as s (crit), and the diagonal element on the P matrix associated with this satellite, referred to as P (crit). This is a common step in satellite navigation RAIM methods.

8. For the specified Probability of False Alarm (P FA), and for the number of observations available, n, evaluate the standard Least Squares Residuals RAIM method test statistic T (n, P FA). This may be approximated by treating the statistic as a chi-square distributed variable with n-4 degrees of freedom, or more precisely by numerical integration. This is a common step in satellite navigation RAIM methods.

9. For Sat (crit), evaluate the amount of bias which, when applied to this satellite with no noise or bias on any other satellite, will generate the LSR test statistic. This is referred to henceforth as the "Critical Bias" :

10. Produce an n x 1 vector q with all elements equal to zero except the element corresponding to Sat (crit), which is set to the value of the Critical Bias.

11. Produce a new matrix [G q], by appending the q matrix to the G matrix to form an n x 5 matrix.

12. Produce an n x n matrix C that presents the weight or significance of each element in the set of observations. Depending upon the specific application this may be an identity matrix, the inverse of the W matrix, or some other expression appropriate to the specific application. For this satellite navigation example, the preferred, but not

mandated, expression for the C matrix IS a diagonal matnx with diagonal elements equal to 1/s (i).

13. Produce a 5 x 5 matrix D whose elements are inversely proportional to the relative weight or significance of each column of the augmented matrix, as described in steps 38 to 67. In the preferred embodiment of this invention the D matrix IS pre-set into the integrity monitoring equipment. However if sufficient processing capability is available the D matrix may be optimised in real-time to suit the available observation matrix and standard errors of the available measurements.

14. Multiply together the matrices described in Steps 11,12 and 13 to produce a new n x 5 matrix C [G q] D.

15. Perform a Singular Value Decomposition (SVD) of the matrix C [G q] D. This is ; j well known technique in matrix algebra. The output from the SVD process are three new matrices, U, S and V such that

In which U is an n x 5 orthogonal matrix, V is a 5 x 5 orthogonal matrix and S is a 5 x 5 diagonal matrix, whose non-zero elements are henceforth referred to as the "singular values"of the decomposition. The smallest of these singular values (normally the fifth diagonal element) is henceforth referred to as svs.

16. Partition the U, V, S and D matrices such that:

17. Produce a new augmented matrix [Htls etis], in which Hus is an n x 4 matrix and ells is an n x 1 vector. The augmented matrix is formed from the matrix partitions described at step 16:

The matrix Hus represents the errors in the observation matrix G, and eus represents the errors in the q matrix, when the linear matrix equation q = (G + H). x + e is solved using the Total Least Squares method, hence the use of the"tls"subscripts, to signify Total Least Squares.

18. Evaluate the Frobenius norm (i. e. the square root of the sum of the squares of all the elements) of Hus. This is referred to henceforth as H, IIF- 19. Evaluate the vector norm (i. e. the square root of the sum of the squares of all the elements) of ets. This is referred to henceforth as #etls#.

20. Evaluate the ratio Htis p/ ! ! es i !. This is the Test Statistic Threshold, for the observation geometry defined by G, for a false alarm probability of P FA. This threshold is henceforth referred to as TSTos.

------------------------------------------------------------------------------------------------------------------ 21. Produce the vector of measurement residues, y. y is an n x 1 vector. As is well known, the y matrix contains the raw range measurements minus the expected ranging values based on the location of the satellites and the presumed user position.

This is a common step in satellite navigation.

22. Produce the least squares solution for the four-dimensional position vector (north, east, up and clock). This is a common step in satellite navigation.

23. Output the vertical position of the receiver, Vpos. This IS a common step in satellite navigation 24. Produce a new matrix [G y], by appending the y matrix to the G matrix to form an n x 5 matrix.

25. Multiply together the matnces described in Steps 24, 12 and 13 to produce a new n x 5 matrix C [G y] D.

26. Perform a Singular Value Decomposition (SVD) of the matnx C [G y] D, and partition the output matrices as at Step 16.

27. Repeat Steps 17,18 and 19 above.

28. Evaluate the ratio Htis !) p/ !) e, !. This is the Test Statistic for this set of observations, henceforth referred to as TS.

29. If TS is greater than TSTs, an unacceptable inconsistency in the ranging information has been detected, for which the navigation system should generate a "fault detected"alarm.

30. Produce a new 4 x n matrix L, analogous to the K matrix produced at Step 22, but without using the weighting matrix W: L = (GTG) "1GT 31. Produce a new n x 1 vector b, in which each element is the product of the corresponding element in the third row of L (i. e. L (3i), for the vertical protection case being described) and some arbitrary positive number, m: b (i) = L (3i). m Note that the sign of L (3i) is relevant in this step. The absolute value (as used in earlier steps) must not be used.

32. Produce a new matrix [G b], by appending the b matrix to the G matrix to form an n x 5 matrix.

33. Perform a Singular Value Decomposition (SVD) of the matrix C [G b] D, and partition the output matrices as at Step 16.

34. Repeat Steps 17,18 and 19 above.

35. Evaluate the ratio H ! ! p/ ! ! es ! !.

36. Iterate Steps 31 to 35, changing the value of m until the value of TS found at Step 35 is equal to the test statistic threshold TSTus, to within specified limits, using any appropriate numerical method to solve this iterative problem.

37. Using the final (worst case) value of b, evaluate the worst case least squares solution for the four-dimensional position vector (north, east, up and clock).

38. Evaluate the worst case vertical position error, Vpos (wc). This is the third element of the vector xwc, Xwc (3)' 39. Evaluate Vp, the difference between the nominal vertical position error and the worst

case vertical position error, for a test statistic equal to TSTtls :

Vp = Xwc (3) - [T (n, P FA) x Vslope (max)]

40. Produce a new vector a, in which each element a (i) is the absolute value of the corresponding element b (i), divided by the standard error for that observation, s (i).

The elements of a are equal to the number of standard deviations on each observation required to form the vector b.

41. Produce a new vector ap in which each element equals the cumulative probability associated with the corresponding element of a, assuming a normal distribution, as found on a table of normal distribution curves. Using standard statistical functions for

example on Microsoft Excel, this step becomes :

ap (i) = 1-normsdist (a (i))

42. Replace the element in ap corresponding to Sat (crit) with the value 1.

43. Evaluate the product of all elements in ap and multiply the total by two. This corresponds to the probability of occurrence of the worst case undetected vertical position error, given the condition of a bias on Sat (crit) (the factor of two represents the case of vector b with all signs reversed).

44. Using any appropriate standard technique, evaluate the inverse of the probability found at Step 43. This value, henceforth referred to as k (MD max), is the number of standard deviations corresponding to the probability of occurrence of the worst case undetected vertical position error.

45. Divide Vp evaluated at Step 39 by k (MD max). This value, henceforth referred to as s (tis), is an estimate of the undetected vertical position error, with a probability equivalent to one standard deviation.

46. Evaluate k (P MD), the number of standard deviations associated with the conditional missed detection probability P MD. This is a common step in satellite navigation RAIM methods.

47. Evaluate the available Vertical Protection Limit, VPL, as follows :

VPL = [T (N, P FA) X Vs1ope (maxi + [s (tls) x k (P MD)]

48. If VPL is greater than the specified Vertical Alert Limit requirement, the geometry of the ranging sources relative to the user position is unacceptable, given the specified integrity requirements, in which case the navigation system should generate a"RAIM VAL Unavailable"alarm.

49. The following steps define the process required to set the D matrix as appropriate to the specific application. This process may be undertaken in real-time or near real- time, to modify the D matrix as appropriate to the current observation geometry; it may be performed off-line once only to set a fixed D-matrix to be used for all geometries; or it may be performed with any other frequency as appropriate to the specific application.

50. Create a provisional D matrix as a 5 x 5 Identity matrix (i. e. all elements equal zero, except the leading diagonal elements which are equal to one).

51. Identify the order of magnitude of the Critical Bias (as defined in Step 9) associated with the specific application. If the D matrix is being defined in real-time then the Critical Bias is as calculated at Step 9; otherwise the value to be used may be found by simulation or other form of analysis.

52. Produce a G matrix, as at Step 1. If the D matrix is being defined in real-time then the G matrix is as calculated at Step 1; otherwise representative matrices may be denved by simulation or other form of analysis.

53. Produce a C matrix, as at Step 12. If the D matrix is being defined in real-time then the C matrix is as calculated at Step 12 ; otherwise representative matrices may be derived by simulation or other form of analysis.

54. Set the fifth diagonal element of the D matrix to be the reciprocal of the order of magnitude found in Step 51. For example, if the Critical Bias is expected to take values typically around 1 Om, and the values in the y vector are given in metres, set D (5,5) equal to 0.1. For satellite navigation applications this value will typically be between 0.1 and 0.01, although other values are possible.

55. The remaining four elements of the D matrix allow the process to be modified to emphasise horizontal positioning error, vertical positioning error or clock bias. For the example of vertical position protection, only element D (3, 3) needs to be modified.

56. Set D (3,3) to some arbitrary value between zero and one. Typically, 0.5 will be a reasonable first estimate.

57. Produce an n x 1 vector q with all elements equal to zero except the first element, which is set to some small bias value.

58. Multiply the matrices together to form the matrix C [G q] D.

59. Perform the steps previously described at Steps 15 to 19.

60. Evaluate the ratio) ! H !) p/ en, !).

61. Using the q vector, evaluate the notional least squares solution for the fourdimensional position vector (north, east, up and clock).

Xnotional = (GTWG)-W. q = K. q

62. Evaluate the notional vertical position error, Vpos (notional). This is the third element of the vector Xnotiona),XnotionatO) 63. Plot the point ) Hs !) p/ ! t ens , Vpos (notional)] on a chart of notional vertical position error against test statistic.

64. Iterate Steps 57 to 63, with Increasing values of bias to cover the expected range of possible bias values, to produce a line graph on the chart of notional vertical position error against test statistic.

65. Repeat Steps 57 to 64, applying the bias sequentially to observations 2 to n in turn, to produce a family of lines on the chart of notional vertical position error against test statistic.

66. The objective is to find a value for D (3,3) which shows the smallest dispersion in the slope of the lines shown in the chart of notional vertical position error against test statistic. The closer together are these lines, the smaller will be the value of Vp as found at Step 39, and hence the lower the overall Vertical Protection Limit. This optimisation may be performed by iteration of Steps 56 to 65 to find a value for D (3,3) which meets the requirements of the specific application, or by any other appropriate method. Typically 0 (3,3) may take values from about 0.3 to about 0.7.

67. If, over the range of bias values appropriate to the application, the lines on a chart of notional vertical position error against test statistic are not approximately linear, then D matrix is not correctly formed. Where appropriate the fifth element of the D matrix may again be modified, to further linearise and reduce the dispersion in the slope of the lines shown in the chart of notional vertical position error against test statistic.

68. In practice, this process may be automated without actually producing a chart of notional vertical position error against test statistic. The objective of finding values for the elements of the D matrix that generate linear characteristics of vertical error against test statistic, and minimal dispersion in the value of the slopes of these lines, remains valid.

Claims

CLAIMS 1. A process including a weighted Total Least Squares (TLS) analysis of an over- determined positioning problem, used to generate a test statistic indicating the self consistency of the linear model used in the TLS analysis, formed from the ratio of the norms of mismatches of the observation model and measurement data.
2. A process as claimed in Claim 1 wherein for a given over-determined positioning problem a fault detection threshold is determined for the test statistic, as a function of a specified Probability of False Alarm (P FA), by determination and application of a residual vector comprising a critical bias on one element only, with all other elements set to zero.
3. A process as claimed in Claim 1 wherein appropriate weightings may be determined to be used in Total Least Squares solutions of over-determined positioning problems, which reduces the dispersion of the slopes on a chart of notional position error against test statistic, such that emphasis can be given to required dimensions to be protected, such as Vertical Protection Limits (VPL) or Horizontal Protection Limits (HPL).
4. A process as claimed in Claim 1 or Claim 2 wherein for a given over-determined positioning problem a limit on position errors such as a Vertical Protection Limit or Horizontal Protection Limit can be determined, as a function of a specified Probability of False Alarm (P FA) and Probability of Missed Detection (P MD), by determination and application of a residual vector that causes the largest position error for a given test statistic threshold.