Topics in Multisensor Maneuvering Target Tracking
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee. This
dissertation does not include proprietary or classi ed information.
Soonho Jeong
Certi cate of Approval:
Stanley J. Reeves
Professor
Department of Electrical and Computer
Engineering
Jitendra K. Tugnait, Chair
James B. Davis Professor
Department of Electrical and Computer
Engineering
Thomas S. Denney Jr.
Associate Professor
Department of Electrical and Computer
Engineering
Stephen L. McFarland
Dean
Graduate School
Topics in Multisensor Maneuvering Target Tracking
Soonho Jeong
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Ful llment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
August 8, 2005
Topics in Multisensor Maneuvering Target Tracking
Soonho Jeong
Permission is granted to Auburn University to make copies of this thesis at its discretion,
upon the request of individuals or institutions and at their expense. The author reserves
all publication rights.
Signature of Author
Date
iii
Vita
Soonho Jeong, a Korean Army Major, was born on December 7, 1967 in Seoul,
Republic of Korea (ROK), the third of four children of Wansuk Jeong and Jungja Hwang.
He graduated from Daewon High School, in Seoul, Korea, in 1986. He received the
B.S. and M.S. degrees in electrical engineering from Korea Military Academy, Seoul,
Korea, in 1990 and Naval Postgraduate School, Monterey, CA, in 1996, respectively.
From 1999 to 2002 he served in the Weapons Systems Management Group (WSMG),
Korea Army Headquarters as an Intelligence Technology O cer (ITO). He entered the
graduate program in Electrical Engineering at Auburn University in August 2002.
He married Jeesoo Song, daughter of Kibok Song and Sunmoon Park, on June 13,
1993. He has two daughters, Seoin Jeong and Jahyun Jeong, who were born in November
16, 1994 and in August 26, 1997, respectively.
iv
Dissertation Abstract
Topics in Multisensor Maneuvering Target Tracking
Soonho Jeong
Doctor of Philosophy, August 8, 2005
(M.S., Naval Postgraduate School, Monterey, CA, March 1996)
(B.S., Korea Military Academy, Republic of Korea, March 1990)
207 Typed Pages
Directed by Jitendra K. Tugnait
Tracking uses models of the real environment to estimate the past and present and
even predict the future state of a moving object from noisy observations of uncertain
origin. In a tracking scenario the most critical problem is that of data-association. This
topic has received considerable attention in the literature and a number of solutions have
been proposed. This dissertation considers the problem of tracking highly maneuvering
target(s) using multiple sensors in the presence of clutter. A set of noble algorithms are
developed to handle this problem.
First, the basic interacting multiple model (IMM) approach has been combined
with probabilistic data association (PDA) to develop an IMMPDA (interacting multiple
model probabilistic data association) algorithm with simultaneous measurement update
(SMU) for tracking a maneuvering target in clutter with multiple sensors.
Second, we extend our noble SMU algorithm to a more practical tracking scenario,
that of tracking a maneuvering target with asynchronous (in-sequence but time delayed)
measurements. A state-augmented approach is developed to estimate the time delay
v
between a local sensor (assumed to be collocated and synchronized with a central pro-
cessor) and a remote sensor (assumed to be separately located and not synchronized
with a central processor).
Third, we address one of the most important issues for target tracking in a multisen-
sor fusion network: out-of-sequence measurements (OOSM). However, this dissertation
is not concerned with di erent sampling rate among sensors. Instead, we focus on a
suboptimal ltering algorithm dealing with possibly time delayed, out-of-sequence mea-
surements (OOSM) with a xed relative time-delay (we assume that sampling rate are
all the same for all sensors) among sensor measurements. A state-augmented approach
is also developed to improve tracking performance with the possible presence of OOSM.
The ltering algorithm is developed by OOSM updating with IMMPDA for the target.
Finally, we consider tracking of multiple highly maneuvering targets using multiple
sensors with possibly unresolved measurement. When multiple targets move temporarily
in close formation, it possibly gives rise to a single detection due to the resolution limita-
tions of the sensor. Assuming that there are possibly unresolved measurements from at
least two targets (i.e., measurement association with more than two targets simultane-
ously), any measurement therefore is either associated with a target, a group of merged
targets, or caused by clutter. The ltering algorithm is developed by applying the basic
IMM approach and the joint probabilistic data association with merged measurements
(JPDAM) technique and coupled target state estimation.
vi
Acknowledgments
I would like to express my deep appreciation to Dr. Jitendra K. Tugnait for his
guidance and encouragement throughout this work. During the course of this research,
he has provided all the possible support for conducting my research with lasting patience.
Thanks are due to my committee members, Dr. Stanley J. Reeves and Dr. Thomas S.
Denney Jr., for their guidance and help. I would like to thank Dr. Dean G. Ho man, Dr.
Soo-Young Lee, and Dr. Scottedward A. Hodel for their advice and valuable discussions.
I also would like to thank the Korea Army Headquarters and the O ce of Naval Research
for their support during this research. I also would like to thank Dr. Sang-Bum Lee,
Dr. Junhak Song and Gen. Chi-Kyu Yang for their encouragement on my study. I
also would like to acknowledge the immense support I received from my colleagues at
Auburn University. I will thank them here under the socially acceptable convention of
alphabetical order: Dr. Shuangchi He, Dr. Jung Hyun Jo, Dr. Moon Seong Kang, Dr.
Hyungwook Kim, Dr. Xiaohong Meng, Dr. Sumedh P. Puranik, Mr. Shuying Qi, Dr.
Junhak Song and Dr. Daesub Yoon.
Finally, I would never have been able to nish writing this dissertation without
the support and understanding of my family. Many thanks to my children, Seoin and
Jahyun, for giving Daddy time to work on the research. My hearty thanks to my wife,
Jeesoo, whose love, support, encouragement and patience during the past three years
have made it possible for me to complete my education.
vii
Style manual or journal used Journal of Approximation Theory (together with the
style known as \auphd"). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package TEX (speci cally
LATEX) together with the departmental style- le auphd.sty.
viii
Table of Contents
List of Tables xii
List of Figures xiii
1 Introduction 1
1.1 Overview of Target Tracking Problems . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Data Association Techniques . . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Estimation Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.3 Combination of IMM and PDA . . . . . . . . . . . . . . . . . . . . 7
1.2 The Main Issues of Target Tracking Problem . . . . . . . . . . . . . . . . 9
1.2.1 Tracking a Maneuvering Target . . . . . . . . . . . . . . . . . . . . 9
1.2.2 Tracking a Target in Presence of Clutter . . . . . . . . . . . . . . . 10
1.2.3 Tracking a Target using Multiple Sensors . . . . . . . . . . . . . . 12
1.2.4 Tracking Multiple Targets . . . . . . . . . . . . . . . . . . . . . . . 15
1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . 20
2 Basic Target Tracking Algorithms 22
2.1 Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.1.1 Discrete-Time Kalman Filter . . . . . . . . . . . . . . . . . . . . . 23
2.1.2 Extended Kalman lter . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2 Interactive Multiple Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.1 Optimal Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.2 Basic IMM Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.2.3 Target Dynamic Models . . . . . . . . . . . . . . . . . . . . . . . . 42
2.3 Probabilistic Data Association Filter . . . . . . . . . . . . . . . . . . . . . 48
2.3.1 Measurement Validation . . . . . . . . . . . . . . . . . . . . . . . . 49
2.3.2 PDA approach combined with Kalman lter . . . . . . . . . . . . . 51
2.3.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3 Multisensor Tracking of a Maneuvering Target in Clutter using
IMMPDA Filtering with Simultaneous Measurement Update 54
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2 Simultaneous measurement update . . . . . . . . . . . . . . . . . . . . . . 57
3.3 Problem Formulation for the Multiple Model System . . . . . . . . . . . . 64
3.4 IMM/MSPDAF Algorithm for Simultaneous Measurement Update . . . . 66
3.5 Simulation Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
ix
4 Multisensor Tracking of a Maneuvering Target in Clutter with
Asynchronous Measurements using Augmented State IMMPDA Fil-
tering and Simultaneous Measurement Update 83
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.2 Problem Formulation for Asynchronous measurements . . . . . . . . . . . 85
4.3 State-Augmented System . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.4 AS-IMM/MSPDAF Algorithm for Asynchronous Measurements . . . . . . 88
4.5 Simulation Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5 Multisensor Tracking of a Maneuvering Target in Clutter with
Asynchronous and Possibly Out-of-Sequence Measurements using
Augmented State IMMPDA Filtering and Simultaneous Measurement
Update 106
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.2 Modeling Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.3 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.3.1 Target Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.3.2 Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.4 State-Augmented System . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.5 IMM/MSPDAF Algorithm for Asynchronous and Possibly Out-of-
Sequence Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.6 Simulation Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
6 Tracking of Multiple Maneuvering Targets in Clutter with Possibly
Merged Measurements using IMM and JPDAM Coupled Filtering 136
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
6.2 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
6.2.1 Target Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
6.2.2 Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
6.3 Modeling for the Merged Measurements . . . . . . . . . . . . . . . . . . . 144
6.3.1 Modeling Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . 144
6.3.2 Measurement Model . . . . . . . . . . . . . . . . . . . . . . . . . . 145
6.3.3 Sensor Resolution Model . . . . . . . . . . . . . . . . . . . . . . . . 147
6.4 IMM/JPDAM Coupled Filtering Algorithm . . . . . . . . . . . . . . . . . 148
6.5 Simulation Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
7 Conclusions and Future Work 174
7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
7.2 Suggested Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
x
7.2.1 Unresolved Measurements with Simultaneous Measurement Up-
date . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
7.2.2 Track Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
Bibliography 178
Appendices 183
A Derivation of the mode-conditioned association event probability
Eqn. (3.46) 184
B Derivation of Eqn. (6.62) 188
xi
List of Tables
4.1 Simulation Results: No. of lost tracks obtained from 100 Monte Carlo
runs for xed-but-unknown dk . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.1 Simulation results summery based on 1000 runs . . . . . . . . . . . . . . . 172
xii
List of Figures
1.1 An outline of the components of a tracking system [12, page 5] . . . . . . 2
1.2 A mathematical view of the state estimation problem [12, page 2] . . . . . 2
1.3 An example of origin uncertainties and validation regions: ^Z1k and ^Z2k are
predicted measurements from target 1 and target 2 at time k, respectively.
The parameters of the ellipses are determined by the covariance matrices
of ^Z1k and ^Z2k. Measurement Z1k originates from target 1 or clutter. Z3k
originates from target 2 or clutter. Z2k originates from target 1 or target
2 or both (in case of merged measurement) or clutter. . . . . . . . . . . . 16
2.1 The trajectory of the target (true vs. estimated) . . . . . . . . . . . . . . 32
2.2 Performance of the Kalman lter (read top-to-bottom): (a) RMSE in
velocity. (b) RMSE in position . . . . . . . . . . . . . . . . . . . . . . . . 33
2.3 Optimal solution (circle) vs suboptimal solution (diamond: interacting
multiple model (IMM) algorithm) in term of the number of lters required 38
2.4 An example of a set of kinematic models describing a maneuvering target
using two models: 1) constant velocity model and 2) coordinate turn with
constant acceleration model. . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.1 Target tracking example using synchronized measurements observed from
multiple sensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.2 Trajectory of the maneuvering target (read left to right, top to bottom).
(a) Position in xy plane. (b) x and y velocities. (c) x and y accelerations.
(d) magnitude of accelerations . . . . . . . . . . . . . . . . . . . . . . . . 80
3.3 Performance of the simultaneous measurement updating IMM/MSPDAF
in terms of the RMSE (root mean square error) in position (read
top to bottom). (a) proposed IMM/MSPDAF vs standard sequential
IMM/MSPDAF [21] with collocated sensors (Case 2). (b) proposed
IMM/MSPDAF vs standard sequential IMM/MSPDAF [21] with sepa-
rated sensors (Case 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
xiii
4.1 Target tracking with a xed-but-unknown relative time delay dk between
the remote sensor clock and the central processor clock at sample time tk 86
4.2 Trajectory of maneuvering target (read left to right, top to bottom). (a)
Position in xy plane. (b) x and y velocities. (c) x and y accelerations. (d)
magnitude of accelerations . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.3 Estimation of delay (given unknown but xed timing mismatch between
two separated sensors) based on 100 Monte Carlo runs (read left to right,
top to bottom). (a) d = 0. (b) d = 0.1T. (c) d = 0.3T. (d) d = 0.5T. (e)
d = 0.7T. (f) d = 0.9T. (T = sampling rate). Solid: estimated delay ^d;
dashed: xed delay d. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.4 Comparison of ltered and smoothed (lag = 1) estimate for various delay
values (acceleration, velocity, and position RMS errors (3 rows each), read
left to right, top to bottom). (a) d = 0. (b) d = 0.1T. (c) d = 0.3T.
(d) d = 0.5T. (e) d = 0.7T. (f) d = 0.9T. (T = sampling rate). In the
gure legends, estimation refers to ltering, and smoothing is with lag =
1. Solid: ltered estimate; dash-dot: smoothed estimate. . . . . . . . . . . 103
4.5 RMSE in position using IMM/MSPDAF for various delay values (read left
to right, top to bottom). (a) d = 0. (b) d = 0.1T. (c) d = 0.3T. (d) d =
0.5T. (e) d = 0.7T. (f) d = 0.9T. Unless otherwise stated, the results are
for ltering. Solid: proposed AS-IMM/MSPDAF; dotted: proposed AS-
IMM/MSPDAF (smoothing); dash-dot: standard IMM/MSPDAF [38];
dashed: proposed AS-IMM/MSPDAF with the knowledge of dk. . . . . . 104
5.1 Asynchronous measurements: In-sequence but delayed measurements and
out of sequence measurements (OOSM) in multisensor tracking system . . 112
5.2 Trajectory of maneuvering target (read left to right, top to bottom). (a)
Position in xy plane. (b) x and y velocities. (c) x and y accelerations. (d)
magnitude of accelerations. . . . . . . . . . . . . . . . . . . . . . . . . . . 133
5.3 AS-IMM/MSPDA comparison (RMSE in position, read top to bottom)
for various probabilities of delayed measurement, Pd = 0:25 and 0.4: (a)
Pd=0.25. (b) Pd=0.4. Solid: proposed AS-IMM/MSPDAF algorithm
dealing with OOSM; dashed: AS-IMM/MSPDAF algorithm [50] with
OOSM discarding; dotted: AS-IMM/MSPDAF algorithm applied to the
hypothetical case of Pd = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . 134
6.1 The true trajectories of the maneuvering targets (read left to right, top
to bottom): (a) Position in xy plane, (b) x and y velocities, (c) x and y
accelerations, (d) distance between the targets. . . . . . . . . . . . . . . . 170
xiv
6.2 Performance (RMSE in position) of the proposed IMM/JPDAMCF and
the IMM/JPDACF of [25] based on successful runs (read left to right):
(a) proposed IMM/JPDAMCF, (b) standard IMM/JPDACF [25]. . . . . 171
xv
Chapter 1
Introduction
Target tracking is the estimation (the process of selecting the value of interest from
indirect, inaccurate and uncertain observations) of the states (position, velocity, acceler-
ation, etc.) of moving objects (plane, missile, submarine, etc.) both at the current time
( ltering) and at any point in the future (prediction) based on remote measurements
obtained from sensor(s). The objective of target tracking is to collect sensor data from
a eld of view (FOV) containing potential target information and then to partition the
sensor data into sets of observations or tracks that are originated from the same sources.
Di erent types of sensors can be employed to obtain remote measurements. Active
sensors observe objects by illuminating them with an energy source and measuring the
re ected energy. On the other hand, passive sensors observe objects by measuring char-
acteristic emissions of the object. The sensors may be located at a xed location or on
moving platforms. The system block diagram of target tracking and the corresponding
mathematical view of state estimation are shown in Fig. 1.1 and Fig. 1.2, respectively.
Note that, in Figs. 1.1 and 1.2, the measurements are the only variables to which the
estimator has access and are a ected by the error sources in the form of measurement
noise [12].
Solution of the target tracking problem requires the simultaneous completion of
two tasks: estimation and data association. Estimation can be viewed as the process of
nding the best model parameters to describe the observed data. The process of assigning
observations to each target track is referred to as data association. There are many
1
Figure 1.1: An outline of the components of a tracking system [12, page 5]
Figure 1.2: A mathematical view of the state estimation problem [12, page 2]
2
di erent approaches to both estimation and data association, and these are generally
the distinguishing features that give rise to di erent tracking algorithms. Details of
these tasks are discussed in Chapter 2. In this chapter we present an overview of general
target tracking problems, some speci c issues of target tracking problem, and main
contributions of our works.
1.1 Overview of Target Tracking Problems
1.1.1 Data Association Techniques
Target tracking algorithms have continued to receive signi cant attention over the
past fty years. The problem of target tracking dealing with noise corrupted measure-
ments has attracted signi cant attention since World War II. In 1955, Wax [1] formulated
the problem of track formation, track maintenance, and track rejection based on the prob-
lem of detecting the path of a particle in a bubble chamber, but this paper was written
before the adoption of Kalman ltering techniques for recursive target state estimation.
In 1960, Kalman published his famous paper [2] describing a recursive solution to the
discrete-data linear ltering problem. Since then, the Kalman lter has been the subject
of extensive research and application, particularly in the area of autonomous or assisted
navigation. It was recognized by Sittler [3] in 1964 that there can be an uncertainty
associated with the origin of the measurements in target tracking. This is called the
data association problem; a measurement may not have originated from the target of
interest.
In the early 70?s papers by Bar-Shalom [4, 5] and Singer [6, 7] began the develop-
ment of the modern tracking systems which combined the data association and Kalman
3
ltering theory. In 1975, the paper by Bar-Shalom and Tse [5] introduced a suboptimal,
computationally-bounded extension of the Kalman lter to cases where measurements
were not always available or multiple measurements were simultaneously returned from
the sensor. This extension was termed the probabilistic data association lter (PDAF).
The PDAF is a Bayesian technique that incorporates clusters of measurements that could
have originated from a single target into an updated state estimate for the target. It led
to fewer lost tracks in the presence of clutter when compared to the nearest neighbor
(NN) technique of data association where only the closest measurement to the predicted
measurement was considered to have originated from the target and used to update the
state estimate of the target [5]. Two years later from [5], the multiple hypothesis testing
(MHT) technique was introduced by Reid [8, 9]. This approach splits an existing track
into multiple tracks when more than one measurement is available for the target. How-
ever, in this dissertation we are not concerned with the MHT technique, but instead we
focus on the PDAF technique because it has been recognized as one of the best subopti-
mal data association approaches in terms of combined accuracy and computational cost
[27]. The theory behind the PDAF is discussed in Chapter 2.
In the early 1980?s, as an extension of PDAF, the joint probabilistic data association
lter (JPDAF) was introduced by Fortmann et al. [10] to deal with the multiple target
tracking situation. The JPDAF is identical to the PDAF except for the extra computa-
tion of the association probabilities for all observations and all tracks. The JPDAF was
fairly successful in simulations at tracking multiple targets travelling in close proximity.
However, when two targets are \closely" spaced, they may give rise to a single detection
due to the resolution limitations of the sensor. For instance, in radar ranging, returns
4
from multiple targets could fall in the same range cell, resulting in one unresolved de-
tection only. This violates a principal assumption of the JPDAF: a measurement can
originate from at most one target. Standard tracking algorithms that ignore such a phe-
nomenon can lead to poor performance in multiple target tracking. As a consequence
of the study of this problem, Chang and Bar-Shalom [11] introduced the joint proba-
bilistic data association with merged measurements (JPDAM) technique in 1984. Under
the Bayesian framework, there are two basic methods of track estimation with merged
measurements in multiple target environments [12]: JPDAM and MHT. The MHT l-
ter associates feasible measurements to track and forms multiple hypotheses for track
extension. It is a measurement-oriented approach whereas JPDAM is a target-oriented
approach. Moreover, MHT is a multiscan approach utilizing several scans of measure-
ments to make data association decisions. MHT makes hard decisions where highly
improbable hypotheses are pruned to reduce the computational burden. On the other
hand, the JPDAM lter is a single scan approach which does not make hard decisions;
rather it updates a track with a weighed sum of the measurements which could have
originated from the target in track. The basic theory and simulation of the JPDAM are
developed in Chapter 6. Looking ahead in this dissertation, the JPDAM technique has
better performance than JPDAF especially in terms of the track estimation accuracy
and the loss of tracks when two targets are \closely" spaced and may give rise to a single
detection.
1.1.2 Estimation Algorithms
Optimal techniques, such as the Kalman lter, provide accurate results during cer-
tain motions but also cause large errors during target motions which are not modeled in
5
the Kalman lter. Moreover, the discrete-time Kalman lter requires sensor measure-
ments at every time step, which is not always possible in target tracking scenario. After
the data association problem had been rst addressed by the PDAF, several multiple
model tracking methods appeared in the early 1980?s. The idea behind a multiple-model
approach versus a traditional adaptive lter such as the Kalman lter was that traditional
techniques would respond too slowly in abruptly changing environments [12, 13]. Opti-
mal multiple model techniques were replaced by suboptimal techniques as memory re-
quirements and computational time grew exponentially for the optimal approaches. The
suboptimal estimation techniques introduced included the generalized pseudo-Bayesian
(GPB1), second-order generalized pseudo-Bayesian (GPB2), detection-estimation algo-
rithm (DEA), random sampling algorithm (RSA), and the interacting multiple model
(IMM). The di erences in these estimation approaches were how and when the multiple
models were combined (also known as hypothesis pruning).
The GPB1 and GPB2 algorithms were introduced to combine the history of target
models by Ackerson and Fu [14] in 1970. In the GPB methods, the hypothesis pruning
was performed after the ltering step. The GPB2 di ered from the GPB1 by including
knowledge of the previous possible mode transitions, as modeled by a Markov chain. In
the late 1970?s, both DEA and RSA approaches were proposed for the switching environ-
ment problem and well detailed in [15] and [16], respectively. Tugnait [17] provided good
overviews where suboptimal hypothesis pruning techniques, including the GPB, DEA,
and RSA methods, are compared. These pruning techniques involved approximating
sums of Gaussian distributions by a single Gaussian distribution, but were presented
before the widespread adoption of IMM techniques.
6
The IMM algorithm was introduced by Blom [18] in 1984. He and Bar-Shalom
[19, 20] further developed the idea of IMM in a couple of papers in the late 1980?s.
Di ering from the GPB techniques, the hypothesis pruning was performed before each
ltering step after information was shared among the various lters. Thus, like the GPB2
method, the IMM considers the previous possible mode transitions, but unlike the GPB2
method, the IMM requires not n2 but n lters, where n is the number of models being
considered. It was shown in [19] that the GPB2 and IMM produced slightly smaller
tracking errors than the GPB1 during non-maneuvering target motion, but achieved
signi cantly smaller tracking errors during maneuvering motion. In view of the fact that
the IMM is conceptually similar to the GPB2 but its computational load is similar to that
of the GPB1 [19], the IMM is superior to both the GPB1 and GPB2 when considering
computational complexity and tracking accuracy. The basic theory behind the IMM is
discussed in Sec. 2.2.
1.1.3 Combination of IMM and PDA
In the process of tracking a maneuvering target, for cases when measurements are
not always available for the Kalman lter, or multiple measurements for a target exist,
Houles and Bar-Shalom [21] introduced the combination of the IMM and the PDAF in
the late 1980?s. In [21] the IMM algorithm is combined with the PDA lter in a multiple
sensor scenario to propose a combined IMM/MSPDAF (interacting multiple model/
multisensor probabilistic data association lter) algorithm. This algorithm considers an
estimation of the target state at present time k given measurements up to time k (on
state ltering). This combination is detailed in [21].
7
In [10, 22] multiple targets in clutter (but without using switching multiple mod-
els) have been considered using JPDA lter which, unlike the PDA lter, accounts for
the interference from other targets. Most recently, various versions of IMM/JPDA l-
ters for multiple target tracking using switching multiple models have been presented
in [23, 24, 25, 26]. While [24, 26] present uncoupled lters (i.e., assumes that di erent
target states are mutually independent conditioned on the past measurements), [23, 25]
present coupled lters (i.e., assume that there exists \shared" measurements, yielding
cross-covariances which re ect the correlation between the targets? state estimation er-
rors). In addition, while Blom and Bloem [23] presented an \exact" JPDA coupled lter
for non-switching models using the framework of a linear descriptor system, for switch-
ing models, they also presented IMM/JPDA uncoupled lter approximations [24]. In
[25], Tugnait presented an IMM/JPDA coupled ltering algorithm where a simulation
example resulted in fewer target swapping compared with uncoupled IMM/JPDA. As
noted in [27], IMMPDA lter is in general superior to IMM/MHT lter when the asso-
ciated computational cost and performance are considered. Therefore, our emphasis will
be on IMM/JPDA techniques. Neither [11] nor [28] consider multiple switching kine-
matic models for maneuvering targets; rather they are limited to single (non-switching)
kinematic models per target to achieve much enhanced performance.
When two targets are \closely" spaced, they may give rise to a single detection due
to the resolution limitations of the sensor. For instance, in radar ranging, returns from
multiple targets could fall in the same range cell, resulting in one unresolved detection
only [11, 28]. Standard tracking algorithms that ignore such a phenomenon can lead
to poor performance in multiple target tracking [11, 28]. Despite its importance, prior
work on tracking with unresolved measurements and modeling of resolution capability of
8
a sensor in particular is sparse. Prior work includes [11] and [28] and references therein.
In [11] the resolution phenomena related to tracking have been treated on the basis of
a grid of resolution cells \frozen" in space. In [28] the resolution capability of the event
that two targets are unresolved is conditioned on the relative distance between the two
targets in terms of the measured variables (range, azimuth, etc.). A simple Gaussian
shape is assumed which captures the sensor behavior in a mathematically tractable way.
While [11] considers JPDA for data association, [28] exploits MHT.
In Chapter 6, we propose to use sensor resolution modeling of [28] in conjunction
with JPDAM coupled ltering and IMM approach (see [25] for tracking with resolved
measurements scenario).
1.2 The Main Issues of Target Tracking Problem
The main issues of target tracking problem discussed in this dissertation can be
summarized as follows.
1.2.1 Tracking a Maneuvering Target
Traditional target tracking assumes that the states of the target of interest satisfy
a certain kinematic model. However target model uncertainty typically exists because
generally the target of interest is not a cooperating target in that it does not follow a
prede ned trajectory [27]. Basically the uncertainty in the model can be modeled by
additive noise which compensates for the modeling inaccuracy. This approach works
fairly well when the target kinematics can be closely approximated by a single model.
9
Therefore, if a certain model truly describes the kinematics of the target, a Kalman
lter gives the optimal solution for tracking in the sense of minimizing the mean-square
error in state estimation. The optimal Kalman lter yields large errors when the target
is moving with nearly-constant speed and the noise variance is very high. However, in
some cases, a target can maneuver and can exhibit di erent kinematic characteristics
from time to time. For example, a military aircraft can perform many kinds of abrupt
maneuvers. In that case, a single model cannot describe completely the behavior of the
target, and the estimates based on a single model often lead to poor performance of the
state estimates or loss of track.
The idea of using multiple models for describing the di erent motion phases of the
target comes from the above impracticality of the lter. For instance, if a constant
velocity model is used to describe the behavior of the target all the time, the selection
of additive noise variance which models the uncertainty in the model is di cult. If the
variance of the noise is small and the target undergoes rapid acceleration or sharp turn,
this model su ers the loss of track of the target because the small noise variance can not
cover properly those maneuvers. On the other hand, if the variance of the noise is large,
it would be able to track those maneuvers but estimation accuracy would be degraded.
The interacting multiple model (IMM) estimator has become well accepted in the
literature [12, 22, 29, 30, 33] as the best approach for tracking a maneuvering target.
1.2.2 Tracking a Target in Presence of Clutter
In target tracking, there exist various kind of \undesired" measurements (clutter)
obtained from the sensor which are not generated from the actual target of interest.
10
This clutter (e.g., interfering radar echoes caused by objects such as clouds, sea waves,
etc.) can degrade the tracking accuracy severely. Ideally, this measurement uncertainty
can be described by a random model. However, in practice, there exists clutter that is
persistent or somewhere between random and persistent [12]. We will only consider ran-
dom clutter in the sequel. In the worst case, there exists no true measurement detected
by the sensor, and only clutter may be present as measurement data. This cluttered
environment complicates the target tracking problem by introducing uncertainty in the
measurement origin. As a consequence of this measurement uncertainty, errors can be
made in association of measurements to existing tracks. Hence one has to consider data
association techniques to determine whether the true measurement originated from the
target (Does a given measurement come from the target of interest, or is it a clutter or
interference from the nearby target?).
The simplest way to deal with multiple measurements is to use NN (nearest neigh-
bor) lters which, at any given moment, use the NN measurement as if it were the mea-
surement coming from the target of interest. Thus it discards all other measurements for
the nal state estimation of the target of interest. But sometimes NN measurement may
not originate from the target of interest and could be clutter. This happens when clutter
density is high. So it is not a robust approach, as it often leads to loss of track [12]. The
other way to deal with multiple measurements is to associate each measurement with a
weight according to the probability that a given measurement could have originated from
the target of interest. This is a Bayesian approach. More speci cally, a validation region
(hypersphere) centered at the predicted measurement (based on previous state estimate)
is established. A validation process of all the measurements in the current scan is carried
out. Only those measurements falling within the validation region are considered so that
11
this process eliminates highly improbable measurements. The validation process needs
an appropriate model to describe the clutter. Once the clutter model is established, a
data association process can be developed accordingly. An e ective data association ap-
proach in a Bayesian framework is that of probabilistic data association (PDA) [22, 33]
which is discussed in Sec. 2.3.
1.2.3 Tracking a Target using Multiple Sensors
The objective of using multiple sensors for target tracking is to improve the accuracy
of tracking (i.e., infrared sensors only measures the direction of arrival, hence radar can
be used together to measure the distance of the target, or at least two infrared sensors
at di erent locations can be used to obtain the position of a target). In this dissertation,
we restrict our attention to radar and infrared sensors.
Radar Sensors
Radar (radio detection and ranging) sensor is mostly used as an active sensor since
it initiates the energy and detects the target re ections. It is well known that
radar?s all-weather performance and excellent kinematic measurement capabilities
allow it to play a dominant role in any multisensor tracking system. This active
nature of the radar is the primary factor that makes it such an e ective sensor.
Since the source of the energy is the radar itself, the nature of the returned target
signal is at least partially under control of the radar operator. The range res-
olution and aperture size of the radar is primarily determined by the frequency
band used by the radar. It is also possible to use the sensor in a passive mode,
searching for targets emitting radar signals which is often referred as radar warning
12
receivers (RWR). The measurement model for a 2-D radar sensor considered in this
dissertation is
z = h(x) =
2
66
4
R
3
77
5+w
where x is the system state, z is the (true) measurement vector, w is the measure-
ment noise, R is the range, and is the azimuth angle. A common model for the
measurement noise is to assume an additive white Gaussian noise model.
Infrared Sensors
The passive nature of the infrared (IR) sensors leads to distinct advantages and
disadvantages. The excellent measurement accuracy of the IR sensors makes them
an important element of any multiple sensor tracking system. Emphasis in military
applications on reduced target signatures places new emphasis on the development
of IR system in mainly two ways. First, use of active sensors must be minimized
since the use of active sensors potentially uncovers one?s existence (e.g., the en-
emy can easily detect radar tracking by using RWR). Second, the enemy?s e orts
at reduced RCS (radar cross section) signature and advanced ECM (electronic
countermeasures) increase the need to exploit other target signatures such as its
infrared emissions. In addition, the use of IR sensors enhances resolution capability
for closely spaced targets, accurate angle measurements, and resistance to ECM.
The measurement model for this sensor considered in this dissertation is
z = h(x) =
2
66
4
3
77
5+w
13
where is the azimuth angle and is the elevation angle.
Unlike single sensor case, multiple sensors introduce the problems of temporal and
spatial characteristics of the measurement data. Sensors may not always be located
at the same place. Hence measurements available at the sensors could be at di erent
time instants. So sensors may be either synchronous (ideal case), which means that
the sampling times of all sensors are the same and there is no delay (timing mismatch)
among sensors, or asynchronous, which means that the sampling times may vary or
there is delay (timing mismatch) from each other. This gives rise to architectural issues
of the sensor network, and one may carry out the estimation process in a centralized
way (estimation at one central node, called the master node) or in a decentralized way
(estimation at each sensor locally and then obtain global estimate by track fusion) [27].
If sensors are distant from each other, then there could be delay in the estimation
due to data transmission time and channel bandwidth issues. So one has to deal with
synchronous as well as asynchronous measurement data to obtain the estimate. Usually
if synchronous measurement data is available from all the sensors then either ?parallel
updating? (updating state by processing data from all the sensors simultaneously at the
same time) or ?sequential updating? (updating state based on the available measurement
data from that particular sensor) can be carried out. In [34], it has been shown that
those two schemes are equivalent (for linear models) if provided sensors are synchronized
and measurement noise across the sensors are uncorrelated (a realistic assumption).
But parallel updating is computationally more expensive than sequential updating. On
the other hand, if the sensors are not synchronized, then one has to deal with possible
OOSM (out of sequence measurements). This asynchronous measurements case has been
14
considered in [12, 35, 36, 37] and a noble state-augmented algorithm to deal with OOSM
is presented in Chapter 5.
1.2.4 Tracking Multiple Targets
The presence of multiple targets in the detection range of sensors (also called the
scope) makes the target tracking very complicated by introducing the additional uncer-
tainty in the origin of the measurement. One has to solve the data association problem
to determine whether the available measurement is generated from one of the targets
present in the scope and from which target it originated. Data association in a multiple
target scenario is carried out by considering all the targets simultaneously.
Since each target can constitute a \persistent" interference to other targets, it is
a great challenge to track multiple targets moving close to each other. The merged
measurement case arise when two (or more) targets are so \close" that the sensor detects
them as one target due to a lack of resolution (see Fig. 1.3). The objective of the
multitarget tracking is to partition the sensor data into sets of observations, or tracks,
produced by the same source (target). Once tracks are formed and con rmed, the
task is to maintain each track. All of the above target tracking issues have already
been studied by numerous researchers and, as a consequence, many e ective solutions
have been provided for these problems. However, more e ective and cheap algorithms
are still required to improve the estimation accuracy and to reduce the computational
cost. Compared to single target tracking, multitarget tracking is a very complicated
and computationally intensive process. When there are n targets in the detection range
(scope) of the sensors, the measurement origin uncertainty problem is the worst one as the
15
1?
kZ
2?
kZ
2
kZ
1
kZ
3
kZ
Figure 1.3: An example of origin uncertainties and validation regions: ^Z1k and ^Z2k are pre-
dicted measurements from target 1 and target 2 at time k, respectively. The parameters
of the ellipses are determined by the covariance matrices of ^Z1k and ^Z2k. Measurement Z1k
originates from target 1 or clutter. Z3k originates from target 2 or clutter. Z2k originates
from target 1 or target 2 or both (in case of merged measurement) or clutter.
number of combinations of association of measurement and target increases dramatically.
An objective of our research is to develop advanced signal processing algorithms for
multiple target tracking by exploiting ltering theory and data association techniques.
There are mainly two approaches to deal with the merged measurement problem.
Multiple hypothesis tracker (MHT)
No prior knowledge is needed about the number of targets. One has to consider
simultaneously all the targets and also consider simultaneously all the scans of
measurements. The main theme of MHT is that it evaluates the probabilities that
there is a target from which a sequence of measurements have originated. It is
16
called a measurement-oriented approach. The more complex MHT provides im-
proved performance, but it is di cult to implement and in cluttered environments
a large number of hypotheses may have to be maintained, which requires extensive
computational resources.
Joint probabilistic data association(JPDA)
The JPDA techniques are based on PDA (probabilistic data association), which
uses a weighted average of all the measurements falling inside a track?s validation
region (i.e., gating) at the current time to update the track state. In this approach,
it is assumed that the number of targets is known already. The tracks of these
targets have already been formed. The task is to evaluate the measurement-to-
target association probabilities for the latest set of measurements and then obtain
the state estimates by combining those measurements. It is called a target-oriented
approach. These techniques give much better performance than the simpler data
association techniques (i.e., NN lter) and requires less computational resources
than MHT approach.
JPDA with merged measurements (JPDAM)
When multiple targets are close enough in the measurement range, they will give
merged (unresolved) measurement due to the resolution limitation of sensor. Since
the JPDA algorithm can track targets with measurements originating from one of
the targets in track or clutter, additional data association for the merged measure-
ments needs to be done in JPDA technique. To deal with the possibility that a
17
measurement originated from more than one target, JPDAM can be applied. JP-
DAM uses a probabilistic model and corresponding data association for the merged
measurements.
1.3 Contributions
Our contributions in this dissertation are:
IMMPDA simultaneous measurement update for a maneuvering target track-
ing using multiple sensors
We present a noble suboptimal ltering algorithm for tracking a highly maneuver-
ing target in a cluttered environment using multiple sensors. The ltering algorithm is
developed by applying the basic IMM approach and the PDA technique to a two sensor
(radar and infrared) problem for state estimation for the target. A simultaneous mea-
surement update approach is followed where the raw sensor measurements are passed
to a fusion node and fed directly to the target tracker. A multisensor (MS) PDA is
developed for parallel sensor processing for target tracking under clutter. The algorithm
is illustrated via a highly maneuvering target tracking simulation example where two
sensors, a radar and an infrared sensor, are used. Compared with an existing IMMPDA
ltering algorithm with sequential sensor processing, the proposed algorithm achieves
signi cant improvement in the accuracy of track estimation.
IMMPDA simultaneous measurement update for a maneuvering target track-
ing with delayed but in-sequence measurements
We extend IMMPDA simultaneous measurement update algorithm to possibly asyn-
chronous (in-sequence but time delayed) measurements. A state augmented approach is
18
developed to estimate the time delay between local and remote sensors. A multisensor
probabilistic data association lter (PDAF) is developed for parallel sensor processing for
target tracking under clutter. Compared with an existing IMMPDA ltering algorithm
with the assumption of synchronous measurements sensor processing, the proposed algo-
rithm achieves considerable improvement (especially in the case of larger delays) in the
accuracy of track estimation.
IMMPDA simultaneous measurement update with possibly out-of-sequence
measurements
We propose a suboptimal ltering algorithm dealing with possibly time delayed,
out-of-sequence measurements (OOSM) with a xed relative time-delay among sensor
measurements. The ltering algorithm is developed by applying the basic IMM approach,
the PDA technique, and OOSM updating for the target. A state augmented approach is
developed to improve tracking performance with the presence of OOSM. A multisensor
PDA lter is developed for parallel sensor processing for target tracking under clutter.
The algorithm is illustrated via a highly maneuvering target tracking simulation example
where two sensors, a radar and an infrared sensor, are used. Compared with an exist-
ing IMMPDA ltering algorithm with in-sequence only sensor processing, the proposed
algorithm achieves considerable improvement in the accuracy of track estimation.
Multitarget tracking with possibly merged measurements
We propose a noble suboptimal ltering algorithm for tracking multiple maneuvering
targets in a cluttered environment using multiple sensors. This algorithm is an extension
of IMM/JPDA ltering algorithm [25] to deal with possibly merged measurements. We
19
concentrate on the case of two targets which temporarily move in close formation. The
ltering algorithm is developed by applying the basic IMM approach and the joint prob-
abilistic data association with merged measurements (JPDAM) technique and coupled
target state estimation to a Markovian switching system. The algorithm is illustrated
via a simulation example involving tracking of two highly maneuvering, at times closely
spaced, targets with possibly unresolved measurements. Compared with an existing
IMM/JPDA ltering algorithm developed without allowing for merged measurements
[25], the proposed algorithm achieves signi cant improvement in the accuracy of track
estimation during target merging period.
1.4 Organization of the Dissertation
The remainder of this dissertation is organized as follows:
In Chapter 2, we elaborate on some basic target tracking algorithms, mainly the
Kalman lter, IMM lter, and PDA lter. Theory associated with those algorithms
serves as the foundation for our contributions. Chapter 3 presents a noble suboptimal l-
tering algorithm for tracking a highly maneuvering target in a cluttered environment us-
ing multiple sensors. We develop a noble algorithm \simultaneous measurement update
using IMMPDA ltering". The proposed algorithm provides signi cant improvement
over an existing standard IMMPDA ltering algorithm with sequential sensor processing
especially tracking a maneuvering target. In Chapter 4, we extend the proposed ltering
algorithm to the asynchronous measurement case. Compared with an existing IMM-
PDA ltering algorithm with the assumption of synchronous (no delay) measurements
sensor processing, the proposed algorithm achieves considerable improvement (especially
20
in the case of larger delays) in the accuracy of track estimation. A suboptimal ltering
algorithm dealing with possibly time delayed, out-of-sequence measurements (OOSM)
with a xed relative time-delay among sensor measurements is presented in Chapter 5.
In Chapter 6, we focus on multitarget tracking problem with possibly merged measure-
ments. The simulation example shows signi cant improvement in the position estimate
compared to the multisensor tracking by IMM/JPDA coupled ltering algorithm of [25].
Finally, Chapter 7 presents the future work that we would like to continue in the area.
21
Chapter 2
Basic Target Tracking Algorithms
In this chapter we review some basic algorithms. We divide the problem of track-
ing into state estimation and data association. We brie y describe some existing state
estimation algorithms: Kalman lter, extended Kalman lter, and interacting multiple
model (IMM) algorithms. For data association, the probabilistic data association lter
(PDAF) is brie y described in this chapter. This is background work and is quite useful
in later chapters where we discuss our contributions.
2.1 Kalman Filter
Kalman lter theory is fundamental for the development of target tracking algo-
rithms. Given a collection of sensor measurements, it is often possible to determine
several optimal state estimates based on di erent optimality conditions which can be
stated in terms of the probability density function (pdf) of the states given the obser-
vation. Hence, the problem of state estimation can be formulated by calculating the
probability density of the state. Under general conditions, it is not possible to obtain
a closed-form solution to the problem of state estimation. However, when the system
is linear and the random variables are assumed to be Gaussian, a closed-form solution
can be obtained. Under this conditions, the optimal estimator is the Kalman lter and
it is usually speci ed as a nite length algorithm based on a recursion of the mean and
covariance of the state. It can be applied to both continuous and discrete linear sys-
tems. The target of interest, which is being tracked, is continuous by nature, but the
22
processor, which runs the tracking algorithm, is discrete. The original Kalman lter [2]
was de ned in continuous-time, but a discrete version also derived after [2]. Since the
discrete Kalman lter is only used in this dissertation, we brie y describe the steps for
implementing a Kalman lter in discrete time domain.
2.1.1 Discrete-Time Kalman Filter
Consider a discrete time linear dynamic system with additive white Gaussian noise
that models the disturbances. The discrete time index is denoted by k. A general
time-varying state space model for the Kalman lter is
xk+1 = Fkxk +Gkuk +vk (2.1)
where
xk = the state vector of dimension nx at time k
Fk = nx nx transition matrix
Gk = nx nu input distribution matrix
uk = a known input vector of dimension nu
vk = the random process noise vector of dimension nx.
It is assumed that vk is the sequence of zero-mean white Gaussian process noise with
covariance
Efvjv0kg= Qjk = Q jk (2.2)
23
where Ef:g denotes the expectation, 0 denotes the transpose operation and jk is the
Dirac delta function. The statistical model of the measurement is described by
zk = Hkxk +wk (2.3)
where zk is the mx-dimensional measurement vector, Hk is the mx nx observation
matrix, and wk is zero mean white Gaussian measurement noise with covariance
Efwjw0kg= Rjk = R jk: (2.4)
Since a complete knowledge of the statistical model constitutes the knowledge of Fk, Gk,
Qk, Hk, and Rk, these system matrices are assumed known and can be time varying. The
initial condition of the state x0 is a random variable with known mean x0 and covariance
P0. Process noise vk, measurement noise wk, and initial state x0 are assumed to be
mutually independent. De ne ^xkjj (:= EfxkjZjg) as an estimate of the state vector xk
based upon the knowledge of the measurements up to time j where Zj :=fz1;z2; ;zjg.
Specially, k >j denotes a predicted estimate, k M2
model switch
M1-> M2
model switch
M2-> M1
model switch
M2-> M1
Model M1: constant velocity
Model M2: coordinate turn w/ constant acceleration
Figure 2.4: An example of a set of kinematic models describing a maneuvering target
using two models: 1) constant velocity model and 2) coordinate turn with constant
acceleration model.
kinematic models (discretized model is generated via a continuous-time white process
noise) and direct discrete-time kinematic models (it is directly de ned in discrete-time
domain and the noise is a discrete-time process noise) [13]. We mainly focus on direct
discrete-time kinematic models as follows:
Piecewise Constant Wiener Process Acceleration Model
Assume that there is an object moving in a 3-dimensional Cartesian space plane
with non-zero acceleration. Let the state vector at time kT (T = sample period)
xk = [ k _ k k k _ k k k _ k k]0 (2.61)
where k, k and k are the Cartesian position coordinates, and where the velocity
components are denoted by _ k, _ k and _ k, and where the acceleration components are
denoted by k, k and k. The state dynamics for the piecewise constant Wiener process
43
acceleration model can be approximated using Taylor series expansion:
f (t) = f (t0) + _f (t0) (t t0) + f (t0) (t t0)
2
2! + : (2.62)
Let t = (k + 1)T and t0 = kT. Then
k+1 = k + _ kT + kT22 +p k; (2.63)
_ k+1 = _ k + kT +v k; (2.64)
k+1 = k +a k: (2.65)
De ne
xk ( ) :=
2
66
66
66
4
k
_ k
k
3
77
77
77
5
; xk :=
2
66
66
66
4
xk ( )
xk ( )
xk ( )
3
77
77
77
5
; wk ( ) :=
2
66
66
66
4
p k
v k
a k
3
77
77
77
5
; and wk :=
2
66
66
66
4
wk ( )
wk ( )
wk ( )
3
77
77
77
5
:
(2.66)
From Eqns. (2.62)-(2.66), the state dynamics for the piecewise constant Wiener process
acceleration model can be written as in (2.67) [13]. (Note that, only in this section,
we denote \wk" (instead of \vk") as the process noise at time k to distinguish it from
velocity components in (2.64) and (2.66).)
xk+1 = Fkxk +Gkwk (2.67)
44
where
Fk =
2
66
66
66
4
Fk ( ) 0 0
0 Fk ( ) 0
0 0 Fk ( )
3
77
77
77
5
; Fk ( ) = Fk ( ) = Fk ( ) =
2
66
66
66
4
1 T 12T2
0 1 T
0 0 1
3
77
77
77
5
(2.68)
and
Gk =
2
66
66
66
4
Gk ( ) 0 0
0 Gk ( ) 0
0 0 Gk ( )
3
77
77
77
5
; Gk ( ) = Gk ( ) = Gk ( ) =
2
66
66
66
4
1
2T
2
T
1
3
77
77
77
5
: (2.69)
Piecewise Constant White Acceleration Model
Assume that there is an object moving in a 3-dimensional Cartesian space plane
with constant speed (zero acceleration). Let the state vector at time kT (T = sample
period)
xk = [ k _ k k _ k k _ k]0: (2.70)
Then similar to Eqns. (2.62)-(2.66),
k+1 = k + _ kT +p k; (2.71)
_ k+1 = _ k +v k: (2.72)
45
De ne
xk ( ) :=
2
66
4
k
_ k
3
77
5; xk :=
2
66
66
66
4
xk ( )
xk ( )
xk ( )
3
77
77
77
5
; wk ( ) :=
2
66
4
p k
v k
3
77
5 and wk :=
2
66
66
66
4
wk ( )
wk ( )
wk ( )
3
77
77
77
5
:
(2.73)
From Eqns. (2.70)-(2.73), the state dynamics for the piecewise constant Wiener process
acceleration model can be written as [13]
xk+1 = Fkxk +Gkwk (2.74)
where
Fk =
2
66
66
66
4
Fk ( ) 0 0
0 Fk ( ) 0
0 0 Fk ( )
3
77
77
77
5
; Fk ( ) = Fk ( ) = Fk ( ) =
2
66
4
1 T
0 1
3
77
5 (2.75)
and
Gk =
2
66
66
66
4
Gk ( ) 0 0
0 Gk ( ) 0
0 0 Gk ( )
3
77
77
77
5
; Gk ( ) = Gk ( ) = Gk ( ) =
2
66
4
T
1
3
77
5: (2.76)
Coordinated Turn Model
Assume that there is an object moving in a xy plane with constant speed and turning
with a constant angular rate (i..e., coordinated turn in aviation language) [13, 44]. Let
46
the state vector
xk = [ _ _ ]0 (2.77)
where and are the Cartesian position coordinates, and where the velocity components
are denoted _ and _ . If we assume a constant speed =
q
_ 2 + _ 2 and a constant angular
rate = _ ( > 0 implies a counterclockwise turn), where we have
_ = cos( ); _ = sin( ): (2.78)
Now the acceleration (with _ = 0) components are
= d
dt
_ = sin ( ) = _ (2.79)
= ddt _ = cos ( ) = _ : (2.80)
The state equation can be expressed as
x(t) = ddt
2
66
66
66
66
66
4
_
_
3
77
77
77
77
77
5
=
2
66
66
66
66
66
4
_
_
_
_
3
77
77
77
77
77
5
: (2.81)
47
By discretizing the system, it gives
xk+1 =
2
66
66
66
66
66
4
1 sin( T) 0 1 cos( T)
0 cos ( T) 0 sin ( T)
0 1 cos( T) 1 sin( T)
0 sin ( T) 0 cos ( T)
3
77
77
77
77
77
5
xk +Gkwk: (2.82)
This allows us to generate the state trajectories for such turns. These coordinated turns
are very common for any ying objects.
2.3 Probabilistic Data Association Filter
The multiple model approach discussed above section is e ective when we track a
single target under the assumption that measurements of the target are available at all
time steps from the sensor(s). This assumption means that we must be able to associate
measurements from sensor(s) to each target of interest before we do any ltering. In
addition, in the case of multiple radar returns received at the radar sensor, it is possible
for sensor(s) to miss the radar return from the target of interest for several consecutive
scans. Therefore, we must not only associate the available measurements with existing
target tracks but also consider the cases when there is no measurement available for the
target of interest.
The basic algorithms discussed so far in the previous section do not take into ac-
count the problem of measurement origin uncertainty. In this section we discuss the
basic probabilistic data association lter (PDAF) that was proposed by Bar-Shalom in
[5]. PDAF is a Bayesian technique that handles these two issues: measurement origin
48
uncertainty and no measurement available for the target. It is mainly used for tracking
a single target in the presence of clutter. For tracking multiple targets, the same ap-
proach of PDA has been extended to JPDA (joint probabilistic data association) [10].
In JPDA one has to do data association jointly by considering simultaneously all the
targets present in the surveillance region. Details about JPDA are given in Chapter 6.
Here we give the basic PDA algorithm of [5]. PDA algorithm calculates the association
probabilities of each validated measurement to the target at the current time. It is usu-
ally combined with a proper tracking lter to track a single target in clutter. Let us
discuss now the validation of measurements.
2.3.1 Measurement Validation
Assume that there is a target whose track has been initialized. One can then set up
a validation gate [12], based on predicted measurement and residual covariance obtained
in the appropriate lter prediction stage. A key assumption is that the past informa-
tion about the true target state conditioned on the past measurements is summarized
approximately by
p(xkjZk 1) =N(xk; ^xkjk 1;Pkjk 1): (2.83)
In addition, the true measurement conditioned on the past measurements is also Gaussian
distributed with probability density function given by
p(zkjZk 1) =N(zk; ^zkjk 1;Sk): (2.84)
49
The validation region is given by
V(k; ) :=
n
z : [z ^zkjk 1]0S 1k [z ^zkjk 1]
o
(2.85)
with gate threshold . The gate probability
PG := P[zk2V(k; )] (2.86)
gives the probability that the true measurement will fall within the validation region.
The probability PG is a function of and the dimension of the measurement nz. The
volume V of the validation region V corresponding to the threshold is given by [12]
Vk = cnz nz2 jSkj12 (2.87)
where cnz is the volume of the unit hypersphere of dimension nz (c1 = 2;c2 = ;c3 =
4 =3, etc.). Therefore, an increase in the design parameter, cnz or Sk, will increase the
size of the validation area. Practically we take = 16 and the corresponding PG = 0:9989
(or 0.9997) when nz = 3 (or 2). The following assumptions are made [12]:
(AS1) Among the possibly several validated measurements, at most one of them can be
target originated provided the target is detected and the corresponding measurement
falls within the validation gate.
(AS2) All the remaining measurements are assumed to be due to clutter or false alarm
and modeled as independent and identically distributed (i.i.d.) with uniform spatial
distribution.
(AS3) Target detection occurs independently over time with known probability PD.
50
2.3.2 PDA approach combined with Kalman lter
The following section discusses how to apply the PDA to a typical Kalman lter
which tracks a single target in clutter. From among all the raw measurements at time
k, i.e.,
Zk :=fz(1)k ;z(2)k ; ;z(mk)k g; (2.88)
de ne the set of validated measurement at time k as
Yk :=fy(1)k ;y(2)k ; ;y( mk)k g (2.89)
where mk is total number of validated measurement at time k. And
y(i)k := z(j)k 1 i j mk (2.90)
where 1 mk mk when mk 6=0. The cumulative set of validated measurements is
denoted by
Zk :=fY1;Y2; ;Ykg: (2.91)
We also de ne set of association events [12]
ik = yik is originated from the target: i = 1; ; mk
0k = none of the measurements is originated from the target: (2.92)
Note that the eventsf ikg mki=0 are mutually exclusive and exhaustive based on (AS1).
Now the estimate of the state xk conditioned on the above events is obtained by the
51
total probability theorem as
^xkjk := EfxkjZkg
=
mkX
i=0
Efxkj ik;ZkgP[ ikjZk] (2.93)
=
mkX
i=0
^xikjk ik
where xikjk := Efxkj ik;Zkg is the updated state conditioned on the event that the i-th
validated measurement is correct and ik := P[ ikjZk] is the association probability with
this event. The conditional state estimate based on measurement i being correct is
^xikjk = ^xkjk 1 +Wk ik i = 1; ; mk (2.94)
where the innovations
ik := zik ^zkjk 1 (2.95)
and gain
Wk := Pkjk 1H0kS 1k : (2.96)
The nal state update equation is
^xkjk = ^xkjk 1 +Wk k (2.97)
where the combined innovations
k :=
mkX
i=1
ik ik: (2.98)
52
For i = 0 (when no measurement is originated from the target of interest), we have
^x0kjk = ^xkjk 1, where the updated state estimate is identical to the predicted state. The
covariance associated with the updated state is [12]
Pkjk = Pkjk 1 (1 0k)WkSkW0k +Wk(
mkX
i=1
ik ik i0k k 0k)W0k: (2.99)
2.3.3 Summary
In this chapter we brie y discussed the background work needed to understand the
fundamentals of target tracking. The Kalman lter provides the optimum solution in
the MMSE sense, but its application is limited to a single model. For the multiple
model approach, suboptimal algorithms such as IMM lter are applicable. To tackle the
measurement origin uncertainty, one has to combine the IMM lter with data association
techniques like PDA. Tracking a target using multiple sensors can improve the tracking
accuracy. With this background work we are ready to discuss our noble algorithms in
the following chapters.
53
Chapter 3
Multisensor Tracking of a Maneuvering Target in Clutter using
IMMPDA Filtering with Simultaneous Measurement Update
In this chapter, we present a suboptimal ltering algorithm for tracking a highly
maneuvering target in a cluttered environment using multiple sensors. The ltering algo-
rithm is developed by applying the basic interacting multiple model (IMM) approach and
the probabilistic data association (PDA) technique to a two sensor (radar and infrared,
for instance) problem for state estimation for the target. A simultaneous measurement
update approach is followed where the raw sensor measurements are passed to a central
processor and fed directly to the target tracker. A multisensor probabilistic data associ-
ation lter is developed for simultaneous measurement update (SMU) for target tracking
under clutter. A past approach using SMU has ignored certain data association proba-
bilities leading to an inaccurate implementation. Another existing approach applies only
to non-maneuvering targets. The algorithm is illustrated via a highly maneuvering tar-
get tracking simulation example where two sensors, a radar and an infrared sensor, are
used. Compared with an existing IMMPDA ltering algorithm with sequential sensor
processing, the proposed algorithm achieves signi cant improvement in the accuracy of
track estimation.
3.1 Introduction
We consider the problem of tracking a maneuvering target in clutter. This class of
problem has received considerable attention in the literature [5, 12, 21, 22, 44, 45]. We
54
develop a simultaneous measurement update technique by applying the basic interacting
multiple model (IMM) algorithm and probabilistic data association (PDA) technique.
The switching multiple model approach has been found to be quite e ective in modeling
highly maneuvering targets [12, 14, 21, 29, 30, 33, 45]. In this approach various \modes"
of target motion are represented by distinct kinematic models, and in a Bayesian frame-
work, the target maneuvers are modeled by switching among these models controlled by
a Markov chain. To accommodate the fact that the target can be highly maneuvering,
we will follow a switching multiple model formulation as in [12, 21, 29, 30, 33] and refer-
ences therein. It is assumed that a track has been formed (initiated) and our objective
is that of track maintenance. In [21] such a problem has been considered using multiple
sensors, PDA and switching multiple models. The optimal solution (in the minimum
mean-square error sense) to target state estimation, given sensor measurements and ab-
sence of clutter, requires exponentially increasing (with time) computational complexity;
therefore, one has to resort to suboptimal approximations. For the switching multiple
model approach, the IMM algorithm of [30] has been found to o er a good compromise
between the computational and storage requirements and estimation accuracy [29]. In
the presence of clutter, the measurements at the sensors may not all have originated
from the target of interest. Therefore, one has to account for measurements of uncertain
origin (target or clutter?). In this case one has to solve the problem of data association.
An e ective approach in a Bayesian framework is that of PDA [12, 22]. Here too, in a
Bayesian framework, one has to resort to approximations to reduce the computational
complexity, resulting in the PDA lter [5, 12, 21, 22, 45].
In [21] the IMM algorithm has been combined with the PDA lter in a multiple
sensor scenario to propose a combined IMM/MSPDAF (interacting multiple model/
55
multisensor probabilistic data association lter) algorithm. While [21] uses sequential
updating of the state estimates with measurements (i.e., updating the state estimates
sequentially with measurements from di erent sensors), the multisensor approach of [45]
falls in the category of simultaneous measurement update, or parallel sensor processing,
where the raw measurements from all sensors are passed to a central processor to be
processed simultaneously (i.e., updating the state estimates with all the measurements
at the same time as if they were from a single sensor). Sequential updating results in
computational savings but this approach is not necessarily the best. For linear systems,
both sequential and parallel updating methods are algebraically equivalent but the par-
allel updating is computationally more expensive [12]. Ref. [45] uses SMU but has some
errors: during data association, all measurements at the same time from di erent sensors
are assumed to be either from clutter or from the target. The possibility that a measure-
ment from sensor 1 may be from target while the measurement from sensor 2 may be
clutter-induced (and vice-versa) is implicitly not allowed in [45] - this is clearly incorrect.
Ref. [46] allows for such distinctions (hypotheses), but it is limited to non-maneuvering
targets.
In this chapter, we extend the multisensor approach of [46] to maneuvering targets.
\Standard" assumptions are used for PDA ltering [12, 22]: a measurement can have
only one source; among the possibly several validated measurements, at most one of
them can be target-originated and the remaining validated measurements are assumed
to be due to false alarms or clutter, and are modeled as independently and identically
distributed (i.i.d) with uniform spatial distribution over the entire validation region.
The remainder of this chapter is organized as follows. Sec. 3.2 presents the ba-
sic implementation of simultaneous measurement update technique. Sec. 3.3 presents
56
the problem formulation for multiple model system. Sec. 3.4 describes the proposed
IMM/MSPDAF algorithm with simultaneous measurement update to multiple model
system. Simulation results using the proposed algorithm for a realistic problem are
given in Sec. 3.5. Finally, concluding remarks may be found in Sec. 3.6.
3.2 Simultaneous measurement update
In the simultaneous measurement update procedure, the state can be updated simul-
taneously with all the synchronized measurements observed from multiple sensors. We
consider a system of a 2-D radar and an infrared sensor located separately at [x1;y1;z1]
and [x2;y2;z2], respectively, and covering a common cluttered surveillance region that
is being traversed by a single non-maneuvering target. Assume that the target dynamic
equation can be described as
xk = Fk 1xk 1 +Gk 1vk 1: (3.1)
Assume that there are q synchronized sensors. At a given sampling time k, there are ml
measurements from each sensor l. The measurement from sensor l at time k is
zlk = Hlkxk +wlk for l = 1; ;q; (3.2)
where xk is the system state at tk and of dimension nx, zlk is the (true) measurement
vector (i.e., due to the target) from sensor l at tk and of dimension nzl, where Hlk is the
Jacobian matrix of hl evaluated at some value of the estimate of state xk. The process
noise vk 1 and the measurement noise wlk are mutually uncorrelated zero-mean white
57
Gaussian processes with covariance matrices Qk 1 and Rlk, respectively. Note that, in
general, at any time k, some measurements may be due to clutter and some due to the
target, i.e., there can be more than a single measurement at time k at sensor l. The
measurement set (not yet validated) generated by sensor l at time k is denoted as
Zlk :=fzl(1)k ;zl(2)k ; ;zl(ml)k g (3.3)
where ml is the number of measurements generated by sensor l at time k. Variable
zl(i)k (i = 1; ;ml) is the ith measurement within the set. The cumulative set of
measurements (not yet validated) from sensor l up to time k is denoted as Zk(l) :=
fZl1;Zl2; ;Zlkg. In heavily cluttered environment, validation gates can be applied to
reduce the number of unwanted measurements for further processing. Following [12, 21],
one sets up a validation gate for sensor l centered at the predicted measurement, ^zlk.
Then measurement zl(i)k (i=1,2, ,ml) is validated if and only if
[zl(i)k ^zlk]0[Slk] 1[zl(i)k ^zlk] < (3.4)
where is an appropriate threshold. The volume of the validation region with the
threshold is
Vlk := cnzl nzl=2jSlkj1=2 (3.5)
where nzl is the dimension of the measurement and cnzl is the volume of the unit hyper-
sphere of this dimension (c1=2, c2= , c3=4 =3, etc.). Choice of is discussed in more
detail in [12, Sec. 2.3.2]. From among all the raw measurements from sensor l at time
k, i.e., Zlk :=fzl(1)k ;zl(2)k ; ;zl(ml)k g, de ne the set of validated measurement for sensor
58
l at time k as
Ylk :=fyl(1)k ;yl(2)k ; ;yl( ml)k g (3.6)
where ml is total number of validated measurement for sensor l at time k, and
yl(i)k := zl(li)k (3.7)
where 1 l1
><
>>:
2
66
4
1(i)k
2(i)k
3
77
5
1(i)0k 1(i)0k
jZk 1
9>
>=
>>;
=
2
66
64
H1k
H2k
3
77
75Pkjk 1
H10k H20k
+
2
66
64
R1k 0
0 R2k
3
77
75: (3.13)
De ne the association event probabilities as
a;bk := P[ a;bk jY 1k;Y 2k;Zk 1]: (3.14)
Exploiting the di use model for clutter in [12, 21], it turns out that
0;0k = C(1 PD1PG1)(1 PD2PG2)(V1
k )
m1(V2
k )
m2 ; a = 0;b = 0
0;bk = C
PD2(1 PD1PG1)N
h
2(b)k ;0;S2k
i
(V2k ) m2 1 m2 ; a = 0; b = 1; ; m2
a;0k = C
PD1(1 PD2PG2)N
h
1(a)k ;0;S1k
i
(V1k ) m1 1 m1 ; a = 1; ; m1; b = 0
a;bk = C PD1PD2N
a;bk ;0;Sk
m1 m2(V1k ) m1 1(V2k ) m2 1; a = 1; ; m1; b = 1; ; m2
(3.15)
61
where PD1 and PD2 are the detection probabilities that the sensors 1 and 2 detect the
target, respectively, PG1 and PG2 are the probabilities that the target is in the validation
region observed from sensors 1 and 2, respectively, C is a normalization constant such
that P m1a=0P m2b=0 a;bk = 1. The related likelihood function is
k := p
Y 1k;Y 2kjZk 1
=
m1X
a=0
m2X
b=0
p
Y 1k;Y 2k; a;bk jZk 1
(3.16)
where using the Bayes? formula, we obtain
p
Y 1k;Y 2k; a;bk jZk 1
= p
Y 1k;Y 2kj a;bk ;Zk 1
P[ a;bk ]
=
8
>>>>
>>>>
>>>>
>>>
<
>>>
>>>>
>>>>
>>>
>:
(1 PD1PG1)(1 PD2PG2)
[V1k ] m1[V2k ] m2 ; a = 0; b = 0
(1 PD1PG1)(PD2PG2)= m2
PG2[V2k ] m2 1 N
h
2(b)k ; 0;S2k
i
; a = 0; b = 1; ; m2
(PD1PG1)(1 PD2PG2)= m1
PG1[V1k ] m1 1 N
h
1(a)k ; 0;S1k
i
; a = 1; ; m1; b = 0
(PD1PG1)(PD2PG2)=( m1 m2)
PG1[V1k ] m1 1PG2[V2k ] m2 1 N
h
a;bk ; 0;Sk
i
; a = 1; ; m1; b = 1; ; m2:
(3.17)
Using state prediction ^xkjk 1 and its covariance Pkjk 1 at time k 1 (see Eqns. (2.12)
and (2.12), respectively), one computes the partial update ^xkjk and its covariance Pkjk
according to the standard PDAF [21] except that the augmented state is conditioned on
a;bk with data combination from sensors 1 and 2.
De ne the combined mode-conditioned innovations
k :=
m1X
a=0
m2X
b=0
(a;b)6=(0;0)
a;bk a;bk : (3.18)
62
Therefore, partial update of the state estimate is given by
^xa;bkjk := E
n
xkj a;bk ;Zk 1;Y 1k;Y 2k
o
= ^xkjk 1 +Wa;bk a;bk (3.19)
where Kalman gains, Wa;bk , are computed as
Wa;bk =
8>
>>>>
>>>
>>>>
>><
>>>
>>>>
>>>>
>>>:
W0;0k = 0; for a = 0; b = 0
Wa;0k = Pkjk 1[H10k [S1k] 1 0]; for a6= 0; b = 0
W0;bk = Pkjk 1[0 H20k [S2k] 1]; for a = 0; b6= 0
Wa;bk = Pkjk 1H0k[Sk] 1; for a6= 0; b6= 0;
(3.20)
and H0k =
h
H10k H20k
i
. Despite the fact that here we follow [46] for the association
events (hypotheses) a;bk to deal with existing uncertain measurements? origins, there
are \oversights" in [46, Sec. 2.1]. When both sensor measurements are associated with
the target, [46] states (p. 62, 2nd para., [46, Sec. 2.1]) that W[k] = [W1[k] W2[k]]
where W[k] = Kalman lter gain for the \overall" lter and Wi[k] = Kalman lter gain
corresponding to the measurements from sensor i (i = 1;2). This is NOT true. Since
the two sensors observe the same state, there is \cross-coupling" as described in Eqn.
(3.20).
The mode-conditioned update of the state estimate is
^xkjk := E
n
xkjMk;Zk 1;Y 1k;Y 2k
o
=
m1X
a=0
m2X
b=0
a;bk ^xa;bkjk 1 (3.21)
63
and the covariance of ^xkjk is (follow the steps in [12, 22] for the standard PDAF)
Pkjk := Pkjk 1
m1P
a=0
m2P
b=0
(a;b)6=(0;0)
a;bk Wa;bk Sa;bk Wa;b0k +
m1P
a=0
m2P
b=0
a;bk Wa;bk a;bk a;b0k Wa;b0k
m
1P
a=0
m2P
b=0
a;bk Wa;bk a;b
m
1P
a=0
m2P
b=0
a;bk Wa;bk a;b
0
:
(3.22)
Now we are ready to extend the above simultaneous measurement update technique to
a multiple model scenario.
3.3 Problem Formulation for the Multiple Model System
Assume that the target dynamics can be modeled by one of n hypothesis models.
The model set is denoted asMn :=f1; ;ngand there are total q sensors from which q,
or fewer (if the probability of target detection is less than one) or more (due to clutter),
measurement vectors are generated at a time. The event that model j is in e ect during
the sampling period (tk 1;tk] is denoted by Mjk. For the j-th hypothesized model(mode),
the state dynamics and measurements, respectively, are modeled as
xk = Fjk 1xk 1 +Gjk 1vjk 1 (3.23)
and
zlk = hl(xk) +wj;lk for l = 1; ;q; (3.24)
where xk is the system state at tk and of dimension nx, zlk is the (true) measurement
vector (i.e., due to the target) from sensor l at tk and of dimension nzl, Fjk 1 and Gjk 1
are the system matrices when model j is in e ect over the sampling period (tk 1;tk],
and hl is the nonlinear transformation of xk to zlk (l = 1; ;q) for model j. Henceforth,
64
time tk will be denoted by k. A rst-order linearized version of (3.24) is given by
zlk = Hj;lk xk +wj;lk for l = 1; ;q; (3.25)
whereHj;lk is the Jacobian matrix ofhl evaluated at some value of the estimate of statexk.
The process noise vjk 1 and the measurement noise wj;lk are mutually uncorrelated zero-
mean white Gaussian processes with covariance matrices Qjk 1 and Rj;lk , respectively.
At the initial time t0, the initial conditions for the system state under each model j
are assumed to be Gaussian random variables with the known mean xj0 and the known
covariance Pj0 . The probability of model j at t0, j0 = P[Mj0], is also assumed to be
known. The switching from model Mik 1 to model Mjk is governed by a nite-state
stationary Markov chain with known transition probabilities pij = P[MjkjMik 1].
The following notations and de nitions are used regarding the measurements at
sensor l. Note that, in general, at any time k, some measurements may be due to clutter
and some due to the target, i.e., there can be more than a single measurement at time
k at sensor l. The measurement set (not yet validated) generated by sensor l at time k
is denoted as
Zlk :=fzl(1)k ;zl(2)k ; ;zl(ml)k g (3.26)
where ml is the number of measurements generated by sensor l at time k. Variable
zl(i)k (i = 1; ;ml) is the ith measurement within the set. The cumulative set of
measurements (not yet validated) from sensor l up to time k is denoted as Zk(l) :=
fZl1;Zl2; ;Zlkg and the cumulative set of measurements (not yet validated) from all
sensors up to time k is denoted as Zk :=fZk(1);Zk(2); ;Zk(q)g where q is the number
of sensors. The validated set of measurements of sensor l at time k will be denoted
65
by Ylk, containing ml ( ml) measurement vectors. The cumulative set of validated
measurements from sensor l up to time k is denoted as
Yk(l) :=fYl1;Yl2; ;Ylkg: (3.27)
The cumulative set of validated measurements from all sensors up to time k is denoted
as
Zk :=fYk(1);Yk(2); ;Yk(q)g: (3.28)
Our goal is to nd the state estimate
^xkjk := EfxkjZkg (3.29)
and the associated error covariance matrix
Pkjk := Ef[xk ^xkjk][xk ^xkjk]0jZkg (3.30)
where x0k denotes the transpose of xk.
3.4 IMM/MSPDAF Algorithm for Simultaneous Measurement Update
We now modify the IMMPDA algorithms of [21] and [47] to derive the proposed
IMM/MSPDAF with simultaneous measurement update system. We con ne our atten-
tion to the case of 2 sensors; however, the algorithm can be easily adapted to the case of
arbitrary q sensors. We assume that all measurements observed from sensors are synchro-
nized with the same sampling rate (see Fig. 3.1). We will only brie y outline the basic
66
Local
sensor
Remote
sensor
Central
processor
j107
j107
j107j45j49
j107j45j49
j107j43j49
j107j43j49
j107j107j45j49 j107j43j49
Figure 3.1: Target tracking example using synchronized measurements observed from
multiple sensors
steps in ?one cycle? (i.e., processing needed to update for a new set of measurements) of
the IMM/MSPDA lter.
Assumed available : Given the state estimate ^xjk 1jk 1 := Efxk 1jMjk 1;Zk 1g,
the associated covariancePjk 1jk 1, and the conditional mode probability jk 1 :=P[Mjk 1jZk 1]
at time k 1 for each mode j2Mn.
Step 1. Interaction - mixing of the estimate from previous time (8j2Mn):
predicted mode probability:
j k := P[MjkjZk 1] =
nX
i=1
pij ik 1: (3.31)
mixing probability:
ijj := P[Mik 1jMjk;Zk 1] = pij ik 1= j k : (3.32)
67
mixed estimate:
^x0jk 1jk 1 := Efxk 1jMjk;Zk 1g=
nX
i=1
^xik 1jk 1 ijj: (3.33)
covariance of the mixed estimate:
P0jk 1jk 1 := Ef[xk 1 ^x0jk 1jk 1][xk 1 ^x0jk 1jk 1]0jMjk;Zk 1g (3.34)
=
nX
i=1
fPik 1jk 1 + [^xik 1jk 1 ^x0jk 1jk 1][^xik 1jk 1 ^x0jk 1jk 1]0g ijj:
Step 2. Predicted state and measurements for sensors 1 and 2 (8j2Mn) :
state prediction:
^xjkjk 1 := EfxkjMjk;Zk 1g= Fjk 1^x0jk 1jk 1: (3.35)
state prediction error covariance:
Pjkjk 1 := Ef[xk ^xjkjk 1][xk ^xjkjk 1]0jMjk;Zk 1g
= Fjk 1P0jk 1jk 1Fj0k 1 +Gjk 1Qjk 1Gj0k 1: (3.36)
The mode-conditioned predicted measurement for sensor l is
^zj;lk := hl(^xjkjk 1): (3.37)
68
Using the linearized version (3.25), the covariance of the mode-conditioned residual
j;l(i)k := zl(i)k ^zj;lk
is given by (assume q=2, the case of 2 sensors)
Sj;1k := Ef j;1(i)k j;1(i)0k jMjk;Zk 1g= Hj;1k Pjkjk 1Hj;10k +Rj;1k (3.38)
Sj;2k := Ef j;2(i)k j;2(i)0k jMjk;Zk 1g= Hj;2k Pjkjk 1Hj;20k +Rj;2k (3.39)
where Hj;lk is the rst order derivative (Jacobian matrix) of hl(:) evaluated at the state
prediction ^xjkjk 1 (see (3.37)). Note that (3.38) and (3.39) assume that zl(i)k originates
from the target.
As we mentioned earlier, since our approach to the problem deals with the multiple
simultaneous measurements arising from two separate sensors that are tracking a single
target through a common surveillance region, a method for combination of multiple
measurements has to be devised. In order to do this, the combined covariance Sjk of the
mode-conditioned residual obtained from (3.38) and (3.39) also needs to be considered
as follows:
Sjk := E
8>
><
>>:
2
66
4
j;1(i)k
j;2(i)k
3
77
5
j;1(i)0k j;1(i)0k
jMjk;Zk 1
9>
>=
>>;
=
2
66
64
Hj;1k
Hj;2k
3
77
75Pjkjk 1
Hj;10k Hj;20k
+
2
66
64
Rj;1k 0
0 Rj;2k
3
77
75: (3.40)
69
Step 3. Measurement validation for sensors 1 and 2 (8j2Mn):
Since the measurement validation process was well explained in Sec. 3.2 for non-
maneuvering target tracking scenario, we brie y describe the di erence between non-
maneuvering target tracking and maneuvering target tracking in this step. Following
[12, 21], one sets up a validation gate for sensor l centered at the mode-conditioned
predicted measurement, ^zj;lk . Let (jAj = det(A))
ja := arg
maxj2M
n
Sj;lk
: (3.41)
Then measurement zl(i)k (i=1,2, ,ml) is validated if and only if
[zl(i)k ^zja;lk ]0[Sja;lk ] 1[zl(i)k ^zja;lk ] < : (3.42)
The volume of the validation region with the threshold is
Vlk := cnzl nzl=2jSja;lk j1=2 (3.43)
where nzl is the dimension of the measurement and cnzl is the volume of the unit hyper-
sphere of this dimension (c1=2, c2= , c3=4 =3, etc.). Choice of is discussed in more
detail in [12, Sec. 2.3.2].
Step 4. State estimation with validated measurement from sensors 1 and 2
(8j2Mn) :
70
Using the de nition of the association events (hypotheses) a;bk (see Sec. 3.2), de ne
the mode-conditioned association event probabilities as
j;a;bk := P[ a;bk jMjk;Y 1k;Y 2k;Zk 1]: (3.44)
The mode-conditioned innovations j;a;bk can be de ned as
j;0;0k =
2
66
64
0nz1 1
0nz2 1
3
77
75; a = 0;b = 0
j;0;bk =
2
66
64
0nz1 1
j;2(b)k
3
77
75; a = 0; b = 1; ; m2
j;a;0k =
2
66
64
j;1(a)k
0nz2 1
3
77
75; a = 1; ; m1; b = 0
j;a;bk =
2
66
64
j;1(a)k
j;2(b)k
3
77
75; a = 1; ; m1; b = 1; ; m2:
(3.45)
71
Exploiting the di use model for clutter in [12, 21], it turns out that (see Appendix 1
for details)
j;0;0k = C(1 PD1PG1)(1 PD2PG2)(V1
k )
m1(V2
k )
m2 ; a = 0;b = 0
j;0;bk = C
PD2(1 PD1PG1)N
h
j;2(b)k ;0;Sj;2k
i
(V2k ) m2 1 m2 ; a = 0; b = 1; ; m2
j;a;0k = C
PD1(1 PD2PG2)N
h
j;1(a)k ;0;Sj;1k
i
(V1k ) m1 1 m1 ; a = 1; ; m1; b = 0
j;a;bk = C PD1PD2N
j;a;bk ;0;Sjk
m1 m2(V1k ) m1 1(V2k ) m2 1; a = 1; ; m1; b = 1; ; m2:
(3.46)
The likelihood function for each mode j is
jk := p
Y 1k;Y 2kjMjk;Zk 1
=
m1X
a=0
m2X
b=0
p
Y 1k;Y 2k; a;bk jMjk;Zk 1
(3.47)
where using the Bayes? formula, we obtain
p
Y 1k;Y 2k; a;bk jMjk;Zk 1
= p
Y 1k;Y 2kjMjk; a;bk ;Zk 1
P[ a;bk ]
=
8
>>>>
>>>>
>>>>
>>>
<
>>>
>>>>
>>>>
>>>
>:
(1 PD1PG1)(1 PD2PG2)
[V1k ] m1[V2k ] m2 ; a = 0; b = 0
(1 PD1PG1)(PD2PG2)= m2
PG2[V2k ] m2 1 N
h
j;2(b)k ; 0;Sj;2k
i
; a = 0; b = 1; ; m2
(PD1PG1)(1 PD2PG2)= m1
PG1[V1k ] m1 1 N
h
j;1(a)k ; 0;Sj;1k
i
; a = 1; ; m1; b = 0
(PD1PG1)(PD2PG2)=( m1 m2)
PG1[V1k ] m1 1PG2[V2k ] m2 1 N
h
j;a;bk ; 0;Sjk
i
; a = 1; ; m1; b = 1; ; m2:
(3.48)
72
Using ^xjkjk 1 (from (3.35)) and its covariance Pjkjk 1 (from (3.36)), one computes the
partial update ^xjkjk and its covariance Pjkjk according to the standard PDAF [21], except
that the augmented state is conditioned on a;bk with data combination from sensors 1
and 2.
De ne the combined mode-conditioned innovations
jk :=
m1X
a=0
m2X
b=0
(a;b)6=(0;0)
j;a;bk j;a;bk : (3.49)
The partial update of the state estimate is given by
^xj;a;bkjk := E
n
xkj a;bk ;Mjk;Zk 1;Y 1k;Y 2k
o
= ^xjkjk 1 +Wj;a;bk j;a;bk (3.50)
where the mode-conditioned Kalman gains, Wj;a;bk , are computed as
Wj;a;bk =
8
>>>>
>>>>
>>>>
>><
>>>>
>>>>
>>>>
>>:
Wj;0;0k = 0; for a = 0; b = 0
Wj;a;0k = Pjkjk 1
h
Hj;10k [Sj;1k ] 1 0
i
; for a6= 0; b = 0
Wj;0;bk = Pjkjk 1
h
0 Hj;20k [Sj;2k ] 1
i
; for a = 0; b6= 0
Wj;a;bk = Pjkjk 1Hj0k [Sjk] 1; for a6= 0; b6= 0
(3.51)
and Hjk0 =
h
Hj;10k Hj;20k
i
. Therefore, the mode-conditioned update of the state estimate
is
^xjkjk := E
n
xkjMjk;Zk 1;Y 1k;Y 2k
o
=
m1X
a=0
m2X
b=0
j;a;bk ^xj;a;bkjk 1 (3.52)
73
and the covariance of ^xjkjk is
Pjkjk := Pjkjk 1
m1P
a=0
m2P
b=0
(a;b)6=(0;0)
j;a;bk Wj;a;bk Sj;a;bk Wj;a;b0k +
m1P
a=0
m2P
b=0
j;a;bk Wj;a;bk j;a;bk j;a;b0k Wj;a;b0k
m
1P
a=0
m2P
b=0
j;a;bk Wj;a;bk j;a;b
m
1P
a=0
m2P
b=0
j;a;bk Wj;a;bk j;a;b
0
:
(3.53)
Step 5. Update of mode probabilities(8j2Mn) :
jk := P
h
MjkjYk
i
= 1C j k jk (3.54)
where C is a normalization constant such that P
j
jk = 1.
Step 6. Combination of the mode-conditioned estimates(8j2Mn) :
The nal state estimate update at time k is given by
^xkjk = Xj ^xjkjk jk (3.55)
and its covariance is given by
Pkjk = Xj
Pjkjk +
h
^xjkjk ^xkjk
ih
^xjkjk ^xkjk
i0
jk: (3.56)
3.5 Simulation Example
The following example of tracking a highly maneuvering target in clutter is consid-
ered.
The True Trajectory: The target starts at location [21689 10840 40] in Cartesian
74
coordinates in meters. The initial velocity (in m/s) is [-8.3 -399.9 0] and the target stays
at constant altitude with a constant speed of 400 m/s. Its trajectory is:
a straight line with constant velocity between 0 and 20s,
a coordinated turn (0.15 rad/s) with constant acceleration of 60 m=s2 between 20
and 35s,
a straight line with constant velocity between 35 and 55s,
a coordinated turn (0.1 rad/s) with constant acceleration of 40 m=s2 between 55
and 70s,
a straight line with constant velocity between 70 and 90s.
Fig. 3.2 shows the true trajectory of the target.
The Target Motion Models: The target motion models are patterned after [21]. In
each mode the target dynamics are modelled in Cartesian coordinates as
xk = Fxk 1 +Gvk 1
where the state of the target is position, velocity, acceleration in each of the 3 Cartesian
coordinates (x;y; and z). Thus xk is of dimension 9 (nx=9). Three maneuver models
are considered in the following discussion. The system matrices F, G, are de ned as
F =
2
66
66
66
66
4
Fb 0 0
0 Fb 0
0 0 Fb
3
77
77
77
77
5
and G =
2
66
66
66
66
4
Gb 0 0
0 Gb 0
0 0 Gb
3
77
77
77
77
5
:
75
Model 1. Nearly constant velocity model with zero mean perturbation in accel-
eration
F1b =
2
66
66
66
66
4
1 T 0
0 1 0
0 0 0
3
77
77
77
77
5
and G1b =
2
66
66
66
66
4
T2
2
T
0
3
77
77
77
77
5
;
where T is the sampling period. The standard deviation of the process noise of
M1 is 5m=s2 (as in [21]).
Model 2. Wiener process acceleration (nearly constant acceleration motion)
F2b =
2
66
66
66
66
4
1 T T22
0 1 T
0 0 1
3
77
77
77
77
5
and G2b =
2
66
66
66
66
4
T2
2
T
1
3
77
77
77
77
5
The standard deviation of the process noise of M2 is 7.5m=s2 (as in [21]).
Model 3. Wiener process acceleration (model with large acceleration increments,
for the onset and termination of maneuvers), with F3b = F2b and G3b = G2b. The
standard deviation of the process noise of M3 is 40m=s2 (as in [21]).
76
The initial model probabilities are 10 = 0:8, 20 = 0:1 and 30 = 0:1. The mode
switching probability matrix is given by (as in [21])
2
66
66
66
66
4
p11 p12 p13
p21 p22 p23
p31 p32 p33
3
77
77
77
77
5
=
2
66
66
66
66
4
0:8 0:0 0:2
0:0 0:8 0:2
0:3 0:3 0:4
3
77
77
77
77
5
:
The Sensors: Two sensors are used to obtain the measurements:
Case 1. Sensor 1 and Sensor 2 are located at [x1;y1;z1]=[-4000 4000 0] m and
[x2;y2;z2]=[5000 0 0] m, respectively, and the central processor is collocated with
sensor 1 platform.
Case 2. Sensor 1 and Sensor 2 are collocated at [x1;y1;z1]=[-4000 4000 0] m
together with the central processor.
The measurements from sensor l for model j are zlk = hl(xk) + wj;lk for l = 1 and
2, re ecting range and azimuth angle for sensor 1 (radar) and azimuth and elevation
angles for sensor 2 (infrared). The range, azimuth, and elevation angle transformations,
respectively, are given by
rl = f(x xl)2 + (y yl)2 + (z zl)2g1=2 (3.57)
al = tan 1[(y yl)=(x xl)] (3.58)
el = tan 1[(z zl)=f(x xl)2 + (y yl)2g1=2]: (3.59)
77
The Jacobian matrices of hl(:) for sensors 1 and 2 are
Hj;1 =
2
66
64
cos (e1) cos (a1) 0 0 cos (e1) sin (a1) 0 0 sin (e1) 0 0
sin (a1)r1 cos (e1) 0 0 cos (a1)r1 cos (e1) 0 0 0 0 0
3
77
758j (3.60)
Hj;2 =
2
66
64
sin (a2)r2 cos (e2) 0 0 cos (a2)r2 cos (e2) 0 0 0 0 0
sin (e2) cos (a2)r2 0 0 sin (e2) sin (a2)r2 0 0 cos (e2)r2 0 0
3
77
758j; (3.61)
respectively. The measurement noise wj;lk for sensor l is assumed to be zero-mean white
Gaussian with known covariances, R1=diag[qr;qa1]=diag[400m2,49 mrad2] with qr and
qa1 denoting the variances for the radar range and azimuth measurement noises, re-
spectively, and R2 = diag[qa2;qe] = diag[4mrad2;4mrad2] with qa2 and qe denoting the
variances for the infrared sensor azimuth and elevation measurement noises, respectively.
The sampling interval was T=1s and it was assumed that the probability of detection
PD=0.997 for both sensors.
The Clutter: For generating false measurements in simulations, the clutter was as-
sumed to be Poisson distributed with expected number of 1 = 50 10 6=m mrad for
sensor 1 and 2 = 3:5 10 4=m2 mrad for sensor 2. These statistics were used for
generating the clutter in all simulations. However, a nonparametric clutter model was
used for implementing all the algorithms for target tracking.
Other Parameters: The gates for setting up the validation regions for both the sensors
were based on the threshold =16 corresponding to a gate probability PG=0.9997.
78
Simulation Results: Fig. 3.3 shows the RMSE (root mean-square error) in posi-
tion for the proposed IMM/MSPDAF and the standard sequential IMM/MSPDAF [21]
based on 200 Monte Carlo runs. It is seen that the proposed simultaneous measurement
updating can signi cantly improve the accuracy of track estimation during the periods
following the onset of the target maneuvers. The rst maneuver starts at 20 sec and in
Fig. 3.3 one can see a signi cant improvement from 22 sec through 26 sec. The second
maneuver starts at 55 sec and in Fig. 3.3 one can see a signi cant improvement from 58
sec through 59 (in (a) or 60 in (b)) sec. That is, simultaneous measurement updating
responds faster to maneuvers. Once the target is \settled" in a particular mode, there
is insigni cant di erences between the two approaches. It is also seen that tracking
with separated sensors can improve the accuracy of track estimation compared to using
collocated sensors.
To assess the computational requirements of the two approaches, we computed the
CPU time needed to execute 90 time steps in each run (averaged 100 Monte Carlo runs
excluding data/clutter generation) in MATLAB 6.5 on a 2.8 GHz (Mobile) Pentium 4
operating under Windows XP (professional). The standard sequential IMM/MSPDAF
needs 0.4862 secs (for all 90 time steps) compared to 0.5806 secs required by proposed
IMM/MSPDAF. Thus there is, on the average, a 20.3% increase in computational cost.
3.6 Conclusions
We investigated an IMM/MSPDAF algorithm with simultaneous measurement up-
date for tracking a highly maneuvering target in clutter. A past approach [45] us-
ing parallel sensor processing has ignored certain data association probabilities leading
79
0.8 1 1.2 1.4 1.6 1.8 2 2.2x 1040
2000
4000
6000
8000
10000
12000
14000
16000
18000
x coord (m)
y coord (m)
t=0
0 10 20 30 40 50 60 70 80 90?400
?300
?200
?100
0
100
200
300
400
time (sec)
velocity (m/sec)
xy
0 10 20 30 40 50 60 70 80 90?60
?40
?20
0
20
40
60
time (sec)
acc. (m/sec
2 )
xy
0 10 20 30 40 50 60 70 80 900
10
20
30
40
50
60
70
time (sec)
acceleration magnitude (m/sec
2 )
Figure 3.2: Trajectory of the maneuvering target (read left to right, top to bottom). (a)
Position in xy plane. (b) x and y velocities. (c) x and y accelerations. (d) magnitude of
accelerations
80
0 10 20 30 40 50 60 70 80 900
10
20
30
40
50
60
70
80
90 Sequential vs Simultaneous Update: collocated sensors
time (sec)
Sequential update [21]Simultaneous update (proposed)
RMS error in position (m)
0 10 20 30 40 50 60 70 80 900
10
20
30
40
50
60
70
80
90 Sequential vs Simultaneous Update: separated sensors
time (sec)
Sequential update [21]Simultaneous update (proposed)
RMS error in position (m)
Figure 3.3: Performance of the simultaneous measurement updating IMM/MSPDAF in
terms of the RMSE (root mean square error) in position (read top to bottom). (a)
proposed IMM/MSPDAF vs standard sequential IMM/MSPDAF [21] with collocated
sensors (Case 2). (b) proposed IMM/MSPDAF vs standard sequential IMM/MSPDAF
[21] with separated sensors (Case 1).
81
to an inaccurate implementation. Another existing approach [46] applies only to non-
maneuvering targets. Our proposed approach has extended the multisensor approach of
[46] to maneuvering targets by employing a switching multiple model approach.
The proposed algorithm was illustrated via a simulation example where it outper-
formed a standard IMM/MSPDAF algorithm with sequential updating [21] during the
periods following the onset of the target maneuvers. The simultaneous updating is ex-
pected to be more accurate [12] since it considers all association hypotheses coupled
across multisensor while the sequential updating considers the separate hypothesis for
each sensor. This improvement in accuracy is seen in our simulation example only during
the periods following the onset of the target maneuvers. Once the target is settled in
a particular mode, there is insigni cant di erences between the two approaches. The
improvement in accuracy comes at the expense of a slight increase (20%) in the compu-
tational cost.
82
Chapter 4
Multisensor Tracking of a Maneuvering Target in Clutter with
Asynchronous Measurements using Augmented State IMMPDA Filtering
and Simultaneous Measurement Update
In this chapter, we discuss a suboptimal ltering algorithm for tracking a highly
maneuvering target in a cluttered environment using multiple sensors dealing with pos-
sibly asynchronous (time delayed) measurements. The ltering algorithm is developed
by applying the basic IMM approach, the PDA technique, and asynchronous measure-
ment updating for state-augmented system estimation for the target. A state augmented
approach is developed to estimate the time delay between local and remote sensors. A
multisensor probabilistic data association lter is developed for parallel sensor processing
for target tracking under clutter. The algorithm is illustrated via a highly maneuvering
target tracking simulation example where two sensors, a radar and an infrared sensor,
are used. Compared with an existing IMMPDA ltering algorithm with the assumption
of synchronous (no delay) measurements sensor processing which is presented in Chapter
3, the proposed algorithm achieves considerable improvement (especially in the case of
larger delays) in the accuracy of track estimation.
83
4.1 Introduction
We extend our simultaneous measurement update technique presented in Chapter 3
to asynchronous (delayed) measurements problem. In target tracking systems measure-
ments are typically collected in \scans" or \frames" and then transmitted to a process-
ing center [27, 35]. Asynchronous (delayed) measurements arise in a multisensor central
tracking system due to communication network delays, varying preprocessing times at
the sensor platforms and possibly lack of sampling time synchronization among sensor
platforms (see Fig. 4.1).
One of the asynchronous measurement problems is that of out-of-sequence measure-
ments (OOSM) where measurements at various sensors may arrive out-of-sequence (not
in correct time order) at the central processor. OOSM has been considered using IMM
[35, 36, 37]. In this chapter we do not consider OOSM (OOSM scenario is presented in
Chapter 5.) but instead consider \in-sequence" measurements with a xed-but-unknown
relative time-delay among sensor measurements. Various sensor measurements are as-
sumed to be at the same rate but not necessarily time synchronized. All measurements
over one sampling interval (based on the local clock of the central processor) are collected
at the central processor, attributed to one time instant and processed simultaneously.
We exploit IMM and PDA techniques. It is assumed that a track has been formed
(initiated) and the objective of this work is to investigate xed-but-unknown relative
time-delay (measurement timing mismatch) arising in a multisensor central tracking
system.
In [26], xed-lag smoothing techniques have been investigated using IMM algo-
rithm combined with PDA lter in a multiple sensor scenario to propose a combined
84
IMM/MSPDAF (interacting multiple model multiple sensor probabilistic data associa-
tion lter). We exploit the basic structure of [21] in combination with a state-augmented
approach to deal with the xed-but-unknown relative time-delay. In [21] and [45] it is
assumed that the sensors are collocated and (time) synchronized with the sampling rate.
In contrast, the sensor collocation and (time) synchronization are no longer assumed in
this chapter. Also, we use simultaneous measurement updating of the state estimates.
In this chapter, we also extend the multisensor approach of [46] to maneuvering targets
(see Step 4 in Sec.4.4).
This chapter is organized as follows. Sec. 4.2 presents the problem formulation.
Sec. 4.3 describes the state-augmented system approach. Sec. 4.4 describes the proposed
augmented state IMM/MSPDAF (AS-IMM/MSPDAF) algorithm for asynchronous mea-
surements. Simulation results using the proposed algorithm for a realistic problem are
given in Sec. 4.5. Finally, Sec. 4.6 presents a discussion of the results and some conclu-
sions.
4.2 Problem Formulation for Asynchronous measurements
The system dynamics and measurement equation which is synchronized with central
processor are the same as in the problem formulation for the multiple model system in
Chapter 3 (see (3.23)-(3.25)). Hence we do not reiterate them in this chapter. The basic
scenario of multisensor target tracking system dealing with asynchronous measurements
can be seen from Fig. 4.1. Assume that there is a xed-but-unknown relative time
delay dk between the remote sensor clock and the central processor clock at sample
time tk as shown in Fig. 4.1. This time delay could be due to unsynchronized clocks
at the two locations or due to inherent delay due to congestion, insu cient bandwidth
85
Local
sensor
Remote
sensor
Central
processor
In-sequence
w/ delay j100
j107j108
j107
j100j108
j107
j100j108
j107
j107
j107j45j49
j100j108
j107j45j49
j100j108
j107j45j49
j107j45j49
j100
j107j108
In-sequence
w/ delay j100
j107j108
Figure 4.1: Target tracking with a xed-but-unknown relative time delay dk between the
remote sensor clock and the central processor clock at sample time tk
etc. in the communication link between the remote sensor platform and the central
processor. The measurements from sensor l are sent to the central processor where all
measurements collected between local sampling interval (tk 1;tk] are attributed to time
tk. The state dynamics and measurements reported from the remote sensor platform at
time tkdl (henceforth will be denoted by kdl) to the central processor at time tk can be
modeled as (see Sec. 4.5 for some details)
xkdl = Fjkdl;k 1xk 1 +Gjkdl;k 1vjk 1 (4.1)
and
zlk = hl(xkdl) +wj;lk (4.2)
where tkdl = tk dkl and dkl is the time di erence between the sampling time at the
central processor and the measurement time at the local sensor (assume that 0 dkl >>>
>>>>
>>>>
>>>
>>><
>>>>
>>>>
>>>
>>>>
>>>:
2
66
4
z2k
z2k 1
3
77
5; for 11
z2k; for 10
z2k 1; for 01
no measurement; for 00:
(5.6)
5.4 State-Augmented System
De ne the augmented state ~xk from xk as
~x0k = [x0k;v0k;x0k 1;v0k 1;x0k 2;v0k 2] (5.7)
where x0k denotes the transpose of xk. Assume that there is a xed-and-known delay,
dkl, between the central processor and the remote sensor l platform. Using the above
de nitions and (5.1), the augmented state equation may be written more compactly as
~xk = ~Fjk;k 1~xk 1 + ~Gjk;k 1vjk; (5.8)
~xkdl = ~Fjkdl;k 1~xk 1 + ~Gjkdl;k 1vjkdl (5.9)
where ~Fjk;k 1, ~Gjk;k 1, ~Fjkdl;k 1, and ~Gjkdl;k 1 are de ned in Sec. 5.6 (see (5.53)-(5.60)).
Note that the process noise in (5.8) is vjk (at time k, not at time k 1). Using the
augmented state (5.7) the counterparts to (5.2) and (5.4), respectively, are
z1k = h1(~xk) +wj;lk = h1([I;0;0;0;0;0]~xk) +wj;lk (5.10)
114
and
~z2k =
8
>>>>
>>>>
>>>>
>>>
>>><
>>>>
>>>>
>>>
>>>>
>>>:
2
66
4
z2k
z2k 1
3
77
5 = h2
0
BB
@
2
66
4
0 0 FjdlGjdl 0 0
0 0 0 0 FjdlGjdl
3
77
5~xk
1
CC
A+
2
66
4
wj;2k
wj;2k 1
3
77
5; for 11
z2k = h2
h
0 0 FjdlGjdl 0 0
i
~xk +wj;2k ; for 10
z2k 1 = h2
h
0 0 0 0 FjdlGjdl
i
~xk +wj;2k 1; for 01
no measurement; for 00:
(5.11)
The following notations and de nitions are used regarding the measurements at sensor l.
Note that, in general, at any time some measurements may be due to clutter and some
due to the target, i.e., there can be more than a single measurement at time k at sensor
l (l = 1;2). The measurement set (not yet validated) generated by sensor l at time k is
denoted as
Zlk :=fzl(1)k ;zl(2)k ; ;zl(ml)k g (5.12)
where ml is the number of measurements generated by sensor l at time k. Variable zl(i)k
(i = 1, ,ml) is the ith measurement within this set. The validated set of measurements
of sensor l at time k will be denoted by Ylk, containing ml ( ml) measurement vectors.
The cumulative set of validated measurements from sensor l up to time k is denoted as
Yk(l) :=fYl1;Yl2; ;Ylkg: (5.13)
The cumulative set of validated measurements from all sensors up to time k is denoted
as
Zk :=fYk(1);Yk(2); ;Yk(q)g (5.14)
115
where q is the number of sensors.
Our goal is to nd the state estimate
^~xkjk := Ef~xkjZkg (5.15)
and the associated error covariance matrix
~Pkjk := Ef[~xk ^~xkjk][~xk ^~xkjk]0jZkg (5.16)
where x0k denotes the transpose of xk.
5.5 IMM/MSPDAF Algorithm for Asynchronous and Possibly Out-of-Sequence
Measurements
We now modify the IMMPDA algorithms of [47] and [25] to apply to the multi-
sensor asynchronous measurements system. We con ne our attention to the case of 2
sensors; however, the algorithm can be adapted to the case of arbitrary q sensors. We
will only brie y outline the basic steps in ?one cycle? (i.e., processing needed to update
for a new set of measurements) of the IMM/MSPDA lter.
Assumed available: Given the state estimate ^~xjk 1jk 1 := Ef~xk 1jMjk 1;Zk 1g, the
associated covariance ~Pjk 1jk 1, and the conditional mode probability jk 1 :=P[Mjk 1jZk 1]
at time k 1 for each mode j2Mn.
Step 1. Interaction { mixing of the estimate from the previous time (8j 2
Mn) :
116
predicted mode probability:
j k := P[MjkjZk 1] = X
i
pij ik 1: (5.17)
mixing probability:
ijj := P[Mik 1jMjk;Zk 1] = pij ik 1= j k : (5.18)
mixed estimate:
^~x0jk 1jk 1 := Ef~xk 1jMjk;Zk 1g= X
i
^~xik 1jk 1 ijj: (5.19)
covariance of the mixed estimate:
~P0jk 1jk 1 := En[~xk 1 ^~x0jk 1jk 1][~xk 1 ^~x0jk 1jk 1]0jMjk;Zk 1o (5.20)
= X
i
n~
Pik 1jk 1 + [^~xik 1jk 1 ^~x0jk 1jk 1][^~xik 1jk 1 ^~x0jk 1jk 1]0
o
ijj:
Step 2. Predicted state and measurements for sensors 1 and 2 (8j2Mn) :
state prediction:
^~xjkjk 1 := Ef~xkjMjk;Zk 1g= ~Fjk 1^~x0jk 1jk 1: (5.21)
state prediction error covariance:
~Pjkjk 1 := Ef[~xk ^~xjkjk 1][~xk ^~xjkjk 1]0jMjk;Zk 1g
117
= ~Fjk 1 ~P0jk 1jk 1 ~Fj0k 1 + ~Gjk 1Qjk 1 ~Gj0k 1: (5.22)
The mode-conditioned predicted measurements for sensors 1 and 2 are
^zj;1k := h1(^~xjkjk 1); (5.23)
^~zj;2k := h2 ^~xjkjk 1 =
8>
>>>>
>>>>
>>>>
>>>
>><
>>>>
>>>
>>>>
>>>>
>>>:
2
66
4
^zj;2kjk 1
^zj;2k 1jk 1
3
77
5 for 11
^zj;2kjk 1 for 10
^zj;2k 1jk 1 for 01
no measurement for 00:
(5.24)
Using the linearized version (5.3) and (5.4), the covariance of the mode-conditioned
residual
j;1(i)k := z1(i)k ^zj;1k (5.25)
j;2(i;r)k := z2(i;r)k ^~zj;2k =
8
>>>>
>>>>
>>>
>>>>
>>><
>>>>
>>>
>>>>
>>>>
>>>:
2
66
4
z2(i)k ^zj;2kjk 1
z2(r)k 1 ^zj;2k 1jk 1
3
77
5 for 11
z2(i)k ^zj;2kjk 1; for 10
z2(r)k 1 ^zj;2k 1jk 1; for 01
0; for 00
(5.26)
are given by
Sj;1k := Ef j;1(i)k j;1(i)0k jMjk;Zk 1g= ~Hj;1k ~Pjkjk 1 ~Hj;10k +Rj;1k (5.27)
118
Sj;2k :=
8>
>>>
>>>>
>>>>
>><
>>>
>>>>
>>>>
>>>
:
2
66
64
~Hj;2k
~Hj;2k 1
3
77
75 ~Pjkjk 1
~Hj;20k ~Hj;20k 1
+
2
66
64
Rj;2k 0
0 Rj;2k 1
3
77
75; for 11
~Hj;2k ~Pjkjk 1 ~Hj;20k +Rj;2k ; for 10
~Hj;2k 1 ~Pjkjk 1 ~Hj;20k 1 +Rj;2k 1; for 01
(5.28)
where ~Hj;lk ( ~Hj;lk 1) is the rst order derivative (Jacobian matrix) of hl(:) evaluated at the
state prediction ^~xjkjk 1 (^~xjk 1jk 1). Note that (5.27) and (5.28) assume that z1(i)k , z2(i)k ,
and z2(r)k 1 originate from the target. The results (5.27) and (5.28) do not depend upon
the actual measurements.
As mentioned earlier, since our approach to the problem deals with the OOSM as
well as multiple simultaneous measurements [46, 48] arising from two separate sensors
that are tracking a single target through a common surveillance region, a method for
fusion of multiple measurements has to be devised. In order to do this, now the combined
covariance of the mode-conditioned residual obtained from (5.27) and (5.28) also needs
to be considered as follows
Sjk := E
8>
><
>>:
2
66
4
j;1(i)k
j;2(i;r)k
3
77
5
j;1(i)0k j;2(i;r)0k
jMjk;Zk 1
9>
>=
>>; (5.29)
119
=
8
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>>>>
>><
>>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>>>:
2
66
66
66
66
4
~Hj;1k
~Hj;2k
~Hj;2k 1
3
77
77
77
77
5
~Pjkjk 1
~Hj;10k ~Hj;20k ~Hj;20k 1
+
2
66
66
66
66
4
Rj;1k 0 0
0 Rj;2k 0
0 0 Rj;2k 1
3
77
77
77
77
5
; for 11
2
66
64
~Hj;1k
~Hj;2k
3
77
75 ~Pjkjk 1
~Hj;10k ~Hj;20k
+
2
66
64
Rj;1k 0
0 Rj;2k
3
77
75; for 10
2
66
64
~Hj;1k
~Hj;2k 1
3
77
75 ~Pjkjk 1
~Hj;10k ~Hj;20k 1
+
2
66
64
Rj;1k 0
0 Rj;2k 1
3
77
75; for 01
~Hj;1k ~Pjkjk 1 ~Hj;10k +Rj;1k ; for 00:
(5.30)
Step 3. Measurement validation for sensors 1 and 2 (8j2Mn) :
There is uncertainty regarding the measurements? origins. Therefore, we perform
validation for each target separately. One sets up a validation gate for sensor l centered
at the mode-conditioned predicted measurement, ^zj;lk . Let (jAj= det(A))
ja := arg
maxj2M
n
Sj;lk
and ja := arg
maxj2M
n
Sj;2k 1
: (5.31)
Then measurement zl(i)k (z2(r)k 1) is validated if and only if
[zl(i)k ^zja;lk ]0[Sja;lk ] 1[zl(i)k ^zja;lk ] <
[zl(i)k ^z ja;2k 1]0[S ja;2k 1] 1[zl(i)k ^z ja;2k 1] <
(5.32)
120
where is an appropriate threshold. The volume of the validation region with the
threshold is [12, Sec. 2.3.2]
Vlk := cnzl nzl=2
Sja;lk
1=2
V 2k 1 := cnz2 nz2=2
S ja;2k 1
1=2
: (5.33)
After performing the validation for each target separately, we deal with all the validated
data for measurement fusion.
Step 4. State estimation with validated measurement from sensors 1 and 2
(8j2Mn) :
From among all the raw measurements from sensorl at timek, i.e., Zlk :=fzl(1)k ;zl(2)k ; ;zl(ml)k g,
de ne the set of validated measurements for sensor l at time k as
Ylk :=fyl(1)k ;yl(2)k ; ;yl( ml)k g (5.34)
where ml is total number of validated measurements for sensor l at time k and
yl(i)k := zl(li)k (5.35)
where 1 l1 >>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
<
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>:
(1 PD1PG1)(1 PD2PG2)2
(V1k ) m1(V2k ) m2(V2k 1) mo2 ; a = b = c = 0
PD1(1 PD2PG2)2N
h
j;1(a)k ;0;Sj;1k
i
(V1k ) m1 1 m1 ; a = 1; ; m1; b = c = 0
PD2(1 PD1PG1)(1 PD2PG2)N
h
j;2(b)k ;0;Sj;2k
i
(V2k ) m2 1 m2 ; a = 0; b = 1; ; m2; c = 0
N
h
j(a;b;0)k ;0;Sj;a;bk
i
PD1PD2(1 PD2PG2)
m1 m2(V1k ) m1 1(V2k ) m2 1 ; a = 1; ; m1; b = 1; ; m2; c = 0
(1 PD1PG1)(1 PD2PG2)PD2N
h
j;2(c)k 1 ;0;Sj;2k 1
i
(V2k 1) mo2 1 mo2 ; a = b = 0; c = 1; ; m
o2
PD1PD2(1 PD2PG2)N
h
~ j(a;0;c)k ;0;~Sj;a;ck
i
(V1k ) m1 1(V2k 1) mo2 1 m1 ; a = 1; ; m1; b = 0; c = 1; ; m
o2
P2D2(1 PD1PG1)N
h
~ j(0;b;c)k ;0;~Sj;b;ck
i
(V2k ) m2 1(V2k 1) mo2 1 m2 mo2 ; a = 0; b = 1; ; m2; c = 1; ; m
o2
PD1P2D2N
h
~ j(a;b;c)k ;0;~Sj;a;b;ck
i
(V1k ) m1 1(V2k ) m2 1(V2k 1) mo2 1 m1 m2 mo2; a = 1; ; m1; b = 1; ; m2; c = 1; ; m
o2
(5.42)
Using ^~xjkjk 1 (from (5.21)) and its covariance ~Pjkjk 1 (from (5.22)), one computes the
partial update ^~xjkjk and its covariance ~Pjkjk according to the standard PDAF [21], except
that the augmented state is conditioned on a;b;ck with data fusion from sensors 1 and 2.
De ne the combined mode-conditioned innovations
jk =
m1X
a=0
m2X
b=0
m12X
c=0
j;a;b;ck j;a;b;ck : (5.43)
Therefore, partial update of the state estimate
^~xj;a;b;ckjk := Enxkj a;b;ck ;Mjk;Zk 1;Y 1k;Y 2k;Y 2k 1o = ^~xjkjk 1 +Wj;a;b;ck j;a;b;ck (5.44)
125
where Kalman gains, Wj;a;b;ck , are computed as
8>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>>>>
<
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>>>>
>:
Wj;0;0;0k = 0; for a = b = c = 0
Wj;a;0;0k = ~Pjkjk 1
h ~
Hj;10k [Sj;1k ] 1 0 0
i
; for a6= 0; b = c = 0
Wj;0;b;0k = ~Pjkjk 1
h
0 ~Hj;20k [Sj;2k ] 1 0
i
; for a = 0; b6= 0; c = 0
Wj;a;b;0k = ~Pjkjk 1
h
[ ~Hj;10k ~Hj;20k ][Sjk] 1 0
i
; for a6= 0; b6= 0; c = 0
Wj;0;0;ck = ~Pjkjk 1
h
0 0 ~Hj;20k 1[Sj;2k 1] 1
i
; for a = b = 0; c6= 0
Wj;a;0;ck = ~Pjkjk 1
h ~
Hj;10k [Sj;1k ] 1 0 ~Hj;20k 1[Sj;2k 1] 1
i
; for a6= 0; b = 0; c6= 0
Wj;0;b;ck = ~Pjkjk 1
h
0 [ ~Hj;20k ~Hj;20k 1][Sj;2k ] 1
i
; for a = 0; b6= 0; c6= 0
Wj;a;b;ck = ~Pjkjk 1
h
[ ~Hj;10k ~Hj;20k ~Hj;20k 1][Sjk] 1
i
; for a6= 0; b6= 0; c6= 0:
(5.45)
Therefore, mode-conditioned update of the state estimate
^~xjkjk := EnxkjMjk;Zk 1;Y 1k;Y 2k;Y 2k 1o =
m1X
a=0
m2X
b=0
m12X
c=0
j;a;b;ck ^~xj;a;b;ckjk 1 (5.46)
and covariance of ^~xjkjk
~Pjkjk = ~Pjkjk 1 m1P
a=0
m2P
b=0
mo2P
c=0
(a;b;c)6=(0;0;0)
j;a;b;ck Wj;a;b;ck Sj;a;b;ck Wj;a;b;c0k
+
m1P
a=0
m2P
b=0
m12P
c=0
j;a;b;ck Wj;a;b;ck j;a;b;ck j;a;b;c0k Wj;a;b;c0k
"
m1P
a=0
m2P
b=0
m12P
c=0
j;a;b;ck Wj;a;b;ck j;a;b;c
#"
m1P
a=0
m2P
b=0
m12P
c=0
j;a;b;ck Wj;a;b;ck j;a;b;c
#0
:
(5.47)
126
Step 5. Update of mode probabilities (8j2Mn) :
jk := P
h
MjkjZk
i
= 1C j k jk (5.48)
where C is a normalization constant such that P
j
jk = 1.
Step 6. Combination of the mode-conditioned estimates (8j2Mn) :
The nal augmented state estimate update at time k is given by
^~xkjk = X
j
^~xjkjk jk (5.49)
and its covariance is given by
~Pkjk = X
j
~Pjkjk +h^~xjkjk ^~xkjkih^~xjkjk ^~xkjki0
jk: (5.50)
From the nal augmented state (see (5.49)), the state ltered vector ^xkjk and the state
smoothing vector ^xk 1jk can be easily obtained.
5.6 Simulation Example
The following example of tracking a highly maneuvering target in clutter is consid-
ered. The target starts at location [21689 10840 40] in Cartesian coordinates in meters.
The initial velocity (in m/s) is [-8.3 -399.9 0] and the target stays at constant altitude
with a constant speed of 400 m/s. Its trajectory is a straight line with constant velocity
between 0 and 20s, a coordinated turn (0.15 rad/s) with constant acceleration of 60
m/s2 between 20 and 35s, a straight line with constant velocity between 35 and 55s, a
127
coordinated turn (0.1 rad/s) with constant acceleration of 40 m/s2 between 55 and 70s,
and a straight line with constant velocity between 70 and 90s. The target motion models
are patterned and modi ed after [21]. In each mode the target dynamics are modeled in
Cartesian coordinates as
~xk = ~Fjk;k 1~xk 1 + ~Gjk;k 1vjk; (5.51)
~xkdl = ~Fjkdl;k 1~xk 1 + ~Gjkdl;k 1vjkdl (5.52)
where the augmented state of the target consists of position, velocity, acceleration, and
the process noise in each of the three Cartesian coordinates (x;y; and z) at tk, tk 1 and
tk 2. Thus both ~xk and ~xkdl are of dimension 36 (nx = 36). Three maneuver models
are considered in the following discussion. From (5.8) and (5.9), ~Fjk;k 1, ~Gjk;k 1, ~Fjkdl;k 1
and ~Gjkdl;k 1 are de ned as
~Fjk;k 1 =
2
66
66
66
66
66
66
66
66
66
4
Fjk;k 1 Gjk;k 1 0 0 0 0
0 0 0 0 0 0
I 0 0 0 0 0
0 I 0 0 0 0
0 0 I 0 0 0
0 0 0 I 0 0
3
77
77
77
77
77
77
77
77
77
5
; ~Gjk;k 1 =
2
66
66
66
66
66
66
66
66
66
4
0
I
0
0
0
0
3
77
77
77
77
77
77
77
77
77
5
(5.53)
128
~Fjk
dl;k 1 =
2
66
66
66
66
66
66
66
66
66
4
Fjkdl;k 1 Gjkdl;k 1 0 0 0 0
0 0 0 0 0 0
I 0 0 0 0 0
0 I 0 0 0 0
0 0 I 0 0 0
0 0 0 I 0 0
3
77
77
77
77
77
77
77
77
77
5
; ~Gjkdl;k 1 =
2
66
66
66
66
66
66
66
66
66
4
0
I
0
0
0
0
3
77
77
77
77
77
77
77
77
77
5
(5.54)
where the system matrices Fk;k 1, Gk;k 1, Fkdl;k 1 and Gkdl;k 1 are de ned as
Fjk;k 1 =
2
66
66
66
66
4
Fj 0 0
0 Fj 0
0 0 Fj
3
77
77
77
77
5
; Gjk;k 1 =
2
66
66
66
66
4
Gj 0 0
0 Gj 0
0 0 Gj
3
77
77
77
77
5
; (5.55)
Fjkdl;k 1 =
2
66
66
66
66
4
Fjd 0 0
0 Fjd 0
0 0 Fjd
3
77
77
77
77
5
; Gjkdl;k 1 =
2
66
66
66
66
4
Gjd 0 0
0 Gjd 0
0 0 Gjd
3
77
77
77
77
5
: (5.56)
Model 1. Nearly constant velocity model with zero mean perturbation in acceleration
F1 =
2
66
66
66
66
4
1 T 0
0 1 0
0 0 0
3
77
77
77
77
5
; G1 =
2
66
66
66
66
4
T2
2
T
0
3
77
77
77
77
5
; (5.57)
129
F1d =
2
66
66
66
66
4
1 (T dkl) 0
0 1 0
0 0 0
3
77
77
77
77
5
; G1d =
2
66
66
66
66
4
(T dkl)2
2
(T dkl)
0
3
77
77
77
77
5
; (5.58)
where T is the sampling period. The standard deviation of the process noise of M1 is 5
m/s2 (as in [21]).
Model 2. Wiener process acceleration (nearly constant acceleration motion)
F2 =
2
66
66
66
66
4
1 T T22
0 1 T
0 0 1
3
77
77
77
77
5
; G2 =
2
66
66
66
66
4
T2
2
T
1
3
77
77
77
77
5
; (5.59)
F2d =
2
66
66
66
66
4
1 (T dkl) (T dkl)22
0 1 (T dkl)
0 0 1
3
77
77
77
77
5
; G2d =
2
66
66
66
66
4
(T dkl)2
2
(T dkl)
1
3
77
77
77
77
5
: (5.60)
The standard deviation of the process noise of M2 is 7.5 m/s2 (as in [21]).
Model 3. Wiener process acceleration (model with large acceleration increments, for the
onset and termination of maneuvers), with F3 = F2, G3 = G2, F3d = F2d and G3d = G2d.
The standard deviation of the process noise of M3 is 40 m/s2 (as in [21]).
130
The initial model probabilities are 10 = 0:8, 20 = 0:1 and 30 = 0:1. The mode
switching probability matrix is given by (as in [21])
2
66
66
66
66
4
p11 p12 p13
p21 p22 p23
p31 p32 p33
3
77
77
77
77
5
=
2
66
66
66
66
4
0:8 0:0 0:2
0:0 0:8 0:2
0:3 0:3 0:4
3
77
77
77
77
5
: (5.61)
The Sensors: Two sensors are used to obtain the measurements. Sensor 1 and Sensor
2 are located at [x1;y1;z1]=[-4000 4000 0] m and [x2;y2;z2]=[5000 0 0] m, respectively,
and the central processor is collocated with sensor 1 platform (we assume that there is no
time delay between sensor 1 and central processor and, on the other hand, there is xed-
and-known time delay between sensor 2 and central processor). The measurements from
sensor l for model j are zlk = hl(xk) +wj;lk for l = 1 and 2, re ecting range and azimuth
angle for sensor 1 (infrared) and azimuth and elevation angles for sensor 2 (radar) . The
range, azimuth, and elevation angle transformations, respectively, are given by
rl = [(x xl)2 + (y yl)2 + (z zl)2]1=2 (5.62)
al = tan 1[(y yl)=(x xl)] (5.63)
el = tan 1
(z zl)=[(x xl)2 + (y yl)2]1=2
: (5.64)
The measurements obtained from sensors 1 and 2 can be expressed as we see from (5.1),
(5.2), (5.4) and (5.11). The measurement noise wj;lk for sensor l is assumed to be zero-
mean white Gaussian with known covariances, R1 = diag[qa1;qe] = diag[4mrad2;4mrad2]
131
with qa1 and qe denoting the variances for the infrared sensor azimuth and elevation mea-
surement noises, and R2 = diag[qr;qa2] = diag[400m2;49mrad2] with qr and qa2 denoting
the variances for the radar range and azimuth measurement noises, respectively. The
sampling interval was T=1s and it was assumed that the probability of detection PD=1
for both sensors. The time di erence between the sampling time at the central processor
(local sensor) and the measurement time at the remote sensor dkl is xed-and-known to
be 0:5T. The data at any time k has a probability Pd = 0:4 (Pd = 0:25 also applied for
comparison) of being delayed where the delay is uniformly distributed with a maximum
time delay less than or equal to T where T is sampling rate.
The Clutter: For generating false measurements in simulations, the clutter was as-
sumed to be Poisson distributed with expected number of 1 = 2 10 4= mrad2 for
sensor 1 (infrared) and 2 = 20 10 6/m-mrad for sensor 2 (radar) [21, case 1]. These
statistics were used for generating the clutter in all simulations. However, a nonpara-
metric clutter model was used for implementing all the algorithms for target tracking.
Other Parameters: The gates for setting up the validation regions for both the sensors
were based on the threshold =16. With the measurement vector of dimension 2, this
leads to a gate probability PG=0.9997 (see [12, pages 95-96]).
Simulation Results: The results were obtained from 100 Monte Carlo runs. Fig. 5.2
shows the true trajectory of the target. Fig. 5.3 (a) and (b) show RMS error comparison
in position among proposed AS-IMM/MSPDAF algorithm dealing with OOSM, standard
AS-IMM/MSPDAF algorithm [50] with OOSM discarding, and AS-IMM/MSPDAF al-
gorithm applied to the hypothetical case of Pd = 0. The rst maneuver starts at 20 sec
and in Fig. 5.3 (b) one can see a signi cant improvement from 25 sec through 36 sec. The
132
0.8 1 1.2 1.4 1.6 1.8 2 2.2x 1040
2000
4000
6000
8000
10000
12000
14000
16000
18000
x coord (m)
y coord (m)
t=0
0 10 20 30 40 50 60 70 80 90?400
?300
?200
?100
0
100
200
300
400
time (sec)
velocity (m/sec)
xy
0 10 20 30 40 50 60 70 80 90?60
?40
?20
0
20
40
60
time (sec)
acc. (m/sec
2 )
xy
0 10 20 30 40 50 60 70 80 900
10
20
30
40
50
60
70
time (sec)
acceleration magnitude (m/sec
2 )
Figure 5.2: Trajectory of maneuvering target (read left to right, top to bottom). (a)
Position in xy plane. (b) x and y velocities. (c) x and y accelerations. (d) magnitude of
accelerations.
second maneuver starts at 55 sec and in Fig. 5.3 (b) one can see a signi cant improve-
ment from 60 sec through 68 sec. That is, the proposed AS-IMM/MSPDAF algorithm
responds faster to maneuvers. Once the target is \settled" in a particular mode, there is
insigni cant di erences between the two approaches. It is seen from Fig. 5.3 (a) and (b)
that the higher delay probability Pd, the more signi cant performance improvement can
be obtained compared with the standard AS-IMM/MSPDAF algorithm [50] that ignores
and discards possibly existing OOSM.
133
0 10 20 30 40 50 60 70 80 900
50
100
150
200
250
300
time (sec)
x?y position RMS error (m)
AS?IMM/MSPDAF w/ OOSM, delay = 0.5T, Pd = 0.25, (100 MCRs)
proposed AS?IMM/MSPDAFstandard AS?IMM/MSPDAF [50]
In?sequence measurements only
0 10 20 30 40 50 60 70 80 900
50
100
150
200
250
300
AS?IMM/MSPDAF w/ OOSM, delay = 0.5T, Pd = 0.4, (100 MCRs)
time (sec)
x?y position RMS error (m)
proposed AS?IMM/MSPDAFstandard AS?IMM/MSPDAF [50]
In?sequence measurements only
Figure 5.3: AS-IMM/MSPDA comparison (RMSE in position, read top to bottom) for
various probabilities of delayed measurement, Pd = 0:25 and 0.4: (a) Pd=0.25. (b)
Pd=0.4. Solid: proposed AS-IMM/MSPDAF algorithm dealing with OOSM; dashed:
AS-IMM/MSPDAF algorithm [50] with OOSM discarding; dotted: AS-IMM/MSPDAF
algorithm applied to the hypothetical case of Pd = 0.
134
5.7 Conclusions
We investigated an augmented state IMM/MSPDAF algorithm with asynchronous
measurement (there is xed-and-known timing mismatch between sensor platforms with
possible OOSM) for tracking a highly maneuvering target in clutter. Simultaneous mea-
surement update technique is applied for better data association and is expected to be
more accurate [12] since it considers all association hypotheses coupled across multisen-
sor while in the sequential updating considers the separate hypothesis for each sensor.
Our proposed approach has extended the multisensor approach of Chapter 4 (or see [50])
to OOSM by employing additional data association. The proposed algorithm was illus-
trated via a simulation example where it outperformed a standard AS-IMM/MSPDAF
algorithm with OOSM discarding [50] especially during the periods following the onset
of the target maneuvers. This improvement in accuracy is seen in our simulation ex-
ample only during the periods following the onset of the target maneuvers. Once the
target is settled in a particular mode, there is insigni cant di erences between the two
approaches. As one can easily notice, the higher delay probability Pd (with more OOSM
appearances), the more signi cant performance improvement can be obtained compared
with the standard AS-IMM/MSPDAF algorithm [50] that ignores and discards possibly
existing OOSM.
135
Chapter 6
Tracking of Multiple Maneuvering Targets in Clutter with Possibly
Merged Measurements using IMM and JPDAM Coupled Filtering
In this chapter, we present a suboptimal ltering algorithm for tracking multiple
highly maneuvering targets in a cluttered environment using multiple sensors. We con-
centrate on two targets which temporarily move in close formation, giving rise to a single
detection due to the resolution limitations of the sensor. The ltering algorithm is de-
veloped by applying the basic IMM approach and the joint probabilistic data association
with merged measurements (JPDAM) technique and coupled target state estimation to
a Markovian switching system. The algorithm is illustrated via a simulation example in-
volving tracking of two highly maneuvering, at times closely spaced, targets with possibly
unresolved measurements. Compared with an existing IMM/JPDA ltering algorithm
developed without allowing for merged measurements, the proposed algorithm achieves
signi cant improvement in the accuracy of track estimation during target merging period.
6.1 Introduction
In this chapter, we consider the problem of tracking multiple maneuvering targets
which temporarily operate in close formation in clutter. This class of problem has re-
ceived considerable attention in the literature [5, 12, 21, 22, 44, 45]. The switching
multiple model approach has been found to be quite e ective in modeling highly ma-
neuvering targets [12, 14, 21, 29, 30, 33, 45]. In this approach various \modes" of target
motion are represented by distinct kinematic models, and in a Bayesian framework, the
136
target maneuvers are modeled by switching among these models controlled by a Markov
chain. In the presence of clutter, the measurements at the sensors may not all have
originated from the target-of-interest. In this case one has to solve the problem of data
association. An e ective approach in a Bayesian framework is that of probabilistic data
association (PDA) [12, 22, 33] for a single target in clutter and that of joint probabilistic
data association (JPDA) [10, 12, 22, 27, 33] for multiple targets in clutter.
Typically it is assumed that the number of targets is known and for each target, a
tracks have been formed (initiated), so that the objective is that of track maintenance.
In [21] such a problem has been considered using multiple sensors, PDA, and switching
multiple models. The optimal solution (in the minimum mean-square error sense) to
target state estimation given sensor measurements and absence of clutter, requires expo-
nentially increasing (with time) computational complexity; therefore, one has to resort
to suboptimal approximations. For the switching multiple model approach, the interact-
ing multiple model (IMM) algorithm of [30] has been found to o er a good compromise
between the computational and storage requirements and estimation accuracy [29]. In
the presence of clutter, one has to account for measurements of uncertain origin (target
or clutter?). Here too, in a Bayesian framework, one has to resort to approximations
to reduce the computational complexity, resulting in the PDA lter [5, 12, 21, 22, 45].
In [21] the IMM algorithm has been combined with a PDA lter in a multiple sensor
scenario to propose a combined IMM/MSPDAF (interacting multiple model/ multisen-
sor probabilistic data association lter) algorithm. In [10, 22] multiple targets in clutter
(but without using switching multiple models) have been considered using JPDA lter
which, unlike the PDA lter, accounts for the interference from other targets. Various
137
versions of IMM/JPDA (interacting multiple model/ joint probabilistic data associa-
tion) lters for multiple target tracking using switching multiple models may be found
in [12, 23, 24, 25, 26]. While [24, 26] present uncoupled lters (i.e., assume that dif-
ferent target states are mutually independent conditioned on the past measurements),
[12, 23, 25] present coupled lters (i.e., assume that there exists \share" measurements,
yielding cross-covariances which re ect the correlation between the targets? state esti-
mation errors). [23] presents an \exact" JPDA coupled lter for non-switching models
using the framework of a linear descriptor system. For switching models, [24] presents
IMM/JPDA uncoupled lter approximations. In [25], an IMM/JPDA coupled lter-
ing algorithm has been presented where a simulation example resulted in fewer target
swapping compared with uncoupled IMM/JPDA.
When two targets are \closely" spaced, they may give rise to a single detection due
to the resolution limitations of the sensor. For instance, in radar ranging, returns from
multiple targets could fall in the same range cell, resulting in one unresolved detection
only [11, 28]. Standard tracking algorithms that ignore such a phenomenon can lead
to poor performance in multiple target tracking [11, 28]. Despite its importance, prior
work on tracking with unresolved measurements in general and modeling of resolution
capability of a sensor in particular, is sparse. Prior work includes [11] and [28] and
references therein. In [11] the resolution phenomena related to tracking have been treated
on the basis of a grid of resolution cells \frozen" ([28]) in space. In [28] the resolution
capability of a sensor is described in terms of a conditional probability of the event
that two targets are unresolved, conditioned on the relative distance between the two
targets in terms of the measured variables (range, azimuth, etc.). A simple Gaussian
shape is assumed which captures the sensor behavior in a mathematically trackable way.
138
While [11] considers JPDA for data association, [28] exploits multiple hypothesis tracking
(MHT). Under the Bayesian framework, there are two basic methods of measurement-
to-track association in multiple target environments [12]: JPDA and multiple hypothesis
tracking (MHT). The MHT lter associates feasible measurements to track and form
multiple hypotheses for track extension. It is measurement-oriented approach whereas
JPDA is a target-oriented approach. Moreover, MHT is a multiscan approach utilizing
several scans of measurements to make data association decisions. MHT makes hard
decisions where highly improbable hypotheses are pruned to reduce the computational
burden. The JPDA lter is a single scan approach which does not make hard decisions;
rather it updates a track with a weighed sum of the measurements which could have
(reasonably) originated from the target in track.
In this chapter we propose to use sensor resolution modeling of [28] in conjunction
with JPDA coupled ltering and interacting multiple model (IMM) approach (see e.g.
[25] for tracking with resolved measurements). As noted in [27], IMMPDA lter is in
general superior to IMM/MHT lter when the associated computational cost and per-
formance are considered. Therefore, our emphasis will be on IMM/JPDA techniques.
Neither [11] nor [28] consider multiple switching kinematic models for maneuvering tar-
gets; rather they are limited to single (nonswitching) kinematic models per target to
achieve much enhanced performance.
This chapter is organized as follows. The problem formulation is presented in Sec.
6.2. The modeling scheme for the merged measurements is focused in Sec. 6.3. Sec.
6.4 describes the proposed IMM/JPDAM algorithm with coupled ltering. Simulation
results using the proposed algorithm for a realistic problem are given in Sec. 6.5. Finally,
Sec. 6.6 presents a discussion of the results and some conclusions.
139
6.2 Problem Formulation
Assume that there are a total of two targets with the target set denoted as T2.
Assume that the target dynamics can be modeled by one of n hypothesized models.
The model set is denoted as Mn :=f1; ;ng and there are total q sensors from which
q 2, or less (if probability of target detection is less than one) or more (due to clutter),
measurement vectors are generated at a time. For target r (r2T2), the event that model
j is in e ect during the sampling period (tk 1;tk] will be denoted by Mjk(r). Although
two targets share a common model set, they may be in a di erent motion status from
time to time.
6.2.1 Target Dynamics
For the j-th hypothesized model (mode), the state dynamics of target r (r2T2),
are modeled as
xk(r) = Fjk 1(r)xk 1(r) +Gjk 1(r)vjk 1(r) (6.1)
where xk(r) is the system state of target r at tk and of dimension nx (assuming all targets
share a common state space), Fjk 1(r) and Gjk 1(r) are the system matrices when model
j is in e ect over the sampling period (tk 1;tk] for target r. The process noise vjk 1(r) is
a zero-mean white Gaussian process with covariance matrix Qjk 1 (same for all targets).
At the initial time t0, the initial conditions for the system state of target r under each
model j are assumed to be Gaussian random variables with the known mean xj0(r) and
the known covariance Pj0 (r). The probability of model j at t0, j0(r) = P[Mj0(r)], is also
assumed to be known. The switching from model Mik 1(r) to model Mjk(r) is governed
140
by a nite-state stationary Markov chain (same for all targets) with known transition
probabilities pij = P[Mjk(r)jMik 1(r)]. Henceforth, time tk will be denoted by k.
In coupled state estimation the states of two targets are estimated jointly [12]. To
this end de ne the \global coupled" state
xk := colfxk(1);xk(2)g (6.2)
and the corresponding matrices/vectors J := colfj1;j2g where jr 2Mn is model j for
target r,
FJk := block diagfFj1k (1);Fj2k (2)g (6.3)
GJk := block diagfGj1k (1);Gj2k (2)g (6.4)
vJk := colfvj1k (1);vj2k (2)g: (6.5)
Then we have the state equation for two targets as
xk = FJk 1xk 1 +GJk 1vJk 1 (6.6)
where
EfvJkvJk0g= QJk := block diagfQj1k ;Qj2kg: (6.7)
De ne the global mode
MJk :=fMj1k (1);Mj2k (2)g: (6.8)
141
The two targets are assumed to evolve independently of each other. Therefore, the
transition probability for the global modes are given by
pIJ := P[Mj1k (1);Mj2k (2)jMi1k 1(1);Mi2k 1(2)] =
2Y
r=1
pirjr: (6.9)
Similarly we have
J0 := P[Mj10 (1);Mj20 (2)] =
2Y
r=1
jr0 (r): (6.10)
6.2.2 Measurements
For the j-th hypothesized model (mode), measurements of target r (r 2T2), are
modeled, when resolved, as
zlk(r) = hl(xk(r)) +wj;lk (r) for l = 1; ;q (6.11)
where zlk(r) is the (true) measurement vector (i.e., due to target r) from sensor l at tk and
of dimension nzl, and hl is the nonlinear transformation of xk(r) to zlk(r) (l = 1; ;q).
The measurement noise wj;lk (r) is a zero-mean white Gaussian process with covariance
matrix Rlk (same for all targets) and is mutually uncorrelated with the process noise
vjk 1(r). Similarly de ne the global measurement vector at sensor l as
zlk := colfzlk(1);zlk(2)g (6.12)
and related vectors
hl(xk) = colfhl(xk(1));hl(xk(2))g; (6.13)
142
wJ;lk := colfwj1;lk (1);wj2;lk (2)g (6.14)
where
EfwJ;lk wJ;lk 0g= Rlk := block diagfRlk;Rlkg: (6.15)
Then the measurement equation for two targets at sensor l (assuming no clutter and
perfect detections) is given by
zlk = hl(xk) +wJ;lk for l = 1; ;q: (6.16)
Regarding the measurements at sensor l, we follow the notations and de nitions used in
[26]. Note that, in general, at any time k, some measurements may be due to clutter
and some due to the target(s). The measurement set (not yet validated) generated by
sensor l at time k is denoted as Zlk :=fzl(1)k ;zl(2)k ; ;zl(ml)k g where ml is the number of
measurements generated by sensor l at time k. Variable zl(i)k (i = 1; ;ml) is the ith
measurement within the set. The validated set of measurements of sensor l at time k
will be denoted by Ylk := fyl(1)k ;yl(2)k ; ;yl( ml)k g where ml is total number of validated
measurement for sensor l at time k. And yl(i)k := zl(li)k where 1 l1 0 is a \di use" prior (for nonparametric modeling of clutter) whose exact value is
irrelevant. De ne
PJ;1;ak =
Z
P(Ajxk)N
xk; ^xJkjk 1;PJkjk 1
dxk (6.62)
=
2 R1;d
1=2N
0
BB
@0;
h
Hj1;1k Hj2;1k
i
^xJkjk 1;
Hj1;1k Hj2;1k
PJ;1kjk 1
2
66
4
Hj1;10k
Hj2;10k
3
77
5+R1;d
1
CC
A:
The last line of (6.62) can be obtained by substituting (6.34) into the rst line of (6.62).
For details, see [28] and Appendix B. Then we have
P[ jMJk;Zk 1] =
8
>>><
>>>:
D( )(1 PJ;1;ak ) for resolved target(s)
D( )PJ;1;ak for merged targets:
(6.63)
Unlike [26], we do not assume that the states of the targets (including the modes)
conditioned on the past observations are mutually independent. Then we have
p
Y 1kj ;MJk;Zk 1
= V ( )1 p
~
Y 1k ( )jMJk;Zk 1
(6.64)
156
where ~Y 1k ( ) Y 1k is a subset of the validated measurements Y 1k , consisting of the mea-
surements associated with the targets as speci ed by . The number of measurements
in ~Y 1k ( ) is equal m1 ( ) where ( ) is the number of false alarms.
De ne a m1 [ m1 ( )] matrix ^ ( ) as a submatrix of ^ ( ) obtained by deleting
the rst column and all null columns of ^ ( ). Then for a given , we have a measurement
vector ~Y 1k ( ) of dimension ( m1i=1min[1; i( )])nz1 given by
~Y 1k ( ) = (In
z1
^ 0( ))colfy1(i)k ;i = 1;2; ; m1g (6.65)
where we stack up all target-associate validated measurements in in ascending order of
targets, In is the n n identity matrix, and the symbol denotes the Kronecker product.
De ne a [( m1 ( ))zz1] [2nx] matrix HJ;1k ( ) as a submatrix of HJ;1k obtained by
deleting all i-th block rows of HJ;1k for which i( ) = 0. That is, we have modi ed HJ;1k
to keep only the block elements associated with target-associated measurements in .
To further simplify the equation for ~Y 1k ( ), one has to consider all the possible joint
association events . De ne a related set of mutually exclusive and exhaustive data
interpretations as follows (here we follow [28])
11: Both targets were resolved and detected ( ( ) = m1 2),
10: Both targets were resolved and only target 1 was detected ( ( ) = m1
1; 1( ) = 1; 2( ) = 0),
01: Both targets were resolved and only target 2 was detected ( ( ) = m1
1; 1( ) = 0; 2( ) = 1),
157
1: Both targets were detected but merged as a single measurement ( ( ) =
m1 1; 1( ) = 1; 2( ) = 1),
0: No target was detected ( ( ) = m1).
It then follows that the linearized measurement equation for ~Y 1k ( ) is given by
~Y 1k ( ) =
8
>>>>
>>>>
>>>
>>>>
>>>>
<
>>>>
>>>>
>>>
>>>>
>>>>
:
HJ;1k ( )xk +wJ;1k ; for 2 11
Hj1;1k ( )xk +wj1;1k (1); for 2 10
Hj2;1k ( )xk +wj2;1k (2); for 2 01
2
66
4
HJ;1;ak ( )
HJ;1;dk ( )
3
77
5xk +
2
66
4
wJ;1;ak
wJ;1;dk
3
77
5; for 2 1:
(6.66)
Conditioned on the joint association event and mode J, the \coupled" innovation
is given by
J;1k ( ) =
8
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
<
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>:
~Y 1k ( ) ^zJ;1k ( ); for 2 11
~Y 1k ( ) ^zj1;1k ( ); for 2 10
~Y 1k ( ) ^zj2;1k ( ); for 2 01
~Y 1k ( )
2
66
4
^zJ;1;ak ( )
^zJ;1;dk ( )
3
77
5; for 2 1;
0; for 2 0
(6.67)
where ^zJ;1k ( ) (^zJ;1;ak ( ), ^zJ;1;dk ( )) are subvector(s) of ^zJ;1k (^zJ;1;ak , ^zJ;1;dk ) obtained by
deleting all i-th block rows (nz1 1) of ^zJ;1k (^zJ;1;ak , ^zJ;1;dk ) for which i( ) = 0. The
covariance of mode-conditioned residual conditioned on the joint association event is
158
given by
SJ;1k ( ) = HJ;1k ( )PJkjk 1HJ;1k 0( ) +RJ;1k ; for 2 11
Sj1;1k ( ) = Hj1;1k ( )Pj1kjk 1(1)Hj1;1k 0( ) +Rj1;1k ; for 2 10
Sj2;1k ( ) = Hj2;1k ( )Pj2kjk 1(2)Hj2;1k 0( ) +Rj2;1k ; for 2 01
SJ;1;ak ( ) =
2
66
4
HJ;1;ak ( )
HJ;1;dk ( )
3
77
5PJkjk 1
HJ;1;a0k ( ) HJ;1;d0k ( )
+
2
66
4
RJ;1;ak 0
0 R1;d
3
77
5; for 2 1;
(6.68)
where Pj1kjk 1(1) and Pj2kjk 1(2) are diagonal submatrices of PJkjk 1.
There are a total of ( m1 + 1) ( m1 + 1) possible association hypotheses, each of
which has an association probability. Then we have
p
Y 1kj ;MJk;Zk 1
8>
>>>>
>>>>
>>>
>>>>
>>><
>>>>
>>>
>>>>
>>>>
>>>>
:
V 2 m11 p( ~Y 1k ( )jMJk;Zk 1); for 2 11
V 1 m11 p( ~Y 1k ( )jMJk;Zk 1); for 2 10
V 1 m11 p( ~Y 1k ( )jMJk;Zk 1); for 2 01
V 1 m11 p( ~Y 1k ( )jMJk;Zk 1); for 2 1
V m11 ; for 2 0
(6.69)
159
where the conditional pdf (probability density function) of the validated measurements
~Y 1k ( ) given their origins (speci ed by ) and the global mode J, is given by
p
~
Y 1k ( )jMJk;Zk 1
=
8
>>>>
>>>>
>>>
>>>>
>>>>
<
>>>>
>>>>
>>>
>>>>
>>>>
:
N
~
Y 1k ( ); ^zJ;1k ( );SJ;1k ( )
; for 2 11
N
~
Y 1k ( ); ^zj1;1k ( );Sj1;1k ( )
; for 2 10
N
~
Y 1k ( ); ^zj2;1k ( );Sj2;1k ( )
; for 2 01
N
0
BB
@~Y 1k ( );
2
66
4
^zJ;1;ak ( )
^zJ;1;dk ( )
3
77
5;SJ;1;ak ( )
1
CC
A; for 2 1
(6.70)
The probability of the joint association event given that global mode J is e ective
from time k-1 through k is
J;1k ( ) := P[ jMJk;Zk 1;Y 1k ]
= 1cp(Y 1kj ;MJk;Zk 1)P[ jMJk;Zk 1] (6.71)
where the rst term can be calculated from (6.64)-(6.68), the second term from (6.61)-
(6.63), and c is a normalization constant such that
P[ jMJk;Zk 1;Y 1k ] = 1: (6.72)
Using ^xJkjk 1 (from 6.34) and its covariance PJkjk 1 (from 6.35), one computes the partial
update ^xJkjk and its covariance PJkjk following the standard PDAF [12], except that the
global state is conditioned on , not the marginal events ir; details follow.
160
Kalman gain:
WJ;1k ( ) =
8>
>>>
>>>>
>>>>
>>>>
>>>
<
>>>
>>>>
>>>>
>>>
>>>>
>:
PJkjk 1HJ;1k ( )0
h
SJ;1k ( )
i 1
; for 2 11
Pj1kjk 1Hj1;1k ( )0
h
Sj1;1k ( )
i 1
; for 2 10
Pj2kjk 1Hj2;1k ( )0
h
Sj2;1k ( )
i 1
; for 2 01
PJkjk 1
2
66
4
HJ;1;ak ( )
HJ;1;dk ( )
3
77
5
0
h
SJ;1;ak ( )
i 1
; for 2 1
(6.73)
Partial update of the state estimate:
^xJ;1kjk( ) := E
n
xkj ;MJk;Zk 1;Y 1k
o
=
8
>>>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
<
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>:
^xJkjk 1 +WJ;1k ( ) J;1k ( ); for 2 11
^xJkjk 1 +
2
66
4
WJ;1k ( ) J;1k ( )
0nx 1
3
77
5; for 2 10
^xJkjk 1 +
2
66
4
0nx 1
WJ;1k ( ) J;1k ( )
3
77
5; for 2 01
^xJkjk 1 +WJ;1k ( ) J;1k ( ); for 2 1
^xJkjk 1; for 2 0:
(6.74)
Mode-conditioned update of the state estimate:
^xJ;1kjk := E
n
xkjMJk;Zk 1;Y 1k
o
= X
J;1k ( )^xJ;1kjk( ) (6.75)
161
Covariance of ^xJ;1kjk( ):
PJ;1kjk( ) =
8>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
<
>>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>:
PJkjk 1 WJ;1k ( )SJ;1k ( )WJ;1k ( )0; for 2 11
PJkjk 1
2
66
4
WJ;1k ( )Sj1;1k ( )WJ;1k ( )0 0nx nx
0nx nx 0nx nx
3
77
5;for 2 10
PJkjk 1
2
66
4
0nx nx 0nx nx
0nx nx WJ;1k ( )Sj2;1k ( )WJ;1k ( )0
3
77
5; for 2 01
PJkjk 1 WJ;1k ( )SJ;1;ak ( )WJ;1k ( )0; for 2 1
PJkjk 1; for 2 0:
(6.76)
Covariance of ^xJkjk: PJ;1kjk := E
n
(xk ^xJ;1kjk)(xk ^xJ;1kjk)0jY 1k;Zk 1;MJk
o
= PJkjk 1 P
2 10
J;1k ( )WJ1k ( )SJ1;1k ( )WJ1k ( )0 P
2 01
J;1k ( )WJ2k ( )SJ2;1k ( )WJ2k ( )0
P
2 11; 1
J;1k ( )WJ;1k ( )SJ;1k ( )WJ;1k ( )0 + P
2 11; 1
J;1k ( )WJ;1k ( ) J;1k ( ) J;1k ( )0WJ;1k ( )0
+ P
2 10
J;1k ( )WJ1k ( ) J1;1k ( ) J1;1k ( )0WJ1k ( )0
+ P
2 01
J;1k ( )WJ2k ( ) J2;1k ( ) J2;1k ( )0WJ2k ( )0
"
P
2 10
J;1k ( )WJ1k ( ) J1;1k ( )0
#"
P
2 10
J;1k ( )WJ1k ( ) J1;1k ( )0
#0
"
P
2 01
J;1k ( )WJ2k ( ) J2;1k ( )0
#"
P
2 01
J;1k ( )WJ2k ( ) J2;1k ( )0
#0
"
P
2 11; 1
J;1k ( )WJ;1k ( ) J;1k ( )0
#"
P
2 11; 1
J;1k ( )WJ;1k ( ) J;1k ( )0
#0
(6.77)
162
where
WJ1k ( ) =
2
66
4
WJ;1k ( ) 0nx nzl
0nx nzl 0nx nzl
3
77
5; WJ2k ( ) =
2
66
4
0nx nzl 0nx nzl
0nx nzl WJ;1k ( )
3
77
5;
SJ1k ( ) =
2
66
4
Sj1k ( ) 0nzl nzl
0nzl nzl 0nzl nzl
3
77
5; SJ2k ( ) =
2
66
4
0nzl nzl 0nzl nzl
0nzl nzlSj2;1k ( )
3
77
5;
J1k ( ) =
2
66
4
J;1k ( )
0nzl 1
3
77
5; J2k ( ) =
2
66
4
0nzl 1
J;1k ( )
3
77
5:
(6.78)
Eqn. (6.77) follows in a manner similar to eqn. (3.4.2-10) in [12]. For details, see [12].
Step 5. The mode-conditioned predicted measurements for sensor 2 (8J 2
Mn):
For two resolved targets: The \predicted" measurement for sensor 2 is given by
^zJ;2k := h2(^xJ;1kjk): (6.79)
The covariance of the global mode-conditioned residual J;2(I)k := z2(I)k ^zJ;2k , where
z2(I)k := colfz2(i1)k ;z2(i2)k g, is given by
SJ;2k := Ef J;2(I)k J;2(I)0k jMJk;Zk 1;Y 1kg:= HJ;2k PJ;1kjkHJ;20k +R2k: (6.80)
where HJ;2k := block diagfHj1;2k ;Hj2;2k g is the rst order derivative (Jacobian matrix)
of h2(:) evaluated at the state prediction ^xJ;1kjk.
163
For unresolved targets: The \predicted" measurement for sensor 2 is given by
^zJ;2;ak := h2;a(^xJ;1kjk): (6.81)
De ne
^zJ;2;dk := h2;d(^xJ;1kjk): (6.82)
Using the linearized (6.19) and (6.24), the covariance of the mode-conditioned residual
J;2;ak :=
2
66
4
zJ;2;ak
0
3
77
5
2
66
4
^zJ;2;ak
^zJ;2;dk
3
77
5 (6.83)
is given by
SJ;2;ak := E
n
J;2;ak J;2;a0k jMJk;Zk 1;Y 1k
o
=
2
66
4
HJ;2;ak
HJ;2;dk
3
77
5PJ;1kjk
HJ;2;a0k HJ;2;d0k
+
2
66
4
R2;ak 0
0 R2;d
3
77
5 (6.84)
where HJ;2;ak = [ kHj1;2 (1 k)Hj2;2] (HJ;2;dk = [Hj1;2 Hj2;2]) is the rst order
derivative (Jacobian matrix) of h2;a(:) (h2;d(:)) evaluated at the state prediction ^xJ;1kjk.
Step 6. Measurement validation measurements for sensor 2 (8J2 Mn):
This is similar to Step 3 where we replace SJ;1k with SJ;2k , z1(i)k with z2(i)k , m1 with
m2, V 1k (r) with V 2k (r), and V 1k with V 2k . Details are similar to that in Step 3, hence
omitted.
Step 7. Update with validated measurements for sensor 2 (8J2 Mn):
164
This is similar to Step 4. Using the validated measurements obtained from Step 6
and starting from ^xJ;1kjk and PJ;1kjk, one computes the nal updates ^xJkjk and PJkjk, and the
likelihood
J;2k := p
h
Y 2k;jMJk;Zk 1;Y 1k
i
= X
p
h
Y 2kj ;MJk;Zk 1;Y 1k
i
Pf jMJk;Zk 1;Y 1kg: (6.85)
Details are similar to that in Step 4, hence omitted.
Step 8. Update of mode probabilities (8j2Mn, 8r2T2):
Jk := P[MJkjZk] = 1c J k J;1k J;2k (6.86)
where c is a normalization constant such that J Jk = 1. For individual targets we have
j1k (1) := P[Mj1k (1)jZk] =
nX
j2=1
j1;j2k ; j2k (2) =
nX
j1=1
j1;j2k (6.87)
with J = (j1;j2) in (6.86).
Step 9. Combination of the mode-conditioned estimates (8r2T2):
The nal global state estimate update at time k is given by
^xkjk = XJ ^xJkjk Jk (6.88)
and its covariance is given by
Pkjk = XJ
PJkjk +
h
^xJkjk ^xkjk
ih
^xJkjk ^xkjk
i0
Jk: (6.89)
165
The state estimate ^xkjk(r) for target r is the nx-subvector of ^xkjk consisting of elements
(r 1)nx +m, m = 1;2.
Remark 1. In the above algorithm we used sequential updating of the state es-
timates with measurements (one sensor at a time - see Steps 4 and 7) as in [21] and
[26]. This approach is suboptimal but leads to computational savings as one does not
have to simultaneously associate measurements across sensors (as in [46], [48], [53]).
In Step 4 we are interested in EfxkjMJk;Zk 1k ;Y 1kg which is decomposed as in (6.79)
conditioned on ?s; Measurements Y 2k are not considered in this step. If one were to
seek EfxkjMJk;Zk 1k ;Y 1k;Y 2kg, then we would have to follow the approach of [46]-[53] by
picking all possible association events across sensors also.
Remark 2. Compared to the uncoupled ltering of [26] where the equations are
formulated conditioned on marginal association events ir, here we have conditioning on
joint association events for couples ltering. Eqn. (6.70) does not decompose into the
product of marginal probabilities as in [26].
Remark 3. Partition the set of all s into disjoint sets is such that
i := f j r( ) = r( ~ ) 8r; ~ 2 ig where i = 1;2; ;K. For instance, for N=2,
we have K=5 with 1=all s in which there are two validated measurements associ-
ated with two targets, 2=all s in which one validated measurements associated with
target 1 and none with target 2, 3=all s in which one validated measurements associ-
ated with target 2 and none with target 1, 4=all s in which one validated unresolved
measurement associated with two targets, 5=all s in which none of the validated mea-
surements are associated with any target. It is then easily seen that WJ;1k ( ), HJ;1k ( ),
SJ;1k ( ) and J;1k ( ) in Step 4, all are invariant for 2 i. This fact can be used to
simplify computations in (6.77). Similar comments apply to Step 7.
166
Remark 4. If one substitutes (6.65) into (6.66), then one obtains a linear descriptor
system type of equation such as (12) in [23]. Therefore, the standard state-space system
framework used in this dissertation and the linear descriptor system framework used in
[23] are equivalent (except that [23] uses non-switching models whereas we use Markovian
switching models).
6.5 Simulation Example
The following example of tracking two highly maneuvering targets in clutter is pat-
terned after [25].
The True Trajectory: We consider the same scenario as that in [25] except for a lin-
ear shift in the x and y-directions in the trajectory of target 1 to achieve merged target
with a long-term duration. Target 1 starts at location [21689+dx 10840+dy 40] with
dx=-3040m and dy=5500m in Cartesian coordinates in meters. The initial velocity (in
m/s) is [-8.3 -399.9 0] and the target stays at constant altitude with a constant speed of
400 m/s. Its trajectory is:
a straight line with constant velocity between 0 and 27s,
a coordinated turn (0.15 rad/s) with constant acceleration of 60 m/s2 between 27
and 42s,
a straight line with constant velocity between 42 and 47s,
a coordinated turn (0.1 rad/s) with constant acceleration of 40 m/s2 between 47
and 65s,
a straight line with constant velocity between 65 and 87s.
167
Target 2 starts at location [30000 -3040 40] in Cartesian coordinates in meters. The
initial velocity (in m/s) is [-382 157 0] and the target stays at constant altitude with a
constant speed of 413 m/s. Its trajectory is:
a straight line with constant velocity between 0 and 44s,
a coordinated turn (0.075 rad/s) with constant acceleration of 30 m/s2 between 44
and 59s,
a straight line with constant velocity between 59 and 87s.
The Target Motion Models: These are exactly as in [25]. In each mode the target
dynamics are modelled in Cartesian coordinates as xk(r) = F(r)xk 1(r) + G(r)vk 1(r)
where the state of the target is position, velocity, and acceleration in each of the 3
Cartesian coordinates (x;y; and z). Model 1 for Nearly constant velocity model with
zero mean perturbation in acceleration; Model 2 for Wiener process acceleration (nearly
constant acceleration motion); Model 3 for Wiener process acceleration (model with large
acceleration increments, for the onset and termination of maneuvers). The details re-
garding these models have been described in previous chapters and identical. The initial
model probabilities for two targets are identical: 10 = 0:8, 20 = 0:1 and 30 = 0:1. The
mode switching probability matrix for two targets is also identical and is as in previous
chapters.
The Sensors: Two sensors (we assume collocation, and time synchronization of obser-
vations, etc.) are used to obtain the measurements. The measurements from sensor l for
model j are zlk = hl(xk) + wlk, l = 1;2, re ecting range and azimuth angle for sensor 1
(radar) and azimuth and elevation angles for sensor 2 (infrared). The range, azimuth, and
elevation angle transformations would be given by r = (x2 +y2 +z2)1=2, a = tan 1[y=x],
168
e = tan 1[z=(x2 +y2)1=2], respectively. The measurement noise wlk for sensor l is as-
sumed to be zero-mean white Gaussian with known covariances R1 = diag[qr;qa1] =
diag[400m2;49mrad2] with qr and qa1 denoting the variances for the radar range and az-
imuth measurement noises, respectively, and R2 = diag[qa2;qe] = diag[4mrad2;4mrad2]
with qa2 and qe denoting the variances for the infrared sensor azimuth and elevation
measurement noises, respectively. Resolutions of both sensors are selected after from
[28] (twice of the standard deviations for the corresponding sensor measurement noise):
a range resolution of sensor 1 R = 2 pqr = 40m, a angular resolution of sensor 1 1 =
2 pqa1 = 14 10 3rad, a angular resolution of sensor 2 2 = 2 pqa2 = 4 10 3rad
and a elevation angle resolution of sensor 2 = 2 pqe = 4 10 3rad. The noise
for merged measurements wJ;l;ak for sensor l is assumed to be the same with resolved
measurement noise wlk for sensor l. The measurement distance noise wl;dk for sensor l is
assumed to be zero-mean white Gaussian with known covariances R1;d = diag[ R; 1],
and R2;d = diag[ 2; ]. To generate the true target trajectories, following (6.28),
given the measurements of the two targets at a given sensor, the targets are unresolved
with the conditional probability Pl;ak (6.29) and they are merged into one by the linear
combination model (6.19) with the signal strength ratio k = 0:5 for all time k. (The
tracking algorithm does not have this knowledge of how Pl;ak is used to generate data.)
Both sensors are assumed to be located at the coordinate system origin. The sampling
interval was T=1s and it was assumed that the probability of detection PD=0.997 for
both sensors.
The Clutter: For generating false measurements in simulations, the clutter was as-
sumed to be Poisson distributed with expected number of 1 = 20 10 6=m-mrad for
sensor 1 and 2 = 2 10 4=mrad2 for sensor 2. These statistics were used for generating
169
0 0.5 1 1.5 2 2.5 3x 104?5000
0
5000
10000
15000
20000
x coord (m)
y coord (m)
target 1target 2
t=0
t=0
0 10 20 30 40 50 60 70 80 90?500
?400
?300
?200
?100
0
100
200
300
400
500
time (sec)
velocity (m/sec)
target 1?xtarget 1?ytarget 2?x
target 2?y
0 10 20 30 40 50 60 70 80 90
?60
?40
?20
0
20
40
60
time (sec)
acceleration (m/sec
2 )
target 1?xtarget 1?y
target 2?xtarget 2?y
0 10 20 30 40 50 60 70 80 900
0.5
1
1.5
2
2.5x 104
time (sec)
distance between targets (m)
Figure 6.1: The true trajectories of the maneuvering targets (read left to right, top to
bottom): (a) Position in xy plane, (b) x and y velocities, (c) x and y accelerations, (d)
distance between the targets.
the clutter in all simulations. However, a nonparametric clutter model was used for
implementing all the algorithms for target tracking.
Other Parameters: The gates for setting up the validation regions for both the sensors
were based on the threshold =16 corresponding to a gate probability PG=0.9997.
Simulation Results: The results were obtained from 1000 Monte Carlo runs. Fig.
6.1(a)-(d) shows the true trajectory of the two targets and the distance between the two
170
0 10 20 30 40 50 60 70 800
50
100
150
200
250
300
time (sec)
x?y position rms error (m)
JPDAM SM num of run/ err/ swap = 953/44/3 (1000 MCRs)
target 1
target 2
0 10 20 30 40 50 60 70 800
50
100
150
200
250
300
time (sec)
x?y position rms error (m)
JPDA SM num of run/ err/ swap = 902/94/4 (1000 MCRs)
target 1
target 2
Figure 6.2: Performance (RMSE in position) of the proposed IMM/JPDAMCF and
the IMM/JPDACF of [25] based on successful runs (read left to right): (a) proposed
IMM/JPDAMCF, (b) standard IMM/JPDACF [25].
171
index proposed IMM/JPDAMCF standard IMM/JPDACF [25]
No. of lost tracks 44/1000 94/1000
No. of swapped tracks 3/1000 4/1000
No. of successful tracks 953/1000 902/1000
Table 6.1: Simulation results summery based on 1000 runs
targets as a function of time. The two targets start out far apart, move close to each
other from 38 to 42 seconds, and then move apart again. Fig. 6.2(a) shows the results
of the proposed IMM/JPDAMCF based on 953 successful runs (target swap occurred
in 3 runs with 44 track failures). Fig. 6.2(b) shows the standard IMM/JPDACF [25]
(the merged target case was not considered) based on 902 successful runs (target swap
occurred in 4 runs with 94 track failures) in terms of the RMSE in position. Table 1
shows the number of successful runs (excluding target swapping) for the two approaches
IMM/JPDAMCF and IMM/JPDACF. It is seen from Table 1 and Figs. 6.1 and 6.2
that the proposed IMM/JPDAMCF is better than the standard IMM/JPDACF [25] in
performance especially in terms of the track estimation accuracy and the loss of tracks.
To assess the computational requirements of the two approaches, we computed the
CPU time needed to execute 87 time steps in each run (averaged 100 Monte Carlo runs ex-
cluding data/clutter generation) in MATLAB 6.5 on a 2.8 GHz (Mobile) Pentium 4 oper-
ating under Windows XP (professional). The standard IMM/JPDACF needs 7.6178 secs
(for all 87 time steps) compared to 14.3848 secs required by proposed IMM/JPDAMCF.
Thus with a 88.8% increase in computational cost, it is seen that the proposed algorithm
can signi cantly improve the accuracy of track estimation during the periods when tar-
gets are temporarily move in close formation.
172
6.6 Conclusions
We investigated on two targets which temporarily move in close formation, giving
rise to a single detection due to the resolution limitations of the sensor. We proposed
a noble IMM/JPDAMCF algorithm for state estimation for two highly maneuvering
targets in clutter. While a past IMM/JPDACF approach [25] applies only to resolved
measurements, our proposed approach has extended the multisensor approach of [25]
to deal with possibly existing merged measurements arrived from sensors to the central
processor.
The proposed algorithm was illustrated via a simulation example and it outper-
formed a standard IMM/JPDACF algorithm [25] which does not consider the possibility
of the data association to the merged measurements, in terms of RMS position error, the
number of track losses, and target swappings. The proposed IMM/JPDAMCF resulted
in a 25% fewer target swappings and a 54% track losses compared with the standard
MM/JPDACF [25]. The improvement in accuracy and track loss comes at the expense
of an increase (88.8%) in the computational cost. This improvement in accuracy is seen
in our simulation example only during the periods when targets are temporarily move
in close formation and gives rise to a single detection due to the resolution limitations
of the sensor.
173
Chapter 7
Conclusions and Future Work
7.1 Conclusions
In this dissertation, we developed several noble algorithms for maneuvering tar-
get tracking in the presence of clutter using multiple sensors. Our main concern was
to investigate solutions to the target tracking problems which require the simultaneous
completion of two tasks: estimation and data association. There are many di erent
approaches to both estimation and data association and these are generally the distin-
guishing features that give rise to di erent tracking algorithms.
In Chapter 3, a suboptimal ltering algorithm for tracking a highly maneuvering
target in a cluttered environment using multiple sensors was proposed. We developed
a noble algorithm \simultaneous measurement update using IMMPDA ltering". The
existing IMM ltering algorithm and PDA technique were combined with simultane-
ous measurement update algorithm in our proposed algorithm. This algorithm provides
signi cant improvement over an existing IMMPDA ltering algorithm with sequential
sensor processing, because it considers all association hypotheses coupled across multi-
sensor while in the sequential updating considers the separate hypothesis for each sensor.
This improvement in accuracy was remarkable especially during the periods following
the onset of the target maneuvers. The improvement in accuracy comes at the expense
of a slight increase in the computational cost.
174
In a multisensor central tracking system, asynchronous (delayed) measurements arise
due to communication network delays, varying preprocessing times at the sensor plat-
forms and possibly lack of sampling time synchronization among sensor platforms. A
state augmented approach was combined with simultaneous measurement update tech-
nique to deal with asynchronous (delayed) measurements problem in Chapter 4. In
this chapter, we investigated \ xed-but-unknown" delay between sensor network and
developed a noble algorithm \simultaneous measurement update using AS-IMMPDA
ltering". Compared with an existing IMMPDA ltering algorithm with the assump-
tion of synchronous (no delay) measurements sensor processing, the proposed algorithm
achieves considerable improvement (especially in the case of larger delays) in the accu-
racy of track estimation. Using the proposed AS-IMMPDA algorithm, the smoothed
estimate (lag = 1) can be easily obtained from the augmented state estimate and it
shows better performance (in terms of RMS error) than the ltered estimate.
In practice, in a multisensor central target tracking system, one can possibly have
situations where measurements from the same target do not arrive in correct time order.
This OOSM can occur due to communication network delays, varying preprocessing
times at the sensor platforms, and lack of sampling time synchronization among sensor
platforms. In Chapter 5, the AS-IMMPDA algorithm proposed in Chapter 4 is then
extended to deal with \ xed-and-known" delay and possible OOSM between sensor
network. Compared with an existing IMMPDA ltering algorithm without dealing with
possibly existing OOSM (which assumes only \ xed-and-known" delay), the proposed
algorithm achieves considerable improvement in the accuracy of track estimation as the
probability of OOSM arrival increases.
175
In Chapter 6, we considered the problem of tracking multiple maneuvering targets
which temporarily operate in close formation in clutter. When multiple targets are
close, there might be a single detection originating from more than one target, yielding
unresolved (merged) measurements due to a sensor?s inherent resolution limitation. The
existing IMM/JPDA can not handle data association on the merged measurements.
Combining the existing IMM algorithm with JPDAM technique, we presented a noble
IMM/JPDAM algorithm for multisensor tracking of multiple maneuvering targets in
clutter with possibly merged measurements. The algorithm is illustrated via a simulation
example compared with an existing IMM/JPDA ltering algorithm which does not deal
with the merged measurements data association hypotheses. The proposed algorithm
achieved signi cant improvement in the accuracy of track estimation especially during
target merging period, and caused fewer track losses.
7.2 Suggested Future Work
7.2.1 Unresolved Measurements with Simultaneous Measurement Update
In Chapters 3, it has been shown that the IMMPDA ltering algorithm with simul-
taneous measurement update achieves signi cant improvement in the accuracy of track
estimation compared with an existing IMMPDA ltering algorithm with sequential sen-
sor processing. As we have presented proposed IMM/JPDAMCF algorithm to deal with
unresolved measurements by sequential updating approach in Chapter 5, we can ap-
ply the extension of our IMM/JPDAMCF algorithm by combining with simultaneous
measurement update algorithm. This integration of existing IMM/JPDAM algorithms
176
and simultaneous measurement update algorithms will lead to better performance than
earlier algorithms designed for tracking multiple maneuvering targets.
7.2.2 Track Initialization
Until now we have limited our investigation to track maintenance in this disserta-
tion (i.e., the initial state of each target is assumed to have been already obtained). In
practice, the track formation needs to be done before all the algorithms in this disserta-
tion can be applied. Measurement-to-measurement association is the core requirement
for track formation in the presence of measurement uncertainty. Some existing target
track formation approaches can be found in [12, 33]. MHT has been recognized as the
theoretically best approach to the multitarget tracking problem, yet it requires a con-
siderable amount of computation and memory. Therefore, in the past it was considered
una ordable because of the high computational requirements. Now interest is renewed
because of the continuous advances in computer technology and algorithm design. Since
MHT approach, unlike IMM/JPDA, is a measurement-oriented approach and considers
the track initialization as well, integration of existing track formation algorithms and
our track maintenance algorithms is a good challenge for future research topics.
177
Bibliography
[1] N. Wax,\Signal-to-noise improvement and the statistics of track populations," J.
of Applied Physics, Vol. 26, pp. 586-595, May 1955.
[2] R.E. Kalman,\A new approach to linear ltering and prediction problems," Journal
of Basic Engineering, pp. 35-46, March 1960.
[3] R.W. Sittler,\An optimal data association problem in surveillance theory," IEEE
Trans. Military Electronics, Vol. MIL-8, pp. 125-139, April 1964.
[4] A.J. Ja er and Y. Bar-Shalom, \On optimal tracking in multiple target environ-
ments," Proceedings of the Third Symposium on Non-Linear Estimation Theory
and Its Applications,San Diego, CA, pp. 112-117, Sept. 11-13, 1972.
[5] Y. Bar-Shalom and E. Tse, \Tracking in a cluttered environment with probabilistic
data association," Automatica, Vol. 11, pp. 451-460, Sept. 1975.
[6] R.A. Singer and J.J. Stein, \An optimal tracking lter for processing sensor data
of imprecisely determined origin in surveillance system,"Proceedings of the 1971
IEEE Conference on Decision and Control, Miami Beach, FL, pp. 171-175, Dec.
1971.
[7] R.A. Singer, R.G. Sea, and K.B. Housewright, \Derivation and evaluation of im-
proved tracking lters for use in multi-target environments,"IEEE Transactions on
Information Theory, IT-20, pp. 423-432, July 1974.
[8] D.B. Reid, \A multiple hyphothesis lter for tracking multiple targes in a clut-
tered environment," Lockheed Missiles and Space Company Report, No. LMSC,
D-560254, Spet. 1977.
[9] D.B. Reid, \An algorithm for tracking multiple targets," IEEE Transaction on
Automatic Control, Vol. AC-24, pp. 843-854, Dec. 1979.
[10] T.E. Fortmann, Y. Bar-Shalom, and M. Sche e \Sonar tracking of multiple targets
using joint probabilistic data association", IEEE Journal of Oceanic Engineering,
Vol. OE-8, pp. 173-184, July 1983.
[11] Kuo-Chu Chang and Y. Bar-Shalom, \Joint probabilistic data association for mul-
titarget tracking with possibly unresolved measurements and maneuvers", IEEE
Transaction on Automatic Control, Vol. 29, No.7, pp. 585-594, July 1984.
178
[12] Y. Bar-Shalom and X.R. Li, Multitarget-Multisensor Tracking: Principles and
Techniques, Storrs, CT: YBS Publishing, 1995.
[13] Y. Bar-Shalom, X.R. Li, and T. Kirubarajan, Estimation with Applications to
Tracking and Navigation, New York: John Wiley & Sons, 2001.
[14] G.A. Ackerson and K.S. Fu, \On state estimation in swithching environments",
IEEE Trans. on Automatic Control, Vol. AC-15, pp. 10-17, Feb. 1970.
[15] J.K. Tugnait and A.H. Haddad, \A detection-estimation scheme for state estima-
tion in switching environment", Automatica, Vol. 15, pp 477-481, 1979.
[16] H. Akashi and H. Kumamoto, \Random sampling approach to state estimation in
switching environments", Automatica, Vol. 13, pp 429, 1977.
[17] J.K. Tugnait, \Detection and estimation for abruptly changing systems", Auto-
matica, Vol. 18, pp. 607-615, Sept. 1982.
[18] H.A.P. Blom, \An e cient lter for abruptly changing systems," Proceedings of
the 23rd IEEE Conf. on Decision and Control, pp. 656-658, Las Vegas, NV, 1984.
[19] H.A.P. Blom and Y. Bar-Shalom, \The interacting multiple model algorithm for
systems with Markovian switching coe cients," IEEE Transaction on Automatic
Control, Vol. 33, No 8, pp. 780-783, Aug. 1988.
[20] Y. Bar-Shalom, K.C. Chang, and H.A.P. Blom, \Tracking a maneuvering target
using input estimation versus the interacting multiple model algorithm," IEEE
Transaction on Aerospace and Electronic Systems, Vol. AES-25(2), pp. 296-300,
March 1989.
[21] A. Houles and Y. Bar-Shalom, \Multisensor tracking of a maneuvering target in
clutter", IEEE Trans. Aerospace and Electronic Systems, Vol. AES-25, pp. 176-188,
March 1989.
[22] Y. Bar-Shalom and T. E. Fortmann, Tracking and Data Association, New York:
Academic Press, 1988.
[23] H.A.P. Blom and E.A. Bloem, \Probabilistic data association avoiding track coa-
lescence," IEEE Trans. Automatic Control, Vol. AC-45, pp.247-259, Feb. 2000.
[24] H.A.P. Blom and E.A. Bloem, \Combining IMM and JPDA for tracking multiple
manuvering targets in clutter," Proc. 5th Intern. Conf. on Information Fusion,
Annapolis, MD, Vol. 1, pp.705-712, July 8-11. 2002.
[25] J.K. Tugnait, \Tracking of multiple maneuvering targets in clutter using multiple
sensors, IMM and JPDA coupled ltering", IEEE Trans. Aerospace & Electronic
Systems, Vol. AES-40, pp. 320-330, Jan. 2004.
179
[26] B. Chen and J.K. Tugnait, \Tracking of multiple maneuvering targets in clutter
using IMM/JPDA ltering and xed-lag smoothing", Automatica, Vol. 37, pp.
239-249, Feb. 2001.
[27] Y. Bar-Shalom and W.D. Blair, Multitarget-Multisensor Tracking: Application and
Advances, Vol. III, Artech House, Norwood, MA, 2000.
[28] W. Koch and G.V. Keuk, \Multiple hypothesis Track Maintenance with Possi-
bly Unresolved Measurements", IEEE Transaction on Aerospace and Electronic
Systems, Vol. 33, No.3, pp. 883-892, July 1997.
[29] X.R. Li and Y. Bar-Shalom, \Design of an interacting multiple algorithm for air
tra c control tracking", IEEE Transactions on Control System Technology, pp.
186-194, Sept. 1993.
[30] H.A.P. Blom and Y. Bar-Shalom, \The interacting model multiple algorithm for
systems with Markovian switching coe cients", IEEE Trans. Automatic Control,
Vol. AC-33, pp. 780-783, Aug. 1988.
[31] X.R. Li and Y. Bar-Shalom, \Multiple model estimation with variable structure",
IEEE Transaction on Automatic Control, Vol. 41, No.4, pp. 478-493, April 1996.
[32] X.R. Li, X. Zhi, and Y. Zhang, \Multiple model estimation with variable structure
Part III: model group switching algorithm.", IEEE Transaction on Aerospace and
Electronic Systems, Vol. 351, No.1, pp. 225-240, Jan. 1999.
[33] Y. Bar-Shalom and X.R. Li, Estimation and Tracking: Principles, Techniques and
Software, Norwood, MA: Artech House, 1993.
[34] D. Willner, C.B. Chang, and K.P. Dunn, \Kalman lter algorithms for a multisen-
sor system.", Proceedings of 1976 IEEE Conf.on Decision and Control, Clearwater
Beach, FL, Dec. 1976.
[35] Y. Bar-Shalom, \Update with out-of-sequence measurements in tracking: exact
solution", Proc. SPIE Conf. Signal and Data Processing of Small Targets, Vol.
4048, pp. 541-556, Orlando, FL, April 2000.
[36] Y. Bar-Shalom, M. Mallick, H. Chen, and R. Washburn, \One-step solution for the
general out-of-sequence measurement problem in tracking", Proc. IEEE Aerospace
Conference, Big Sky, MT, Vol. 70, pp. 4.1551-4.1559, March 2002.
[37] M. Mallick, J. Krant, and Y. Bar-Shalom, \Multi-sensor multi-target tracking
using out-of-sequence measurements", Proc. Fifth International Conference on In-
formation Fusion(Fusion 2002), Annapolis, MD, USA.
180
[38] S. Jeong and J.K. Tugnait, \Parallel detection fusion for multisensor tracking of
a maneuvering target in clutter using IMMPDA ltering," in Proc. 37th Annual
Asilomar Conf. Signals Systems Computers, Paci c Grove, CA, pp. 1213-1217,
Nov. 2003.
[39] B.D.O. Anderson and J.B. Moore, Optimal Filtering, Eaglewood Cli s, NJ: Pren-
tice Hall, 1979.
[40] R.G. Brown and P.Y.C. Hwang, \Introduction to random signals and applied
Kalman ltering: with matlab exercises and solutions", John wiley and sons, 1997.
[41] M.S. Grewal and A.P. Andrews, \Kalman ltering: Theory and Practice" , Eagle-
wood Cli s, NJ: Prentice Hall, 1993.
[42] K.V. Ramachandra, Kalman Filtering Techniques for Radar Tracking, New York,
NY: Marcel Dekker, 2000.
[43] X.R. Li, \Hybrid estimation Techniques.",Control and dynamic systems: advances
in theory and application, pp. 1-76, Vol 76, C.T. Leondes(Ed), New York Academic,
1996.
[44] S. Blackman and R. Popoli, Design and Analysis of Modern Tracking Systems,
Norwood, MA: Artech House, 1999.
[45] F. Dufour and M. Mariton, \Tracking a 3D maneuvering target with passive sen-
sors," IEEE Trans. Aerospace and Electronic Systems, Vol. AES-27, pp. 725-738,
July 1991.
[46] N.N. Okello and G.W. Pulford, \Simultaneous registration and tracking for multi-
ple radars with cluttered measurements", Proc. 8th IEEE Workshop on Statistical
Signal and Array Processing, Corfu, Greece, pp. 60-63, June 1996.
[47] B. Chen and J.K. Tugnait, \Multisensor tracking of a maneuvering target in clutter
using immpda xed-lag smoothing", IEEE Transactions on Aerospace and Elec-
tronic Systems, Vol. 36, no.3, pp. 983-992, July 2000.
[48] G.W. Pulford and R.J. Evans, \Probabilistic sata association for systems with mul-
tiple simultaneous measurements", Automatica, Vol. 32(9), pp. 1311-1316, Sept.
1996.
[49] Y. Bar-Shalom, H. Chen, and M. Mallick, \One-step solution for the multistep out-
of-sequence measurement problem in tracking," IEEE Trans. Aerospace & Elec-
tronic Systems, Vol. 40, No.1, pp. 27-35, Jan. 2004.
[50] S. Jeong and J.K. Tugnait, \Multisensor tracking of a maneuvering target in clutter
with asynchronous measurements using IMMPDA ltering and parallel detection
fusion," in Proc. 2004 American Control Conf., Boston, MA, pp. 5350-5355, June
30-July 2, 2004.
181
[51] S. Challa, R. Evans, and X. Wang \A Bayesian solution and its approximations to
out-of-sequence measurement problems," in Journal of Information Fusion, Vol. 4,
No.3, pp. 185-199, Sept. 2003.
[52] X. Wang and S. Challa, \Augmented state IMM-PDA for OOSM solution to ma-
neuvering target tracking in clutter," in Radar Conference, 2003. Proceedings of
the International, pp. 479-485, Sept. 2003.
[53] G.W. Pulford and R. J. Evans, \A multipath data association tracker for over-
the-horizon radar", IEEE Transactions on Aerospace and Electronic Systems, Vol.
AES-34, pp. 1165-1183, Oct. 1998.
[54] L.C. Ng and Y. Bar-Shalom, \Modeling of unresolved measurements for multitarget
tracking", Proc. OCEANS ?81 Conf., Boston, MA, pp. 982-987, Sept. 1981.
[55] G.V. Trunk, \Range resolution of targets using automatic detectors", IEEE Trans.
Aerospace and Electronic Systems, Vol. AES-14, pp. 750-755, Sept. 1978.
[56] G.V. Trunk and J.D. Wilson, \Track initiation of occasionally unresolved radar
targets", IEEE Trans. Aerospace and Electronic Systems, Vol. AES-17, pp. 122-
130, Jan. 1981.
182
Appendices
183
Appendix A
Derivation of the mode-conditioned association event probability Eqn.
(3.46)
To evaluate the mode-conditioned association event probabilities, the conditioning is
broken down into the past data Zk 1 as well as the latest data Y 1k and Y 2k . A probabilistic
interference can be made on both the number of measurements in the validation region
(from the clutter density, if known) and on their location. This can be written out as
[21]
j;a;bk : = P[ a;bk jMjk;Y 1k;Y 2k;Zk 1] = p(
a;b
k ;M
j
k;Y
1k;Y 2k;Zk 1)
P[Mmk ;Y 1k;Y 2k;Zk 1]
= p(Y
1k;Y 2kj a;b
k ;M
j
k;Z
k 1)P[ a;b
k jM
j
k;Z
k 1]P[Mj
k;Z
k 1]
P
p(Y 1k;Y 2kj a;bk ;Mjk;Zk 1)P[ a;bk jMjk;Zk 1]P[Mjk;Zk 1]
= p(Y
1k;Y 2kj a;b
k ;M
j
k;Z
k 1)P[ a;b
k jM
j
k;Z
k 1]p[Zk 1jMj
k]P[M
j
k]P
p(Y 1k;Y 2kj a;bk ;Mjk;Zk 1)P[ a;bk jMjk;Zk 1]p[Zk 1jMjk]P[Mjk]
= 1cp(Y 1k;Y 2kj a;bk ;Mjk;Zk 1)P[ a;bk jMjk;Zk 1] (A.1)
where c is a normalization constant such that Pm1a=0Pm2b=0 j;a;bk = 1. The rst term of
the right side of (A.1) is [12, 21]
pfY 1k;Y 2kj a;bk ;Mjk;Zk 1g
184
=
8>
>>>
>>>>
>>>>
>><
>>>
>>>>
>>>>
>>>
:
V 1
k
m1 V 2
k
m2; a = 0; b = 0
P 1G2 V 2k m2+1 N
h
j;2(b)k ; 0;Sj;2k
i
; a = 0; b = 1; ; m2
P 1G1 V 1k m1+1 N
h
j;1(a)k ; 0;Sj;1k
i
; a = 1; ; m1; b = 0
P 1G1 V 1k m1+1P 1G2 V 2k m2+1 N
h
j;a;bk ; 0;Sjk
i
; a = 1; ; m1; b = 1; ; m2:
(A.2)
As in [12, 21], we have
P[ a;bjMjk;Zk 1] = P[ a;bj m1 = m1; m2 = m2] (A.3)
where we have conditioning on the total number of validated measurements obtained from
sensor l, ml = ml; in this notation ml is the random variable and ml is its realization.
Denoting by l the number of false measurements, one has
P[ a;bj m1 = m1; m2 = m2] = P[ a;bj 1 = m1; 2 = m2]P[ 1 = m1; 2 = m2]
+ P[ a;bj 1 = m1; 2 = m2 1]P[ 1 = m1; 2 = m2 1]
+ P[ a;bj 1 = m1 1; 2 = m2]P[ 1 = m1 1; 2 = m2]
+ P[ a;bj 1 = m1 1; 2 = m2 1]P[ 1 = m1 1; 2 = m2 1]
=
8
>>>>
>>>>
>>>>
>><
>>>>
>>>>
>>>>
>>:
1 P [ 1 = m1; 2 = m2] a = 0; b = 0
1
m2 P [ 1 = m1; 2 = m2 1] a = 0; b = 1; ; m2
1
m1 P [ 1 = m1 1; 2 = m2] a = 1; ; m1; b = 0
1
m1 m2 P [ 1 = m1 1; 2 = m2 1] a = 1; ; m1; b = 1; ; m2
(A.4)
185
since l must be either ml 1 (if the target has been detected from sensor l and its
measurement fell in the validation gate) or ml. Using Bayes? formula one has
P[ 1 = m1; 2 = m2jm1 = m1;m2 = m2] = (1 PD1PG1)(1 PD1PG1) F( m1) F( m2)Pfm1= m1;m2= m2g
P[ 1 = m1; 2 = m2 1jm1 = m1;m2 = m2] = (1 PD1PG1)(PD2PG2) F( m1) F( m2 1)Pfm1= m1;m2= m2g
P[ 1 = m1 1; 2 = m2jm1 = m1;m2 = m2] = (PD1PG1)(1 PD2PG2) F( m1 1) F( m2)Pfm1= m1;m2= m2g
P[ 1 = m1 1; 2 = m2 1jm1 = m1;m2 = m2] = (PD1PG1)(PD2PG2) F( m1 1) F( m2 1)Pfm1= m1;m2= m2g
:
(A.5)
where F(ml) is the probability mass function (pmf) of the number of false measurements
and PDlPGl is the probability that the target has been detected and its measurements
fell in the gate for sensor l. The common denominator in (A.5) is
P[m1 = m1;m2 = m2] = (1 PD1PG1) (1 PD1PG1) F( m1) F( m2)
+ (1 PD1PG1) (PD2PG2) F( m1) F( m2 1)
+ (PD1PG1) (1 PD2PG2) F( m1 1) F( m2)
+ (PD1PG1) (PD2PG2) F( m1 1) F( m2 1): (A.6)
From (A.5) and (A.6) we have
P[ a;bjMjk;Zk 1] = P[ a;bj m1; m2]
186
=
8>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
<
>>>
>>>>
>>>>
>>>
>>>>
>>>>
>>>>
>>>
>>>>
>:
(1 PD1PG1) F( m1) F( m1 1)
h
PD1PG1 + (1 PD1PG1) F( m1) F( m1 1)
i 1
(1 PD2PG2) F( m2) F( m2 1)
h
PD2PG2 + (1 PD2PG2) F( m2) F( m2 1)
i 1
; a = 0; b = 0
(1 PD1PG1) F( m1) F( m1 1)
h
PD1PG1 + (1 PD1PG1) F( m1) F( m1 1)
i 1
1 m2PD2PG2
h
PD2PG2 + (1 PD2PG2) F( m2) F( m2 1)
i 1
; a = 0; b = 1; ; m2
1
m1PD1PG1
h
PD1PG1 + (1 PD1PG1) F( m1) F( m1 1)
i 1
(1 PD2PG2) F( m2) F( m2 1)
h
PD2PG2 + (1 PD2PG2) F( m2) F( m2 1)
i 1
; a = 1; ; m1; b = 0
1
m1PD1PG1
h
PD1PG1 + (1 PD1PG1) F( m1) F( m1 1)
i 1
1 m2PD2PG2
h
PD2PG2 + (1 PD2PG2) F( m2) F( m2 1)
i 1
; a = 1; ; m1; b = 1; ; m2:
(A.7)
Using the (nonparametric) di use prior model [12]
F( ml) = F( ml 1) = (A.8)
where is irrelevant since it cancels out. Exploiting the di use model for clutter, we
obtain Eqn. (3.46).
187
Appendix B
Derivation of Eqn. (6.62)
For the coupled state vector xk of a group of two targets de ned in (6.2), the
conditional probability of unresolved targets given xk is de ned as (6.28)
P(Ajxk) = exp
12
h
hl(xk(1)) hl(xk(2))
iT
Rl;d 1
h
hl(xk(1)) hl(xk(2))
i
=
2 Rl;d
1=2N
0;hl;d(xk);Rl;d
: (B.1)
From the de nition of Eqns. (6.62) and (6.28), PJ;1;ak can be expressed as
PJ;1;ak =
Z
P(Ajxk)N(xk; ^xJkjk 1;PJkjk 1)dxk
=
Z
2 Rl;d
1=2N
0;hl;d(xk);Rl;d
N
xk; ^xJkjk 1;PJkjk 1
dxk: (B.2)
= 1
j2 Pkj1=2
Z
exp
12
h
( xk)T A 1 xk + (xk ^xk)T P 1k (xk ^xk)
i
dxk
where we rede ne
^xk := ^xJkjk 1; Pk := PJkjk 1; := HJ;1;dk ; A := R1;d: (B.3)
Let
=
h
( xk)T A 1 xk + (xk ^xk)T P 1k (xk ^xk)
i
: (B.4)
Recall that the trace of an n n matrix is the sum of its diagonal elements, and, from
the elementary properties of trace (?tr? stands for ?trace?):
1. tr(A) := Pni=1aii = tr(A0) where aii is the elements of A
188
2. tr(AB) = tr(BA) where A is m n and B is n m
3. tr(ABC) = tr(BCA)= tr(CAB) where A, B, and C are n n
4. xTAx = tr(xTAx) = tr(xTAx), where x is a n 1 column vector and A is n n
matrix,
it is clear that is equal to tr( ). Thus we can manipulate (B.4) as follows;
= tr
h
( xk)T A 1 xk + (xk ^xk)T P 1k (xk ^xk)
i
1 1
= tr
h
TA 1 +P 1k
xkxTk P 1k ^xkxTk ^xTkP 1k xk + ^xTkP 1k ^xk
i
= tr
h
R 1xkxTk QxTk QTxk + ^xTkP 1k ^xk
i
= tr
h
xTkR 1xk R 1RQxTk xkQTRR 1 + ^xTkP 1k ^xk
i
= tr
h
xTkR 1xk R 1UxTk xkUTR 1 +UTR 1U UTR 1U + ^xTkP 1k ^xk
i
=
h
xTkR 1xk R 1UxTk xkUTR 1 +UTR 1U UTR 1U + ^xTkP 1k ^xk
i
= (xk U)T R 1 (xk U) + ^xTkP 1k ^xk UTR 1U (B.5)
where, for simplicity,
Q := P 1k ^xk; R :=
TA 1 +P 1k
1
;
U := RQ =
TA 1 +P 1k
1
P 1k ^xk: (B.6)
Next from (B.3) and (B.5), Eqn. (6.62) can be rewritten as
PJ;1;ak = j2 Rj
1=2
j2 Pkj1=2 exp
12
h
^xTkP 1k ^xk UTR 1U
i
189
Z 1
j2 Rj1=2 exp
12 (xk U)T R 1 (xk U)
dxk: (B.7)
From the second term (integration part) in (B.7), since xk is a Gaussian random vector
with parameters U and R, the integration part of (B.7) turns out to be
Z 1
j2 Rj1=2 exp
12 (xk U)T R 1 (xk U)
dxk = 1: (B.8)
Therefore, now we get
PJ;1;ak = j2 Rj
1=2
j2 Pkj1=2 exp
12
h
^xTkP 1k ^xk UTR 1U
i
: (B.9)
Then from (B.6), the exponential part of (B.9) can be obtained as
^xTkP 1k ^xk UTR 1U
= ^xTkP 1k ^xk QTRQ
= ^xTkP 1k ^xk ^xTkP 1k RP 1k ^xk
= ^xTk
P 1k P 1k
h
TA 1 +P 1k
i 1
P 1k
^xk
= ^xTk
P 1k P 1k
Pk Pk T
h
Pk T +A
i 1
Pk
P 1k
^xk
= ^xTk
T
h
Pk T +A
i 1
^xk
= ( ^xk)T
h
Pk T +A
i 1
^xk:
(B.10)
Substituting ( ^xk)T
h
Pk T +A
i 1
^xk for ^xTkP 1k ^xk UTR 1U (B.10) thus yields
190
PJ;1;ak = jRj
1=2
jPkj1=2 exp
12 ( ^xk)T
h
Pk T +A
i 1
^xk
=
jRj1=2 (2 )nz=2
Pk T +A
1=2
jPkj1=2 N
0; ^xk; Pk T +A
: (B.11)
Recall the elementary identities of determinants (?det? stands for ?determinant?):
1. det(AB) = det(A)det(B) = det(BA) where A and B are square matrices
2. det(In AB) = det(Im BA) where A is n m and B is m n.
The determinant operation part of (B.11) can by simpli ed as (jAj = det(A))
jRj1=2j Pk T+Aj1=2
jPkj1=2 =
det
h
Pk T +A
RP 1k
i 1=2 : (B.12)
Then we also rewrite (B.12) as
det
h
Pk T +A
RP 1k
i
= det
Pk T +A
Pk Pk T
h
Pk T +A
i 1
P
P 1k
= det
Pk T +A
I Pk T
h
Pk T +A
i 1
= det
Pk T +A
det
I Pk T
h
Pk T +A
i 1
= det
Pk T +A
det
I
h
Pk T +A
i 1
Pk T
= det
Pk T +A Pk T
= det (A) =jAj:
(B.13)
191
From (B.3), (B.11), (B.12), and (B.13), it turns out that
PJ;1;ak := (2 )nz=2jAj1=2N
0; ^xk; Pk T +A
=
2 R1;d
1=2N
0
BB
@0;
h
Hj1;1k Hj2;1k
i
xJkjk 1;
Hj1;1k Hj2;1k
PJ;1kjk 1
2
66
4
Hj1;10k
Hj2;10k
3
77
5+R1;d
1
CC
A:
(B.14)
192