Model Parameter Adaption-Based Multi-Model Algorithm for Extended Object Tracking Using a Random Matrix

Traditional object tracking technology usually regards the target as a point source object. However, this approximation is no longer appropriate for tracking extended objects such as large targets and closely spaced group objects. Bayesian extended object tracking (EOT) using a random symmetrical positive definite (SPD) matrix is a very effective method to jointly estimate the kinematic state and physical extension of the target. The key issue in the application of this random matrix-based EOT approach is to model the physical extension and measurement noise accurately. Model parameter adaptive approaches for both extension dynamic and measurement noise are proposed in this study based on the properties of the SPD matrix to improve the performance of extension estimation. An interacting multi-model algorithm based on model parameter adaptive filter using random matrix is also presented. Simulation results demonstrate the effectiveness of the proposed adaptive approaches and multi-model algorithm. The estimation performance of physical extension is better than the other algorithms, especially when the target maneuvers. The kinematic state estimation error is lower than the others as well.


Introduction
During the past several decades, the resolution of radar and other sensors continues to increase, thereby generating complex requirements for object tracking (OT) technology [1]. Traditional OT technology characterizes the target as a point source object whose physical extension is ignored [2]. The centroid is usually utilized to represent the target, and its kinematic state such as position, velocity, and acceleration is estimated in OT. This simplification is meaningful only when physical extension can be neglected compared with the sensors' measurement error. However, considering the increasing sensor resolution, physical extension state should be estimated jointly in OT approaches, especially when tracking relatively large targets such as aircraft carriers and large transport planes. This method is referred to as extended object tracking (EOT) [1,[3][4][5]. When targets form a closely spaced formation (e.g., an aircraft formation), the individual target becomes indistinguishable from the formation because of limited sensor capabilities and sensor-to-target geometry. Thus, closely spaced group objects can also be considered extended objects [1,3]. Data association in EOT is difficult because extended objects produce a highly varying number of measurements corresponding to a single extended object [3].
A number of techniques and approaches have been proposed to solve the EOT problem, which requires estimating physical extension and the centroid's kinematic state jointly. Baum et al. [6] proposed a novel modeling approach to describe physical extension. The proposed model was called random hypersurface model. A recursive Bayesian estimator was derived based on this model. Mahler [7,8] applied probability hypothesis density (PHD) filter and cardinalized PHD (CPHD) filter to EOT, including both extended and unresolved targets. Vo et al. [9] presented a mathematically strict Bayesian filter for EOT, which could be reduced to a CPHD filter under some assumptions when tracking a single target. Carmi et al. [10] utilized the Gaussian mixture model to describe the time-varying dynamics of group objects. A Markov chain Monte Carlo particle filter for multi-target tracking was then obtained. Multiple hypotheses tracking was also a feasible option for EOT [11]. Zhang and Bar-Shalom [12] proposed a new framework based on multiple kernel centers for visual object tracking.
Koch [1] proposed the random matrix-based approach for extended object tracking. A recursive form filter is derived within the Bayesian framework, and the final form of this random matrix-based filter (RMF) is more or less similar to that of the standard Kalman filter. RMF describes the physical extension of the target with an ellipse or ellipsoid. Thus, RMF is reasonable and can be applied in many practical applications as shown in Figure 1. Given that symmetrical positive definite (SPD) matrices can represent all the ellipses/ellipsoids centering at the origin, RMF utilizes the random SPD matrix to characterize the physical extension of the target. The random SPD matrix is assumed to obey Wishart-related distributions. Then the Bayesian formula can be applied to obtain a very simple filter according to the properties of the Kronecker product and Gaussian and Wishart-related distributions. Feldmann [3] introduced a more accurate distribution of measurement noise with respect to the target's physical extension and derived a new RMF approach. This new approach can be integrated into an interacting multi-model (IMM) framework to improve estimation performance if the object state switches between maneuvering and non-maneuvering. However, the derivation process cannot be maintained within the Bayesian framework, and several approximations are applied to obtain a recursive form filter. Lan and Li [4,5] further improved extension dynamic model and measurement noise model by describing extension evolution and observation distortion with real matrices A k and B k . This improvement is logical given the relationship between the SPD matrix and ellipse/ellipsoid. It can be applied within the Bayesian framework, which requires only an insignificant approximation. RMF approach is much simpler than the other EOT approaches mentioned above, thereby making it more promising. The simulation results demonstrate the effectiveness of the three abovementioned RMF approaches. However, the physical extension estimation error obtained by these approaches significantly increases when the target maneuvers (e.g., turning motion). The extension dynamic model in [4] is improved in this study based on RMF approach to simplify the description of extension evolution. Model parameter adaptive approaches for both extension dynamic and measurement noise are proposed based on properties of SPD matrix and ellipse/ellipsoid. The proposed adaptive RMF is called ARMF. The preliminary results of single-model ARMF have been published in a previous conference paper [13]. A multi-model algorithm is proposed by integrating ARMF into the widely utilized IMM framework to further improve the performance of RMF when the target maneuvers. Simulation results show that the ARMF approach and ARMF-based multi-model (ARMF-MM) algorithm are effective.

Bayesian EOT Using Random Matrix
The dynamic state of an extended object at scan k is described by both the centroid's kinematic state x k and physical extension X k , where x k is an sd × 1 random vector and X k is a d × d SPD random matrix. Spatial dimension d is usually set to 2 or 3 in EOT. s is the dimension of the state in one spatial dimension. s = 2 means that position and velocity are considered. If acceleration is also considered, then s = 3.
The basic idea of random matrix-based EOT is to characterize joint density p[x k ,X k Z k ] as a product of a Gaussian density and an inverted Wishart density based on Bayes' formula: where denotes the set of n k measurements at scan k, n k  1 is assumed to be independent of x k and X k , , and c always denotes the normalization factor in this paper.
This section summarizes the key results and steps. More detailed assumptions and derivations can be found in [1,3,4,14].

Dynamic and Measurement Models
The dynamic model of x k for EOT is very similar to the models utilized in many Kalman-related filters: where k k d  F F I , "" is the right Kronecker product [15], k F represents the state matrix in one spatial dimension, and I d is the d × d identity matrix. Independent process noise  k follows a normal distribution [1]: where N(μ,Σ) stands for normal distribution with mean μ and covariance matrix Σ, and k D denotes the covariance matrix of  k in one spatial dimension. k F and k D can be selected according to specific conditions of EOT. Singer's model [1], constant velocity model, and constant acceleration model [16] are all feasible alternatives.
The following assumptions are introduced for physical extension [1,4]: where δ k > d − 1 stands for the degrees of freedom,  k−1k−1 > 2d is a scalar parameter, subscript k−1k−1 indicates the estimated value of corresponding variable at scan k − 1, and A k is a d × d nonsingular real matrix describing physical extension evolution. W(X; a,A), the Wishart density of SPD matrix X, is defined by: where etr[•] is short for exp[trace(•)], a  d, and expectation E[X] = aA. IW(X; a,A), the inverted Wishart density of SPD matrix X, is defined by: with a > 2d + 2 and expectation E The assumption in Equation (4) is different from that in Equation (10) in [4]. It is simpler and more plausible, which will be explained in Section 3.1.
Measurement z r k is modeled as follows: denotes the meas-urement matrix in one spatial dimension [1].
Gaussian measurement noise  r k is independent of  j k (j = 1,2,…,n k , j  r). Its distribution is assumed to be [4]: where B k is a d × d nonsingular real matrix describing the measurement distortion of extension. Equation (9) indicates that the measurements are affected by physical extension X k .

Prediction
The prediction density of x k and X k can be factorized as: The first factor on the right of the previous equation is assumed to have the following structure [1]: where subscript kk−1 indicates the predicted value of corresponding variable at scan k .
We then obtain the prediction equations of kinematic state: can be calculated based on Equations (4) and (5): where GB II d (•) is the probability density function (pdf) of "generalized beta type II" distribution. An inverted Wishart distribution is utilized to obtain a recursive-form estimator and approximate GB II d (•) via first and second moment matching [1,4]: where: and: with The substitution of Equations (11) and (14) into Equation (10) yields:

Update
In consideration of Equation (1) and based on Equations (8) and (9), likelihood function p[Z k x k ,X k ] is obtained as follows [4]: where: The substitution of Equations (17) and (18) into Equation (1) yields [1,4]: The product of the two Gaussian probability density functions in the previous equation can be converted to [1,4]: with: and: where X kk−1 is used to approximate X k . Based on Equations (21) and (20) can be rewritten as [1,4]: ; Subsequently, the first three pdfs on the right side of the previous equation can be combined as IW(X k ; ν kk ,X kk ) with: S , which can be obtained by Cholesky factorization. They are utilized to keep N kk−1 and X kk in positive definite or positive semi-definite structure [3]. The derivation of Equation (25) is similar to the corresponding process in [1].
Then, p[x k ,X k Z k ] is converted to the product of a Gaussian density and an inverted Wishart density as follows: The expectation of physical extension X k is employed as the final extension estimation: The Bayesian recursive estimator for EOT with random matrix (i.e., RMF) is composed of Equations (12)

Model Parameter Adaptive Approaches
For a certain model used in multi-model (MM) algorithm, the model parameters include k D in Equation (3), δ k and A k in Equation (4), B k in Equation (9), etc. The model parameter adaptive approaches act on extension dynamic model parameter A k and measurement noise model parameter B k based on the relationship between the SPD matrix and ellipse/ellipsoid. The basic principle of parameter adaptive approach when d = 2 and d = 3 is the same. d = 2 is regarded as an example to introduce the parameter adaptive approaches.
The following assignments are adopted in the sequel for the purpose of convenience: d = 2, s = 3, and , which means that the kinematic state in one spatial dimension is [position, velocity, acceleration] T and only the position coordinates of the target's centroid are measured. Therefore, the ellipse describing physical extension can be represented by [4]: with x denoting the coordinates of the points on the ellipse. The following factors must be identified when determining an ellipse on the 2D plane.
(1) Location: the coordinates of the ellipse's center point, which are determined by kinematic state x k ; (2) Size: the lengths of the semi-axes, which are equal to the positive square roots of X k 's eigenvalues; (3) Orientation: the direction of the semi-axes, which is determined by the angle between either of the semi-axes and either of the coordinate axes when d = 2. Several definitions of orientation angles are available. For example, the angle between the major semi-axis and x axis is suitable when the counterclockwise direction is considered the positive direction.
Size and orientation are determined entirely by physical extension X k . As long as the size or orientation of an ellipse changes, the SPD matrix representing the ellipse would change. This occurrence means that the SPD matrices represent only the ellipses centering at the origin. However, all the ellipses centering at any position on the 2D plane will be covered if X k is combined with kinematic state x k .
The following lemma is necessary to decompose SPD matrix X k into a form wherein the relationship between X k and the ellipse can be easily analyzed. ( 1) 2 dd "Givens rotations" [17]. When d = 2, X can be transformed as follows:

Q.E.D.
Decomposing SPD matrix X into the product of orthogonal matrices and a diagonal matrix is relatively easy in mathematics. Except for "Givens rotations," singular value and eigenvalue decomposition are also possible. However, the orthogonal matrix contains not only rotation transformation but also symmetry transformation (rotation matrix is the orthogonal matrix whose determinant is +1). The form and sequence of symmetry transformation are uncertain. Thus, neither of the two decomposition methods is applicable for model parameter adaption.
According to Lemma 1, SPD matrix X k representing a 2D ellipse could be factorized as follows: where ,1 k  and ,2 k  are X k 's eigenvalues describing the ellipse's size, and , 22 is the orientation angle describing the ellipse's orientation.

Extension Dynamic Model Parameter Adaption
From the extension dynamic model in Equation (4), we obtain: No δ k exists on the right side of Equation (35), which is different from Equation (11) in [4]. However, employing Equation (35) allows A k to cover all the possible changes in the size and orientation of the ellipses represented by SPD matrices X k and X k-1 . According to Equation (15), δ k is still incorporated in the Bayesian RMF estimator, describing the uncertainty of physical extension evolution [1]. Therefore, the improvement of Equation (4) where φ k and φ k−1 are orientation angles described in Equation (34). Let φ G (X k ) and φ G (X k−1 ) denote the angles calculated by Lemma 1. If the difference between φ G (•) and k  is a constant . However, no matter "+" or "−" in Equation (30) is selected, the difference between φ G (•) and φ k is not a constant, i.e., φ G (X k ) − φ G (X k−1 )  ∆φ k . Thus, the following two restrictions are introduced. First of all, at the very beginning of EOT, the initial value of extension state should be selected by: According to the previous discussion and Equation (34), 0 X should be able to cover all the possible SPD matrices representing the ellipses centering at the origin. Secondly, we select "+" or "−" in Equation (30) to satisfy the following condition: where  1 and  2 are the same as those in Equation (29), and sign(•) is the sign function.
If these two restrictions are met, then φ G (•) differs from φ k by a constant during the whole EOT process, i.e., φ G (X k ) − φ G (X k−1 ) = ∆φ k ). This property is verified by the simulation results.
However, X k and X k−1 are unknown and k X is also unavailable at the beginning of period k, thus we use to approximate X k and X k−1 . k  A is then calculated according to the change of orientation angles of the ellipse represented by From the above conditions, k  A can be calculated by: The next task is to determine . can be determined by the ratios of the corresponding eigenvalues of where tr(•) is short for trace(•). Therefore, A k can be adaptively determined according to Equations (36), (39) and (40).

Measurement Noise Model Parameter Adaption
Let  v k denote the real covariance matrix of measurement noise ν r k . According to Equation (9), k as accurately as possible. Given that X k and  v k are unknown and k X is also unavailable at the beginning of period k, B k is calculated according to the difference between 1 k  X (k > 1) and  * k , where  * k is an approximation of  v k . Let Ψ k denote the noise covariance matrix of measurement sensor, which is usually known as prior information. Then  * k is calculated by [3]: where η is a scaling coefficient,

Model Parameter Adaption When d = 3
For 3D EOT (d = 3), X k can be transformed to a diagonal matrix by up to three "Givens rotations", which means that at most three angles are required to determine a 3 × 3 SPD matrix [17]. Thus, the following equation is introduced to calculate A k : where is the 3 × 3 rotation matrix. The calculation of and in Equation (45) is similar to that in 2D EOT. However, it should be noted that the three rotation matrices in Equation (45) must be maintained in the same order during the whole EOT process. In addition, when utilizing Lemma 1 to decompose X k , the order of rotation matrices should be the same as that of Equation (45). The adaptive approach of B k when d = 3 is similar to that of A k .
where and j k m indicates that model j is in effect at scan k.

Simulation Results and Discussions
Two typical EOT scenarios are simulated in this section to compare the proposed ARMF approach and other RMF approaches.

Scenario 1
The route of an aircraft carrier, which is approximately 300 m long and 80 m wide, is shown in Figure 2. The route starts from origin (0, 0) T , moves along the trajectory, and then maintains a uniform circular motion on the 2D plane. Constant velocity is set to 27 knots (approximately 50 km/h). The aircraft carrier's extension is characterized by an ellipse with radiuses 170 m and 40 m in RMF approach. Scan period T is set to 10 s. The measurement scattering centers are assumed to be distributed uniformly over the physical extension (i.e., the ellipse) [3], and the number of measurements in every scan obeys Poisson distribution with mean 20. The sensor measurement noise follows a zero-mean Gaussian distribution with covariance matrix  The single-model ARMF is compared with Koch's approach [1] and Feldmann's approach [3] through N M = 500 Monte Carlo runs.
Average k  A , which describes the change in the ellipse's orientation angle, is shown in Figure 3.
The curve shows that the model parameter adaptive approach is effective and that ARMF is able to adjust k  A according to the extension's changes.
The root mean square error (RMSE) of physical extension X k is defined by [3]: where l k X stands for the expectation of physical extension in lth Monte Carlo run.  The RMSE of physical extension of the three RMF approaches is shown in Figure 4. When the object makes a turn (uniform circular motion after scan 10), the extension RMSE of ARMF becomes lower than that of the other two approaches.

Scenario 2
The extended object scenario in [3] is simulated below. The aircraft carrier whose size is the same as the one in scenario 1 starts from the origin and moves along the trajectory on the 2D plane as shown in Figure 5. Three turns are implemented.
The following three approaches are compared through N M = 500 Monte Carlo runs.
(1) Koch's approach: parameters are the same as those in scenario 1;  for all three models. The simulation results reconfirmed the validity of the model parameter adaptive approach. ARMF-MM identified the angle changes correctly when the object maneuvers as shown in Figure 6. The position RMSEs of three RMFs are close to one another. The error curve is shown in Figure 7. The position RMSE of ARMF-MM is lower than that of the others. The RMSEs of velocity and physical extension are shown in Figures 8 and 9. ARMF-MM is able to maintain a low level of velocity RMSE in the entire process and performs better than the others in physical extension estimation. The physical extension RMSE of ARMF-MM is significantly reduced when the target maneuvers.

Conclusions
Bayesian EOT using a random matrix is a very effective and simple approach to incorporate physical extension in the OT framework. The Bayesian RMF approach was reviewed first in this paper, and certain improvements were made. Afterward, model parameter adaptive approaches were derived for both extension dynamic evolution and measurement noise based on the analysis of SPD matrix decomposition. The proposed adaptive RMF can be integrated into an IMM framework to further improve estimation performance. The validity of ARMF and ARMF-MM was verified by the simulation results. ARMF is able to change model parameters adaptively when the extended object maneuvers, thereby resulting in low physical extension estimation error.