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