# Tracking Maneuvering Group Target with Extension Predicted and Best Model Augmentation Method Adapted

Authors: Linhai Gan; Gang Wang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The random matrix (RM) method is widely applied for group target tracking. The assumption that the group extension keeps invariant in conventional RM method is not yet valid, as the orientation of the group varies rapidly while it is maneuvering; thus, a new approach with group extension predicted is derived here. To match the group maneuvering, a best model augmentation (BMA) method is introduced. The existing BMA method uses a fixed basic model set, which may lead to a poor performance when it could not ensure basic coverage of true motion modes. Here, a maneuvering group target tracking algorithm is proposed, where the group extension prediction and the BMA adaption are exploited. The performance of the proposed algorithm will be illustrated by simulation.


---

## Body

## 1. Introduction

Groups are structured objects and formations of entities moving in a coordinated manner [1]. Limited by poor sensor resolution and little requirement in application, the objects to be tracked are considered as point sources in conventional sense [2]. With ever-increasing sensor resolution capability and requirements for tracking a crowd, a herd, or an air fleet, group target tracking has been a hot topic in recent years. An increasing capability of sensor technology makes it possible to learn more feature information about a group. However, there will be a fluctuation number of detections for a group of closely spaced targets, accounting for limited sensor resolution, which will prevent a successful tracking of the individual targets [2]. Thus, there are some new challenges in tracking a group target compared with tracking point one.The methods for group target tracking mainly include tracking via Poisson likelihoods [3, 4], group’s evolution modeled by a dynamic Bayesian network [5], the approach with random matrices [2, 6–11], group tracking using random finite sets [12, 13], and so on. A random matrix method was proposed by Koch [8] in 2008, which described the extension of the group by a symmetrical positive definite (SPD) random matrix, and characterized by simple filter equations, a small number in preset parameters, strong robustness, and so on. On the basis of [8], an interacting multiple model (IMM) structure was introduced by Feldmann [2, 7], and the performance for maneuvering extended objects tracking was improved; moreover, with the measurement noise being considered in the innovation covariance, the application area of the algorithm was expanded. A suitable measurement model to address the radar’s measurement noise and its conversion from polar coordinates to Cartesian coordinates was proposed in [14]. A Bayesian framework was developed in [15], within which the probability density function (PDF) of the object state and extension and the probability mass function of the object class were jointly obtained; a joint tracking and classification algorithm for an extended object was proposed using random matrix consequently. Within the random matrix framework, a nonellipsoidal extended object (NEO) was approximated by multiple ellipsoidal subobjects, with each subobject described by a random matrix in [11, 16]. In [17], the orientation of the extension in the next time was predicted by the extension estimation of the last two moments; the assumption that the extension of the adjacent two moments keeps invariant was no longer valid, but there was no analysis in theory about this. A multiple group target tracking algorithm under clutter environment was proposed based on random matrix method in [18].IMM approach is one of the most wildly used tracking algorithms for maneuvering target; with the interacting of models, fusion estimation for maneuvering target was obtained. With a fixed structure, many more models are needed for an efficient tracking performance, as the target maneuverings in a more complex mode, which will lead to an increase in computation burden and model competition consequently. One of the variable structure multiple model (VSMM) methods called best model augmentation (BMA) was proposed in [19], in which the state estimation model set was adjusted adaptively at each step, which will alleviate the problem of the IMM approach referred above.The paper is organized as follows. Section2 gives a derivation of a RM method for group target tracking under Bayesian framework, where the group extension prediction is available. Section 3 describes the BMA approach and its modified form. Section 4 outlines the procedure of the proposed algorithm. The performance of the algorithm is presented via simulations in Section 5. Finally, conclusions are given in Section 6.


## 2. Bayesian Group Target Tracking

### 2.1. Background

Considering a target group with an ellipsoidal shape, we use a SPD random matrixXk to describe the extension of scan k. Consider a 2-dimensional extension, Xk, can be expressed as Xk=Ckdiagr12,r22CkT, with r1 and r2 denoting the main axis and minor axis of the ellipse, respectively, and Ck is a rotation matrix.It is assumed that in each scank there is a random number of nk independent position measurements:(1)zkj=Hxk+wkj,where xk and zkjj=1,2,…,nk denote the state to be estimated and the jth sensor measurement of scan k, respectively. Moreover, H denotes the measurement matrix, wkj is the jth measurement noise and is assumed to be a zero mean normally distributed random vector with covariance λXk∣k-1+R, where λ is a scaling factor, Xk∣k-1 is the predicted extended state, and R is the covariance matrix of the sensor measurement error.In a Bayesian view, a group target tracking algorithm is an iterative updating scheme for conditional probability densitiespxk,Xk∣Zk at each time tk [8]. Its product representation can be written as (2)pxk,Xk∣Zk=pxk∣Xk,ZkpXk∣Zkor(3)pxk,Xk∣Zk=pXk∣xk,Zkpxk∣Zkwith its Bayes’ formula as(4)pxk,Xk∣Zk=pxk,Xk∣Zk-1pZk,nk∣xk,Xk∫pxk,Xk∣Zk-1pZk,nk∣xk,XkdxkdXk,where Zk=zkjj=1nk and Zk=Zl,nll=1k are used to denote the set of the nk measurements in a particular scan and for the accumulated measurements set among the time. The conditional probability densities pxk∣Xk,Zk and pXk∣xk,Zk describe the kinematical object properties and the extended properties, respectively. And pZk,nk∣xk,Xk denotes the likelihood function and pxk,Xk∣Zk-1 is a predicted joint probability density function.With the assumption that the extension does not tend to change over time in [8], that is, Xk=Xk-1, moreover Xk∣k-1=Xk-1∣k-1, where Xk-1∣k-1 represents the estimation of the extended state. The joint prediction of kinematic state and extension can be derived from (2):(5)pxk,Xk∣Zk-1=∫pxk,Xk∣xk-1,Xk-1,Zk-1pxk-1,Xk-1∣Zk-1dxk-1dXk-1,pxk,Xk∣Zk-1=pxk∣Xk,Zk-1pXk∣Zk-1.The prediction of kinematic state and extension can be expressed, respectively, as follows.(A) Kinematic state is as follows:(6)pxk∣Xk,Zk-1=Nxk;xk∣k-1,Pk∣k-1⊗Xk,xk∣k-1=Fk∣k-1⊗Idxk-1∣k-1,Pk∣k-1=Fk∣k-1Pk-1∣k-1Fk∣k-1T+Qk-1, where xk∣k-1 and xk-1∣k-1 denote prediction and estimation of kinematic state, respectively; Id is a d-dimension identity matrix; Pk∣k-1 and Pk-1∣k-1 denote the predicted covariance and estimated covariance, respectively; Fk∣k-1 denotes the state transform matrix, that is, for a constant velocity model, Fk∣k-1=1T01, where T is the sample period; Qk-1 denotes the covariance of process noise; ⊗ stands for Kronecker products.(B) Extension is as follows:(7)pXk∣xk,Xk-1=WXk;δk∣k-1,Xk-1δk∣k-1,pXk∣Zk-1=GBdIIXk;12δk∣k-1,12νk-1∣k-1-d-1,Xk-1∣k-1δk∣k-1,where W· denotes a Wishart distribution function; δk∣k-1=δe-T/τ denotes freedom degree, with δ and τ denoting the evolution parameters of extension; T denotes the sensor sample period; GBdIIXk;a,b,A,B denotes “Generalized Beta Type II” density.In reality, the assumption that extension does not tend to change over time is not so valid in many situations, for example, when the group target is moving in a turning mode, the orientation of the group extension will change timely; a proper prediction to the extension may improve the tracking performance.


### 2.2. Innovation Prediction

It can be learned heuristically from [2] that the joint prediction pxk,Xk∣Zk-1 can be derived as (8)pxk,Xk∣Zk-1=∫pXk∣xk,Xk-1pxk∣xk-1×pXk-1∣xk-1,Zk-1pxk-1∣Zk-1dxk-1dXk-1.While the group target has a relatively large extension and a dense scattering, the kinematic state may have little influence on the extension. So, posterior probability density of the extension is nearly independent from the kinematic state, but mainly dependent on the measurement. Thus we make the following assumption.Assumption 1.
The posterior probability density of extension can be approximated by marginal probability(9)pXk∣xk,Zk≈pXk∣Zk.Thus(10)pXk-1∣xk-1,Zk-1≈pXk-1∣Zk-1=IWXk-1;νk-1∣k-1,αk-1∣k-1Xk-1∣k-1,where IWXk-1;α,A denotes the inverse Wishart density ofd-dimension SPD random matrix Xk-1, α denotes freedom degree; αk-1∣k-1=νk-1∣k-1-d-1, and νk-1∣k-1 is the estimation of the freedom degree. Thus, (8) can be approximated as (11)pxk,Xk∣Zk-1=∫pxk∣xk-1pxk-1∣Zk-1dxk-1×∫pXk∣xk,Xk-1pXk-1∣Zk-1dXk-1=pxk∣Zk-1pXk∣xk,Zk-1.


#### 2.2.1. Kinematic State

It can be learned from (3) and (11) that the prediction and estimation of extension and kinematic state can be calculated, respectively. Assume that the process noise is a zero mean normally distributed random vector with covariance Qk-1. Thus, the prediction of kinematic state is the same as the standard Kalman filter’s; xk∣k-1 and Pk∣k-1 can be given by(12)xk∣k-1=Fk∣k-1xk-1∣k-1,Pk∣k-1=Fk∣k-1Pk-1∣k-1Fk∣k-1T+Qk-1.


#### 2.2.2. Extension

Naval vessels, Submarines, or ground moving convoys show a clear orientation [9]. As for airplane formation, the angle between its moving orientation and the extension appearance will keep invariant (or a relatively small change) in general. Thus, the group’s extension can be predicted according to its moving orientation in many applications. Assume that(13)Xk∣k-1=Ck∣k-1Xk-1∣k-1Ck∣k-1T,where Ck∣k-1 denotes a predicted rotation matrix, while in a 2-dimensional situation, it satisfies: (14)Ck∣k-1=cos⁡θk∣k-1-sin⁡θk∣k-1sin⁡θk∣k-1cos⁡θk∣k-1,θk∣k-1=atan⁡vk∣k-1,yvk∣k-1,x-atan⁡vk-1∣k-1,yvk-1∣k-1,x,where θk∣k-1 denotes the prediction of the angle that the group target rotates in two adjacent moments. vx and vy denote the velocity in x- and y-axis, respectively, and are parts of kinematic state x. According to Assumption 1, Xk-1∣k-1 and Ck∣k-1 are approximately independent. Thus, pXk∣xk,Xk-1 can be expressed as (15), according to the property of Wishart distribution: (15)pXk∣xk,Xk-1=WXk;δk∣k-1,Ck∣k-1Xk-1Ck∣k-1Tδk∣k-1.Thus the distribution of pXk∣xk,Zk-1 can be derived as [8] (16)pXk∣xk,Zk-1=GBdIIXk;12δk∣k-1,12νk-1∣k-1-d-1,αk-1∣k-1Ck∣k-1Xk-1Ck∣k-1Tδk∣k-1≈IWXk-1;νk-1∣k-1∗,Xk-1∣k-1∗.


### 2.3. State Update

With the regularities of distribution being invariant, the state update procedure is similar to [2].(17)xk∣k=xk∣k-1+Kk∣k-1z-k-Hxk∣k-1,Pk∣k=Pk∣k-1-Kk∣k-1Sk∣k-1Kk∣k-1T,where(18)Kk∣k-1=Pk∣k-1HTSk∣k-1-1,Sk∣k-1=HPk∣k-1HT+Yk∣k-1nk,Yk∣k-1=λXk∣k-1+R,with xk∣k and Pk∣k denoting the estimation of the kinematic state and its covariance, respectively.Extension Update
(19)
X
k
∣
k
=
1
α
k
∣
k
α
k
∣
k
-
1
X
k
∣
k
-
1
+
N
^
k
∣
k
-
1
+
Y
^
k
∣
k
-
1
,where Xk∣k denotes the estimation of the extension. And(20)Xk∣k-1=Xk∣k-11/2Xk∣k-11/2T,N^k∣k-1=Xk∣k-11/2Sk∣k-1-1/2Nk∣k-1Sk∣k-1-1/2TXk∣k-11/2T,Y^k∣k-1=Xk∣k-11/2Yk∣k-1-1/2Y-k∣k-1Yk∣k-1-1/2TXk∣k-11/2T,Nk∣k-1=z-k-Hxk∣k-1z-k-Hxk∣k-1T,Y-k∣k-1=∑j=1nkzkj-z-kzkj-z-kT,z-k=1nk∑j=1nkzkj,αk∣k=αk∣k-1+nk,αk∣k-1=2+exp⁡-Tταk-1∣k-1-2.


## 2.1. Background

Considering a target group with an ellipsoidal shape, we use a SPD random matrixXk to describe the extension of scan k. Consider a 2-dimensional extension, Xk, can be expressed as Xk=Ckdiagr12,r22CkT, with r1 and r2 denoting the main axis and minor axis of the ellipse, respectively, and Ck is a rotation matrix.It is assumed that in each scank there is a random number of nk independent position measurements:(1)zkj=Hxk+wkj,where xk and zkjj=1,2,…,nk denote the state to be estimated and the jth sensor measurement of scan k, respectively. Moreover, H denotes the measurement matrix, wkj is the jth measurement noise and is assumed to be a zero mean normally distributed random vector with covariance λXk∣k-1+R, where λ is a scaling factor, Xk∣k-1 is the predicted extended state, and R is the covariance matrix of the sensor measurement error.In a Bayesian view, a group target tracking algorithm is an iterative updating scheme for conditional probability densitiespxk,Xk∣Zk at each time tk [8]. Its product representation can be written as (2)pxk,Xk∣Zk=pxk∣Xk,ZkpXk∣Zkor(3)pxk,Xk∣Zk=pXk∣xk,Zkpxk∣Zkwith its Bayes’ formula as(4)pxk,Xk∣Zk=pxk,Xk∣Zk-1pZk,nk∣xk,Xk∫pxk,Xk∣Zk-1pZk,nk∣xk,XkdxkdXk,where Zk=zkjj=1nk and Zk=Zl,nll=1k are used to denote the set of the nk measurements in a particular scan and for the accumulated measurements set among the time. The conditional probability densities pxk∣Xk,Zk and pXk∣xk,Zk describe the kinematical object properties and the extended properties, respectively. And pZk,nk∣xk,Xk denotes the likelihood function and pxk,Xk∣Zk-1 is a predicted joint probability density function.With the assumption that the extension does not tend to change over time in [8], that is, Xk=Xk-1, moreover Xk∣k-1=Xk-1∣k-1, where Xk-1∣k-1 represents the estimation of the extended state. The joint prediction of kinematic state and extension can be derived from (2):(5)pxk,Xk∣Zk-1=∫pxk,Xk∣xk-1,Xk-1,Zk-1pxk-1,Xk-1∣Zk-1dxk-1dXk-1,pxk,Xk∣Zk-1=pxk∣Xk,Zk-1pXk∣Zk-1.The prediction of kinematic state and extension can be expressed, respectively, as follows.(A) Kinematic state is as follows:(6)pxk∣Xk,Zk-1=Nxk;xk∣k-1,Pk∣k-1⊗Xk,xk∣k-1=Fk∣k-1⊗Idxk-1∣k-1,Pk∣k-1=Fk∣k-1Pk-1∣k-1Fk∣k-1T+Qk-1, where xk∣k-1 and xk-1∣k-1 denote prediction and estimation of kinematic state, respectively; Id is a d-dimension identity matrix; Pk∣k-1 and Pk-1∣k-1 denote the predicted covariance and estimated covariance, respectively; Fk∣k-1 denotes the state transform matrix, that is, for a constant velocity model, Fk∣k-1=1T01, where T is the sample period; Qk-1 denotes the covariance of process noise; ⊗ stands for Kronecker products.(B) Extension is as follows:(7)pXk∣xk,Xk-1=WXk;δk∣k-1,Xk-1δk∣k-1,pXk∣Zk-1=GBdIIXk;12δk∣k-1,12νk-1∣k-1-d-1,Xk-1∣k-1δk∣k-1,where W· denotes a Wishart distribution function; δk∣k-1=δe-T/τ denotes freedom degree, with δ and τ denoting the evolution parameters of extension; T denotes the sensor sample period; GBdIIXk;a,b,A,B denotes “Generalized Beta Type II” density.In reality, the assumption that extension does not tend to change over time is not so valid in many situations, for example, when the group target is moving in a turning mode, the orientation of the group extension will change timely; a proper prediction to the extension may improve the tracking performance.


## 2.2. Innovation Prediction

It can be learned heuristically from [2] that the joint prediction pxk,Xk∣Zk-1 can be derived as (8)pxk,Xk∣Zk-1=∫pXk∣xk,Xk-1pxk∣xk-1×pXk-1∣xk-1,Zk-1pxk-1∣Zk-1dxk-1dXk-1.While the group target has a relatively large extension and a dense scattering, the kinematic state may have little influence on the extension. So, posterior probability density of the extension is nearly independent from the kinematic state, but mainly dependent on the measurement. Thus we make the following assumption.Assumption 1.
The posterior probability density of extension can be approximated by marginal probability(9)pXk∣xk,Zk≈pXk∣Zk.Thus(10)pXk-1∣xk-1,Zk-1≈pXk-1∣Zk-1=IWXk-1;νk-1∣k-1,αk-1∣k-1Xk-1∣k-1,where IWXk-1;α,A denotes the inverse Wishart density ofd-dimension SPD random matrix Xk-1, α denotes freedom degree; αk-1∣k-1=νk-1∣k-1-d-1, and νk-1∣k-1 is the estimation of the freedom degree. Thus, (8) can be approximated as (11)pxk,Xk∣Zk-1=∫pxk∣xk-1pxk-1∣Zk-1dxk-1×∫pXk∣xk,Xk-1pXk-1∣Zk-1dXk-1=pxk∣Zk-1pXk∣xk,Zk-1.


### 2.2.1. Kinematic State

It can be learned from (3) and (11) that the prediction and estimation of extension and kinematic state can be calculated, respectively. Assume that the process noise is a zero mean normally distributed random vector with covariance Qk-1. Thus, the prediction of kinematic state is the same as the standard Kalman filter’s; xk∣k-1 and Pk∣k-1 can be given by(12)xk∣k-1=Fk∣k-1xk-1∣k-1,Pk∣k-1=Fk∣k-1Pk-1∣k-1Fk∣k-1T+Qk-1.


### 2.2.2. Extension

Naval vessels, Submarines, or ground moving convoys show a clear orientation [9]. As for airplane formation, the angle between its moving orientation and the extension appearance will keep invariant (or a relatively small change) in general. Thus, the group’s extension can be predicted according to its moving orientation in many applications. Assume that(13)Xk∣k-1=Ck∣k-1Xk-1∣k-1Ck∣k-1T,where Ck∣k-1 denotes a predicted rotation matrix, while in a 2-dimensional situation, it satisfies: (14)Ck∣k-1=cos⁡θk∣k-1-sin⁡θk∣k-1sin⁡θk∣k-1cos⁡θk∣k-1,θk∣k-1=atan⁡vk∣k-1,yvk∣k-1,x-atan⁡vk-1∣k-1,yvk-1∣k-1,x,where θk∣k-1 denotes the prediction of the angle that the group target rotates in two adjacent moments. vx and vy denote the velocity in x- and y-axis, respectively, and are parts of kinematic state x. According to Assumption 1, Xk-1∣k-1 and Ck∣k-1 are approximately independent. Thus, pXk∣xk,Xk-1 can be expressed as (15), according to the property of Wishart distribution: (15)pXk∣xk,Xk-1=WXk;δk∣k-1,Ck∣k-1Xk-1Ck∣k-1Tδk∣k-1.Thus the distribution of pXk∣xk,Zk-1 can be derived as [8] (16)pXk∣xk,Zk-1=GBdIIXk;12δk∣k-1,12νk-1∣k-1-d-1,αk-1∣k-1Ck∣k-1Xk-1Ck∣k-1Tδk∣k-1≈IWXk-1;νk-1∣k-1∗,Xk-1∣k-1∗.


## 2.2.1. Kinematic State

It can be learned from (3) and (11) that the prediction and estimation of extension and kinematic state can be calculated, respectively. Assume that the process noise is a zero mean normally distributed random vector with covariance Qk-1. Thus, the prediction of kinematic state is the same as the standard Kalman filter’s; xk∣k-1 and Pk∣k-1 can be given by(12)xk∣k-1=Fk∣k-1xk-1∣k-1,Pk∣k-1=Fk∣k-1Pk-1∣k-1Fk∣k-1T+Qk-1.


## 2.2.2. Extension

Naval vessels, Submarines, or ground moving convoys show a clear orientation [9]. As for airplane formation, the angle between its moving orientation and the extension appearance will keep invariant (or a relatively small change) in general. Thus, the group’s extension can be predicted according to its moving orientation in many applications. Assume that(13)Xk∣k-1=Ck∣k-1Xk-1∣k-1Ck∣k-1T,where Ck∣k-1 denotes a predicted rotation matrix, while in a 2-dimensional situation, it satisfies: (14)Ck∣k-1=cos⁡θk∣k-1-sin⁡θk∣k-1sin⁡θk∣k-1cos⁡θk∣k-1,θk∣k-1=atan⁡vk∣k-1,yvk∣k-1,x-atan⁡vk-1∣k-1,yvk-1∣k-1,x,where θk∣k-1 denotes the prediction of the angle that the group target rotates in two adjacent moments. vx and vy denote the velocity in x- and y-axis, respectively, and are parts of kinematic state x. According to Assumption 1, Xk-1∣k-1 and Ck∣k-1 are approximately independent. Thus, pXk∣xk,Xk-1 can be expressed as (15), according to the property of Wishart distribution: (15)pXk∣xk,Xk-1=WXk;δk∣k-1,Ck∣k-1Xk-1Ck∣k-1Tδk∣k-1.Thus the distribution of pXk∣xk,Zk-1 can be derived as [8] (16)pXk∣xk,Zk-1=GBdIIXk;12δk∣k-1,12νk-1∣k-1-d-1,αk-1∣k-1Ck∣k-1Xk-1Ck∣k-1Tδk∣k-1≈IWXk-1;νk-1∣k-1∗,Xk-1∣k-1∗.


## 2.3. State Update

With the regularities of distribution being invariant, the state update procedure is similar to [2].(17)xk∣k=xk∣k-1+Kk∣k-1z-k-Hxk∣k-1,Pk∣k=Pk∣k-1-Kk∣k-1Sk∣k-1Kk∣k-1T,where(18)Kk∣k-1=Pk∣k-1HTSk∣k-1-1,Sk∣k-1=HPk∣k-1HT+Yk∣k-1nk,Yk∣k-1=λXk∣k-1+R,with xk∣k and Pk∣k denoting the estimation of the kinematic state and its covariance, respectively.Extension Update
(19)
X
k
∣
k
=
1
α
k
∣
k
α
k
∣
k
-
1
X
k
∣
k
-
1
+
N
^
k
∣
k
-
1
+
Y
^
k
∣
k
-
1
,where Xk∣k denotes the estimation of the extension. And(20)Xk∣k-1=Xk∣k-11/2Xk∣k-11/2T,N^k∣k-1=Xk∣k-11/2Sk∣k-1-1/2Nk∣k-1Sk∣k-1-1/2TXk∣k-11/2T,Y^k∣k-1=Xk∣k-11/2Yk∣k-1-1/2Y-k∣k-1Yk∣k-1-1/2TXk∣k-11/2T,Nk∣k-1=z-k-Hxk∣k-1z-k-Hxk∣k-1T,Y-k∣k-1=∑j=1nkzkj-z-kzkj-z-kT,z-k=1nk∑j=1nkzkj,αk∣k=αk∣k-1+nk,αk∣k-1=2+exp⁡-Tταk-1∣k-1-2.


## 3. Improved BMA Approach

### 3.1. Rationale of BMA

Models with fixed parameter are used to form the candidate model set of BMA, whose structure is different from that of the basic model set. The model who can best match the true motion mode of the target in the candidate model set will be activated. Moreover, the activated model will be selected to form a state estimation model set (SEMS) of a particular moment combination with the basic model set. All the estimations of the models in SEMS are fused to output the final estimation [19].Assume thatZk denotes the measurement sequence; the SEMS Mk is(21)Mk=Mb+Mka,where Mb denotes a small set of necessary models and it is adopted to ensure a basic coverage of true modes; Mka denotes the activated model set at k, which is activated from Mkc to improve the estimation performance; Mkc denotes the candidate model set of k.In order to select the model to be activated, a Kullback-Leibler (KL) criterion will be adopted to quantify the difference between a candidate model and the true motion mode; the one who has minimum divergence in KL with true mode will be selected to replace the activated model set of the previous. The SEMS is updated consequently.The discrepancy of modelmc,kj∈Mkc from true mode sk can be provided by a metric Dysk,mc,kj:(22)Dysk,mc,kj=12ln⁡ΣkjΣksk-n+trΣkj-1Σksk+y-ksk-y-kj·T,where (23)y-ksk=Ey∣Mk-1,sk,Z,Σksk=Ey-y-ksk·T∣Mk-1,sk,Z,y-kj=Ey∣Mk-1,mc,kj,Z,Σkj=Ey-y-kj·T∣Mk-1,mc,kj,Z,with tr· denoting the matrix trace; Mk-1=M1,M2,…,Mk-1 denotes the sequence of model set through k-1, and n denotes the dimension of y, with y acting as a common variable of sk and mc,kj. Particularly, y≜zk denotes the measurement vector; Z≜Zk-1 denotes the measurement sequence. As the true model is unknown in practice, sk is approximated by M^k∣k-1(24)M^k∣k-1=mkj∣pmkj∣mk-1j>0,mk-1j∈Mk-1,where pmkj∣mk-1j denotes the model transition probability.Ifmkj∈M^k∣k-1(25)y-ksk≈Ezk∣Mk-1,M^k∣k-1,Zk-1=EEzk∣mkj∈M^k∣k-1,Mk-1,Zk-1∣Mk-1,M^k∣k-1,Zk-1=∑mkj∈M^k∣k-1zk∣k-1jμk∣k-1j,Σksk≈Ezk-y-ksk·T∣Mk-1,M^k∣k-1,Zk-1=∑mkj∈M^k∣k-1Pk∣k-1j+zk∣k-1j-y-ksk·Tμk∣k-1j,where (26)zk∣k-1j=Ezk∣mkj,Mk-1,Zk-1,Pk∣k-1j=Ezk-zk∣k-1j·T∣mkj,Mk-1,Zk-1,μk∣k-1j≜Pmkj∣Mk-1,M^k∣k-1,Zk-1.Otherwise(27)y-kj=Ezk∣Mk-1,mc,kj,Zk-1≈Ezk∣mc,kj,xk-1∣k-1,Pk-1∣k-1,Σkj=Ezk-y-kj·T∣Mk-1,mc,kj,Zk-1≈Ezk-y-kj·T∣mc,kj,xk-1∣k-1,Pk-1∣k-1.The model to be activated from Mkc is(28)m^k=argminmc,kj∈MkcDysk,mc,kj.


### 3.2. BMA Adaptation

As shown in Section3.1, the SEMS of BMA are mainly composed of the basic model set, which keeps invariant all over the tracking period. This will bring the following two problems:①
With the basic model set selected offline, more prior knowledge is needed to ensure a basic coverage of true modes, which is a difficult task in practice.②
To cover a complicate maneuvering mode, the basic model set should be augmented, which will make only a part of model efficiently match the maneuvering at each moment; the other models just do little help in improving the tracking performance, moreover, will bring a heavy burden in computation and model competition.To overcome the drawbacks of BMA method mentioned above, an adaptive structure is considered for a better performance in maneuvering tracking. Suppose that(29)Mk=Mk-1-Mkd+Mka,Mkc=Mk-1c+Mkd-Mka,where Mkd∈Mkc∩Mk-1 denotes the model set eliminated from Mk-1, which is then added to form Mkc at k; Mka∈Mk-1c∩Mk denotes the optimal model set selected from Mk-1c, which is activated at k. Equation (29) shows the update of the estimation model set and the candidate model set; moreover, the basic model set of k step can be updated as Mkb=Mk-1-Mkd.The metricDysk,m~kj represents the discrepancy between model m~kj and true motion mode sk, which can be calculated just as (22). Moreover, m~kj∈Mk-1c+Mk-1 will be selected from either Mk-1c or Mk-1.If(30)minm~kj∈Mk-1cDysk,m~kj≥maxm~kj∈Mk-1Dysk,m~kj,then(31)Mk=Mk-1,Mkc=Mk-1c.Otherwise, the activated model is(32)m^k=argminm~kj∈Mk-1cDysk,m~kj.The model eliminated from SEMSMk-1 is(33)m^kd=argmaxm~kj∈Mk-1Dysk,m~kj.


## 3.1. Rationale of BMA

Models with fixed parameter are used to form the candidate model set of BMA, whose structure is different from that of the basic model set. The model who can best match the true motion mode of the target in the candidate model set will be activated. Moreover, the activated model will be selected to form a state estimation model set (SEMS) of a particular moment combination with the basic model set. All the estimations of the models in SEMS are fused to output the final estimation [19].Assume thatZk denotes the measurement sequence; the SEMS Mk is(21)Mk=Mb+Mka,where Mb denotes a small set of necessary models and it is adopted to ensure a basic coverage of true modes; Mka denotes the activated model set at k, which is activated from Mkc to improve the estimation performance; Mkc denotes the candidate model set of k.In order to select the model to be activated, a Kullback-Leibler (KL) criterion will be adopted to quantify the difference between a candidate model and the true motion mode; the one who has minimum divergence in KL with true mode will be selected to replace the activated model set of the previous. The SEMS is updated consequently.The discrepancy of modelmc,kj∈Mkc from true mode sk can be provided by a metric Dysk,mc,kj:(22)Dysk,mc,kj=12ln⁡ΣkjΣksk-n+trΣkj-1Σksk+y-ksk-y-kj·T,where (23)y-ksk=Ey∣Mk-1,sk,Z,Σksk=Ey-y-ksk·T∣Mk-1,sk,Z,y-kj=Ey∣Mk-1,mc,kj,Z,Σkj=Ey-y-kj·T∣Mk-1,mc,kj,Z,with tr· denoting the matrix trace; Mk-1=M1,M2,…,Mk-1 denotes the sequence of model set through k-1, and n denotes the dimension of y, with y acting as a common variable of sk and mc,kj. Particularly, y≜zk denotes the measurement vector; Z≜Zk-1 denotes the measurement sequence. As the true model is unknown in practice, sk is approximated by M^k∣k-1(24)M^k∣k-1=mkj∣pmkj∣mk-1j>0,mk-1j∈Mk-1,where pmkj∣mk-1j denotes the model transition probability.Ifmkj∈M^k∣k-1(25)y-ksk≈Ezk∣Mk-1,M^k∣k-1,Zk-1=EEzk∣mkj∈M^k∣k-1,Mk-1,Zk-1∣Mk-1,M^k∣k-1,Zk-1=∑mkj∈M^k∣k-1zk∣k-1jμk∣k-1j,Σksk≈Ezk-y-ksk·T∣Mk-1,M^k∣k-1,Zk-1=∑mkj∈M^k∣k-1Pk∣k-1j+zk∣k-1j-y-ksk·Tμk∣k-1j,where (26)zk∣k-1j=Ezk∣mkj,Mk-1,Zk-1,Pk∣k-1j=Ezk-zk∣k-1j·T∣mkj,Mk-1,Zk-1,μk∣k-1j≜Pmkj∣Mk-1,M^k∣k-1,Zk-1.Otherwise(27)y-kj=Ezk∣Mk-1,mc,kj,Zk-1≈Ezk∣mc,kj,xk-1∣k-1,Pk-1∣k-1,Σkj=Ezk-y-kj·T∣Mk-1,mc,kj,Zk-1≈Ezk-y-kj·T∣mc,kj,xk-1∣k-1,Pk-1∣k-1.The model to be activated from Mkc is(28)m^k=argminmc,kj∈MkcDysk,mc,kj.


## 3.2. BMA Adaptation

As shown in Section3.1, the SEMS of BMA are mainly composed of the basic model set, which keeps invariant all over the tracking period. This will bring the following two problems:①
With the basic model set selected offline, more prior knowledge is needed to ensure a basic coverage of true modes, which is a difficult task in practice.②
To cover a complicate maneuvering mode, the basic model set should be augmented, which will make only a part of model efficiently match the maneuvering at each moment; the other models just do little help in improving the tracking performance, moreover, will bring a heavy burden in computation and model competition.To overcome the drawbacks of BMA method mentioned above, an adaptive structure is considered for a better performance in maneuvering tracking. Suppose that(29)Mk=Mk-1-Mkd+Mka,Mkc=Mk-1c+Mkd-Mka,where Mkd∈Mkc∩Mk-1 denotes the model set eliminated from Mk-1, which is then added to form Mkc at k; Mka∈Mk-1c∩Mk denotes the optimal model set selected from Mk-1c, which is activated at k. Equation (29) shows the update of the estimation model set and the candidate model set; moreover, the basic model set of k step can be updated as Mkb=Mk-1-Mkd.The metricDysk,m~kj represents the discrepancy between model m~kj and true motion mode sk, which can be calculated just as (22). Moreover, m~kj∈Mk-1c+Mk-1 will be selected from either Mk-1c or Mk-1.If(30)minm~kj∈Mk-1cDysk,m~kj≥maxm~kj∈Mk-1Dysk,m~kj,then(31)Mk=Mk-1,Mkc=Mk-1c.Otherwise, the activated model is(32)m^k=argminm~kj∈Mk-1cDysk,m~kj.The model eliminated from SEMSMk-1 is(33)m^kd=argmaxm~kj∈Mk-1Dysk,m~kj.


## 4. Fusion Procedure

(
1) Reinitialization. ① Mixed probability is(34)μk-1∣k-1i∣j=pMkj∣Mk-1i,Zk-1=π-ijμk-1ic-j,where μk-1i denotes the model probability of k-1 step; c-j=∑i=1mπ-ijμk-1i, with π-ij denoting the transition probability and m denoting the number of model in the state estimation model set.② Mixed estimation is(35)xk-1∣k-10j=∑i=1mμk-1∣k-1i∣jxk-1∣k-1i,Pk-1∣k-10j=∑i=1mμk-1∣k-1i∣jPk-1∣k-1i+xk-1∣k-1i-xk-1∣k-10jxk-1∣k-1i-xk-1∣k-10jT,Xk-1∣k-10j=∑i=1mμk-1∣k-1i∣jXk-1∣k-1i,ek-1∣k-10j=∑i=1mμk-1∣k-1i∣jek-1∣k-1i+trXk-1∣k-1i-Xk-1∣k-10j2.(
2) Update. ① The prediction xk∣k-1j, Pk∣k-1j, Xk∣k-1j can be calculated by (12)~(13).② Likelihood function is(36)Λk∣k-1j=Nz-k;Hxk∣k-1j,Sk∣k-1j×2πVarΔk∣k-1i-d+1/4×etr-0.5Δk∣k-1iVarΔk∣k-1i-1Δk∣k-1i,where Δk∣k-1i=z-k-(nk-1)Yk∣k-1j(37)VarΔk∣k-1i≈nk-1trYk∣k-1jYk∣k-1j+Yk∣k-1j2+λ2nk-12Vk∣k-1j,Vk∣k-1j=VarXk∣k-1j∣Zk-1=αk∣k-1jtrXk∣k-1jXk∣k-1j+αk∣k-1j+2Xk∣k-1j2αk∣k-1j+1αk∣k-1j-2.Hereetr· is the abbreviation of exp⁡tr·.③ state and covariance update is as follows.The estimationsxk∣kj, Pk∣kj, Xk∣kj, ek∣kj, and αk∣kj can be calculated according to Section 2.3.Withek∣kj=trVk∣kj,Vk∣kj=VarXk∣kj∣Zk.④ Model probability update is as follows:(38)μkj=Λk∣k-1jc-j∑j=1mΛk∣k-1jc-j.(
3) Fusion Estimation
(39)
x
k
∣
k
=
∑
j
=
1
m
μ
k
j
x
k
j
,
P
k
∣
k
=
∑
j
=
1
m
μ
k
j
P
k
∣
k
i
+
x
k
∣
k
-
x
k
∣
k
j
x
k
∣
k
-
x
k
∣
k
j
T
,
X
k
∣
k
=
∑
j
=
1
m
μ
k
j
X
k
j
,
e
k
∣
k
=
∑
i
=
1
m
μ
k
j
e
k
∣
k
j
+
t
r
X
k
∣
k
-
X
k
∣
k
j
2
.


## 5. Simulation Results

Assume that the initial kinematic state of the group is-2km,0m/s,5km,-80m/sT; initial extended state Xk∣k=r1200r22=2020052m2. The parameter γ=20 denotes measurement ratio. The measure sustains 40 s, with sampling interval Ts=1s and measurement error covariance R=15200152. The group take a constant turning (CT) with Ω1=π/16rad/s, Ω2=π/8rad/s, and Ω3=-π/8rad/s being angular speeds of time periods 9~16 s, 17~24 s, and 25~32 s, respectively, and move with a constant velocity (CV) mode in other time periods. The track of the group is as shown in Figure 1, where the black spots denote the real measurements of the group target; the ellipse shows the confidence region of the group extension with confidence level being 0.9.Figure 1
Measurement of group target.Basic model set is composed of 13 CV models. Candidate model set is composed of 8 CT models. The dynamic model, model parameters, transform matrix, and the process noise matrix are the same as in [19]. Simulation results are shown in Figures 2–6 over 50 Monte Carlo runs.Figure 2
RMSE of extended state.Figure 3
RMSE of centroid position.Figure 4
RMSE of centroid velocity.Figure 5
ANEES of kinematics.Figure 6
ANEES of extension.In Figures2–6, “IMM13” denotes IMM algorithm with 13 basic models. “IMM21” denotes IMM algorithm with all the 21 models. “BMA13+1” denotes BMA algorithm with 13 basic models and 1 activated model. “BMA13+1+A” denotes the BMA algorithm with model set adaptively adjusted. “BMA13+1+New” denotes the proposed algorithm of this paper with the group extension predicted and BMA adapted.Figures2–4 show the root mean square error (RMSE) of centroid position, velocity, and extended state. It can be learned that, comparing IMM with BMA, a varied structure in model set can do a great help in improving the estimation performance in both kinematics and extension. The difference between IMM13 and BMA13+1 in extension estimation may be attributed to IMM’s divergence performance on kinematic state estimation in later stage, which has a terrible influence on innovation covariance prediction, moreover, worsening the estimation of extension.Moreover, it can be seen from Figures2–4 that IMM21 has an evident advantage in kinematic state estimation compared with BMA13+1+A but nearly has the same performance in extension estimation. It indicates that improving the kinematic performance to gain an improvement in extension estimation is limited. On the other hand, it indicates that the group extension is relatively independent of kinematic state in some situation, and this may explain the rationality of Assumption 1 to a certain extent.As to BMA13+1+A, it can be seen from Figures2–4 that an adaptive structure of BMA can improve the estimation performance of kinematic state obviously, but it does not work in extension estimation, which, on the one hand, indicates that the BMA adaption algorithm we derive in Section 3.2 is effective in improving the estimation performance of kinematic state; on the other hand, only BMA adaption does little help in improving the estimation performance of group extension; more needs to be done to improve the overall performance of group target tracking.The proposed algorithm BMA13+1+New, with extension predicted in orientation, can improve the estimation performance of extension obviously, especially when it is maneuvered; moreover, an adaptive structure of BMA improves the performance of its kinematic state estimation; thus the proposed algorithm can improve the overall performance of group target tracking compared with the original BMA13+1.Figures5 and 6 are average normalized estimation error squared (ANEES) of kinematic state and extension; they are credibility measure of the estimation results of these algorithms, while in the values larger than 1, the filter is overly confident about its estimation quality. The ANEES can be calculated as follows.① ANEES of kinematics is(40)ANEESx=1Mdimxk∑u=1Mdk∣ku,where the subscript u indicates tracking results concerning the uth run of a Monte Carlo simulation totaling M runs; dk∣k=xk∣k-xkTPk∣k-1xk∣k-xk.② ANEES of extension is(41)ANEESX=1M∑u=1MtrXk∣k-Xk2ek∣ku.It can be learned from Figures5 and 6 that the proposed algorithm has an excellent performance in both kinematic state and extension estimation with respect to confidence region.A comparison is presented in Table1, which contains average running time of a single step, average RMSE of centroid position, centroid velocity, and extended state. It can be learned that the proposed algorithm can do a good improvement in performance with a small computation cost in group target tracking.Table 1
A comparison on average performance.
IMM13
IMM21
IMM13+1
IMM13+1+A
IMM13+1+New
Average time (s)
0.0217
0.0493
0.0283
0.0324
0.0325
Position RMSE (m)
14.2303
5.5390
8.2438
6.4474
6.4699
Velocity RMSE (m/s)
19.9825
7.3435
12.9626
8.9902
8.8283
Extension RMSE (m2)
551.6187
214.2907
237.3733
214.2780
161.2098


## 6. Conclusion

This article focuses on the maneuvering group target tracking based on random matrix method. As the group target’s motion shows a clear orientation in many occasions, a random matrix approach is derived based on Bayesian theory for maneuvering target tracking, which does have a prediction on target extension, and the assumption that the extension keeps invariant is eliminated. Then, the model set of BMA is adjusted in real time and improves the algorithm’s tracking ability for target maneuvering. Finally, a comparison is presented by simulation, and the results shows that the proposed algorithm can gain an obvious improvement in tracking accuracy in both kinematic state and extension estimation with a relatively small computation cost. As the work of this paper only focuses on the group target with ellipse extension, the future work should pay more attention to group target with a more general extension appearance.



---

*Source: 5752870-2017-09-24.xml*