Evaluation of Tracking Filters for Tracking of Manoeuvring Targets

(1)

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2020

Evaluation of Tracking

Filters for Tracking of

Manoeuvring Targets

Ludvig Junler

(2)

Ludvig Junler LiTH-ISY-EX--20/5322--SE

Supervisor: Robin Forsling

isy_{, Linköpings universitet}

Johan Rasmusson

Saab AB

Examiner: Gustaf Hendeby

isy_{, Linköpings universitet}

Division of Automatic Control Department of Electrical Engineering

(3)

Abstract

This thesis evaluates different solutions to the target tracking problem with the use of airborne radar measurements. The purpose of this report is to present and compare options that can improve the tracking performance when the target is performing various manoeuvres while the radar measurements are noisy. A simulation study is done to evaluate and compare the presented solutions, where the evaluating criteria are the estimation errors and the computational complex-ity. The algorithms investigated are the general pseudo Bayesian of order one (gpb(1)) filter and the interacting multiple model (imm) filter, each using three motion models, along with several single model Kalman filters. Additionally, the impact on the tracking performance by different choices of radar parameters is also examined.

The results show that filters using multiple models are best suited for tracking targets performing different manoeuvres. The tracking performance is improved with both the gpb(1) and imm algorithms compared to the filters using a sin-gle model. Looking at the estimation errors, imm outperforms the other algo-rithms and achieves a better general performance for different kinds of manoeu-vres. However, imm have a much higher computational complexity than the fil-ters with a single model. gpb(1) could therefore be more suited for applications where computational power poses a problem, since it is less computationally de-manding than imm.

Furthermore, it is shown that different radar parameters have an impact on the tracking performance. The choice of pulse repetition frequency (prf) and duty cycle used by the radar affects the accuracy of the measurements. The estimation errors of the tracking filters become larger with poor measurements, which also makes it more difficult for the multiple model algorithms to make good use of the different motion models. In most cases, imm is however less sensitive to the choice of prf in relation to how the models are used in the algorithm, compared to gpb(1). Nevertheless, the study shows that there are cases where some combi-nations of radar parameters drastically reduces the tracking performance and no clear improvement can be seen, not even for imm.

(4)

(5)

Acknowledgments

I would like to thank my supervisor at Saab, Johan, for all the valuable help and support throughout these months. I also want to express my gratitude to Saab for letting me do an interesting project at an exciting workplace.

Finally, a special thanks to my supervisor Robin and examiner Gustaf at the uni-versity, for their guidance and feedback on the report, which helped me improve my work.

Linköping, June 2020 Ludvig Junler

(6)

(7)

Notation

Abbreviations Abbreviation Meaning ca Constant acceleration cv Constant velocity ct Coordinated turn

crlb Cramér-Rao lower bound

ekf Extended Kalman filter

gpb Generalised pseudo Bayesian

imm Interacting multiple model

mc _{Monte Carlo}

prf _{Pulse repetition frequency}

pri _{Pulse repetition interval}

rcs _{Radar cross section}

rgpo Range gate pull-off

rms Root mean square

rmse Root mean square error

snr Signal to noise ratio

(10)

Notations

Notation Meaning

θ Azimuth angle

P Covariance matrix of state vector

dt Duty cycle

φ Elevation angle

ˆ

x Estimated state vector

θ3dB Half power antenna beam width

PD Probability of detection

PFA Probability of false alarm

τ Pulse width

σ Radar cross section

r Radial distance to target

(11)

1

Introduction

Airborne radars are often used in military applications to track targets, such as other aircraft. In many situations it is important to achieve a good tracking per-formance since this can help the pilot to make better decision about how to act in the future, for example, if to engage the target or not.

This thesis investigates the problem of tracking a target performing different ma-noeuvres, using measurements generated by a radar. Different algorithms are compared while also studying how properties of the radar affects the tracking performance. The work was carried out at Saab AB in Linköping.

In this chapter a background to the problem is given together with what this thesis aims to solve. An overview of the report is also presented in Section 1.5.

1.1 Background

A radar is used to retrieve information about the position and velocity of a target. The measurements from a radar are not perfect and the size of the error depends, among other things, on the distance to the target. For many applications it is desirable to get more exact information on the target than the measurements can provide on their own. Improving the tracking performance is therefore of great value.

The basic concept of a radar is that it sends out an electromagnetic signal. If the signal is reflected by a target it is possible for the radar to detect it. This gives a measurement of the target’s position. Taking advantage of the Doppler effect the radar can also produce a measurement of the target’s radial velocity. Since the radar measurements are not perfect they will have an uncertainty. A standard

(12)

solution for target tracking with airborne radar is therefore to use a Kalman filter (kf). Measurements of position and velocity are then used to estimate the states in the kf.

Using a kf may however be an insufficient solution for cases where the target is performing heavy manoeuvres. One filter might trace a target with great accu-racy for one type of motion, for example, when it is moving with constant velocity. However, if the target starts to accelerate or turn, the uncertainty of the model would have to be increased in the filter to follow the manoeuvres. The tracking performance will then get worse for the non-manoeuvring case. In real life ap-plications the target often performs many different manoeuvres which makes the target tracking more difficult. It is therefore of interest to find good solutions to the tracking problem when there is an uncertainty in the radar measurements and the target is performing different types of manoeuvres.

1.2 Problem Formulation

The goal of this master thesis is to improve tracking of a single target, perform-ing various manoeuvres, by evaluatperform-ing and comparperform-ing different solutions to the problem. The tracking will be based on measurements generated by a radar. The aim is to find answers to the following questions:

• Which tracking algorithms are suitable for tracking of manoeuvring tar-gets?

• Which of these algorithms is the best alternative when focusing on tracking performance and computational complexity?

• How is the performance affected by different properties of the radar?

1.3 Method

To answer the questions stated in Section 1.2, different methods for target track-ing are evaluated in a simulation study. This is done ustrack-ing a radar model im-plemented in Matlab®

to generate measurements, which are then used by the tracking algorithms. By implementing different functions found in real radar applications, the goal is to mimic the behaviour of a real airborne radar when producing the measurements. This makes it possible to change parameters in the radar model to investigate how some of them affect the tracking performance. Furthermore, three different trajectories for a target are used to evaluate the track-ing performance for different kinds of manoeuvres.

(13)

1.4 Limitations 3

1.4 Limitations

To match the scope of a master’s thesis, only tracking of a single target is consid-ered. This is due to the fact that tracking multiple targets with good accuracy can get quite complicated in many cases. It allows the focus to be on improving the tracking performance for different manoeuvres and choices of radar parameters instead.

Furthermore, the number of different methods investigated is limited to the three most promising. These methods are widely used in target tracking applications and they are described in a wide range of literature. This makes them good can-didates to the problem studied.

Finally, the radar is assumed to operate in a clutter free environment where the measurements originates from an actual target. Real radars used in aircraft often have functions to reject clutter in the measurements directly, which motivates this choice.

1.5 Report Overview

Chapters 2 and 3 present a theoretical background on radars and different filter algorithms. Based on this theory, Chapter 4 describes how the simulation study is done by implementing a radar model and filter algorithms. The trajectories used to test the filters are also presented here. Chapter 5 presents and discusses the results of the different solutions, followed by conclusions in Chapter 6.

(14)

(15)

2

Radar Theory

This chapter presents a theoretical background to radars. The theory is later used to implement a realistic model in Matlab®

, making it possible to evaluate how different radar parameters affect the tracking performance.

2.1 Measurements and Coordinates

A radar produces measurements of the position and the range rate of a target in spherical coordinates (r, ˙r, θ, φ), where r is the radial distance to the target and ˙r is the range rate between the radar source and the target [1]. The angle

θ, commonly referred to as the azimuth angle, is the angle of the target in the xy-plane from the x-axis. The elevation angle, denoted φ, is the angle towards

the z-axis. An illustration of the situation is shown in Figure 2.1, where the radar source is located in the origin.

Figure 2.1:Radar measurements.

(16)

The transformation from spherical to Cartesian coordinates is given by

x = r cos(θ) cos(φ), y = r sin(θ) cos(φ), z = r sin(φ).

(2.1)

Transforming the other way around, from Cartesian to spherical coordinates, is done by r = q x2_{+ y}2_{+ z}2_, ˙r = xvx+ yvy+ zvz px2_{+ y}2_{+ z}2 , θ = tan−1y x, φ = tan−1 z (x2_{+ y}2₎1/2 ! , (2.2)

where vi are the velocity components in the x−, y− and z−direction.

Taking advantage of the Doppler effect, the radar can measure the range rate, ˙r, directly. The Doppler frequency shift is a result of the relative motion between the radar and the target, as a result of either the radar source or the object moving, or both [17, p. 257]. Therefore the range rate consists of the velocity of the radar source, VR, and the velocity of the target, VT, projected onto the line of sight to

the target [17, p. 263] according to

˙r = −(V1+ V2). (2.3)

where the definitions of V1 and V2 are found in Figure 2.2. The velocity V2 is

defined to be positive if target is approaching the radar. For example, if the radar travels with V1= 300 m/s and the target is approaching the radar with V2= 400

m/s, the measured range rate is ˙r = −(300 + 400) = −700 m/s.

Figure 2.2:Range rate, ˙r, measured by the radar. The dotted line is the line of sight to the target.

(17)

2.2 Radar Waveforms 7

2.2 Radar Waveforms

There are two general types of radar, pulsed and continuous wave [17, p. 149]. The focus here is only radars performing pulsed operations, meaning that radio waves are transmitted in pulses and the radar then listens for the effects of echo between the pulses.

2.2.1 Pulsed Operation

The parameters describing the pulsed operation are: [18, p. 5]

• The carrier frequency fc (Hz), which is related to the wavelength, λ (m),

according to: λ = c/fc.

• The pulse width τ (s).

• The pulse repetition interval pri (s).

• The pulse repetition frequency prf (Hz), which is the inverse of pri. A general waveform with the parameters mentioned is depicted in Figure 2.3.

Figure 2.3:Radar waveform.

The power of each transmitted rectangular pulse is P τ, where P is the peak power of the pulse, which is the output power when the transmitter is on. This means that the average power transmitted between two pulses is [17, p. 154]

Pavg= P τ

pri= P τprf, (2.4)

where the termpriτ often is referred to as the duty cycle, dt.

2.2.2 Measurement Ambiguities

The choice of prf is very important when designing a radar. It determines to what extent measurements of range and Doppler frequencies will be ambiguous. For range measurements to be unambiguous, each pulse sent out by the radar must be returned before the next pulse is transmitted. This implies that a low

(18)

prf_{has to be used to be able to measure the distance to a target far away.}

The ambiguity of Doppler measurements depends on the maximum velocity of an approaching target and the velocity of the aircraft [17, p. 380]. The Doppler frequency shift is measured with a bank of narrowband filters, where each filter has a narrow passband. A target’s Doppler frequency is then determined by ob-serving in which of the filters in the filter bank the target appears. However, the spectrum of the received signal will be periodically repeated, which is always the case for discrete-time filters. The spectrum will reappear in sidebands separated by the prf.

If the prf is smaller than the maximum Doppler shift, there is no way of telling if the observed shift lies in one of the side bands, or even in which one of the side bands. This is illustrated by Figure 2.4.

Figure 2.4: With a prf lower than the incoming Doppler frequencies the placement of the filter bank passband, the green stripe, does not matter. It will not be possible to determine the true Doppler shift [17, p. 322].

This implies that a high prf has to be used to resolve targets moving fast.

2.2.3 Pulse Repetition Frequency

When designing a radar, three basic categories of prf are often used. These cate-gories are low, medium and high [17, p. 383]. While a low prf is good for range measurements at large distances, a high prf is good for range rate measurements when a target is moving fast. As a result of this conflict, a compromise sometimes has to be made. This can be done by choosing a medium prf. Both range- and range rate measurements can however be ambiguous for a medium prf.

The classification of the prf depends on the operating conditions [17, p. 387]. For example, a low prf implies that range measurements are unambiguous for the re-quired operating range. Furthermore, the maximum velocities of the radar and all targets affects the maximum unambiguous Doppler shift. Therefore a high prfimplies that all Doppler measurements are unambiguous for the specific op-erating conditions.

(19)

2.3 Radar Range Equation 9

2.2.4 Resolving Ambiguities

To deal with ambiguities, multiple predetermined prfs can be used [18, p. 132]. Transmitting pulse trains with different prfs will change the ambiguous range. The radar preserves the phase of the waveform during transmission and recep-tion of the pulse train, a property called coherence. This makes it possible to match the received signals, using a coincidence filter, to determine the true range to a target, as shown in Figure 2.5.

Figure 2.5:Three different prfs used to solve range ambiguities [18, p. 132]. The same concept can be used to resolve Doppler ambiguities as well. The side-bands of the the Doppler shift will change position when changing the prf, and thus the true Doppler shift can be determined.

2.3 Radar Range Equation

The radar range equation is of great importance when describing radar perfor-mance. It specifies how the received signal power can be related to the maximum detection range when tracking an object.

2.3.1 Received Signal Power

A certain amount of the power transmitted by a radar will be reflected back by an object. One version of the radar range equation is [17, p. 168]

Pr =

PavgGσ Ae

(4π)2_r4 , (2.5)

where:

• Pr= received signal power (W)

• Pavg= average transmitted power (W)

• G = transmitter antenna gain (unitless) • σ = radar cross section (m2)

(20)

• Ae= effective area of the receiver antenna (m2)

• r = distance between radar and object (m)

The radar cross section (rcs), denoted σ in (2.5), is used to describe how much of the transmitted power that will be reflected by a certain target [17, p. 166]. It de-pends on the geometry of the target and how the transmitted waves are scattered on its surface.

The antenna gain is related to the effective area by [18, p. 9]

G = 4πAe

λ2 =⇒ Ae=

Gλ2

4π , (2.6)

where λ is the wavelength of the transmitted signal. By inserting this into (2.5) another version of the received signal power is obtained

Pr =

PavgG2λ2σ

(4π)3_r4 . (2.7)

2.3.2 Receiver Noise

There is always electrical noise of random amplitude and frequency affecting the received signal [18, p. 9]. A great deal of the noise interfering with the frequen-cies used by most radars is generated within the receiver itself. The noise perfor-mance of a receiver is often described by a quantity called the noise figure F [17, p. 160]. This determines the mean noise power of a receiver, which is also related to the temperature of the resistor connected across the input terminals and the bandwidth of the receiver. The mean noise power of a receiver is therefore given by [17, p. 162]

Pn = kT0BF, (2.8)

where:

• Pn= receiver thermal noise (W)

• k = Boltzmann’s constant (1.38 × 10−23Ws/K) • T0= receiver temperature (K)

• B = receiver filter bandwidth (Hz) • F = receiver noise figure (unitless)

2.3.3 Signal to Noise Ratio

Equations (2.7) and (2.8) can be used to express the ratio between the received signal power and the receiver noise [18, p. 10] according to

Pr Pn = PavgG 2_λ2_σ (4π)3_r4_FkT 0B . (2.9)

(21)

2.4 Volume Search 11

The quantity Pr/Pn is referred to as the signal to noise ratio (snr) and is often

expressed in decibels. To account for other system losses, such as absorption and scattering in the atmosphere, transmission loss and the target not being centered in the path of the antenna beam, the factor Lscan be included in the equation [18,

p. 10]. If all system losses are combined in Ls we get

snr= PavgG 2_λ2_σ (4π)3_r4_FkT 0BLs . (2.10)

2.4 Volume Search

To cover a greater search volume, a number of beams can be transmitted by the radar antenna in different directions. The search pattern can for instance be de-fined by a rectangle with a certain amount of rows and columns as in Figure 2.6.

Figure 2.6:Example of search pattern.

The search pattern will then cover the whole search volume displayed in Figure 2.7.

Figure 2.7: Radar search volume. The dashed lines is the direction of the antenna beam, called the radar boresight.

The beamwidth is defined as the angle between two opposite edges of the beam. Because the strength of the beam decreases as the angle from the centre increases, one must decide what are considered to be the edges [17, p. 112]. It is common to define them as the points where the transmitted signal strength has dropped by half its power, which is about −3dB in a logarithmic scale. The beamwidth is

(22)

then referred to as θ3dB. It is shown in Figure 2.8 and can be calculated by [18,

p. 24]

θ3dB =

0.88λ

D radians, (2.11)

where D is the antenna length.

Figure 2.8:Beamwidth, θ3dB, where the power has dropped by 3 dB.

To cover the whole search volume, the number of beam positions is therefore given by

Number of beam positions = Nb≈

αβ

(θ3dB)2

, (2.12)

where α is the azimuth angle and β is the elevation angle of the full search vol-ume. The time at each beam position, called the dwell time, is then given by

Dwell time = tf

Nb

= tf

(θ3dB)2

αβ , (2.13)

where the frame time, tf, is the time period for one iteration of the volume search

[17, p. 185].

2.5 Probability of Detection

One problem when it comes to tracking is to determine a detection threshold on the snr. It should be high enough so the noise and clutter rarely exceeds the limit, but low enough so the target reflection signals rarely falls below it [18, ch.3]. If the threshold is too high there would be many missed detections and if it is to low there is a great risk of falsely detecting a target.

2.5.1 Binary Hypothesis Test

During the detection process, a hypothesis test is made with a null hypothesis, H₀_{, and an alternative hypothesis, H}₁_{. This is known as a binary hypothesis test} [9, p. 60]. An example is illustrated in Figure 2.9, where we observe a normal distributed random variable which can have one of two different means. Based on received data, the hypothesis test is then to determine which is the correct mean. A reasonable approach would be to select a threshold and decide that H1

(23)

2.5 Probability of Detection 13

is true if the data exceeds the threshold, otherwise H0 is decided to be true. For

a binary detection test, two errors can be made. Deciding H1 when H0is true is

called a Type I error and deciding H0 when H1is true is called a Type I I error.

As seen in Figure 2.9, the probability of making these errors depends on how the threshold is selected.

Figure 2.9:Probability density functions of H0and H1for a binary

hypoth-esis test [9, p. 62]. The threshold is in this case selected to be the value where the two distributions intersect.

For a radar application, the goal is to decide if a target is present or not, based on a measurement. This would mean that the null hypothesis is that only noise is present and the alternative hypothesis is that both a signal and noise is present [9, p. 60]. Some probabilities used when describing the detection process are:

• Probability of detection, PD – The probability that a target is detected

when it is present

• Probability of false alarm, PFA – The probability that a target is detected

when it is in fact not present. This corresponds to a Type I error

• Probability of miss, PM– The probability that a target is not detected when

it is present. This corresponds to a Type I I error

Both types of errors, PFAand PM, can not be reduced at the same time and

there-fore a trade-off has to be made [9, p. 62]. A typical approach when deciding the threshold during the detection process is to keep PFAconstant and then select a

threshold so that PDis maximised.

2.5.2 Signal and Noise Probability Distributions

The thermal noise in a receiver is usually assumed to have a Gaussian distribution with zero mean and mean noise power σ2= kT0BF, as described in Section 2.3.2.

The probability that the magnitude of the noise will have a given value, after being passed through the narrowband filter used by the radar, will then have a

(24)

Rayleigh distribution [17, p. 188].

The signal, without noise, will have different statistics with just an impulse at the given signal strength as probability density function. Adding the signal to the noise will result in a Rayleigh distribution shifted to the right by the signal amplitude, as shown in Figure 2.10 [18, p. 52].

Figure 2.10:Probability density function for signal and noise. The probabil-ity of detection is represented by the grey area under the S + N -curve.

The noise was earlier described as Gaussian with zero mean and the signal-plus-noise has a mean at the signal power. Under the assumption that the signal is much stronger than the noise, meaning that the signal-plus-noise distribution will be far away from zero, it has been shown that the Rayleigh distribution can be treated as Gaussian [18, p. 53]. For a radar to operate effectively the signal must be much stronger than the noise, which motivates the choice of treating the distribution as Gaussian. This is also described in [14, p. 261].

2.5.3 Calculating Probabilities

The average time between two false alarms can be calculated using

Tf a= 1

PFAB

, (2.14)

where B (Hz) is the receiver bandwidth [16, p. 41]. It is desirable to get as few false alarms as possible. This is however in conflict with getting good detection performance and therefore PFA is selected to achieve a certain false alarm rate

when designing the radar.

For a given PFAand snr, the probability of detection can be calculated. Because

the signal-plus-noise probability density function could be treated as Gaussian, the PD is found by integrating the Gaussian distribution [18, p. 54]. The integral,

(25)

2.6 Measurement Accuracy 15

which does not result in a closed-form solution, can be approximated using the complementary error function, given by

erfc(T ) = 1 −√2 π T Z 0 e−x2dx. (2.15)

The probability of detection is then calculated as

PD ≈0.5 erfc

p

−_{ln P}_FA− √

snr+ 0.5. (2.16)

2.5.4 Improving Detection Performance

To increase the snr, and thus increasing the chances to detect a target, the num-ber of pulses transmitted by the radar can be increased [17, p. 172]. The radar makes use of coherence when transmitting the pulses so that the phase is con-stant between pulses, which is needed to measure the Doppler shift. The phase and amplitude of the noise will however vary randomly. This means that the noise, when integrating multiple pulses, will mostly cancel out, while the co-herent signal will add up to a greater signal strength. The snr after a burst of multiple coherent pulses is therefore given by [16, p. 46]

snr_n

p = npsnr1, (2.17)

where snr1is the single pulse snr and snrnp is the snr after integration of np

number of coherent pulses.

One way of increasing the detection performance even further is to use multi-ple individual scans or measurements. A common approach is to use multimulti-ple pulse bursts and then require a certain amount of detections out of the individ-ual bursts, called M-out-of-N detection. The PD for each individual burst is first

computed and then the binomial probability theorem is used to calculate the cu-mulative probability given by [18, p. 63]

P (M, N , p) = N ! M(N − M)!p

M_{(1 − p)}N −M_, _(2.18)

where P (N , M, p) is the probability of crossing the detection threshold M times out of N attempts and p is the probability of crossing the detection threshold for a single attempt. The cumulative probability of false alarm can be calculated in the same way by replacing p with PFA. This makes it possible to get a much lower

PFAwhile maintaining almost the same detection probabilities.

2.6 Measurement Accuracy

In reality, measurements from the radar are not perfect. There is a limited res-olution to the measurements which will affect how accurate they are. The mea-surement errors are random and will have a standard deviation which depends

(26)

on the snr [4, p. 166].

Range is computed based on the measured time delay between transmitting and receiving a pulse. To separate two targets they need to be some distance apart from each other, which depends on the used waveform [18, p. 108]. The range rate on the other hand is acquired by measuring the Doppler shift. For a pulse burst waveform, the Doppler filter bandwidth, ∆fd=

prf

Number of pulses, determines

the accuracy [18, p. 127]. Finally, the resolution of angle measurements depends on the beam width, θ3dB. To resolve two targets in angle they need to be separated

by at least one beam width, meaning that θ3dBwill decide the angle measurement

accuracy.

Based on the resolution, the accuracy is the statistical measure of how well the measurement can be estimated. The standard deviation of the different measure-ments are given by Table 2.1 [18, p. 133].

Table 2.1:Standard deviation of radar measurements.

Parameter Relationship

Range accuracy, δr cτ

2 √

2snr

Range rate accuracy, δ ˙r ∆fdλ

2 √

2snr

Angle accuracy, δθ √θ3dB

(27)

3

Target Tracking Theory

This chapter presents the Kalman filter algorithm, motion models that can be used in the filters and strategies for target tracking using multiple models.

3.1 State-Space Representation

A model describing the motion of a target and the measurements related to it can be represented as a state-space system. The conversion between radar measure-ments and Cartesian coordinates is a nonlinear operation as described in Section 2.1. A general nonlinear system can be represented in discrete time as

xk+1= f (xk, uk, vk), (3.1a)

yk = h(xk, uk, ek), (3.1b)

where xkis the state vector, vkis the state noise, ykis the sensor measurements, ek

is the measurement noise, ukis the input signal and k is the time index [6, p. 125].

A special case is when the noise is additive and the input is unknown, which then gives a system in the form

xk+1= f (xk) + vk, Cov(vk) = Qk, (3.2a)

yk = h(xk) + ek, Cov(ek) = Rk. (3.2b)

In the case studied in this report, yk is the radar measurements which are given

in spherical coordinates. The measurement equation, h(xk), is then the

trans-formation from Cartesian to spherical coordinates given by (2.2) and ek is the

measurement errors described in Section 2.6.

(28)

3.2 Motion Models

In target tracking applications it is important to have a good model describing the motion of the target, to be able to estimate the state of the object. In most cases of target tracking the input uk is usually not known, and many kinematic

models therefore model the input as process noise, e.g., [6, p. 343]

xk+1= Fkxk+ Gkwk, (3.3)

where xk = x(tk). The noise wk is a white noise process with process noise

vari-ance (σ2

w)i, i = x, y, z, in each coordinate [11].

Many different motion models exist, and a lot of research has been done in the area. The following sections present some of the most commonly used models in the area of target tracking.

3.2.1 Constant Velocity

For a non-manoeuvring target, with velocity components v = ( ˙x, ˙y, ˙z), the mo-tion is described by vk+1 = vk+ T wk ≈vk, meaning that the target has a nearly

constant velocity. The small white acceleration noise, wk, accounts for

unpre-dictable modelling errors [12, p. 1335]. With states x =x v_x y v_y z v_zT the discrete cv model, in the form (3.3), is given by1

Fcv= diag[ ˜Fcv, ˜Fcv, ˜Fcv], (3.4a) Gcv= diag[ ˜Gcv, ˜Gcv, ˜Gcv], (3.4b) ˜ Fcv="1 T₀ ₁ # , G˜cv ="T 2_/2 T # , (3.4c) with

Qcv_k = cov(Gcvwk) = diag[(σw2)xQ˜cv, (σw2)yQ˜cv, (σw2)zQ˜cv], (3.5a)

˜ Qcv="T 4_/3 _T3_/2 T3/2 T2 # . (3.5b)

3.2.2 Constant Acceleration

The simplest model describing a manoeuvring target is the constant acceleration model. In fact, the acceleration is nearly constant and modelled as a white-noise

1_{The notation A = diag[A}

1, A2, A2] means that A is a block-diagonal matrix on the form         A1 0n 0n 0n A2 0n 0n 0n A3        

, where Ai are matrices of size n × n, 0nare matrices with zeroes of the same size and n is the number of states in each dimension.

(29)

3.2 Motion Models 19

process. With acceleration components a = ( ˙vx, ˙vy, ˙vz) the motion is described

by ak+1= ak+ T wk≈ak, where wkis small jerk noise [12, p. 1336].

With states x =x v_x a_x y v_y a_y z v_z a_zT the discrete ca model, in the form (3.3), is given by

Fca= diag[ ˜Fca, ˜Fca, ˜Fca], (3.6a) Gca= diag[ ˜Gca, ˜Gca, ˜Gca], (3.6b) ˜ Fca=         1 T T2/2 0 1 T 0 0 1         , G˜ca=         T2/2 T 1         , (3.6c) with

Qca_k = cov(Gcawk) = diag[(σw2)xQ˜ca, (σw2)yQ˜ca, (σw2)zQ˜ca], (3.7a)

˜ Qca=         T4/4 T3/2 T2/2 T3/2 T2 T T2/2 T 1         . (3.7b)

3.2.3 Coordinated Turn

Another way of representing a target manoeuvre is to use a turn model. There are many different versions of turn models and one option is to assume that the target has a nearly constant speed and a nearly constant turn rate, ω. Deriving a turn model in 3D can get quite complicated. However, a simplified model is the planar ct model, presented in [12, p. 1353].

With states x =x vx ax y vy ay z vz az

T

the ct model is given by

Fct = diag[ ˜Fct, ˜Fct, ˜Fct], (3.8a) ˜ Fct(ω) =          1 sin ωT_ω 1−cos ωT_ω2 0 cos ωT sin ωT_ω 0 −_{ω sin ωT} _{cos ωT}          , (3.8b) with Qct_k = diag[(σ_w2)xQ˜ct(ω), (σw2)yQ˜ct(ω), (σw2)zQ˜ct(ω)], (3.9a) ˜ Qct(ω) =             6ωT −8 sin ωT +sin 2ωT 4ω5 2 sin4_{(ωT /2)} ω4 −_{2ωT +4 sin ωT −sin 2ωT} 4ω3 2 sin4(ωT /2) ω4 2ωT −sin 2ωT 4ω3 sin 2_ωT 2ω2 −2ωT +4 sin ωT −sin 2ωT 4ω3 sin 2_ωT 2ω2 2ωT +sin 2ωT4ω             . (3.9b)

(30)

3.3 Hybrid System

For a target performing different manoeuvres, multiple models can be used to describe the motion. For the multiple model (mm) approach to target tracking, the state process is described as a hybrid process, meaning that it has both con-tinuous and discrete components [13, p. 1257].

The base state, x, is often continuous while the system mode, s, is discrete. Fur-thermore, the set of possible modes is called the mode space and is denoted by S. A system mode refers to a pattern of behaviour, or the structure of a system, and has the property that it either jumps or stays unchanged. The problem of hybrid estimation, and the mm approach, is then to estimate x and s.

A simple discrete-time hybrid system is the Markov jump-linear system (MJLS), which is given by

xk+1= Fk(sk+1)xk+ Gk(sk+1)wk(sk+1), (3.10a)

yk = Hk(sk)xk+ vk(sk), (3.10b)

where s is a Markov chain. This means that P {s(j)_k+1|_s(i)

k }= pij,k, where s

(i)

k denotes

that mode i is in effect at time k and pij,kis the probability of transitioning from

mode i to j.

The basic idea of the mm approach is to use a set of models, M, as candidates to the true system mode in affect at the time, and then use a bank of filters with the different models [13, p. 1258]. The models can be seen as an approximation of the true modes. There are often more modes than models, meaning that M has fewer elements than S.

Each model in M obeys the state space representation (3.3) and can for example be a cv, ca or ct model. Additionally, the occurring jumps between the models has the transition probabilities

P {m(j)_k+1|_m(i)

k }= πij, (3.11)

where m(i)_k means that the model at time k matches the system mode.

3.4 The EKF Algorithm

The Extended Kalman filter is a nonlinear version of the Kalman filter, used to estimate the state vector of a nonlinear system. The filter recursion consists of a measurement update and a time update. The time update predicts the state vector and the measurement update corrects the predicted states.

(31)

3.4 The EKF Algorithm 21

Given a nonlinear system on the form (3.1), the filter provides the estimated state vector ˆx together with the covariance matrix of the states, P . The state vector is

initialised using ˆx1|0 with covariance P1|0. The time index k|m should be

inter-preted as at time k, given measurements up to time m [6, p. 153].

When a measurement is available, the filter performs a measurement update re-sulting in ˆxk|k and Pk|k. It is computed based on the prediction error, which is

the difference between the actual measurement and the measurement equation evaluated with the predicted states. This can also be called the innovation. A time update is performed at every time instance, even if no measurements are available, resulting in a prediction of the state vector, ˆxk+1|k, and its covariance,

Pk+1|k.

The measurement and time update are computed with the Jacobians of the state and measurement equations, defined as f_k0 = ∂fk(x)

∂x _x=x k and h0_k = ∂hk(x) ∂x _x=x k . All steps in the ekf algorithm are given by Algorithm 1 [6, p. 197].

(32)

3.5 Multiple Model Filters

There are some choices to be made when implementing an mm estimator. The general structure is based on the following parts: [13, p. 1259]

• Model-set determination – The choice of models is very important in terms of performance. They should be selected so that the system modes are cov-ered, meaning that the behaviour of a target can be captured for different manoeuvres.

• Cooperation strategy – Strategies must be used to determine which model is most likely at different time instants. This also includes merging of esti-mates with similar model sequences and pruning of model sequences that are unlikely.

• Conditional filtering – The conditional filtering is the state estimation of the hybrid process conditioned on a mode sequence. This is done in the same way as for an ordinary discrete system.

• Output processing – Based on measurements and state estimations from all filters, the information has to be fused to create an overall estimate that is as good as possible.

An overview of the structure, using two models, is shown in Figure 3.1. The arrows indicate how the models interact within one recursion of the algorithm.

Figure 3.1:General structure of mm algorithms [13, p. 1259]

As time increases, so will the possible sequences of different models. At time k, there are Mk _{possible model sequences, where M = |M|. This means that there}

will be a rapidly growing tree, as depicted in Figure 3.2 [6, p. 265]. Different cooperation strategies have been developed to deal with this difficulty, by ap-proximating the large tree with a simpler one. This is done by merging branches that are similar to each other and pruning unlikely ones.

(33)

3.5 Multiple Model Filters 23

Figure 3.2:Tree of model sequences for two modes.

3.5.1 Generalised Pseudo Bayesian

The Generalised Pseudo Bayesian algorithm of order n, gpb(n), is one example of how this can be done. The idea is to reduce the tree by merging hypothesis that are the same in the latest n time steps [6]. With two models and a time horizon of n = 1, branches (1, 3, 5, 7) in Figure 3.2 would be merged to one. Remaining branches (2, 4, 6, 8) would then be merged to a second branch. If the time horizon is instead set to n = 2, branches (1, 5), (2, 6), (3, 7) and (4, 8) would be merged together. However, only gpb of order n = 1 will be considered here.

For each iteration of gpb(1), the filters are reinitialised with the previous overall estimate, as depicted in Figure 3.3. The different estimates ˆx(i)_k|k = E[xk|yk, m

(i)

k ]

with mode probabilities µ(i)_k = P {m(i)_k |_y_k}_{, where i denotes the model, are then} lumped together to form a new single estimate [13, p. 1272].

Figure 3.3:Reinitialisations of the gpb(1) algorithm [13, p. 1274].

(34)

Algorithm 2One cycle of the gpb(1) algorithm.

1. Model-conditioned reinitialisation(for i = 1, 2, ..., M): Predicted mode probability: µ(i)_k|k−1, P {m(i)k |yk−1}=

P

jπjiµ

(j)

k−1

2. Model-conditioned filtering using the ekf algorithm(for i = 1, 2, ..., M): Predicted state: ˆx(i)_k|k−1= F_k−1(i) xˆ(i)_k−1|k−1+ G(i)_k−1w¯(i)_k−1

Predicted covariance: P_k|k−1(i) = F_k−1(i) Pˆ_k−1|k−1(i) (F_k−1(i) )0

+ G(i)_k−1Q(i)_k−1(G(i)_k−1)0

Predicted measurement: ˆy_k|k−1(i) = H_k(i)xˆ(i)_k|k−1+ ˆv_k(i)

Innovation: ε(i)_k = yk−yˆ

(i)

k|k−1

Innovation covariance: S_k(i)= H_k(i)P_k|k−1(i) (H_k(i))0+ R(i)_k Kalman gain: K_k(i)= P_k|k−1(i) (H_k(i))0

(S_k(i))−₁

Updated state: ˆx(i)_k|k= ˆx(i)_k|k−1+ K_k(i)ε(i)_k

Updated covariance: P_k|k(i)= P_k|k−1(i) −_K(i)

k S (i) k (K (i) k ) 0

3. Mode probability update(for i = 1, 2, ..., M): Model likelihood:

L(i)_k , p[ε(i)_k |_m(i)

k , yk−1] assume = N_(ε(i) k ; 0, S (i) k ) = exp −1 2(ε (i) k ) 0 (S_k(i))−1_(ε(i) k ) 2πS (i) k 1/2

Mode probability: µ(i)_k|k= µ (i) k|k−1L (i) k P jµ (j) k|k−1L (j) k

3.5.2 Interacting Multiple Model

An improved version of gpb is the Interacting multiple model (imm) algorithm. It is basically the same as gpb(2), but with some calculations rearranged. The filters are not reinitialised with a single overall estimate for the imm and gpb(2) algorithms. Instead, each filter is reinitialised individually by forming a statistic based on all old information [13, p. 1275], referred to as mixing estimates. The difference between gpb(2) and imm is that the mixing estimates are formed after the conditional filtering in the gpb(2) algorithm. This order is reversed in the imm algorithm, which omits a lot of calculations that usually do not

(35)

con-3.5 Multiple Model Filters 25

tribute to the performance [6, p. 269]. Looking at the tree in Figure 3.2 with time horizon n = 2, this would mean that conditional filtering is made for all four merged branches in the gpb(2) algorithm. The imm would however only require two filtering operations, since the mixing estimates are calculated before the mea-surement update. Rearranging the calculations in this way thus makes the imm more efficient. The reinitialisation when using three filters in the imm algorithm is shown in Figure 3.4.

Figure 3.4:Reinitialisations of the imm algorithm [13, p. 1274].

The other steps of the algorithm are similar to gpb(1). An illustration of one cycle of imm, using three models, is shown in Figure 3.5 and the algorithm itself is given by Algorithm 3 [13, p. 1276].

(36)

Algorithm 3One cycle of the imm algorithm.

1. Model-conditioned reinitialisation(for i = 1, 2, ..., M): Predicted mode probability: µ(i)_k|k−1, P {m(i)k |yk−1}=

P

jπjiµ

(j)

k−1

Mixing weight: µj|i_k−1, P {m(j)_k−1|_m(i)

k , yk−1}= πjiµ

(j)

k−1/µ

(i)

k|k−1

Mixing estimate: ¯x(i)_k−1|k−1, E[xk−1|m

Predicted state: ˆx(i)_k|k−1= F_k−1(i) x¯(i)_k−1|k−1+ G(i)_k−1w¯(i)_k−1

Predicted covariance: P_k|k−1(i) = F_k−1(i) P¯_k−1|k−1(i) (F_k−1(i) )0+ G(i)_k−1Q(i)_k−1(G(i)_k−1)0 Predicted measurement: ˆy_k|k−1(i) = H_k(i)xˆ(i)_k|k−1+ ¯v_k(i)

Innovation: ε(i)_k = yk−yˆ

(i)

k|k−1

Innovation covariance: S_k(i)= H_k(i)P_k|k−1(i) (H_k(i))0

+ R(i)_k Kalman gain: K_k(i)= P_k|k−1(i) (H_k(i))0(S_k(i))−1

Updated state: ˆx(i)_k|k= ˆx(i)_k|k−1+ K_k(i)ε(i)_k

Updated covariance: P_k|k(i)= P_k|k−1(i) −_K(i)

k S (i) k (K (i) k ) 0

3. Mode probability update(for i = 1, 2, ..., M): Model likelihood: L(i)_k , p[ε(i)k |m (i) k , yk−1] assume = N_(ε(i) k ; 0, S (i) k ) = exp −1 2(ε (i) k ) 0 (S_k(i))−₁ (ε(i)_k ) 2πS (i) k 1/2

Mode probability: µ(i)_k|k= µ (i) k|k−1L (i) k P jµ (j) k|k−1L (j) k

(37)

4

Simulation Study

This chapter describes the simulation study performed. Three cases of targets per-forming different manoeuvres are presented. The implementation of the radar model and tracking filters are also described, together with the choice of some parameters used. Lastly, different ways of evaluating the performance of the fil-ters are given.

4.1 Benchmark Trajectories

To generate radar measurements and evaluate the performance of different fil-ter algorithms, three benchmark trajectories for manoeuvring targets were taken from [2], which provided both trajectories and rcs for the targets. The cases included different types of targets, performing both heavy and smaller manoeu-vres. This was to evaluate how the filters perform in different situations. The cases included:

• A large aircraft performing simple manoeuvres. • A manoeuvrable small aircraft.

• A heavily manoeuvring fighter aircraft.

Simple trajectories for the radar source were also added to create test cases. These were chosen in such a way that the radar is approaching the target, while keeping the target in the radar’s search volume. This was to make the test cases realistic since an aircraft often follows its target and it allowed the radar to detect more points.

(38)

4.1.1 Large Aircraft Performing Simple Manoeuvres

The first trajectory represents a larger aircraft, such as a military cargo aircraft, with an average rcs of 4 m2 [2, p. 1107]. It begins with a constant course and a speed of 290 m/s, before it makes a 2g turn and continues on a new course. Lastly, it makes a 3g turn and flies away from the radar. The trajectory also in-cludes simulations of range gate pull-off (rgpo), which occur after 15 and 40 seconds. rgpo is a technique used in warfare to break the radar lock-on to. The target produces a signal which is similar to the one that would have been reflected by the target from a radar pulse. This pulse is then delayed in time to make the radar lose track of the real target, making it follow a false target instead.

The radar source has a constant course and a speed of 120 m/s. Both trajectories are shown in Figure 4.1.

Figure 4.1:Benchmark trajectory 1.

4.1.2 Manoeuvrable Small Aircraft

The second trajectory represents a smaller, more manoeuvrable aircraft, with an average rcs of 2 m2 [2, p. 1107]. The target begins with a constant course and a speed of 305 m/s, before it makes a 2.5g turn. It then descends before finally making a 4g turn and continuing with a speed of 305 m/s once again. rgpo starts after 12, 50 and 95 seconds.

The radar source has a constant course and a speed of 200 m/s. Both trajectories are shown in Figure 4.2.

(39)

4.1 Benchmark Trajectories 29

Figure 4.2:Benchmark trajectory 2.

4.1.3 Heavily Manoeuvring Fighter Aircraft

The third trajectory represents a fighter/attack aircraft with an average rcs of 1.2 m2[2, p. 1109]. The target begins with an acceleration before making a 5g turn after 30 seconds, while maintaining full throttle. After another turn of 7g, the target flies with constant speed. Lastly, a 6g turn is performed and a climb is made, before completing the trajectory without acceleration. rgpo starts after 5, 25 and 52 seconds.

The radar source has a speed of 100 m/s while making a small turn. Both trajec-tories are shown in Figure 4.3.

(40)

4.2 Radar Model Implementation

Based on the theory presented in Chapter 2, a radar model was implemented in Matlab®. The aim was to create a realistic model of a radar found in fighter aircraft, to be able to evaluate the tracking performance in this setting.

4.2.1 Volume Search

A rectangular search pattern, defined by azimuth and elevation angles (α, β), is used for the volume search. At each beam position, five tries are made to detect a target. This represents sending out multiple pulse trains from the radar. Requir-ing three out of five detections then lowers the probability of false alarm while maintaining about the same detection range as described in Section 2.5.4. The probability of detection is calculated with (2.16), which is used to decide if the target is detected or not.

4.2.2 Pulse Repetition Frequency

For each detection attempt the prf is increased slightly without changing any other parameter. This is done to represent the ability to resolve ambiguities in the measurements, as described in Section 2.2.2.

There is also an option to choose if a high prf, medium prf or a mixture between the both should be used. In the mixed mode, the prf is changed every time one iteration of the search volume pattern is completed. A mixed mode is often used in real radar applications to get the properties from both prfs.

4.2.3 Measurement Noise

Once a target is detected, a measurement is created by adding a normal dis-tributed error to the true value, with the variance represented by the measure-ment accuracy from Table 4.1. In reality, there is a limit for how accurate a radar can be. Therefore a lower bound on the accuracy was implemented, for range, range rate and angle measurements. The bounds were set to to the values shown in Table 4.1.

Table 4.1:Lower bound on the accuracy of the radar measurements.

Measurement Lower bound

r 30 m

˙r 1 m/s

(41)

4.2 Radar Model Implementation 31

4.2.4 Matlab Algorithm

With given trajectories for the target and the radar source, the radar model gen-erates a set of measured points. The steps in Matlab®

when performing the vol-ume search and generating these measurements are summarised by Algorithm 4.

Algorithm 4Volume search.

1: Determine if target is in the current search cone

2: ifTarget is inside search cone then

3: Calculate the signal to noise ratio

4: Calculate the probability of detection based on the snr

5: Check if target is detected

6: ifTarget is detected then

7: Create measurement by adding calculated measurement noise to the true position and radial velocity

8: end if

9: end if

4.2.5 Choice of Parameters

The radar model uses the parameters in (2.10) to calculate the snr. Realistic val-ues for many of these parameters can be found in, for example, [17]. The PFA

was set to a certain value, which was improved during the detection process by performing multiple dwells. The cumulative PFAwhen requiring three out of five

detections then achieved a satisfactory average time between false alarms. All parameters that were set to a specific value are shown in Table 4.2. Other parameters, such as the receiver gain, G, and the antenna beam width, θ3dBwere

calculated using the formulas described in Chapter 2.

Table 4.2:Fixed parameters used in the radar model.

Parameter Value

Carrier wavelength, λ 0.03 m

Receiver temperature 290 K

Peak antenna power, P 20 kW

Receiver noise figure, F 3 dB Total system losses, Ls 6 dB

Antenna length, D 0.5 m

Duty cycle, dt 0.1

PFA 10−4

P_FAc ∼₁₀−11

Azimuth search angle, α 20◦ Elevation search angle, β 12◦

(42)

To investigate how the prf affects the performance, both a medium, high and mixed prf were tried. The parameters related to the choice of prf are shown in Table 4.3. The choice of prf together with the number of pulses sent out at each beam position determines the dwell time at each beam position during the volume search.

Table 4.3:Parameters for different choices of prf.

Parameter Value for high prf Value for medium prf

prf 200 kHz 12 kHz

Number of pulses 1024 128

Dwell time 0.0256 s 0.0533 s

4.2.6 Detection Performance

Different choices of PFAaffects how big the snr must be for the radar model to

detect a target. In Figure 4.4, the snr is plotted against the probability of detec-tion for different values of PFA. The figure also shows how the PFA is lowered,

while keeping about the same detection performance, by using three out of five detections.

Figure 4.4: snrcompared to probability of detection. The dashed lines rep-resent the cumulative probability when requiring three out of five detec-tions.

The required snr for detection does not depend on which prf is used. The snr is however affected by the distance to the target and the prf, which will give different detection ranges. The PDis plotted against the distance to the target, for

(43)

4.3 Generated Measurements 33

(a)Medium prf. (b)High prf.

Figure 4.5: Range to target compared to probability of detection using the parameters in Table 4.2 with rcs = 1.2 m2. The dashed lines represent the cumulative probability when requiring three out of five detections. In the previous section the PFA was chosen as 0.0001 which is represented by the

red lines.

4.3 Generated Measurements

The benchmark trajectories described in Section 4.1 were used by the radar model to generate test data. The beginning of a volume search is illustrated in Figure 4.6.

Figure 4.6:Graphical illustration of the radar performing its volume search. The blue dot indicates the current position of the target and the blue cone is the current beam position. The full search volume is illustrated by the red rectangle with edges connected to the radar source.

(44)

is shown in figure 4.7. The number of detected points will depend on the snr, which is affected by the choice of prf and the distance to the target.

(a)Detected points after t = 20 seconds. (b)Completed volume search.

Figure 4.7: During the volume search target detections are marked as green points. A green cone means that a high prf is used, while a blue cone indi-cates a medium prf.

After each point is detected, measurement noise is added based on the snr. One volume search, using the third benchmark trajectory, resulted in the measured points displayed in Figure 4.8.

Figure 4.8: Generated measurements for the detected points of a volume search.

The same process was then repeated to generate different sets of points, which was later used for Monte Carlo simulations. This was done for all three trajec-tories, and also for different choices of prf, i.e. a medium, high, and mixed prf. The generated data was then saved to be able to test different filters with the same data.

(45)

4.4 Filter Implementation 35

4.4 Filter Implementation

This section describes which filters are evaluated and how they are initialised, together with the model choices they are built on.

4.4.1 Choice of Filters and Parameters

The cv, ca and ct models cover some of the most common motion patterns that can be expected to be performed by a target. The multiple model algorithms de-scribed in the Chapter 3 are therefore evaluated with this model-set. The three different models are also evaluated separately with individual ekfs.

The process noise used in the filters was tuned ad-hoc by trial and error. Different process noise was used for every individual filter, but also for the mm algorithms. That is because the best values found for the individual filters did not give the best performance when used in the mm algorithms. A considerable amount of different values, and combination of values in the mm algorithms, were evalu-ated to find solutions that performed well for all kinds of manoeuvres studied. A fixed turn rate was used for the ct model, which was also determined by simu-lation results. A turn rate of ω = 9◦/s proved to be a good choice since it captures the behaviour of a turning target better than the cv and ca model in some cases. A too small value proved to make the ct model pointless since it then performed similar to a ca model. If the turn was too large the filter performance became poor in all situations.

4.4.2 Initialisation

The uncertainty of the radar measurements varies depending on the snr. For a target far away, the measurements will generally be less accurate. However, as mentioned earlier the radar model was implemented with a lower bound on the accuracy of the measurements. Therefore the measurement noise covariance matrix, R, was tuned by first selecting it as the lower bound. The values were then increased, which in many cases gave a better filter performance. Sometimes, a slightly larger R than the true measurement noise can give better results [10]. The best value found by trial and error was then used in the algorithms.

The state vector in all filters needs to be initialised with some value, ˆx0. A good

option is to choose the first measurement as the position estimate and select the velocity and acceleration to be zero. The state covariance matrix P will then be initialised with P0. The values corresponding to the position variance were

se-lected to be the measurement variance. For the velocity and acceleration, the cor-responding values in P0were set to large values [10]. These values represented

the maximum velocity and acceleration for a target in realistic operating condi-tions. For example, the velocity of a target being tracked will probably be lower than 800 m/s.

(46)

4.4.3 Transition Probabilities

The state transition probability matrix, Π, given by

Π=         π11 π12 π13 π21 π22 π23 π31 π32 π33         =         0.98 0.01 0.01 0.01 0.98 0.01 0.01 0.01 0.98         , (4.1)

was set to the same values for both the gpb(1) and imm algorithm. Each element in the matrix corresponds to the probability of transitioning from one model to another and every row must sum up to one. For a target being tracked, it is rea-sonable to think that it does not change its motion pattern between every time instant and will remain in the same mode for some time [5]. The three models used in the mm algorithms were also different from each other. Therefore the probability of remaining in the same mode should be high.

Setting the values of transitioning between two different models too high makes it difficult for the algorithms to decide on which model should be trusted and therefore many jumps occur. This makes the general tracking performance worse. A common choice found in literature is thus to set the diagonal elements to values close to one [3]. Different values were first used and tested, before choosing the matrix that gave the best overall performance.

4.5 Performance Measures

The performance of the filters were evaluated by looking at the estimation error. This was also compared to a theoretical lower bound of the filter covariance.

4.5.1 Root Mean Square Error

The root mean square error is often used to evaluate the performance of a fil-ter. This was done by using a number of Monte Carlo (mc) simulations. For M number of mc simulations the rmse is calculated as

rmse_{( ˆ}_x_k|k_{) =} v u t 1 M M X i=1 ˜

x_k|kT (i) ˜xk|k(i), (4.2a)

˜

xk|k(i) = xk(i) − ˆxk|k(i), (4.2b)

where xk is the true state value, ˆxk|k is the estimated state and ˜xk|k is the

Evaluation of Tracking Filters for Tracking of Manoeuvring Targets

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2020