Application of statistical fault detection algorithms to navigation systems monitoring

Automatica, Vol. 29, No. 5, pp. 1275-1290, 1993 Printed in Great Britain.

0005-1098193 $6.00 +0.00 ~) 1993 Pergamon Press Ltd

Application of Statistical Fault Detection Algorithms to Navigation Systems Monitoring* I. NIKIFOROV,t V. VARAVVA~: and V. KIREICHIKOV~t

The quickest detection of soft drifting-type faults in aircraft inertial and inertial~radio navigation systems is proposed and tested by using experimental data from test flights of a new airplane. Key Words--Failure detection; monitoring; navigation; guidance systems; sensor failure; redundancy; stochastic systems.

part of this problem. System degradation should be detected as soon as possible when it leads to an unacceptable growth of the output errors. A system degradation can be due, for instance, to an accumulation of sensor errors, or changes in the system operating conditions. In this paper we concentrate on the sensor fault detection problem. Mathematically this problem can be formulated as that of the quickest detection of abrupt changes in a stochastic process. It is intuitively obvious that the criterion to be used should favor fast detection with few false alarms. Fast detection is necessary because, between the fault onset time and its detection, abnormal measurements are taken in the navigation systems, which is clearly very undesirable. On the other hand, false alarms result in lower accuracy of the estimates because some correct information is not used. The optimal solution involves a tradeoff between these two contradictory requirements. The main goal of this paper is to demonstrate the capability of the modified cumulative sum (CUSUM) algorithm for the detection of soft sensor faults and to illustrate its performance for a significant application. The paper is organized as follows. First, we introduce three types of models. The first basic model is used for the design of the change detection algorithm. We also investigate the statistical properties of the algorithm based on this model. The second (A) and third (B) models are used for the detection of sensor faults. Model (A) is a state-space model. Model (B) is a regression model with redundancy. In these two cases we show how the new change detection problem can be reduced to the basic problem and we discuss some new features which play a key role in models (A) and (B).

Airliner--The fault detection problem for a navigation system is investigated. Mathematically this problem is formulated as that of the quickest detection of abrupt changes in a stochastic process. For solving this problem the modified cumulative sum algorithm (CUSUM) is proposed. The ability of the proposed algorithm to detect soft (drifting-type) faults and its application for the monitoring of the inertial and the inertial/radio navigation systems are demonstrated. Test results describing the application of CUSUM algorithms to experimental data from test flights of a new airplane are reported. CUSUM algorithm INS GLR algorithm GMA algorithm ONS RMS SIRU SPRT

NOMENCLATURE cumulative sum algorithm, inertial navigation system, generalized likelihood ratio algorithm, geometric moving average algorithm, Omega navigation system, root-mean-square (standard deviation), strapdown inertial reference unit, sequential probability ratio test.

1. INTRODUCTION

NAVIGATION SYSTEMSare standard equipment for moving vehicles such as planes and boats. Fault detection is one of the main problems in the design of modern navigation systems. Detection of soft drifting-type faults is the most difficult

* Received 2 December 1991; revised 5 November 1992; received in final form 1 February 1993. The original version of this paper was presented at the IFAC Symposium SAFEPROCESS'91, Baden-Baden, Germany, September 1991. The Published Proceedings of this IFAC Meeting may be ordered from: Pergamon Press Ltd, Headington Hill Hall, Oxford OX3 0BW, U.K. This paper was recommended for publication by Guest Editors Albert Benveniste and Karl Johan /~.str0m. Corresponding author: I. Nikiforov. Tel. (33)99847237; FAX (33)99383832; e-mail [email protected]. t Institute of Control Sciences, Profsoyuznaya St, 65, Moscow, 117342 Russia. I. Nikiforov is now with IRISA-INRIA, Campus de Beaulieu, 35042 Rennes Cedex, France. Institute of Aviation Equipment, Zhukovsky, Moscow region, 140160 Russia. 1275

1276

I. NIKIFOROV et al.

Next, we discuss a simple but representative example of model (A). We compare three change detection algorithms [a modified CUSUM algorithm and two other candidates: the generalized likelihood ratio (GLR) algorithm, and a geometric moving average algorithm (GMA), which is based on Kalman filtering]. Comparison relies on the quickest detection optimality criterion. The main features of these algorithms are the following: the modified CUSUM algorithm is based on some prior information about the magnitude of change and possesses asymptotic optimality according to the quickest detection criterion. The GLR algorithm is based on the assumption that the magnitude of change is unknown; it has a uniformly good performance for a wide range of possible magnitudes. On the other hand this algorithm is very time-consuming. The GMA algorithm relies on the assumption that the magnitude of change is unknown and possesses no optimality property. On the other hand the GMA algorithm is a simple and popular fault detection algorithm, which is frequently used for navigation systems monitoring. Finally, we apply the proposed change detection algorithm to navigation system monitoring. First, we discuss the general methodology, which includes the design of fault detection algorithms, as well as the implementation and tuning of these algorithms. We then discuss three particular cases of the general navigation system monitoring problem: • triple inertial navigation system monitoring, • skewed axes redundant strapdown inertial reference unit monitoring, • joint inertial/radio navigation systems monitoring. Main results of the testing of the proposed triple inertial navigation system monitoring scheme on experimental data from test flights of the IL-96-300 airplane are also reported.

where to is the unknown change time. We assume that the vector 00, the matrix ~ and the Kullback information p(0l, 0o) = ½(0o - 0 I ) T ~ - I ( 0 0 -

2.1. Models This paper is concerned with three types of models. Basic model. Let (Yt),~_l be a m-dimensional independent Gaussian random sequence, where ~(Yt) =-At(0, ~E), 0 is the mean vector and ~ is the covariance matrix of the Gaussian distribution. Under the hypothesis of a change, the model of the observation can be written in the following manner

~(Yt)

[N(00,]E) IN(01, X)

if t < t o if t-->to

(1)

(2)

are known a priori. In other words the 'magnitude' p of the change is known, but not its direction. Model A. We consider a linear dynamic system described by the state-space model for the observed signal (Yt)t_>l Xt+ = ~ ( t + 1, t)Xt + ~t Y, = H(/)Xt + vt

(3)

where Xt e R n is the state vector, Yt E R m is the measurement, (g,)t_~0 and (%),--1 are two independent zero mean Gaussian white noises with covariance matrices Q(t) -> 0 and W(t) > 0 respectively. The initial state Xo is a Gaussian zero mean vector with covariance matrix Po > 0. The matrices ~ , H, Q, W, Po are known. The change is modeled as additional bias Ato in the plant equation Xt+ 1 =

O(t + 1, t)Xt + gt + At,

or in the measurement equation Yt = H(t)X, + vt + At, where A,={O

if t < t o if t-> to.

We assume that the additive change A is unknown, but it is assumed that the Kullback information between the two densities P0 (no change) and Pa (a change occurs) for the signal (Yt),_-I can be approximated by an a priori known value. Kullback information will be defined in subsection 3.2. Model B. We consider a linear static system with redundancy described by the regression representation Y, = HX, + v, + Ato

2. M O D E L S A N D C R I T E R I A

01)

(4)

of the observed signal (Yt)t_>t. In this case X, e R" is the unknown state, (vt),_~ is a white noise with covariance W = O'21m, where Im is the identity matrix m × m. H is a constant full rank matrix of size m × n with m > n. The vector Ato is defined as in model A. The characteristic feature of this model is the existence of redundancy (m > n) in the information contained in the observations (Yt)t->l.

Discussion: The key tool for detecting additive changes in model A consists of the transformation of the observations (Y,),___1 into the innovation sequence of the Kalman filter, which is an independent Gaussian random sequence.

Statistical fault detection algorithms Analogously, for detecting additive changes in model B we transform of the observations (Y,)t~l into the residual sequence of the least square algorithm, which is again an independent Gaussian random sequence. Therefore, the change detection problems for models A and B are reduced to the change detection problem for the basic model. This fact explains a key role of this model. 2.2. Criteria Let (Yt)t>_l be a random sequence with conditional density Po(Y, lYe-l), where 0 is a parameter vector. Until the unknown time to the vector 0 is 0 = 0o and from to becomes 0 = 01. The problem is to detect the change in 0 as soon as possible. The change detection algorithm has to compute the alarm time based on the observations YI . . . . , Y , The alarm time (stopping time in the mathematical literature) t, is the time at which the change is detected. In this case, two situations can occur. If the change is detected after the time to (t, -> to is true) then this detection is correct and the detection delay is

On the contrary, if the change in 0 is detected before the time to (ta < to is true) it corresponds to a false alarm. It is intuitively obvious that the criterion which must be used favors fast detection with few false alarms. In other words t,---to should be stochastically small and t, t, < to should be stochastically large. Therefore, the optimality criterion for the quickest detection algorithm consists in minimizing the mean detection delay (5)

where ~ ( I ) is the conditional expectation, for a given mean time i~ between false alarms 7" = ~(ta I t~ < t0).

according to Po,(Y). Consider the hypothesis (no change occurs) that all observations Y l , . - . , Y t are distributed according to Po.(Y). The log likelihood ratio for testing ~k against ~o is

Stk = lnpk(yl' " " "" Y') 2 lnP°'(Yi) Po(Y~, , Yt) -- i=k ~ "

(6)

The on-line detection problem was formulated in the 1950s and 1960s for a scalar independent random sequence. The optimal Bayesian change detection algotihm was investigated by Shiryaev (1961). The cumulative sum (CUSUM) algorithm was proposed by Page (1954) as a heuristic change detection algorithm. More recently, Lorden proved (Lorden, 1971) its ('first order') asymptotic optimality for the 'worst case' criterion which is discussed in Section 3. The CUSUM algorithm can be described as follows. Assume that (Yt),~I is an independent random sequence with density Po(Y,). Consider for some k : 1 -< k -< t the hypothesis ~k that the observations Yl . . . . , Yk-~ are distributed according to Poo(Y) and Yk. . . . . Yt are distributed

(7)

When the change time k is unknown, the standard statistical approach consists of estimating it by using the maximum likelihood principle which leads to the following decision function g, = max S~.

(8)

l'~--k'~--t

Therefore, the alarm time is set as the first instant ta for which the decision function gt reaches a threshold h which is selected a priori ta = min {t -> 1 :g, -> h}.

(9)

Note that the decision function gt for the CUSUM algorithm (9) can be evaluated recursively through the relation

(

l P°'(Yt)\+-- ,

gt= ~gt-l + n

T = t a - t 0 + 1.

= ~(~ I t~ --- to)

1277

go = 0,

Poo(Y,)) (x) + = max (x, 0).

(10)

A modified CUSUM for dependent random sequences in the vector case has been described previously (Nikiforov, 1978, 1983, 1986). 3. ALGORITHMS

In this section we first discuss the solution of the basic change detection problem for the basic model. Then we continue with change detection problems in more complex situations corresponding to models A and B. In these two cases we first show how the new problem can be reduced to the basic problem and next we discuss new features which play a key role in models A and B. See Basseville and Nikiforov (1993) for more extensive discussions.

3.1. Basic problem Invariant SPRT. In the case of known magnitude but unknown change direction, the change detection algorithm is based upon Wald's idea of weight function for the sequential probability ratio test (SPRT) or upon the theory of the invariant SPRT in the modern interpretation. For some families of distributions, it is useful to use the invariance properties with respect to some transformation of the observations in order to guess this weight function. It is possible to prove that the family ~c(0, |) (where | denotes the identity matrix) remains invariant under the transformation g Y = RY, where R

I. NIKIFOROV et al.

1278

such that 2; = RR x [see details in Ghosh (1970)]. In other words, we consider the problem of testing between the hypotheses

~,

= {0:0 = 0,,}

o j=k

and (11)

ff~a1 = { 0 : ( 0 - - 0 0 ) T ] ~ - - I ( 0 - - 0 0 ) = b 2 } ,

where 0o, I2 and b are known.

Design of the algorithm--scalar case. Let us start with the scalar case (m = 1). Wald's idea consists in averaging the likelihood ratio (7) with respect to all possible values of 0~, using a weighting function dF(0 0 which is some appropriate probability distribution S~ = In (~ Po,(Yk . . . . , Yt) dF(0,). J-~Po,,(Yk,

,

(12)

Y,)

In the case of a scalar Gaussian sequence with known variance tr2, the distribution F(O~) is concentrated on two points 0o - ob and 0o + orb. Therefore, the weighted likelihood ratio (12) becomes S~, = In f_[ exp { 0 ~

Therefore, the decision rule (14) consists in stopping at the first time t such that the 1 cumulative sum - ~ (y~- 00) reaches the upper curved boundary c,-k+l or the lower curved boundary --C,-k+~. When l goes to infinity, the asymptotic behaviour of the function c; is given by

c;=

~._00) 2

202

(15)

ta = min {t - 1 : (g,* -> h + In 2) U (g~ -> h + In 2)}, where

g*=(*g,-1 + (y,- 00o b)) + g;=(gt-'-(Y'-O°+b))

~ (yj-Oo)

(16)

g o = g~, = o.

(t - k + 1)} d F ( 0 , )

=ln{½exp(bj~_k(yj-Oo)-(t-k

~ 1 + O(e-2h).

Thus, in the one-dimensional case, the X 2CUSUM algorithm with threshold h is closely related to the following classical Page's twosided CUSUM algorithm with threshold h + In 2:

j=k

(0,

b

+ l) b---2)

+½exp(-b~(yj-Oo)-(t-k+l)b---2)} Oj=k

= - ( t - k + 1)-~-+ In cosh -~j~_k(yj

Algorithm design--vector case. The general case of the x2-SPRT is described in Ghosh (1970) and the x2-CUSUM is described in Nikiforov (1983, 1986) and Basseville and Nikiforov (1993). The proof of formula (12) in the vector case is complex, and we refer to these references for exact mathematical results and details b2 (m b2tZ(xtk)2] ~ ! + l ) - ~ + In G ~ ,

S~=-(t-k

Finally, we replace S~, in (8) by S~, and get the following z2-CUSUM algorithm: where G(m,x) is geometric function

t, = rain {t -> 1 :g, -> h}

{

gt = lm
(17) hyper-

X

G(m,x) = 1 +--+" • • m (13)

Xn

j=k

m(m + 1 ) . . . (m + n -

where b - [ 0 t -

001 is the signal-to-noise ratio. It o is obvious that equation (13) can be rewritten in the following manner t~ = rain t-> 1 : max

generalized

the

(y; - 00) -> C,-k+l

+..., 1)!

and (xk)2 = ('yk __ 0 0 ) T ' ~ - l ( ' ~ r k __ 00),

1

k

flki = k _ i + l j=i

l~k~t

(14) where ct is the unique positive solution of the following equation: b2 h = - 1 ~ - + In cosh (bct).

Note that S~, is the likelihood ratio for testing that the noncentrality parameter of a X2distribution with m degrees of freedom is either equal to 0 or to (t - k + 1)b z. This fact explains the name of x2-CUSUM. As in scalar case, we replace S~, in (8) by S~, and get the following

1279

Statistical fault detection algorithms

x2-CUSUMalgorithm:

Note here that this 'first-order' optimality was proven in Lorden (1971) for the scalar case, when the parameters 0o and 01 are known. The analogous result for the x2-CUSUM algorithm is based upon the following theorem [see proofs and details in Nikiforov (1993) and Basseville and Nikiforov (1993)]:

to = min { t >- l : max S~k>- h } l~k'~--t

b2

$~, = - ( t - k + 1) -~-

(m +lnG

b 2 ( t - k + 1)2(X~,)2)

~-,

4

.

(18)

First order optimality. Let us pursue our discussion of the optimality criterion. It is obvious that, without knowing a priori the distribution of the change time to, the mean detection delaty ~ defined in (5) is a function of the change time to and of the past values of the random sequence. In some practical cases it is very useful to have an algorithm which is independent from the distribution of the change time to and the sample path of the past observations Y I , - . - , Y,o-~. For this reason we use another optimality criterion (Lorden, 1971), which consists in minimizing the 'worst case' mean time delay

Theorem 3.1. Let to be the stopping time (18) of the x2-CUSUM algorithm. We assume that the vector 0o, the matrix ~, and the Kullbaek information p(01, 0o) are a priori known values, Then the 'first-order' optimality relation (20) still holds for the basic problem, where p(0t, 0o)= b2 2" It is obvious that for on-board navigation system monitoring, a simple recursive algorithm is preferable. The classical CUSUM (8)-(9) can be expressed in the recursive form (10) g, = max S~, = (g,-1 + ln PO,(Y,)~ +

l<-k<-t

~* = sup ess sup ~(z ] y~,-l, to -> to) to> l

(19)

for a given mean time between false alarms (6). Let us motivate the derivation of this 'worst' mean detection delay ~*. For a given change time to, the conditional expectation ~.(y[,,-!) = ~(z [ y~o-,, t~ -> to) is a random value. For this reason we have to use the essential supremum esssup ~(Y2'-~)t in order to have the smallest A(to) such that ~(y~0-,) _
Poo(Y,)/ "

But the x2-CUSUM algorithm (18) cannot be directly expressed in the recursive form because S~,(17) is a nonlinear function of the cumulative sum

~ (Yi-00).

i=k

But

the

effect

of

this

nonlinearity on the behaviour of S~, (17) is negligible from the 'first-order' optimality point of view when h ~ oo. Therefore, this xZ-cUSUM algorithm can be approximated by a repeated version of the SPRT (10) with zero lower threshold, which is a recursive algorithm of the form

to=min {t >- l :St>-h}; S,=-nt-~+lna

2'

~Xt ;

to-->1

detection delay with respect to the unknown change time to. Further details and explanations can be found in Benveniste (1986). Lorden proved (Lorden, 1971) that the 'worst' mean detection delay ~* given by (19) and the mean time T between false alarms given by (6) are asymptotically related via ~*=

in

p( O!, 0o)

(1+o(1))

as

2~---,oo, (20)

where p(O!, 0o) = ~° {" Po,(Yt)'~ the Kullback eo,~m-----'7~.../is information. Po,,t r,) / t Let us a s s u m e that y, 3~, x are the r a n d o m values. We say that the y = e s s s u p x if: (1) p ( x _ < y ) = 1; (2) if P ( x - < ) 0 = 1 then P ( y -<)0 = 1.

AUTO 2g:5-!

Z, = ~2 Vy~c-'V,;

(21)

V, = l{s,_,>0~V,-t + (Y, - 0o) nt = l { s , _ t > o } n t _ ! + 1,

where I{A) is the indicator function of the event A and n, is the counter of the observations since the last nonpositive S,_,. Note that X, = X[-,,+v The initial condition is So = 0.

3.2. Model A We consider the linear dynamic system described by the Gaussian state-space model (3). The likelihood function of this state-space model can be computed using the innovation sequence of the Kalman filter. In other words, we have

I. NIKIFOROV et al.

1280

first to transform the initial data (Yt)t_>l into the innovations sequence (Et)t>_l based upon the nominal (without change) state-space model (3) and, next, we detect a change in the innovation sequence (ct)t_>l. This sequence can be modeled

and, next, we detect the change by using the basic problem statement. Therefore there exists a transformation from the observations (Yt),~ to the orthogonalized residuals (e,)t_>l of the least squares algorithm

as

~.AC(0, Rt) ~(E,) = [JC'01(t, to), Rt)

if t < t0 if t -> to'

where tl(t, to) is the dynamic profile of the innovation sequence after change (Willsky, 1976, 1986), and R, is the covariance matrix of the innovation. Note here that the dynamic profile of the innovation sequence results from applying a step change on the input of the Kalman filter associated with the plant equation (3). Let us assume now that the model (3) is time-invariant, that the Kalman filter corresponding to this model is stable and, moreover, that the steady-state has been reached. In other words, let us assume that the covariance matrix R, = R is constant. Since the innovation before and after the additive change in the model (3) is a Gaussian independent sequence, the basic problem statement can be applied in this case. To do so, we have to approximate the complex profile of the Kullback information p(0, tl(t, to)) = tilT(t, to)RTt~lT( t, to)

et = Tryt.

(22)

(23)

by a constant value ~b 1 2n. Hence we have to choose the parameter b 2. This is the heuristic part of the procedure, which is described in subsections 3.4, 5.2 and 5.4. Finally, the x2-CUSUM algorithm for model A is obtained by replacing (Y, - 0o), X and b 2 by 8, Rt and b 2 respectively in (18) or (21).

In formula (24), T = (tl . . . . . tm-,) is a matrix of size m × ( m - n ) , where tj, j = l , . . . , m - n , are the orthonormal eigenvectors of the projection matrix P p

----I m -- H ( H T H ) - I H

s,

,

X, )

by the generalized likelihood ratio (with respect to the nuisance value X,)

s,

= In supx, pa(Y, IX,) supx, p0(Yt I Xt)

rank (P) = m - n

T,

corresponding to eigenvalue 1. The sequence (e,),>_l can be modeled as

o~f(O, O2lm_n)

Le(e,)

N(TTA, o2I,~_,)

Consequently, information is

in

this

case,

if t < to (25) if t -> 0 " the

p(0, TTA) = ~a 2 ATpA = ~b 1 2a.

Kullback

(26)

Finally, the x2-CUSUM for model B is obtained by replacing ( Y t - 0 0 ) , ~, m and b 2, by e,, o21m_,, m - n, and b~, respectively in (18) or

(21). It is worth noting that the transformation (24) of the measurements Y into a set of ( m - n) linearly independent variables by projection on to the left null space of the matrix H is nothing but the transformation of the vector ¥ into the parity vector of the analytical redundancy approach (Frank, 1990) TTH = 0;

3.3. Model B We now consider the linear static system with redundancy described by the regression model (4). Because the dimension m of the observations is greater than the dimension n of the unknown state, we have to use this redundancy in order to reduce this new change detection problem to the basic one. In this case we follow the minmax approach, which consists of replacing the unknown states (nuisance parameters) X, by their maximum likelihood estimates. In other words, we first replace the likelihood ratio = ,npA(Y, IX,)

(24)

TTT = [m_n .

3.4. Discussion Let us discuss now two features which play a key role in the analysis of the change detection problem for models A and B.

Dynamic profile of the Kullback information after a change. It should be clear that, for each state-space model (3), it is necessary to check a

priori the quality of the approximation of the time-variable Kullback information p(0, q(t, to)) by the constant b 2, because the actual value of the Kullback information plays a crucial role in the mean detection delay ~*. It is obvious that the first order optimality (20) holds for the x2-CUSUM algorithm under the assumption that the assumed value of the KuUback information is equal to the actual one. To check this for a given basis A we have: (i) to compute the dynamic profile *l(t, to) of innovation sequence (22); (ii) to compute the dynamic profile p(0, q(t, to)) of

Statistical fault detection algorithms the Kullback information by using formula (23): 1 t0+ ~*--1 (iii) to compute /)=~-g E p(0,11(t, to))

1281

inertial system heading gyro error:

l=to

where ~* is the assumed value of the mean detection delay; and (iv) to compare ~ with a given value of b 2. If 2~ -> b 2 then the quality of the approximation of the time-variable Kullback information is good. Detectability of a change. It results from formula (20) that the larger the Kullback information is, the faster the detection will be. Therefore, it is natural to use the Kullback information between the two distributions before and after the change as a performance index of statistical detectability of a change. In other words, a change is detectable if the Kullback information is positive. Let us discuss this performance index for model B. It results from formula (26) that the change A: IAI > 0 in the regression model (4) can be detected if the following inequality holds: p(0, TTA) = ~ / ~ T p / ~

>

0.

(27)

From the inequality (27) result two important consequences. First, if we want to detect all changes A such that [A[ > 0, then the matrix P should be strictly positive. But in fact rank(P) = m - n. In other words inf

1-

Q~--"(~ "~)

W~02'

H - (1, 0)

/ ~ (81(~

to))

Here v is the sampling period, and Yg is the gyro error time constant, with v<
p(0, "ETA) = 0.

A: JAI-->•>0

Roughly speaking, it is impossible to detect more than m - n sensor faults which arise simultaneously. For this reason we assume a priori that the vector A has k-< m - n nonzero components only. Furthermore, a necessary and sufficient condition for the statistical detectability of such a change with k nonzero components is that all the principal minors with order 1--
13 '

where

if t < to if t--t0

(28)

181 is the jump magnitude.

4.2. Algorithms We investigate three change detection algorithms based on this model. The first one is the xE-CUSUM algorithm which is practically equivalent to the two-sided CUSUM algorithm (16) as we explained in Section 3 ta = min {t -> 1: (gt* ->h) U (g~ ->h)}

4. A SIMPLEEXAMPLE We consider now a simple but representative example, and use it to compare three change detection algorithms. This preliminary comparison was performed in order to decide which change detection procedure should be used in the actual application and to assess each procedure in terms of performance (~* and 2r), robustness with respect to the magnitude of the change A, and computational cost.

4.1. Model We consider the following state space model (n = 2, m = 1) which is a simplified model of the

t5 +

(29)

g;= (g;_, - e,--~) go = g~ = 0 where t~ > 0 is the assumed change magnitude (tuning parameter). The second one can be described as follows (xKerr, 1982). First, the Kalman filter estimate X(t It) can be computed via the recursion ~ ( t i t ) = ~ X ( t - 1 I t - 1) + K(t)E,, where K(t) is the Kalman gain. In our case, it is

1282

I. NIKIFOROV et al.

relevant to use the second component

.l~2(t [ t):

3~2(t ] t ) = a'3~z(t -- 1 It - 1) + kEet, u where t r - - 1 - ~ g g and k2 is the second component of the Kalman gain K(t). Then, the stopping time of the G M A algorithm is given by t~=min{t>--l:lY2(t]t)]>--hl}.

(30)

The third algorithm is the G L R (Lorden, 1971, 1973):

best from the point of view of the accuracy of the computations of ~* and ~/'. The second possible method is to use simple analytical bounds for ~* and it. The accuracy of these bounds is much lower than that of numerical computations based on Fredholm integral equations, but is sufficient for our purpose. As we want to check superiority of the two-sided CUSUM algorithm, we use an upper bound for . . . . . and a lower bound for ir'c. . . . (Nikiforov, 1983, 1986)

t~ = min {t --- 1 :g, - h2} g, = max - I~k~,r--k + l

-* _< t+ 1 r..... 101 , 6 ,,l,(l_~_~I) \Z/

.= ej .

Discussion. The two-sided CUSUM algorithm (29) has two free values (the tuning parameters): the assumed change magnitude 6 and the threshold h. To reach the quickest detection we need to know a priori the change magnitude. If the 6 is unknown then the optimal properties of the CUSUM algorithm are lost but relative good performances are held. The relevance of the G M A and G L R algorithms for the purpose of comparison is motivated by the following facts. G M A is a very simple and popular algorithm, which is frequently used for navigation systems monitoring. For the design of the G M A and G L R algorithms, an a priori information about the jump magnitude ]6] is not necessary. It results (Lorden, 1971, 1973) that the G L R algorithm has uniformly good performance for a wide range of possible values of 161. On the other hand the G L R algorithm is very time-consuming because of the explicit maximization in (31). Hence, this algorithm may be too timeconsuming for real time on-board implementation, but is a very convenient benchmark. 4.3. Criteria for comparison As was explained in Section 2, the optimality criterion consists in minimizing the mean time detection delay ~ (or 'worst' mean delay ~*) for a given mean time between false alarms T. 4.4. Comparison between two-sided CUSUM and G M A Let us first compare the two-sided CUSUM and the G M A algorithms, by given by (29) and (30) respectively. For this comparison we assume that the change magnitude 161 is known, but not the sign of 6. In other words, we assume that in the formula (29) 6 = 16l. For the two-sided CUSUM algorithm (29) it is possible to compute ~* and T by using two methods. The first one is the numerical solution of some particular Fredholm integral equations. This method is the

e 161h- 1 - ] 6 l h

L.sumZ

(32)

q9( - ~-~)

,~2

1

+~ (33)

a,(x) = f_c~ ~0(x) dx

~(x) =

1 e_(X2/2).

On the other hand, the properties of the G M A algorithm are as follows. It is well known that, for the INS error model where the gyro error time constant Tg is large, the relevant value of the constant ~r in the G M A algorithm (30) is close to 1. However, the computation of 2rgma and ~gma for G M A has been done by Robinson and Ho (1978) for a---0.95. Note here that in this paper the mean delay ~gma (and not ~g*a) is computed by assuming a stationary distribution for the value of the decision function ~2(to1 [ t o - 1) just before the change time to. In order to extend these results to larger values of the constant c~, let us consider the limit case a = 1. In this case the filter equation in (30) is in fact the equation of the cumulative sum. For the comparison of the G M A algorithm with the two-sided CUSUM algorithm, we need to use the lower bound for the 'worst' mean detection delay ~g*a and the upper bound for the mean time between false alarms Tgma- In this case we can use the formula giving the bounds for the mean of the first time at which the cumulative sum ~ 2 ( t l t ) reaches - h i or h~. It results that _. > ~hl + • Tgma- max 0--
-iol 4 101 ~(-161) 0(101)

--

~(-161)e -2(h'+')lol - ~(161) ~(_161)e-2(h,+,)l~l _ ~(]6[)e2(h,-,)lal

4hl Tgma ~ h2 + 1 + W ~ . -

1283

Statistical fault detection algorithms ~-" (÷)

q-*

10 4 10 3

161 =0.25 t / /

.t

~=1 \

J

~

~=o.5\

I~1 =0.25 10 3

102

J

0.5~ ~

0.5 10 2

2

. ,.--

10 ~ 101 f

f

~

4

10 0 10 2

10 3

LO4

lO s

10 ° 0

FIG. 1. Results of the comparison between the two-sided CUSUM algorithm and the G M A algorithm when tr = 0.95. The delays ~*,~,m (solid lines) of the CUSUM and ~sm~ (dotted lines) of the G M A are plotted as functions of the mean time between false alarms T, for 161= 0.25; 0.5; 2; 4.

0.4

0.8

1.2

1.6

2

6

FIG. 3. Results of the comparison between the one-sided GLR algorithm and the one-sided CUSUM algorithms when t~ =0.5 and 6 = 1. The delays ~lr (solid line) and ~*usum (dotted lines) are plotted as functions of the actual change magnitude 6.

i-" 0")

103 j

161=o.25 10 2

J j

f'

/

f

/"

f s

f

. /" ~

J 4" J t

I~1 = 0.25

..

~ 1 / / / I

101

/ J

J

0.5

4

2

__...- ~ I

4

I

10 o 10 2

10 3

10 4

10 5

'r

FIG. 2. Results of the comparison between the two-sided CUSUM algorithm and the GMA algorithm when or= 1. The delays ~*,s,m (solid lines) of the CUSUM and ~m~ (dotted lines) of the GMA are plotted as functions of the mean time between false alarms T, for 161 =0.25; 0.5; 1; 2; 4.

The results of this comparison are represented in Fig. 1 for the case cr =0.95 and in Fig. 2 for --1. In Figs 1 and 2, the functions ~(?') and ~*(T) are depicted for 161=0.25; 0.5; 1;2;4. They show that the two-sided CUSUM algorithm is more efficient in all these cases except when ct = 0.95, 161 = 0.25; 0.5 and T < 103. In this situation the efficiency of the two-sided CUSUM and G M A algorithms are approximately the same. 4.5. Comparison between CUSUM and GLR In the previous comparison, we compared the two-sided CUSUM and G M A algorithms, for

the case where the assumed and actual values of the jump magnitudes are the same: 6 = 16[. Let us compare now the one-sided CUSUM and one-sided G L R algorithms. In other words we assume that the sign of 6 is known a priori, but not its magnitude. We compare now two one-sided CUSUM algorithms designed with different assumed jump magnitudes t~ = 0.5 and = 1 and the G L R algorithm designed without any a prior information about the jump magnitude 6. We make this comparison for different values of the actual jump magnitude 6, for a given value of the mean time between false alarms Tc. . . . = Tg,r=105. For the CUSUM algorithm, the relation between the 'worst' mean delay for detection z . . . . . and the jump magnitude 6 for a given value of the mean time between false alarms Tc. . . . can be computed with the aid of the numerical solution of the Fredholm integral equations, because the precision of this method is maximum. On the other hand, for the G L R algorithm this function can be estimated with the aid of an asymptotic formula (Lorden, 1971, 1973). The results of these comparisons are presented in Fig. 3, which shows that the CUSUM algorithm is slightly more efficient than the G L R algorithm around the optimal jump magnitude 16 - c51-< 0.4 + 0.8 and less efficient in the converse case. 5. NAVIGATION SYSTEMS MONITORING

This section is organized as follows. First, we briefly discuss the general methodology of navigation systems monitoring based upon statistical fault detection algorithms. Then we

1284

I. NIKIFOROV et al.

continue our discussion of the fault detection algorithms for three particular cases: • three inertial navigation systems (INS) monitoring. In this case we use the state-space model A(3). • redundant strapdown inertial reference unit (SIRU) monitoring. In this case we use the regression model B(4). • joint INS-aiding radionavigation systems monitoring. In this case we again use model A.

5.1. General methodology Navigation systems are standard equipment for planes, boats and other moving vehicles. Important examples of such systems are inertial navigation systems and radionavigation systems for planes. Conventional navigation systems use some measurement sources or sensors. For example, an inertial navigation system has two types of sensors: gyros and accelerators. Using information given by these sensors and the equations of the motion of the system, all desired signals (such as the geodetic coordinates and velocities of the plane, etc.) can be estimated. In view of safety and accuracy requirements, redundant fault tolerant navigation systems are employed. Fault detection and isolation of faulty sensors are among the main problems for the design of these navigation systems. We concentrate now on the fault detection problem which can be stated as a statistical change detection problem. The criterion to be used is quickest detection for a given mean time between false alarms. Fast detection is necessary because, between the fault onset time and its detection, we use abnormal measurements in the navigation equations, which is obviously undesirable. On the other hand, false alarms result in lower accuracy of the estimate because some correct information is not used. The optimal solution involves a tradeoff between these two contradictory requirements. Under normal operating conditions, the output signals of the navigation systems contain useful information and normal operation errors. If a sensor fault occurs then the output signals involve and additional error. It is obvious that the additional error is masked by useful signals. For this reason we have to remove useful signals by using the redundant measurements. For example, we use the difference between outputs of two identical INS in order to remove useful signals. Some recent investigations in this topic can be found in Gilmore and McKern (1972), Kerr (1987), Varavva et al. (1988), Jeerage (1990), Kireichikov et al. (1990) and Nikiforov et at. (1991).

5.2. Three inertial navigation systems monitoring Three identical INS are usually installed in big civilian airplanes. Such an INS has two kinds of sensors: laser gyros which measure angular velocities and accelerometers. We consider here the case of a soft fault as a bias in one sensor error, which appears at an unknown time. The problem is to detect the fault and isolate the failed INS as soon as possible for a given value of the mean time i" between false alarms. Design of the algorithm. The INS error model can be reduced to the form (3). This yields (see Appendix for details)

IX,+1 = ~(t + 1, t)X, + gt + A,. Yt H(t)X, + Ut + v,

(34)

where Xt e R t3 is the state vector, Yt e R 7 is the measurement of the INS, U, e R 7 is the useful signal (geodetic coordinates, velocities, etc.), (g,),_>o and (v,)t__1 are two independent zero mean Gaussian white noises with known covariance matrices Q(t)->0 and W ( t ) > 0 , respectively. The covariance matrix W(t) is diagonal. Let us denote the output Yt(i) of the INSi (i = I, 2, 3) and assume that only one INS can fail simultaneously. In this case the difference AY~/ = Y t ( i ) - Yt(j) can be written in the form JAXt+: = . ( t + 1, t)l~kXt + ~ + At, , [aV J H(t)AX, + v;

(35)

where A X t = X t ( i ) - X , ( j ) , (g;),_~o and (v;),_>l have covariance matrices 2Q(t) and 2W(t), respectively. Ato is the bias vector corresponding to one sensor error

{Oift
if t > to '

/~ ~-- (0 . . . . .

O, (~, 0 . . . . .

O)T

Therefore, for fault detection we have to compute three differences y]2, y~3, y23 and Kalman filter innovations c 12, c ]3, E23 by using system (35) and then we need to detect changes in each innovation s e q u e n c e (EI2)t>I, (•)3)t> I, (c~)t~] by using the algorithms described in Section 3. In practice, the third Kalman filter is not needed because: E 23 :

E 13 - - E 12.

(36)

Implementation and tuning. Let us now discuss the main steps of our preliminary investigation of the triple INS monitoring. • We determine the minimal absolute value I~Ol of the additional drift rate t$° in a gyro's error and the minimal absolute value Ic5;I of the additional bias 6 ¢ in an accelerometer's error (see Appendix), which lead to inadmissible

Statistical fault detection algorithms output parameter errors, i.e. to a fault. In the model (35) there are elements of the vector A in the fight-hand side of the plant equation. Simulation showed that for a typical flight duration the 1~t'1--6-8 nominal root-meansquare (RMS) of gyro errors and the I~;I-~ 40 - 80 nominal RMS of accelerometer errors. • We investigate the innovations behaviour under normal operating condition and under a fault, i.e. when one of the sensors has an error that exceeds the minimal fault level ([~e I or IBel). Note here that the simulation model includes many details (additional sensor errors, nonlinearities, etc.) which are omitted in the simple INS error model (35). First, simulation showed that the behaviour of the Kalman filter innovation sequence in stationary flight without fault is in relative good accordance with a white noise zero mean model. Second, because of model mismatches, additional biases in the innovations of the Kalman filter arise during turns of the plane, but the duration of these biases is short. Third, the dynamic profile of the innovation sequence (c~J)t___l after faults is analyzed. It results that the dynamic profile of the Kalman filter innovation sequence can be a highly time-varying function of time, which may even display a change of sign during turns (see Fig. 4). It results from our preliminary investigation and simulation that the greatest additional biases after fault are observed in the velocity innovations (in the case of 'horizontal' gyro or accelerometer fault) and heading innovation (in the case of 'vertical' gyro fault). Thus, for monitoring, the dimension of the measurement vector can be reduced to three. Furthermore, it follows that a scalar processing of these components of the innovation sequence is possible. The additional motivation for a scalar processing is the fact that in this case it is possible to use the Kalman filter with scalar measurement processing which is preferable for on-board implementation. Because the covariance matrix W(t) of the measurement noise vt is diagonal, the optimal properties of this Kalman filter is held. Therefore, the three components of the innovations are processed by the two-sided CUSUM algorithms (29), where el is the innovation component normalized by its RMS under normal conditions. We consider that the fault is detected if the alarm is given by at least one of the two-sided CUSUM. When the alarms are given by two two-sided CUSUM which process the same components of the innovations from different INS the isolation of the failed INS is realized by voting. For example, if the alarms

1285

are detected in the same components of c ~2 and c ~, the following hypothesis: ~2 = {INS2 is the failed INS} is accepted. The tuning parameters of the two-sided CUSUM algorithm (29) are the assumed change magnitude ~ and the threshold h. The parameter t~ is selected for a typical dynamic profile of the Kalman filter innovation for a weak fault. In this case, a sensor error which is lower than its minimal level does not lead to a sustained increase of the cumulative sum, hence alarm does not occur. Consequently, the two-sided CUSUM algorithm is tuned for low level faults and all faults which have change magnitude greater than this level will be detected with nonoptimal mean delay. Note here that it is possible to improve the relation between the mean detection delay ~* and the mean time between false alarms T for this type of faults. To see this, note that the comparison between the CUSUM and GLR algorithms of Section 4 shows that if the actual magnitude of change I~il is significantly different from /~, the GLR algorithm is preferable. As was explained in Section 4, the GLR algorithm is computationally complex. Approximation of this algorithm by a simple CUSUM procedure is thus of interest. A possible approximation of the GLR algorithm can be achieved by a joint use of two (or three) two-sided CUSUM for the same component of the innovation sequence. These two algorithms are designed with large and small values of the parameter tS, respectively. However, in fact we did not use this additional possibility for our real implementation. The threshold h of the two-sided CUSUM algorithm is selected in accordance with the given time between false alarms T and the sampling period v. The inequalities (32)-(33) and the numerical solution of the Fredholm integral equations are used for a first attempt. After this the tuning parameters ~ and h are checked by using the simulation model which includes the normal operating mode of the INS, the faults of different sensors, the model mismatch, etc. Simulation results show that for compensating the model mismatch, it is necessary to increase ~ and (or) h above their original values. For example, assuming that the sampling period v is equal to 30 sec, the simulation shows that for a given time between false alarms l"=5000hr, a gyro fault such that A - 1 4 nominal RMS can be detected with a delay about 8-15 rain. For the isolation of a failed INS, the tuning of

I. NIKIFOROV et al.

1286

f~ I I ! I

/

I

.......

.........

.......

4. Normalized dynamic profiles of the two velocity (rV~, 6Vy) innovations under the jump of the 'horizontal' gyro to 14 nominal RMS. The change time to = 10 min. The innovations of the velocities 6Vx and 6Vy are denoted by a solid line and a dotted line, respectively. Fro.

the two-sided C U S U M is accomplished through an additional simulation. This simulation takes into account the crosscorrelation between different innovations which takes place when the triple INS is operating, and the voting procedure. It results from this simulation that the mean time between false alarms for each two-sided C U S U M algorithm should be about 50000hr. Nevertheless, the detection of the sensor faults and the isolation of the failed INS are achieved in proper time, i.e. without any appreciable increase of the output navigation errors. In particular, if the fault is similar to that of Fig. 4, the detection delay of the sensor fault and isolation of the failed INS is about 15-20 min.

20

alarmtime

isolation time

i0

0

_ .

~

T

0

60

~

I

120 "horizontal"

gyro

min

fault

FIG. 5a. Behaviour of decision functions after the fault in INS2. The flight of the IL-96-300 airplane, 28.09.92. Finally, let us briefly discuss the results of the implementation of the proposed monitoring scheme (a set of two-sided C U S U M + voting) for real data. This monitoring scheme has been tested on experimental data from test flights of the IL-96-300 airplane. i4.

2~

min -|.

The main conclusions of this test are as follows: • the 'horizontal' gyro fault was detected during one flight processing. The behaviour of the decision functions is shown in Fig. 5a; • no false alarms were detected; • the innovations behaviour was in sufficiently good agreement with a white-noise zero mean model. A typical sample of the innovation is shown in Fig. 5b; • the actual innovations RMS were about 0.6-0.9 of the RMS c o m p u t e d during the simulation.

-2

5.3. Redundant strapdown inertial reference unit

-3

monitoring

-4

Fro. 5b. A typical sample of the normalized velocity (SVx) innovation. The flight of the IL-96-300 airplane, 4.12.90.

The use of a strapdown inertial reference unit (SIRU) in tandem with a monitoring (fault detection and isolation) system can significantly

Statistical fault detection algorithms increase its reliability [see e.g. Jeerage (1990)]. Conventional SIRU incorporating four or more sensor pairs (gyro/accelerometer) with skewed axes can provide us with internal integrity monitoring. Integrity monitoring synthesis for redundant combination of similar sensors has been considered by Gilmore and McKern (1972) among others. We assume now that the gyro and accelerometer measurements are considered separately. The measurement model can be represented by model B(4). Design of the fault detection algorithm. The simplified SIRU measurement equation can be written as Yt =

HXt + Vt "4- At,,

(37)

where X t E R 3 is the unknown state (three orthogonal body angular velocities or accelerations), Yt e R m (m = 4, 5, 6) is the measurement vector, H is the known matrix of the orthonormal projection vectors from the sensor axes to the orthogonal body axes, (v,),_>l is a white noise with covariance: W = o21,,. The variance o 2 is known. Therefore, because m - n - 1, there exists information redundancy in this case and for the fault detection we can use the recursive version of the x2-CUSUM algorithm (21) which was described in Section 3: t. = min {t --- 1 : S,>-h}; b2

St = - n , ~- + In Xt = 1

o-n7

- n bZn2xt~

v Vv t;

(38)

Vt = l{s,_,>o}Vt-1 + TTyt

nt

=

l{s~_l>o}nt-i + 1.

Tuning and simulation of the xZ-CUSUM. The tuning parameters of the x2-CUSUM algorithm (38) are the assumed generalized signal-to-noise ratio b and the threshold h. The parameter b is selected for a typical sensor error under the assumption of a weak fault. The threshold h is selected in accordance with some given mean time between false alarms 1" and sampling period v. The asymptotic formula (20) is used in a first stage. The SIRU consisting of six degree-offreedom gyros and six accelerometers is considered (m = 6, n = 3). Sensors are equally spaced on a cone with a half-angle of 54.7. We require that the mean time between false alarms should be 25 000 hr, and the sampling period v be 1 sec. Specified levels of the faults are about 30 nominal RMS for gyro faults and 10 nominal RMS for accelerometer faults. For the simulation, the error model (37) has

1287

been extended in order to include the model mismatch-additional components of gyro/accelerometer error (see Table 1 in Appendix). This model has reproduced the normal operating mode of the SIRU and the faults of gyros/ accelerometers. It results from this simulation that the X e - c u s u M algorithm detects the faults in any gyro with a delay of about 3-4 min and in any accelerometer with a delay of approximately several seconds. 5.4. Joint INS-aiding navigation systems

monitoring As the errors of a triple INS increase with time, its outputs should be corrected. For this purpose, the Omega navigation system (ONS) that supplies latitude and longitude data [see e.g. Pierce (1989)] is widely used. However, the sun's disturbances can result in anomaly errors in the ONS output data. The errors can increase up to 10-20 miles for 15-30 min of flight time (typical ONS errors are 0.5-2 miles). This may be interpreted as a soft fault. The data from INS and ONS are jointly processed by Kalman filtering. Hence, in this case we have an error model of the form A. It is obvious that, if we provide incorrect data from ONS in the INS/ONS processing, then the joint filtering errors (INS/ONS errors) will increase. For this reason, it is of interest to detect a soft fault in the ONS as soon as possible and not use this incorrect measurement for the correction of the triple INS any longer. Design of the algorithm. Let us denote the outputs (latitude and longitude) of the triple INS and of the ONS by Yt(INS) and Yt(ONS) respectively. We assume that a failed INS does not exist in this triplet. In this case the difference AYt = Y t ( I N S ) - Y,(ONS) can be written as

{

Xt+l

= @ X t + ~t

A Y , = H ( t ) X t + v, + At,,

(39)

where Xt • R 2, AY, e R 2, (g,)t>0 and (vt)t>_l have covariance matrices Q=o~12 and W = O2~12 respectively. A , 0 e R 2 is the vector of sun's disturbances. The parameters o~ and O2~ are known. Therefore, for fault detection we need to compute the Kalman filter innovation of the system (39), and then we detect a change in this innovation sequence by using the recursive version of the x2-CUSUM algorithm (21). Tuning and simulation of the z2-CUSUM. In this model we assume that the sampling period v is equal to 30 sec. The matrix @ is diagonal and selected in accordance with the error time

1288

I. NIKIFOROV et al.

constant Yc which is about 60 min:

1 -y--~

0

~= 0

1

v

The RMS oe of the g, is equal to 460 m. The RMS tr~ of the measurement noise vt is equal to 100 m. Let us outline some features which are new with respect to triple INS monitoring. First, in this case, it is not possible to split the sun's disturbances detection problem into two scalar detection problems because physically this type of fault is characterized by the Kullback information which is proportional IAI 2. Hence, the two-dimensional x2-CUSUM algorithm is definitely useful here. Second, it is obvious that a monitoring scheme is practically relevant if ~'* is smaller than the time of ONS errors increasing. For this reason the tuning parameters b and h of the x2-CUSUM algorithm are chosen in accordance with the minimal error level (IAI) detected and with ~*. Then the tuning parameters b and h are checked by using the asymptotic formula (20). If ~P is intolerably small under the given ~* and I,~,1, it is necessary to repeat the tuning parameter selection for greater values of b and (or) ~*. When an acceptable value of T is reached, the tuning parameters are checked by using a simulation model which includes the normal operating mode and the anomalous errors of the ONS under different levels and profiles. For the final investigation of the characteristics of the monitoring algorithm, the full simulation model is used. This simulation model includes the joint operation of the triple INS/ONS, as well as the model mismatch. This simulation shows that, for a mean time between false alarms 2r-~ 30005000 hr, the anomaly errors are detected with a delay of about 5-15min and the INS/ONS errors do not exceed 8-9 miles for ONS errors which become larger than 20 miles for 15-30 min of flight. Thus the accuracy of the joint INS/ONS estimation with the fault detection algorithm is increased by at least a factor of 1.5-2 with respect to a standard joint triple INS-aiding ONS processing system. 6. CONCLUSIONS The on-line change detection problem has been considered for an independent Gaussian vector sequence, where the mean vector changes at some unknown time. The problem is to detect the change as soon as possible. For this purpose, the x2-CUSUM algorithm based upon the theory of the invariant SPRT was proposed. In the

scalar case this algorithm is practically equivalent to the classical Page's two-sided CUSUM algorithm. The generalization of the z2-CUSUM algorithm to other two popular models: the state-space model (A) and the regression model (B), is proposed. It was proved that, for both cases, we need to transform the initial data into another sequence (the innovations for model A and the orthogonalized residuals for model B) and then to apply the xE-cUSUM algorithm. The advantages and disadvantages of this algorithm were studied by performing a comparison of three change detection algorithms for a simple but representative example of gyro fault detection. The proposed methods were applied to navigation system monitoring. Three particular cases of this general problem were discussed: triple INS monitoring, skewed axes SIRU monitoring and joint triple INS-aiding (Omega) navigation system monitoring. The first and third cases are based upon model A and the second one is based upon model B. Our general methodology includes three main steps: the design of the monitoring algorithm, and the implementation and tuning of this algorithm. It is obvious that in practice the models A and B are the result of linearizing and simplifying the significantly more complex 'full-dimensional' (or/and nonlinear) model of the plant. For this reason, after tuning the CUSUM algorithm, we need to validate it by using a simulation model which includes the model mismatch and then to compensate the model inadequacies by increasing the t5 and (or) h parameters for the two-sided CUSUM algorithm and the b and (or) h for the xE-cUSUM.

The main results of testing the proposed triple INS monitoring scheme on experimental data from test flights of the IL-96-300 airplane were reported. Acknowledgements--The authors would like to express their thanks to A. Benveniste and M. Basseville for very fruitful discussions. They are also grateful to the reviewers for their many constructive comments on the paper. REFERENCES Basseville, M. and I. Nikiforov (1993). Detection of Abrupt Changes. Theory and Applications. Information and System Sciences Series, Prentice-Hall, Englewood Cliffs, NJ. Benveniste, A. (1986). Advanced methods of change detection: an overview. In M. Basseville and A. Benveniste (Eds), Detection of Abrupt Changes in Signals and Dynamical Systerrt~, Vol. 77, pp. 216-258. Springer, Berlin. Frank, P. M. (1990). Fault diagnosis in dynamic systems using analytical and knowledge based redundancy--A survey and new results. Automatica, 26, 459-474. Ghosh, B. K. (1970). Sequential Tests of Statistical Hypothesis. Addison-Wesley, Cambridge, MA.

Statistical fault detection algorithms Gilmore, J. P. and R. A. McKern (1972). A redundant strapdown inertial reference unit (SIRU). J. Spacecraft, 9, 29-47. Jeerage, M. K. (1990). Reliability analysis of fault-tolerant IMU architectures with redundant inertial sensors. IEEE AES Mag., 5, 23-27. Kerr, T. H. (1982). False alarm and correct detection probabilities over a time interval for restricted classes of failure detection algorithms. 1EEE Trans. Inform. Theory, 1"1"-24, 619-631. Kerr, T. H. (1987). Decentralized filtering and redundancy management for multisensor navigation. IEEE Trans. Aerosp. Electr. Systs., AES-23, 83-119. Kireichikov, V., V. Mangushev and I. Nikiforov (1990). Investigation and application of CUSUM algorithms to minitoring of sensors. Third All-Union Seminar on Detection of Change in Random Processes. In Statistical Problems of Control, Issue 89, pp. 124-130. Vilnius (in Russian). Lorden, G. (1971). Procedures for reacting to a change in distribution. Annals Math. Statistics, 26, 7-22. Lorden, G. (1973). Open-ended tests for Koopman-Darmois families. Annals Statistics, 1, 633-643. Nikiforov, I. V. (1978). A statistical method for detecting the time at which the sensor properties change, Preprints

IMEKO Symp. on Applic. Stat. Meth. in Measurement, Leningrad, pp. 1-7. IMEKO Publications. Nikiforov, I. V. (1983). Sequential Detection of Abrupt Changes in Time Series Properties. Naouka, Moscow (in Russian). Nikiforov, I. V. (1986). Sequential detection of changes in stochastic systems. In M. Basseviile and A. Benveniste (Eds), Detection of Abrupt Changes in Signals and Dynamical Systems, Vol. 77, pp. 216-258. Springer, Berlin. Nikiforov, I., V. Varavva and V. Kireichikov (1991). Application of statistical fault detection algorithm for navigation system monitoring. Proc. SAFEPROCESS'91, Baden-Baden, lap. 351-356. Nikiforov, I. V. (1993). On first order optimality of change detection algorithm for vector case. To appear in Avtomatica i telemehanika (in Russian). Page, E. S. (1954). Continuous inspection schemes. Biometrika, 41, 100-115. Parusnikov, N. A., V. M. Morozov and V. I. Borzov (1982). Correction Problems in Inertial Navigation. Moscow University Publications, Moscow (in Russian). Pierce, J. A. (1989). Omega. IEEE AES Mag., 4, 4-13. Robinson, P. B. and T. Y. Ho (1978). Average run length of geometric moving average charts by numerical methods. Technometrics, 20, 85-93. Shiryaev, A. N. (1961). The problem of the most rapid detection of a disturbance in a stationary process. Soviet Math. Dokl., 2, 795-799. Varavva, V., V. Kireichikov and I. Nikiforov (1988). Application of change detection algorithms to monitoring of measurement systems of moving objects. Second All-Union Seminar on Detection of Change in Random Processes, In Statistical Problems of Control, Issue 83, pp. 169-174. Viinius (in Russian). Willsky, A. S. (1976). A survey of design methods for failure detection in dynamic systems. Automatica, 12, 601-611. Willsky, A. S. (1986). Detection of abrupt changes in dynamic systems. In M. Basseville and A. Benveniste (Eds), Detection of Abrupt Changes in Signals and Dynamical SystentL Vol. 77, pp. 27-49. Springer, Berlin. APPENDIX--INS E R R O R MODEL This Appendix describes the proposed simple INS error model to be used in the detection algorithm. See Parusnikov et al. (1982) for exact results and details. A.1. Plant equation The dimension n of the state space vector X in the plant equation (34) is equal to 13. To simplify our explanations, we have divided the plant equation into two parts: the dynamic

1289

model of INS errors (8 x 8) and the additive model of the sensor errors (5 x 5). We start from the dynamic model of INS errors and then continue with the additive model of the sensor errors. Dynamic model of INS errors. The plant model in the continuous time has the form x=rx+~ where X = (y~, yy, p~, py, .~, %, flz, y~)v

F=

0 -to~ 0

to~ 0 0

1 0 0

0 1 to~

0 0 -~o~

0 0 0 -tOy

0 0 0 ¢o~

-~o~ 1 0 0

0 0 1 0

0 -to~ 0 to, -to z 0 COy -COx

py

-p~

0

0

0

x~z

0 0 0

toy -to~ 0

01 0 0

0 0 0 0

0 0 0 0

0

0

0

R '

The definitions of the parameters of the plant equation are as follows: • (ax, o% otz)"r is the angular vector of a small rotation from the 'platform' frame to the (ideal) local-level frame, • (Y~, 7y, y~)T is the angular vector of a small rotation from the system-computed frame to the (ideal) local-level frame, • (fix, fly, flz)T is the angular vector of a small rotation from the 'platform' frame to the system-computed frame, ~px, py are the generalized impulse errors. (to~, toy, ~o,)T is the vector of local-level frame angular rates relative to the inertial space, • (p~, py, pz) T is the vector of local-level frame angular rates relative to the earth, • R is the mean earth radius, • too is the Schuler frequency. The sensor errors are described by the following two vectors: • ~ = (~x, ey, ~,)T is the vector of generalized drift rate of the gyros in local-level frame, • e = (P,, Py, O,)T is the vector of generalized accelerometer errors in local-level frame. Additive model of sensors errors. Note here that the component p~ of the vector e is not taken into account in the dynamic model of INS errors because navigation in a horizontal plane is assumed. Each component of the vectors ~ and e has following form

P = Oro + Orb + Prn + 6t,~

where (¢fd, ¢~d, Cw,) are the fixed drift, the random drift and the wideband noise components, respectively, and (Pro, Orb, Or,) are the fixed bias, the random bias and the random noise components, respectively. The assumed RMS values of these error components are summarized in Table 1. We assume that the fault can be modeled as an additional drift rate 6,¢0in one gyro's error or as an additional bias 6,,e m" one accelerometer's error. Because we consider the strapdown INS. the body frame to local-level frame transformation of the sensor errors is used. A.2. Measurement equation The measurement equation (3) Y=HX+v involves the following vector Y and matrix H: Y = (6~0, ~ , ~vx, ,Wy, ,~,, el0, 6~v~)r

I. NIKIFOROV et al.

1290 -cos g

sin g

sin g" cos g" cos q0 cos q0 H=

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

R

0

0

-RaJx

-Rto~

0

0

-R

0

0

0

-R%

-Rt%

0

0

0

0

0

0

0

0

0

0

- c o s ~PG

0

0

0

0

0

0

--sin~pGtanO

1

1

sin ~Pc cos 0

TABLE 1. INS SENSORS ERROR MODEL

Description Ring-laser gyros Fixed drift Rate bias (deg/hr) Random drift Rate bias (deg/hr) Wideband noise (deg/V~r)

RMS value

0.01 0.005

Notes

Not included in the filter error model T = 1 hr; the filter RMS equals 0.01 deg/hr

0.005

Accelerometers Fixed bias (#g)

50

Random bias (#g)

25

Random noise (#g)

5

Not included in the filter error model T = I hr; the filter RMS equals 50 #g

cos ~/,~ cos 0 sin ~Pa - c o s 1pc tan 0

The definitions of the parameters of the measurement equation are as follows: • (if, ;t) are the geodetic coordinates, • V,, Vy are the velocities relative to inertial space, calculated in computed frame, • ¢, 0, lpc are the bank, pitch and gyro heading, • g- is the 'wander azimuth angle' between the x axis of navigation reference-coordinate local-level frame (x, y, z) and the east direction counter-clockwise The parameter Yz is usually ignored because it is small and weakly observable. The covariance matrix W of the measurement noise v is diagonal. The assumed RMS of the measurement errors are 20m [position--(¢, ~.)], 0.04 m/sec (velocities--Vx, Vy), 3 arc min (pitch and b a n k - - 0 , 0) and 1 arc min ( h e a d i n g - - ~ a ) . The measurement errors can be reduced by means of the output synchronization of all INS.