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-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 k and 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.

Figure 1.

Description by ellipses of the target's physical extension. A large ship (left) and an aircraft formation (right).

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 k and physical extension k, where k is an sd × 1 random vector and 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[xk,k|k] as a product of a Gaussian density and an inverted Wishart density based on Bayes' formula: where denotes the set of n measurements at scan k, n ≥ 1 is assumed to be independent of and , , 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 for EOT is very similar to the models utilized in many Kalman-related filters: where , “⊗” is the right Kronecker product [15], represents the state matrix in one spatial dimension, and 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 denotes the covariance matrix of k in one spatial dimension. and 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 is a d × d nonsingular real matrix describing physical extension evolution.

W(; a,), the Wishart density of SPD matrix , is defined by: where etr[·] is short for exp[trace(·)], a ≥ d, and expectation E[] = a.

IW(; a,), the inverted Wishart density of SPD matrix , is defined by: with a > 2d + 2 and expectation E[] = /(a − 2d − 2).

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 is modeled as follows: where and denotes the meas–urement matrix in one spatial dimension [1].

Gaussian measurement noise is independent of (j = 1,2,…,n, j ≠ r). Its distribution is assumed to be [4]: where 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 k.

Prediction

The prediction density of k and 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:

p[k|k−1] can be calculated based on Equations (4) and (5): where 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 via first and second moment matching [1,4]: where: and: with λk−1 = νk−1|k−1 − 2d − 2.

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[k|xk,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 is used to approximate .

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(k; k|k,k|k) with: where and are square roots of SPD matrices and , which can be obtained by Cholesky factorization. They are utilized to keep and 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,|] is converted to the product of a Gaussian density and an inverted Wishart density as follows:

The expectation of physical extension is employed as the final extension estimation: where α = ν − 2d – 2 > 0.

The Bayesian recursive estimator for EOT with random matrix (i.e., RMF) is composed of Equations (12), (15), (16), (22), (23), (25) and (27). A notable similarity exists between the standard Kalman filter and the centroid's kinematic part of RMF.

Model Parameter Adaptive Approaches

For a certain model used in multi-model (MM) algorithm, the model parameters include in Equation (3), δk and k in Equation (4)k in Equation (9), etc. The model parameter adaptive approaches act on extension dynamic model parameter k and measurement noise model parameter 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 denoting the coordinates of the points on the ellipse.

The following factors must be identified when determining an ellipse on the 2D plane.

Location: the coordinates of the ellipse's center point, which are determined by kinematic state ;

Size: the lengths of the semi-axes, which are equal to the positive square roots of 's eigenvalues;

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 . 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 is combined with kinematic state .

The following lemma is necessary to decompose SPD matrix into a form wherein the relationship between and the ellipse can be easily analyzed.

Lemma 1. If × d SPD matrix, then “Givens rotations” [ = 2, where is the rotation matrix, and are the eigenvalues of , , and with .

Therefore,

Proof of Lemma 1. Expand the left portion of Equation (29), then the following equation is obtained by setting the non-diagonal elements to 0: Let ρ denote ; thus: Q.E.D.

Decomposing SPD matrix 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 representing a 2D ellipse could be factorized as follows: where σ,1 and σ,2 are 's eigenvalues describing the ellipse's size, and 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 δ exists on the right side of Equation (35), which is different from Equation (11) in [4]. However, employing Equation (35) allows to cover all the possible changes in the size and orientation of the ellipses represented by SPD matrices and -1. According to Equation (15), δ is still incorporated in the Bayesian RMF estimator, describing the uncertainty of physical extension evolution [1]. Therefore, the improvement of Equation (4) is reasonable and meaningful.

For parameter adaption when d = 2, we propose to calculate by the following equation: where describes the change in ellipse's orientation and are scalars describing the change in size.

Although the range of is reduced, defined by Equation (36) is still able to approximately describe all changes in the ellipses' size and orientation.

According to Equations (35) and (36), describes the change amount of orientation angles of the two ellipses represented by and k−1 respectively. Then could be set to Δφ = φ − φ−1, where φ and φ−1 are orientation angles described in Equation (34). Let φ() and φ(−1) denote the angles calculated by Lemma 1. If the difference between φ(·) and φ is a constant (i.e., φ() − φ(−1) = Δφ), then could be set to φ() − φ(k−1). However, no matter “+” or “−” in Equation (30) is selected, the difference between φ(·) and φ is not a constant, i.e., φ() − φ(k−1) ≠ Δφ. 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 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 φ(·) differs from φ by a constant during the whole EOT process, i.e., φ() − φ(−1) = Δφ). This property is verified by the simulation results.

However, and −1 are unknown and is also unavailable at the beginning of period k, thus we use −1 and −2 to approximate and −1. is then calculated according to the change of orientation angles of the ellipse represented by −1 and −2

From the above conditions, can be calculated by:

The next task is to determine can be determined by the ratios of the corresponding eigenvalues of −1 and −2. When the size of the object does not change rapidly and sharply, the following calculation method is suggested. where tr(·) is short for trace(·).

Therefore, can be adaptively determined according to Equations (36), (39) and (40).

Measurement Noise Model Parameter Adaption

Let denote the real covariance matrix of measurement noise . According to Equation (9), is used to approximate . Thus, the basic idea of adaption is adjusting to meet = as accurately as possible. Given that and are unknown and is also unavailable at the beginning of period k, is calculated according to the difference between −1(k > 1) and , where is an approximation of . Let denote the noise covariance matrix of measurement sensor, which is usually known as prior information. Then is calculated by [3]: where η is a scaling coefficient, |−1 is utilized to approximate , and according to Equation (14).

Similar to Equation (36) is calculated by:

and in Equation (42) can be approximately calculated in a manner similar to and calculation in Equations (39) and (40):

Model Parameter Adaption When d = 3

For 3D EOT (d = 3), 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 : 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 , the order of rotation matrices should be the same as that of Equation (45). The adaptive approach of when d = 3 is similar to that of .

Multi-Model Algorithm

The proposed ARMF approach is integrated into a multi-model algorithm framework [18] (e.g., the widely used interacting multi-model (IMM) algorithm [19]) to further improve the estimation performance.

Moment Matching

Moment matching is applied in step initialization and fusion of IMM. Unlike the standard IMM algorithm [18], moment matching in ARMF-MM algorithm involves both kinematic state and physical extension ; this condition makes the moment matching method more complex.

Step fusion is regarded as an example in this section to introduce moment matching method. In the first place, the substitution of Equation (27) into Equation (26) yields:

Moment matching in step fusion requires that the following approximation be solved: where j = 1,2….,n denotes model j and n is the number of models.

The moment matching method for kinematic state is similar to standard IMM: where and indicates that model j is in effect at scan k.

The problem of moment matching method for physical extension is that the method considers and in Equation (25) simultaneously. The mean square error ê of extension estimation is first defined as follows [3]:

Expand Equation (49), and the following equation is derived [3,20]: where:

Detailed derivation and explanation can be found in [3,20]. Then the moment matching method for physical extension is as follows:

After and ê are determined, is calculated based on Equations (50) and (51). Thus, is obtained as follows:

ARMF-Based Multi-Model Algorithm

The four steps of ARMF-IMM algorithm are summarized below:

Reinitialization

Let π denote the probability of transiting from model i to model j. is the probability that model i is in effect at scan k − 1 given and :

Reinitialization of each model is achieved by replacing corresponding probability variables in Equations (48)–(53):

Model Parameter Adaption and Filtering

For model j, j = 1,2….,n, calculate according to Equations (36), (39) and (40). Then the RMF approach as summarized in Section 2 is performed and is calculated according to Equations (42)–(44) for model j.

Probability Update

According to Equation (54), probability is calculated by:

Likelihood is derived as: where is calculated by Equation (17) and by Equation (18).

The substitution of Equations (17) and (18) into Equation (57) approximately yields: where is used to approximate in , and Γd[·] is the multivariate gamma function. The derivation of is quite similar to that of Equations (11) and (13). When n = 1, .

Fusion

The moment matching formulas of fusion are provided by Equations (48)–(53). The final result of extension estimation is defined in Equation (27).

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 Ψ = diag([502,202]) m2.

Figure 2.

Target trajectory of simulation scenario 1.

The single-model ARMF is compared with Koch's approach [1] and Feldmann's approach [3] through NM = 500 Monte Carlo runs.

Koch's approach: the temporal decay constant τ is set to 8T;

Feldmann's approach: 0.25 + is utilized as the measurement noise covariance matrix; and

ARMF: η in Equation (42) is set to 0.25 via moment matching [3].

Average , 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 according to the extension's changes.

Figure 3.

Average in simulation scenario 1.

The root mean square error (RMSE) of physical extension is defined by [3]: where 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.

Figure 4.

RMSE of physical extension in simulation scenario 1.

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.

Figure 5.

Target trajectory of simulation scenario 2.

The following three approaches are compared through NM = 500 Monte Carlo runs.

Koch's approach: parameters are the same as those in scenario 1;

Feldmann's MM approach: parameters are designed similar to those in [3]; and

ARMF-MM: three models are used. Model 1 with low kinematic process noise and low extension agility (small Q, large δ), model 2 with high kinematic process noise and high extension agility (large Q, small δ), and model 3 with moderate kinematic process noise and high extension agility (moderate Q, small δ) are utilized. η=0.25 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.

Figure 6.

Average of model 3 in simulation scenario 2.

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.

Figure 7.

RMSE of position in simulation scenario 2.

Figure 8.

RMSE of velocity in simulation scenario 2.

Figure 9.

RMSE of physical extension in simulation scenario 2.

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.