A Bayesian Approach to Manoeuvre Tracking and Detection

Copyright © IFAC 12th Triennial World Congress, Sydney. Australia, 1993

A BAYESIAN APPROACH TO MANOEUVRE TRACKING AND DETECTION A..J. Isaksson*'t':t and F. Gustafsson** *DepartfMnl 0/ Signals, SellSors and Systems, Royal/lIStitUle a/Technology, S-1oo 44 Stoclcholm, Sweden **Department a/Electrical Engineering, LinJcijping University, S-58183 Lin/cijping, Sweden

Abstract. This paper describes general modelling of radar target tracking scenarios. Manoeuvres are modelled using a stochastic continuous-time model, and estimated in a Bayesian framework leading to a bank of Kalman filters. To reduce the number of filters a pruning algorithm is used, enabling simultaneous state estimation and manoeuvre detection. Key Words.

Target tracking, extended Kalman filter, multiple models.

1 INTRODUCTION

2 TARGET TRACKING MODELS

The aim of this paper is to study models for radar target tracking. The intended use is civil air traffic controller systems, i.e. typical targets are airplanes and helicopters. Although there is a vast literature on the subject, see e.g. the books BarShalom and Fortmann (1988) and Bar-Shalom (1991), there seems to be some room for improvement. Typical approaches described in the literature include

We will in this section discuss modelling of target trajectories. The presentation is divided into three subsections, where the first and the last basically reviews previously used models for the dynamics and the measurement. The novelty is the way to which we approach modelling of the manoeuvres in Subsection 2.2, and the consequences for the estimation algorithm in the following sections.

• simultaneous estimation of the states and the manoeuvre combined with fault detection, Chang et al. (1977), Thorp (1973),

2.1

The target dynamics

A two-dimensional target tracking situation can be described using the following four states

• using a discretised range of manoeuvres leading to a bank of Kalman filters, Ricker and Williams (1978), Gholson and Moose (1977), Moose (1975).

x(t) y(t) v(t) h(t)

• using a discrete-time model augmented with a manoeuvre model, also leading to a bank of filters, Vacher et al. (1991), Blom and Bar-Shalom (1988).

position along x-axis position along y-axis velocity heading (angle to positive x-axis)

Assuming known acceleration a(t) and turn rate w(t) we get the (nonlinear) state space equation

x(t) iJ( t) ti(t) h(t)

The approach suggested here can be seen as a development of the third one above. We advocate to model the manoeuvres in continuous-time and to discretise at the latest possible stage. The model applied is a nonlinear one which implies use of the Extended Kalman Filter. The manoeuvres have been modelled via acceleration and turn rate that are either changing or not during one sampling interval. This model eventually also leads to a multiple model Kalman filter.

v(t) cos h(t) v(t) sin h(t) a(t) w(t)

(1)

The purpose has not been to compare different nonlinear models as suggested in Vacher et al. (1991). This could perhaps be scope for future research.

2.2

tWork done whilst visiting Dept. of Electr. and Comp. Eng., University of Newcastle, Australia IThe first author gratefully acknowledges financial support from The Swedish Institute and J. Gust. Richert's foundation.

Modelling the manoeuvres

In addition to the four states above, the model needs to be augmented by considering explicit models for the possible manoeuvres. We will do this by introducing two states to capture 841

We finish this subsection by giving the complete augmented state space model. Introducing the state vector

• change of velocity • change of course

x(t) = (x

These two manoeuvres have been chosen since they are almost independent in practise. It is not entirely true since, for example, an aircraft commencing a dive is undergoing a simultaneous change of acceleration and turn rate. But for civil air traffic control it is most likely a valid assumption. Moreover, it is also valid to assume the acceleration and turn rate to be piecewise linear. This is not the case if, for example, accelerations in the x and y directions are used to model the manoeuvres instead. A start of a manoeuvre will be assumed as an almost instantaneous change in the acceleration a(t) or turn rate w(t), that of course can take place at any time during a sampling interval. We will, however, model this as occurring at the onset of a sampling period. A piece-wise constant manoeuvre variable ~(t) is introduced for this purpose. We have, for

the complete nonlinear state equation becomes xa COSX4 X3 smX4 Xs

~(t)

=

~(kT)

~(kT)(2 - ~(kT))VI(t)

~~(kT)(~(kT) - 1)v2(t)

2.3

x(kT) y(kT)

if no manoeuvre if change of acceleration if change of turn rate

= n) = {

for n

J.lI

n= 1 n = 2

J.l2

The transformation (4), however, introduces a cross-correlation between the two noise components. Define the multivariable output vector

.

=

~(kT)(2 - ~(kT))VI(t) ~(kT)/2(~(kT) - 1)v2(t)

y(t)

=(

x(t) ) y(t) .

Then we have the measurement equation

(2)

where

y(kT)

Hx(kT)

H

[

+ e(kT)

100000] o 1 000 0 '

(5)

where we have chosen to denote discrete-time quantities by overbars, and an approximate expression for

is a zero-mean continuous-time white noise with

Ev(t)v(t

(4)

A=(~ ~8)'

At a manoeuvre, we model the change in acceleration and turn rate as driven by random noise, i.e. for t E [kT, kT + T]

iz(t) w(t)

r(kT) sin 6(kT), r(kT) cos 6(kT).

The measurement errors in r(kT) and 6(kT) can be considered as mutually independent discretetime white noise processes, with a covariance matrix

=0

J.lo

The measurement equation

By the nature of radar processing the measurements will of course be discrete-time. Typically they consist of bearing 6 and distance r to the target. A transformation to the global Cartesian coordinates is given by

The case of simultaneous manoeuvres has been omitted since it is a relatively unlikely event, but can of course easily be included. We will model the prior knowledge of different manoeuvres by p(~(kT)

(3)

X6

t E [kT, kT + T] 0 = 1 { 2

v haw)T,

y

+ rf = c5(r)Q.

R( kT)

Here c5(t) denotes a Dirac pulse. Since the manoeuvres can be assumed independent the covariance matrix Q is diagonal

= E e( kT)e( kTf

can be found, for example, in Anderson and Moore (1979) as equation (4.5) page 55.

3 THE STATE ESTIMATION SCHEME

To introduce the model for the manoeuvres like this already in continuous-time is to our knowledge a new idea. The difference may seem to be a very subtle one, but we believe rather important. This will be· commented further in the following sections.

Following the previous section we have by (3) and (5) a nonlinear model that can be described as

X y(kT) 842

f(x) + B(kT)v(t), Hx(kT) + e(kT).

(6)

where T

B(kT) =

0 0 0 e(2 - e) 0 0 0 0

[0 0

0

e/2(e - 1)

]

.

o

An important observation and a special feature with (6) is that conditioned on a particular sequence, denoted et, this is a known time-varying state space model. Since it is nonlinear we have to invoke what is normally referred to as the Extended Kalman Filter (EKF), see e.g. Jazwinski (1970).

Remark 2: The significance in modelling the manoeuvres in continuous-time really lies in the matrix Q. If the manoeuvres are modelled adding a discrete-time noise to, for example, turn rate after first sampling the model takes three samples before there is any effect on the position estimates x and y. This means that M has to be very large to allow for models to live long enough to have an impact on €(t). The matrix Q, however, will have non-zero elements even in the top left-hand corner, leading to immediate adjustment of the gain for x and y. It might be tempting to use an approximation like

Time update:

I

x(t + Tit) = x(tlt) + t

t +T

f(x( T)) dT

P(t + Tit) = F(t)P(tlt)F(tf + Q(t) A(t)

= [Of;(x(tlt))]

,

F(t)

The result of (8) can be computed using any standard program for sampling continuous-time systems, like the routine c2d in the Control System Toolbox for Matlab.

= eA(t)T

OXj

where

Q(t) =

fT

io

T

G

eA(t)T B(t)QB(tf eA(t) T dT

This, however, lead to several false alarms in the simulations. Probably because T is long enough to make the approximation bad, yielding an erroneous distribution of gains over the states.

Measurement update:

x(t + Tit + T) = x(t + Tit) +K(t + T)(y(t P(t + Tit

T) = +

4

+ T)HjT

+K(t + T)R(t K(t

o

+ T) - Hx(t + Tit))

+ T) = [/ - K(t + T)H]P(t + Tit)

x[/ - K(t

T P(t + Tlt)H H P(t + Tlt)HT + R(t + T)

et

et,

filter the data through • For each possible a Kalman filter for the (conditional) known state space model (5).

e,

whose • Choose the particular value of Kalman filter gives the smallest prediction errors. Below we will formalise this in Theorem 1 by introducing the maximum a-posteriori (MAP) estimate of The theorem is formulated for linear systems. For our nonlinear model (6) the theorem does not hold exactly, but inserting the quantities computed by the Extended Kalman Filter into (9) certainly gives a good suboptimal solution. For further details on the theory in this section see Gustafsson (1992).

I.e.

Q= [(A ® /) + (I ® A)] col Q + col BQBT .

et.

(7) Here ® denotes the Kronecker matrix product and col Q stacks the columns of Q on top of each other forming a vector. Hence by (7) colQ =

i

fT .

eAT BdT

ESTIMATING THE MANOEUVRE SEQUENCE

As pointed out in the previous section, for a particular manoeuvre sequence the Extended Kalman Filter can be applied to the model (6). One natural strategy to infer the decision sequence is the following:

+ T)K T (t + T)

Remark 1: The covariance matrix Qfor the corresponding discrete-time noise can be obtained by observing that it is the solution to the Lyapunov equation

col

= iT eAtBdt.

(8)

o

Theorem 1 (MAP estimate of <) Consider the linear model

where

x(t + T) 843

=

F(t)x(t)

+ v(t)

yet)

=

Hx(t)

+ e(t)

the minImum variance (or conditional expectation) estimate as

where we assume the noises to be Gaussian, with

M

Ev(t)v(tf E e(t)e(tf

x(tlt)

where M denotes the number of models. This estimator has the obvious disadvantage that at any time there are 3t possible sequences, i.e. we have to store 3t different filter updates. For this approach to be applicable in practise the number of branches obviously has to be reduced dramatically. There are two main alternatives that have been proposed in the literature: sequence merging ( Blom and Bar-Shalom (1988» or sequence pruning ( Thgnait and Haddad (1979), Andersson (1985». We will here only pursue the latter approach.

e

1

tiT

p y

k=O

ai(t)xi(tlt).

R(t)

is given by maximising The MAP estimate of its a posteriori probability

p(elyt)

=L

= - (t) IT p(~(kT)h([E(kT),

SE(kT».

(9) Here [E(t) and Sdt) are the prediction error and computed its covariance matrix conditioned on recursively from the K alman filter

e,

fEet)

yet) - HXE(tlt - T)

SECt)

H PE(t\t - T)H T

A practical recursive approximation which has turned out to work well in a number of applications is given in the following algorithm.

+ R(t)

and 'Y( x - Jl, P) denotes the Gaussian probability density function, i. e.

Algorithm 1 (Local search) In Theorem 1 only M possible sequences are examined. They are chosen according to the following rules. Proof: By Bayes' rule we have

• At time t, keep only the M most probable sequences ~:, i = 1,2, .. , M. That is, the sequences with the largest weights ai at time t among the considered ones.

p(e) p(yt le) p(yt)

• At time t + 1, let the M considered sequences branch into the 3M branches, ~j+l = (~L~(t+l»forall~: and~(t+1), and update their estimates Xj(t + lit) and corresponding weights aj.

p(e) p(y(t)lyt-T,e)p(yt-Tle-T). p(yt) It is a well-known property of the Kalman filter

that

One advantage with this search algorithm is that there is a minimum of rules for the pruning, making it easy to implement. Another advantage is that since we do pruning, each filter has a manoeuvre history which enables simultaneous state estimation and detection of manoeuvre time instants.

Thus, we get

The desired result now follows iteratively using the expansion

5

o Hence for each manoeuvre sequence ~: a separate Kalman filter has to be run, producing an estimate xi(tlt). At time t there is a weight associated with each filter, viz.

DESIGN PARAMETERS AND INITIAL VALUES

Before we present a simulation example in the next section, let us discuss the choice of design parameters. In addition to M - the number of models - the method has essentially three design parameters, viz. Jl, Q and A. Of those the measurement noise variance A is supposed to be reasonably well known for most types of radars. If this is not the case estimation of A can be included in the algorithm, but we will not discuss that further here. The jump probabilities Jl are only used in the update of ai(t), equation (10).

The one with the largest weight corresponds to the MAP estimate of the manoeuvre sequence at time t. Using the theorem it is also straightforward to obtain an estimate of the state x( t). We get 844

As long as 1J1 and 1J2 are relatively small the performance is basically the same (see also Andersson (1985)). For the choice of Q we have to use insight into the typical tracked objects characteristic. To get some intuition for the impact of Q, assume that we sampled the model for the acceleration separately. At a manoeuvre that would give

aCt + T) = aCt) + Vl(t),

(11)

where we have for the discrete-time noise Vl(t)

= Tql.

EVl(t)2

Assume that we may get an instantaneous jump from zero to the maximum acceleration a max • The standard deviation of Vl(t) should then be chosen accordingly to be able to pick this up in one sample. Choose, for example, two standard deviations equal to a maXl i.e. ql =

a~ax 4T .

Similarly, if we want to design for a maximum turn rate we should pick

The matrix T· Q/100 is introduced in P(O) to facilitate for any minor initial errors in acceleration and bearing. Larger deviations have to be picked up as manoeuvres.

6 SIMULATIONS We will in this section simulate the trajectory marked by circles in Figure 1. The radar is located at the origin, and the target enters from the right approaching the radar with a constant speed of 250 m/so After 50 s it commences a circular 360 degree turn at 1.9° /s (~0.8g). This turn is completed after 240 S. Then finally, back in linear motion, at 290 s the tar~et goes into a longitudinal retardation of 2 m/s . The rotation period of the radar is assumed to be 5 s, i.e. T = 5. The measurements of range and bearing are corrupted by Gaussian noise, with standard deviation 100 m and ±0.1 ° respectively.

Design parameters: The following values have been used for the design parameters IJ

A For the sake of completeness, and to make the simulations reproducible, we will also describe how the initial values x(O) and P(O) could be obtained. For this the first two observations can be used

x(O) = [

y]

reO) sin B(O) ] reO) cos B(O) [

JV'"2+ vy2 ,

tan-le V y /v",)

(r(T) sin B(T) - reO) sin B(O))/T, (r(T)cosB(T) - r(O)cosB(O))/T. The initial variance is computed using the known measurement accuracy given by A and the Jacobian J defined by

(ho/ar(O) ] JT _ axo/M(O) - [ axo/ar(T) . axo/aB(T)

[

Roo

0

TQ/lOO

75 2 0 5 0

0.01 ) ;

(0.l1r~180)2

) ;

0 ); M=5. 0.001

The values of Q according to Section 5 correspond to a maximum acceleration 1 g and turn rate 8° / s.

]

In Figures 2 and 3 we have also marked the detected manoeuvres. Obviously the start and finish of the turn are detected the very next sample, whereas there is a two sample delay for the retardation. This is not surprising since 2 m/s 2 is after all relatively small compared to the measurement noise in range.

7 CONCLUSIONS We have in this paper addressed modelling of target trajectories, and target manoeuvres. The main point has been to model the manoeuvres by the almost independent acceleration and turn rate, and to introduce a stochastic continuoustime model for the manoeuvres. We assumed that in any given sampling interval there are three different alternatives: no manoeuvre, change of speed, or change of course. This lead to a bank of 31 Kalman filters.

We have

P(O)

( (

0.01

In Figure 1 the true trajectory together with the estimated positions are depicted. The estimates of velocity and heading are given in Figure 2, and the estimates of acceleration and turn rate in Figure 3. As can be seen the estimates are very good, in particular the turn rate.

where

Xo

Q

( 0.98

'

845

1104 3.5.-~-~-~-~~-~-~--,

....... . ... ··· ··· . ... ...

2.5

Turn ra&e

::6: :EJ

1 ••••••• ~~ • • • • • • • • • , • • • • • • • • • • • •

0.5

A

o

45

3.5

B

SO

100

150

200

2SO

300

350

110'"

Fig. 3: Estimates of acceleration and turn rate. Manoeuvre detections marked by dashed lines

Fig. 1: True trajectory ( marked" 0") and estimated trajectory (marked "*")

:": : Kl I:GJ:J

Andersson P. (1985). Adaptive forgetting in recursive identification through multiple models. Int. J. Contr., 42:1175-1193. Bar-Shalom Y. and T.E. Fortmann (1988). Tracking and Data Association. Mathematics in Science and Engineering, Vo!. 179, Academic Press, Boston. Blom H.A. and Y. Bar-Shalom (1988). The Interacting Multiple Model Algorithm for Systems with Markovian Switching Coefficients. IEEE Trans. on Automatic Control, 33(8):780783. Chang C.B., R.H. Whiting, and M. Athans (1977). On the state and parameter estimation for maneuvering reentry vehicles. IEEE Trans. on Automatic Control, 22:99-105. Bar-Shalom Y. (Ed) (1991). Multi-target Multi-sensor tracking. Artech House. Gholson N.H. and R.L. Moose (1977). Maneuvering target tracking using adaptive state estimation. IEEE Trans. on Aerospace and Electr. Systems, 13:310-317. Gustafsson F. (1992). Estimation of discrete parameters in linear systems. PhD thesis, Linkoping University, Linkoping, Sweden. Jazwinski A.H. (1970). Stochastic Processes and Filtering Theory. Academic Press, N.Y. Moose R.L. (1975). An adaptive state estimation solution to the maneuvering target problem. IEEE Trans. on Automatic Control, 20:359-362. Ricker G.G. and J.R. Williams (1978). Adaptive Tracking Filter for Maneuvering Targets. IEEE Trans. on Aerospace and Electr. Systems, 14(1):185-193. Thorp J.S. (1973). Optimal tracking of maneuvering targets. IEEE Trans. on Aerospace and Electr. Systems, 9:512-519. Tugnait J.K. and A.H. Haddad (1979). A Detection-Estimation Scheme for State Estimation in Switching Environments. Automatica, 15:477-481. Vacher P., I. Barret, and M. Gauvrit (1991). Design of a tracking algorithm for an advanced ATC system. In Y. Bar-Shalom, editor, Multitarget Multi-sensor tracking. Artech House.

~E 1000

so

100

ISO

200

2SO

300

3SO

300

3SO

.~ o

SO

100

ISO

200

2SO

Fig. 2: Estimates of velocity and heading. Manoeuvre detections marked by dashed lines

To reduce the number of filters we used a pruning algorithm, which worked very well. In fact we are a bit confused about reports that pruning does not work, and suspect that the reason may have been use of a less accurate manoeuvre model. To each filter is associated a particular manoeuvre sequence, enabling simultaneous state estimation and manoeuvre detection. The algorithm was demonstrated in a simulation example giving very good signal estimates, and accurate manoeuvre detection.

Acknowledgements The first author would like to thank Professors Rob Evans and O.K. Kwon, and Drs. Iven Mareels and Rick Middleton for many valuable comments and interesting discussions.

REFERENCES Anderson B.D.O. and J.B. Moore (1979). Optimal Filtering. Prentice-Hall, Englewood Cliffs, NJ.

846