Tracking of Multiple Maneuvering Random Hypersurface Extended Objects Using High Resolution Sensors

Abstract

With the increased resolution capability of modern sensors, an object should be considered as extended if the target extent is larger than the sensor resolution. Multiple maneuvering extended object tracking (MMEOT) uses not only measurements of the target centroid but also high-resolution sensor measurements which may resolve individual features or measurement sources. MMEOT aims to jointly estimate object number, centroid states, and extension states. However, unknown and time-varying maneuvers of multiple objects produce difficulties in terms of accurate estimation. For multiple maneuvering star-convex extended objects using random hypersurface models (RHMs) in particular, their complex maneuvering behaviors are difficult to be described accurately and handled effectively. To deal with these problems, this paper proposes an interacting multiple model Gaussian mixture probability hypothesis density (IMM-GMPHD) filter for multiple maneuvering extended object tracking. In this filter, linear maneuver models derived from RHMs are utilized to describe different turn maneuvers of star-convex extended objects accurately. Based on these, an IMM-GMPHD filtering recursive form is given by deriving new update and merging formulas of model probabilities for extended objects. Gaussian mixture components of different posterior intensities are also pruned and merged accurately. More importantly, the geometrical significance of object extension states is fully considered and exploited in this filter. This contributes to the accurate estimation of object extensions. Simulation results demonstrate the effectiveness of the proposed tracking approach—it can obtain the joint estimation of object number, kinematic states, and object extensions in complex maneuvering scenarios.

1. Introduction

In recent years, extended object tracking has attracted more attention and has been applied widely in many practical situations [1,2,3,4] (e.g., airplane tracking, vehicle tracking, and so on). Conventional multiple target tracking approaches [5,6] usually model targets as point targets. A single point target is assumed to generate at most one measurement. The combinatorial nature between targets and measurements is followed to track multiple targets. However, data associations between measurements and targets usually cause intractable computation loads for these approaches. To deal with this problem, a probability hypothesis density (PHD) filter [7,8] and its variations [9,10,11,12] were proposed based on the finite set statistics (FISST) theory. In these filters, targets’ states and measurements are modeled using random finite sets (RFS), which eliminates the complex calculation of data associations.

With the development of high-resolution sensors, a target may have multiple scattering centers (also called measurement sources). Correspondingly, multiple measurements may be taken from these scattering centers. Such a target is defined as an extended object [13,14], and vivid information (e.g., object shape and size) of an extended object can be obtained from its measurements. For multiple extended object tracking (MEOT) approaches, there are two main categories. The first category simplifies the tracking problem of MEOT by ignoring the geometrical shapes of extended objects. In these approaches, the set of measurements from extended objects is a cluster without any geometric significance. Specifically, the measurement number of an extended object is assumed to obey a Poisson distribution [15] and the distribution of these measurements is determined by a so-called inhomogeneous Poisson point process measurement model [16]. Considering these assumptions, many MEOT approaches were proposed in the PHD framework, such as an extended target probability hypothesis density (ET-PHD) filter [17] and an extended target Gaussian mixture probability hypothesis density (ET-GM-PHD) filter [18]. It is noted that these filters only focus on estimating object number and kinematic states without considering the object extension estimation.

Another category of MEOT approaches models object extensions to estimate extension-involved states (e.g., objects’ shapes, orientations, and sizes). One of the extended object modeling approaches is the random matrix approach [19,20,21,22], in which elliptic shape is described as a random symmetric positive definite matrix. By adopting this approach, a Gaussian-inverse Wishart probability hypothesis density (GIW-PHD) filter and a gamma Gaussian-inverse Wishart probability hypothesis density (GGIW-PHD) filter were proposed for tracking multiple elliptic extended objects. Nevertheless, these filters could not handle the tracking of extended objects with extensions more complex than ellipses. To solve this problem, a so-called random hypersurface model (RHM) [23,24] was proposed, which can track not only simple elliptic objects but also extended objects with complex extensions (e.g., star-convex shapes). Based on RHMs, a variation [25] of the ET-GM-PHD filter was proposed for tracking multiple star-convex extended objects. However, its merging step may not fully consider the geometric significance of extension-involved states, which may lead to the degradation of extension estimate performance. Additionally, the maneuvering behaviors of multiple star-convex extended objects may not be solved effectively in this filter.

Considering RHM’s flexibility to describe complex extensions, this paper focuses on multiple maneuvering star-convex extended object tracking (MMSEOT) using RHMs. For MMSEOT using RHMs, complex maneuvering behaviors of star-convex extended objects are difficult to be described and managed accurately. The maneuvering behaviors of a star-convex extended object using RHM are hard to be modeled by a simple mathematical form due to the tight couple between the dynamic evolution of centroid-involved kinematical parameters and that of extension-involved parameters. What is worse, maneuver models of multiple star-convex extended objects are typically unknown and time-varying, which brings significant challenges for MMSEOT. Though the interacting multiple model (IMM) approach and its variations [26,27,28,29,30,31] seem promising for handling the unknown maneuver models of objects, the model probability update step of these approaches is unfit for MMSEOT. Moreover, these approaches do not consider the estimation of object extensions.

To deal with above problems of MMSEOT using RHMs, an IMM-GMPHD filter is proposed to jointly estimate the object number, centroid dynamics, and extensions of maneuvering star-convex extended objects. To describe different turn maneuvers (with distinct turn rates and directions) of star-convex extended objects, linear maneuver models are integrated to RHMs. These maneuver models can effectively handle the tight couple between the evolution of centroid-involved kinematical states and that of extension-involved states. Thanks to their simple forms, they are easily applied to the IMM-GMPHD filter. To handle the uncertainty of maneuver models, a complete filtering recursion of the proposed filter is formed by adjusting the filtering procedures of both the ET-GM-PHD filter and IMM approach. It strictly derives update formulas of model probabilities for extended objects. More importantly, a new merging and pruning strategy of model probabilities is proposed, which greatly contributes to obtaining the accurate posterior intensity function of extended objects. Furthermore, this filter also proposes a merging and pruning strategy of Gaussian mixture components. The merging of means in Gaussian mixture components is divided into two parts, i.e., the merging of centroid kinematical parameters and the merging of extension-involved parameters. For extension-involved parameters in particular, they cannot be directly merged owing to their geometric meanings. Thus, the merging of extension-involved parameters is turned to that of values of radial functions, which supports achieving more accurate extension estimations to a great extent. Simulation results illustrated that the proposed filter obtains better estimation performance in centroid dynamics and object extensions compared with the traditional multiple extended object tracking approach based on RHMs, especially when objects maneuver simultaneously. The innovations of this work can be further summarized as follows.

(1) The maneuver models of star-convex extended objects can accurately describe turn maneuvers with different turn rates and be easily implemented in the IMM-GMPHD filter due to their concise and linear mathematical forms;

(2) The update and merging formulas of model probabilities are strictly derived in the proposed filter, which facilitates the accurate tracking of multiple extended objects with unknown and time-varying maneuvering behaviors;

(3) The geometrical significance of extension-involved states is fully considered in the merging process of Gaussian mixture components to improve the extension estimation performance;

(4) By propagating the intensity function of multiple extended objects, the proposed IMM-GMPHD filtering recursion achieves the joint estimation of object number, centroid states, and object extensions accurately in complex maneuvering tracking scenarios.

The structure of this paper is as follows: Section 2 primarily analyzes problems in multiple maneuvering extended objects tracking. Section 3 proposes an IMM-GMPHD filter for MMSEOT. This section minutely mentions maneuver models for star-convex extended objects, update formulas of model probabilities for maneuvering extended objects, and specific filtering steps of the proposed filter. In Section 4, simulation results are presented to highlight the effectiveness of the proposed filter. The last section concludes this paper.

2. Problem Formulation

In MMSEOT, the state vector ${\xi }_k$ of a maneuvering extended object contains not only dynamic parameters $x_{k,m}$ (e.g., positions ${[x_k,y_k]}^{\mathrm{T}}$ and velocities ${[{\dot{x}}_k,{\dot{y}}_k]}^{\mathrm{T}}$ in a two-dimensional Cartesian coordinate system) about the centroid but also parameters $x_{k,e}$ involving the object extension (e.g., the shape and orientation). It is crucial to note that objects with different extensions (e.g., ellipses, circles, and star shapes) solely have different $x_{k,e}$, i.e., $x_{k,e}$ totally depends on the specific spatial extension of an extended object. Here, we focus on maneuvering star-convex extended objects using RHMs.

2.1. System Models of A Star-Convex Extended Object

Dynamic Model: As mentioned above, the whole state for a star-convex extended object is denoted by

$${\xi }_k={\left[\begin{array}{cc}{x}_{k,e}& x_{k,m}\end{array}\right]}^{\mathrm{T}}$$

(1)

The dynamic equation of an object state justifies the evolution of centroid-involved kinematical variables $x_{k,m}$ and that of extension-involved variables $x_{k,e}$ simultaneously, which has the following form:

$${\xi }_{k+1}=\left[\begin{array}{cc}{F}_{k,e}& 0\\ 0& F_{k,m}\end{array}\right]\left[\begin{array}{c}{x}_{k,e}\\ x_{k,m}\end{array}\right]+\left[\begin{array}{c}{w}_{k,e}\\ w_{k,m}\end{array}\right]$$

(2)

where $F_{k,e}$ and $F_{k,m}$ denote transition matrices for $x_{k,e}$ and $x_{k,m}$, respectively, and $w_{k,e}$ and $w_{k,m}$ represent process noises for $x_{k,e}$ and $x_{k,m}$, respectively.

It is significant that when a star-shape extended object maneuvers with a constant turn motion, mathematical representations of $F_{k,e}$ and $F_{k,m}$ are complex and difficult to calculate. This is because the rotation of the object extension simultaneously occurs with that of the centroid, resulting in the tight couple between the evolution of $x_{k,e}$ and $x_{k,m}$. Moreover, it may be hard to describe different maneuvering behaviors of the object by a uniform form of $F_{k,e}$ and $F_{k,m}$.

Measurement Model: The measurement model is used to describe the sensor measurement generation process of an object. An extended object with a specific extension may have multiple scattering points (also identified as measurement sources) at per scan of a high-resolution sensor. In this situation, the measurement model of an extended object is usually composed of two parts:

(1) Extension Model: This specifies the locations of measurement sources.

(2) Sensor Model: This produces measurements from measurement sources via a sensor noise, i.e.,

$$z_k^{\left(i\right)}=y_k^{\left(i\right)}+v_k^{\left(i\right)}.$$

(3)

where $y_k^{\left(i\right)}$ is the $i^{th}$ measurement source, $z_k^{\left(i\right)}$ is the measurement obtained by $y_k^{\left(i\right)}$, and $v_k^{\left(i\right)}$ is the sensor noise.

For a star-convex extended object, the entire measurement generation process is illustrated in Figure 1. Here the extension of a star-convex object is modeled by the RHM. The object shape boundary is denoted by $\mathcal{S}\left({\xi }_k\right)$. The measurement source lies on a random scaled contour $\overline{\mathcal{S}}\left({\xi }_k\right)$ of this extended object. The scaled version of $\mathcal{S}\left({\xi }_k\right)$ for a star-convex object can be calculated by a scaling factor $s_k$ (modeled as a one-dimensional random variable). If $s_k=1,$ the scaled boundary is exactly equal to the object shape boundary $\mathcal{S}\left({\xi }_k\right)$. Thus, using scaled contours aids in successfully covering the measurement sources lying in the interior of the boundary as well as on the boundary.

With the help of a radial function $r\left(\theta \right)$ (which gives the distance between the centroid $x_{k,c}$ and a contour point), the polar representation of the scaled contour $\overline{\mathcal{S}}\left({\xi }_k\right)$ is derived by

$$\overline{\mathcal{S}}\left({\xi }_k\right)=\left\{s_{k}r\left(\theta \right)e\left(\theta \right)+x_{k,c}\left|\theta \in \left[0,2\pi \right],s_k\in \left[0,1\right]\right.\right\}$$

(4)

where $e\left(\theta \right)=\left[\begin{array}{c}cos\left(\theta \right)\\ sin\left(\theta \right)\end{array}\right]$, $\theta$ represents the angle between the vector from the centroid to a contour point and the x-axis.

Using Fourier series expansion of degree $N_F$, the radial function $r\left(\theta \right)$ has the following finite dimensional parameterizable form:

$$\begin{array}{cc}\hfill r\left(B_k,\theta \right)& =\frac{a_k^{\left(0\right)}}{2}+\sum _{j=1\cdots N_F}\left(a_k^{\left(j\right)}cos\left(j\theta \right)+b_k^{\left(j\right)}sin\left(j\theta \right)\right)\hfill \\ \hfill & =R\left(\theta \right)B_k\hfill \end{array}$$

(5)

where

$$R\left(\theta \right)=\left[1,cos\left(\theta \right),sin\left(\theta \right),\cdots ,cos\left(N_F\theta \right),sin\left(N_F\theta \right)\right]$$

(6)

$$B_k=\left[a_k^{\left(0\right)},a_k^{\left(1\right)},b_k^{\left(1\right)},\cdots ,a_k^{\left(N_F\right)},b_k^{\left(N_F\right)}\right]$$

(7)

Here, $B_k$ is considered to be extended-involved parameters $x_{k,e}$ in the state ${\xi }_k$, i.e.,

$$\begin{array}{cc}\hfill x_{k,e}& =\left[a_k^{\left(0\right)},a_k^{\left(1\right)},b_k^{\left(1\right)},\cdots ,a_k^{\left(N_F\right)},b_k^{\left(N_F\right)}\right].\hfill \end{array}$$

(8)

Substituting Equation (5) into Equation (4) yields the following representation of $\overline{\mathcal{S}}\left({\xi }_k\right)$:

$$\overline{\mathcal{S}}\left({\xi }_k\right)=\left\{s_{k}R\left(\theta \right)B_{k}e\left(\theta \right)+x_{k,c}\left|\theta \in \left[0,2\pi \right]\right.,s_k\in \left[0,1\right]\right\}$$

(9)

Then, the measurement source $y_k^{\left(i\right)}$ can be drawn from $\overline{\mathcal{S}}\left({\xi }_k\right)$, i.e.,

$$\begin{array}{cc}\hfill y_k^{\left(i\right)}& \in s_{k,i}R\left({\theta }_{k,i}\right)B_{k}e\left({\theta }_{k,i}\right)+x_{k,c}\hfill \end{array}$$

(10)

where $s_{k,i}$ is a random draw from $s_k$, and ${\theta }_{k,i}$ is the angle between the vector from the centroid to $y_k^{\left(i\right)}$ and the x-axis.

Following the sensor model as explained in Equation (3), the measurement equation can be derived by

$$\begin{array}{cc}\hfill z_k^{\left(i\right)}& =y_k^{\left(i\right)}+v_k^{\left(i\right)}=s_{k,i}R\left({\theta }_{k,i}\right)B_{k}e\left({\theta }_{k,i}\right)+x_{k,c}+v_k^{\left(i\right)}\hfill \end{array}$$

(11)

Nevertheless, the angle ${\theta }_{k,i}$ is practically unknown. The angle ${\varphi }_{k,i}$ between the vector from the current centroid estimation to the observed measurement and the x-axis is used to approximate ${\theta }_{k,i}$. To minimize the effect of approximating ${\varphi }_{k,i}$ with ${\theta }_{k,i}$, the measurement equation can be converted into the following form, which maps the state ${\xi }_k$, measurement noise $v_k^{\left(i\right)}$, scaling factor $s_{k,i}$, and observed measurement $z_k^{\left(i\right)}$ to a pseudo-measurement 0:

$$\begin{array}{cc}\hfill 0& =h^{*}({\xi }_k,v_k^{\left(i\right)},s_{k,i},z_k^{\left(i\right)})\hfill \\ \hfill & ={s_{k,i}}^2{∥R\left({\varphi }_{k,i}\right)R\left({\varphi }_{k,i}\right)∥}^2+2s_{k,i}R\left({\varphi }_{k,i}\right)B_{k}e\left({\varphi }_{k,i}\right)v_k^{\left(i\right)}+{∥v_k^{\left(i\right)}∥}^2-{∥z_k^{\left(i\right)}-x_{k,c}∥}^2\hfill \end{array}$$

(12)

2.2. Multiple Extended Object Tracking Framework

In traditional multiple extended object tracking scenarios, object number is mostly unknown, and not all objects are detected by sensors. This implies that some objects in the surveillance area may remain undetected, and their corresponding measurements are not observed. Apart from the observed measurements from objects, there are clutters in the tracking scenarios, which introduces difficulty in jointly tracking multiple objects. Here, the measurements from objects are indistinguishable from clutters, and the order of measurements is listed without any specific significance. This implies that the association between objects, clutters, and measurements is uncertain.

To jointly track multiple extended objects, measurements and states at each time step are modeled as RFSs, respectively. This avoids complex data associations between objects and measurements. Derived from the GM-PHD filter, an ET-GM-PHD filter [18] was proposed for multiple extended object tracking. In the ET-GM-PHD filter, the measurement generation process is modeled using spatial distribution models [16]. It indicates that the state $x_k$ to be estimated only contains kinematical variables (e.g., positions and velocities). Similar to the GM-PHD filter, the ET-GM-PHD filter also propagates the posterior intensity of RFS of extended objects by prediction and update steps. The prediction step here is the same as that of the GM-PHD filter. Unlike the update step in the GM-PHD filter, the partitioning results of measurements are used to update the predicted Gaussian mixture components in the ET-GM-PHD filter. Here, the update step in the ET-GM-PHD filter is to be noticed, which is summarized as follows.

It is assumed that the RFSs of object states and measurements at time step k are denoted by $X_k={\left\{x_k^{\left(i\right)}\right\}}_{i=1}^{N_{x,k}}$ and $Z_k={\left\{z_k^{\left(j\right)}\right\}}_{j=1}^{N_{z,k}}$, respectively. Here, $N_{x,k}$ is the number of objects, and $N_{z,k}$ is the number of observed measurements. Different partitions (each of which contains several subsets) of measurements are obtained by measurement partitioning algorithms (e.g., the distance partitioning algorithm [18]). The predicted intensity $D_{k\left|k-1\right.}\left(x\right)$ is represented as follows:

$$D_{k\left|k-1\right.}\left(x\right)=\sum _{j=1}^{J_{k\left|k-1\right.}}{w}_{{}_{k\left|k-1\right.}}^{\left(j\right)}N\left(x;m_{{}_{k\left|k-1\right.}}^{\left(j\right)},P_{{}_{k\left|k-1\right.}}^{\left(j\right)}\right)$$

(13)

where

• $J_{k\left|k-1\right.}:$ the number of predicted Gaussian mixture components;
• $w_{{}_{k\left|k-1\right.}}^{\left(j\right)}:$ the weight of the $j^{th}$ predicted Gaussian mixture component;
• $m_{{}_{k\left|k-1\right.}}^{\left(j\right)},P_{{}_{k\left|k-1\right.}}^{\left(j\right)}:$ the predicted mean and covariance of the $j^{th}$ predicted Gaussian mixture component, respectively.

The update of $D_{k\left|k-1\right.}\left(x\right)$ consists of the update in an object non-detection case and that in a object detection case. Hence, the posterior intensity $D_{k\left|k\right.}\left(x\left|Z_k\right.\right)$ has the following forms:

$$\begin{array}{cc}\hfill D_{k\left|k\right.}\left(x\left|Z_k\right.\right)& =\left[1-\left(1-e^{-\gamma \left(x\right)}\right)p_D\left(x\right)\right]D_{k\left|k-1\right.}\left(x\right)+\sum _{P\angle Z_k}\sum _{W\in P}{D}_{k\left|k\right.}^D\left(x,W\right)\hfill \\ \hfill & =\left[1-\left(1-e^{-\gamma \left(x\right)}\right)p_D\left(x\right)\right]D_{k\left|k-1\right.}\left(x\right)+e^{-\gamma \left(x\right)}{p}_D\left(x\right)\hfill \\ \hfill & \times \sum _{P\angle Z_k}{\omega }_P\sum _{W\in P}\frac{\gamma {\left(x\right)}^{\left|W\right|}}{d_W}·\prod _{z_k\in W}\frac{{\varphi }_{z_k}\left(x\right)}{{\lambda }_k{c}_k\left(z_k\right)}{D}_{k\left|k-1\right.}\left(x\right)\hfill \end{array}$$

(14)

where

• $D_{k\left|k\right.}^D\left(x,W\right):$ the intensity updated by the measurement subset W;
• $\gamma \left(x\right):$ the expected measurement number of a single extended object, which obeys the Poisson distribution;
• $p_D\left(x\right):$ the probability of object detection;
• ${\lambda }_k:$ the mean number of measurements from clutters;
• $c_k\left(z_k\right):$ the spatial distribution of clutters;
• $P\angle Z_k:$ the $P^{th}$ partition in which the measurement set $Z_k$ is partitioned into several non-empty subsets;
• ${\omega }_P,d_W:$ non-negative coefficients for the partition P and the subset W, respectively;
• ${\varphi }_{z_k}\left(x\right):$ the likelihood function for the measurement $z_k$.

Note that the updated Gaussian mixture components need to be further pruned and merged to obtain the final tracking results because the number of Gaussian mixture components increases significantly after updating. The detailed pruning and merging steps are not mentioned here, but they can refer to [18].

3. An IMM-GMPHD Filter For Tracking Multiple Maneuvering Extended Objects Using RHMs

In this section, an IMM-GMPHD filter is proposed to jointly estimate the object number and states (including dynamic states and extension states) of multiple maneuvering star-convex extended objects using RHMs. For clarifying the presentation of the analytic implementation of this filter, turn maneuver models for star-convex objects are firstly presented based on RHMs. Update formulas of model probabilities for MMSEOT are further derived, followed by proposing an IMM-GMPHD recursion for tracking multiple maneuvering star-convex extended objects using RHMs.

3.1. Maneuver Models for Star-Convex Extended Objects Using RHMs

Suppose that the possible maneuver models for the $jth$ star-convex extended object using RHM are treated as a discrete set $\mathrm{M}={\left\{m^{\left(j\right)}\right\}}_{j=1}^{N_m}$, where $N_m$ is the number of maneuver models, and $m^{\left(j\right)}$ represents the $jth$ unique maneuver model. Assuming that the model in effect at time k is $m_k^{\left(i\right)}$, then the model-matched dynamic equation can be expressed according to Equation (2), i.e.,

$${\xi }_{k+1,i}^{\left(j\right)}=\left[\begin{array}{cc}{F}_{k,e,i}^{\left(j\right)}& 0\\ 0& F_{k,m,i}^{\left(j\right)}\end{array}\right]\left[\begin{array}{c}{x}_{k,e,i}^{\left(j\right)}\\ x_{k,m,i}^{\left(j\right)}\end{array}\right]+\left[\begin{array}{c}{\mathit{w}}_{k,e,i}\\ {\mathit{w}}_{k,m,i}\end{array}\right]$$

(15)

where subscript $\left(i\right)$ stands for the model $m_k^{\left(i\right)}$.

Since objects are treated as rigid bodies, the rotation of the object extension occurs with that of the object centroid. In other words, the rotation angle of the object extension remains exactly equal to that of the object centroid. In view of this, $F_{k,e,i}^{\left(j\right)}$ and $F_{k,m,i}^{\left(j\right)}$ can be explicitly derived.

Suppose that the turn rate of the $j^{th}$ extended object is ${\omega }_{k,i}^{\left(j\right)}$. Then, $F_{k,m,i}^{\left(j\right)}$ for the dynamic state $x_{k,m}^{\left(j\right)}={[x_k,y_k,{\dot{x}}_k,{\dot{y}}_k]}^T$ adopts the following widely applied form [27]:

$$F_{k,m,i}^{\left(j\right)}=\left[\begin{array}{cccc}1& 0& \frac{sin\left({\omega }_{k,i}^{\left(j\right)}T\right)}{{\omega }_{k,i}^{\left(j\right)}}& -\frac{1-cos\left({\omega }_{k,i}^{\left(j\right)}T\right)}{{\omega }_{k,i}^{\left(j\right)}}\\ 0& 1& cos\left({\omega }_{k,i}^{\left(j\right)}T\right)& -sin\left({\omega }_{k,i}^{\left(j\right)}T\right)\\ 0& 0& \frac{1-cos\left({\omega }_{k,i}^{\left(j\right)}T\right)}{{\omega }_{k,i}^{\left(j\right)}}& \frac{sin\left({\omega }_{k,i}^{\left(j\right)}T\right)}{{\omega }_{k,i}^{\left(j\right)}}\\ 0& 0& sin\left({\omega }_{k,i}^{\left(j\right)}T\right)& cos\left({\omega }_{k,i}^{\left(j\right)}T\right)\end{array}\right]$$

(16)

where T is the sampling time period, and ${\omega }_{k,i}^{\left(j\right)}T$ represents the rotation angle at per sampling time period.

For $F_{k,e,i}^{\left(j\right)}$, it can be derived according to the radial function. When the object extension rotates the same angle ${\omega }_{k,i}^{\left(j\right)}T$, the radial function translates the same angle. Hence, the radial function $r(x_{k+1,e,i}^{\left(j\right)},{\theta }_k)$ has the following equation:

$$r(x_{k+1,e,i}^{\left(j\right)},{\theta }_{k+1})=r(x_{k,e,i}^{\left(j\right)},{\theta }_{k+1}-{\omega }_{k,i}^{\left(j\right)}T)+{\mathit{w}}_{k,e,i}$$

(17)

where ${\mathit{w}}_{k,e,i}={[{\mathit{w}}_{k,i}^{\left(0\right)},{\mathit{w}}_{k,i}^{\left(1\right)},{\mathit{w}}_{k,i}^{\left(2\right)},...,{\mathit{w}}_{k,i}^{\left(2N_F\right)},{\mathit{w}}_{k,i}^{(2N_F+1)}]}^T$ indicates the process noise. Substituting the polar representation of the rational function in Equation (5) into Equation (17) results in the following equation:

$$\begin{array}{c}{a}_{k+1}^{\left(0\right)}+{\displaystyle \sum _{j=1,...,N_F}}(a_{k+1}^{\left(j\right)}cos\left(j{\theta }_{k+1}\right)+b_{k+1}^{\left(j\right)}sin\left(j{\theta }_{k+1}\right))\hfill \\ \hfill =a_k^{\left(0\right)}+{\mathit{w}}_{k,i}^{\left(0\right)}+{\displaystyle \sum _{j=1,...,N_F}}\left[\left[a_k^{\left(j\right)}cos\left(j{\omega }_{k,i}^{\left(j\right)}T\right)-b_k^{\left(j\right)}sin\left(j{\omega }_{k,i}^{\left(j\right)}T\right)\right]cos\left(j{\theta }_{k+1}\right)\right.\\ +{\mathit{w}}_{k,i}^{(2j-1)}\left.+\left[a_k^{\left(j\right)}sin\left(j{\omega }_{k,i}^{\left(j\right)}T\right)+b_k^{\left(j\right)}cos\left(j{\omega }_{k,i}^{\left(j\right)}T\right)\right]sin\left(j{\theta }_{k+1}\right)+{\mathit{w}}_{k,i}^{\left(2j\right)}\right]\hfill \end{array}$$

(18)

Equation (18) can be further rewritten in the following matrix form:

$$x_{k+1,e,i}^{\left(j\right)}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& F^{e,\left(1\right)}& 0& 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& F^{e,\left(N_F\right)}\end{array}\right]x_{k,e,i}^{\left(j\right)}+{\mathit{w}}_{k,e,i}$$

(19)

where

$$F^{e,\left(j\right)}=\left[\begin{array}{cc}cos\left(j{\omega }_{k,i}^{\left(j\right)}T\right)& -sin\left(j{\omega }_{k,i}^{\left(j\right)}T\right)\\ sin\left(j{\omega }_{k,i}^{\left(j\right)}T\right)& cos\left(j{\omega }_{k,i}^{\left(j\right)}T\right)\end{array}\right],j=1,\dots ,N_F$$

(20)

According to Equation (19), the transition matrix $F_{k,e,i}^{\left(j\right)}$ for the extension-involved state $x_{k,e,i}^{\left(j\right)}$ can be derived as follows:

$$F_{k,e,i}^{\left(j\right)}=diag\left(1,F^{e,1},\dots ,F^{e,N_F}\right)$$

(21)

With $F_{k,m,i}^{\left(j\right)}$ and $F_{k,e,i}^{\left(j\right)}$, the whole dynamic equation is obtained by Equation (15). It effectively handles the tight couple between the evolution of the dynamic state $x_{k,m,i}^{\left(j\right)}$ and that of the extension-involved state $x_{k,e,i}^{\left(j\right)}$. Moreover, it is a general dynamic equation to describe various turn maneuver models with different turn rates, such as left-turn motions and right-turn motions. It can be noted that when ${\omega }_{k,i}^{\left(j\right)}$ equals 0, the turn maneuver model reduces to a constant velocity motion. In this case, $F_{k,m,i}^{\left(j\right)}$ and $F_{k,e,i}^{\left(j\right)}$ become

$$F_{k,m,i}^{\left(j\right)}=\left[\begin{array}{cccc}1& 0& T& 0\\ 0& 1& 0& T\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(22)

$$F_{k,e,i}^{\left(j\right)}=diag\left(1,1,\dots ,1,1\right).$$

(23)

3.2. Model Probability Update for Maneuvering Extended Object Tracking

Each model in $\mathrm{M}$ has a corresponding model probability. Given the measurement set $W_k={\left\{z_j\right\}}_{j=1}^{\left|W_k\right|}$ at time k, the model probability ${\mu }_k^{\left(i\right)}$ matched to the model $m_k^{\left(i\right)}$ can be calculated as follows:

$$\begin{array}{cc}\hfill {\mu }_k^{\left(i\right)}& =P\left(m_k^{\left(i\right)}\left|W_k={\left\{z_j\right\}}_{j=1}^{\left|W_k\right|}\right.\right)\hfill \\ \hfill & =\frac{P\left(W_k={\left\{z_j\right\}}_{j=1}^{\left|W_k\right|}\left|m_k^{\left(i\right)}\right.\right)P\left(m_k^{\left(i\right)}\right)}{{\displaystyle \sum _{i=1}^{N_m}}P\left(W_k={\left\{z_j\right\}}_{j=1}^{\left|W_k\right|}\left|m_k^{\left(i\right)}\right.\right)P\left(m_k^{\left(i\right)}\right)}\hfill \end{array}$$

(24)

where $\left|W_k\right|$ represents the number of measurements in $W_k$, $P\left(m_k^{\left(i\right)}\right)$ is the probability that the $i^{th}$ model is in effect at time k, and $P\left(W_k={\left\{z_j\right\}}_{j=1}^{\left|W_k\right|}\left|m_k^{\left(i\right)}\right.\right)$ is the mode likelihood function for the measurement subset. Assuming that each measurement is generated independently, $P\left(W_k={\left\{z_j\right\}}_{j=1}^{\left|W_k\right|}\left|m_k^{\left(i\right)}\right.\right)$ can be computed by

$$P\left(W_k={\left\{z_j\right\}}_{j=1}^{\left|W_k\right|}\left|m_k^{\left(i\right)}\right.\right)=\sum _{j=1}^{\left|W_k\right|}P\left(z_j\left|m_k^{\left(i\right)}\right.\right)$$

(25)

where $P\left(z_j\left|m_k^{\left(i\right)}\right.\right)$ represents the probability form of the model likelihood function ${\mathsf{\Lambda }}_{j,k}^{\left(i\right)}$.

By substituting Equation (25) into Equation (24), ${\mu }_k^{\left(i\right)}$ can be derived as

$${\mu }_k^{\left(i\right)}=\frac{P\left(m_k^{\left(i\right)}\right){\displaystyle \sum _{j=1}^{\left|W_k\right|}}P\left(z_j\left|m_k^{\left(i\right)}\right.\right)}{{\displaystyle \sum _{i=1}^{N_m}}\left(P\left(m_k^{\left(i\right)}\right){\displaystyle \sum _{j=1}^{\left|W_k\right|}}P\left(z_j\left|m_k^{\left(i\right)}\right.\right)\right)}$$

(26)

To be specific, the substitution of $c_i$ and of ${\mathsf{\Lambda }}_{j,k}^{\left(i\right)}$ yields the desired representation of Equation (26), i.e.,

$${\mu }_k^{\left(i\right)}=\frac{c_i{\mathsf{\Lambda }}_{W_k}^{\left(i\right)}}{c}$$

(27)

where

$$c=\sum _{i=1}^{N_m}\left(c_i{\mathsf{\Lambda }}_{W_k}^{\left(i\right)}\right)$$

(28)

$${\mathsf{\Lambda }}_{W_k}^{\left(i\right)}=\sum _{j=1}^{\left|W_k\right|}{\mathsf{\Lambda }}_{j,k}^{\left(i\right)}.$$

(29)

3.3. The IMM-GMPHD Filtering Recursion

Using above maneuver models and update formulas of model probabilities, an IMM-GMPHD filter for MMSEOR is obtained. The complete framework of the IMM-GMPHD filter is given in Figure 2. To simultaneously track multiple maneuvering star-convex extended objects, the intensity of object states is propagated in the IMM-GMPHD filtering recursion, instead of a single object state. It is assumed that there are $N_m$ possible turn maneuver models for maneuvering star-convex extended objects. The re-initialized intensity matched to a maneuver model is calculated by the model-conditional re-initialization step. Furthermore, in the model-conditional filtering step, each model-conditional re-initialized intensity is updated by measurement partitioning results to obtain the model-conditional posterior intensity. Then, the model probability update step is performed to compute model probabilities. With model probabilities, the mixed posterior intensity is calculated by the estimate fusion step.

To yield the outputs of each model-conditional filter and the final estimation of the object number and states, the pruning and merging step contains three parts: (1) the pruning and merging of Gaussian mixture components of the mixed posterior intensity; (2) the pruning and merging of Gaussian mixture components of each model-conditional posterior intensity; and (3) the pruning and merging of model probabilities. Here, the merging of means of Gaussian mixture components is divided into the merging of dynamic parameters and that of extension-involved parameters. The merging of dynamic parameters is practically similar to that in the ET-GM-PHD filter [18]. However, extension-involved parameters cannot be directly merged owing to their geometric significance. As described previously, the extension-involved parameters of a star-convex extended object are determined by a radial function. In this manner, the merging of extension-involved parameters can indirectly turn to the merging of values of radial functions. Then, coefficients in the radial function after merging can be obtained by the least-squares method. These coefficients are the desired extension-involved parameters after merging.

The complete filtering procedures of the proposed IMM-GMPHD filter are summarized as follows.

(1) Model-Conditional Re-initialization: In this step, a re-initialized intensity is obtained for each model-matched filter, which is a combination of predicted intensities. The calculations of predicted intensities are not mentioned here, which can refer to the prediction step in [18].

Assuming that the model in effect at time $k-1$ is $m_{k-1}^{\left(i\right)}$ ($\forall i\in \left[1,N_m\right]$), the model-matched predicted posterior intensity $D_{i,k|k-1}\left(\xi \right)$ has the following form:

$$D_{i,k|k-1}\left(\xi \right)=\sum _{j=1}^{J_{k|k-1}}{w}_{i,k|k-1}^{\left(j\right)}N\left(\xi ;m_{i,k|k-1}^{\left(j\right)},P_{i,k|k-1}^{\left(j\right)}\right)$$

(30)

where subscript i denotes the $ith$ model $m_k^{\left(i\right)}$. Then the re-initialized intensity $D_{i,k-1\left|k-1\right.}\left(\xi \right)$ is calculated by

$$\begin{array}{cc}\hfill & D_{i,k-1\left|k-1\right.}\left(\xi \right)=\sum _{j=1}^{J_{k|\mathrm{k}-1}}{w}_{i,k-1\left|k-1\right.}^{\left(j\right)}N\left(\xi ;m_{i,k-1\left|k-1\right.}^{\left(j\right)},P_{i,k-1\left|k-1\right.}^{\left(j\right)}\right)\hfill \end{array}$$

(31)

$$w_{i,k-1\left|k-1\right.}^{\left(j\right)}=\sum _{r=1}^{N_m}{w}_{r,k|k-1}^{\left(j\right)}{\mu }_{(r,i),k-1\left|k-1\right.}^{\left(j\right)}$$

(32)

$$m_{i,k-1\left|k-1\right.}^{\left(j\right)}=\sum _{r=1}^{N_m}{m}_{r,k|k-1}^{\left(j\right)}{\mu }_{(r,i),k-1\left|k-1\right.}^{\left(j\right)}$$

(33)

$$\begin{array}{cc}\hfill P_{i,k-1\left|k-1\right.}^{\left(j\right)}& =\sum _{r=1}^{N_m}\left[P_{r,k|k-1}^{\left(j\right)}+\left(m_{r,k|k-1}^{\left(j\right)}-m_{i,k-1|k-1}^{\left(j\right)}\right){\left(m_{r,k|k-1}^{\left(j\right)}-m_{i,k-1|k-1}^{\left(j\right)}\right)}^{\mathrm{T}}\right]\hfill \\ \hfill & {\mu }_{(r,i),k-1\left|k-1\right.}^{\left(j\right)}\hfill \end{array}$$

(34)

where ${\mu }_{(r,i),k-1\left|k-1\right.}^{\left(j\right)}$ represents the mixing probability of the $jth$ Gaussian mixture component, which is calculated as follows:

$${\mu }_{(r,i),k-1\left|k-1\right.}^{\left(i,j\right)}={\pi }_{r,i}{\mu }_{k-1}^{\left(r\right)}/c_i$$

(35)

$$c_i=\sum _{m=1}^{N_m}{\pi }_{ri}{\mu }_{k-1}^{\left(m\right)}$$

(36)

where ${\pi }_{ri}$ denotes the transition probability from $m_{k-1}^{\left(r\right)}$ to $m_k^{\left(i\right)}$ ($i\ne r$, $\forall r,i\in \left[1,N_m\right]$), and $c_i$ indicates the probability that the $i^{th}$ model is in effect at time k.

(2) Model-Conditional Filtering: Using the above re-initialized intensity $D_{i,k-1\left|k-1\right.}\left(\xi \right)$, the update step is run to compute the posterior intensity $D_{i,k\left|k\right.}\left(\xi \right)$, $i\in \left[1,N_m\right]$. Considering the nonlinear measurement equations of star-convex extended objects, an appropriate update step [25] based on the RHM-based measurement model is introduced, which is a modified version of the update step mentioned in Part C of Section 2.

It is assumed that different partitions of the measurement set $Z_k$ at time k are available. The posterior intensity still has the same form as that in Equation (14). The right-hand side of Equation (14) indicates that each subset in each partition is utilized to update the re-initialized intensity $D_{i,k-1\left|k-1\right.}\left(\xi \right)$. Let $D_{i,k\left|k\right.}^D\left(\xi ,W_k^p\right)$ denote the result that $D_{i,k-1\left|k-1\right.}\left(\xi \right)$ is updated by a subset $W_k^p$. Here, nonlinear Bayes filtering approaches (such as the unscented transformation approach) are used to obtain $D_{i,k\left|k\right.}^D\left(\xi ,W_k^p\right)$.

Take updating the $jth$ Gaussian mixture component of $D_{i,k-1\left|k-1\right.}\left(\xi \right)$ by the $dth$ measurement $z_k^{\left(d\right)}$ in the $pth$ subset $W_k^p$, for example. Suppose that the updated weight, mean, and covariance are denoted by $w_{i,k\left|k\right.}^{\left(j\right)}$, $m_{i,k\left|k\right.}^{\left(j\right)}$, and $P_{i,k\left|k\right.}^{\left(j\right)}$, respectively. By generating L sigma points ${\left\{{\epsilon }_{i,k}^{\left(l\right)}\right\}}_{l=0}^L$ with weights ${\left\{w_{i,k}^{\left(l\right)}\right\}}_{l=0}^L$, $m_{i,k\left|k\right.}^{\left(j\right)}$ and $P_{i,k\left|k\right.}^{\left(j\right)}$ can be calculated as follows:

$$m_{i,k\left|k\right.}^{\left(j\right)}=m_{i,k-1\left|k-1\right.}^{\left(j\right)}-G_{i,k}^{\left(j\right)}{\left(S_{i,k}^{\left(j\right)}\right)}^{-1}\left(z_k^{\left(d\right)}-{\eta }_{i,k\left|k-1\right.}^{\left(j\right)}\right)$$

(37)

$$P_{i,k\left|k\right.}^{\left(j\right)}=P_{i,k-1\left|k-1\right.}^{\left(j\right)}-G_{i,k}^{\left(j\right)}{\left(S_{i,k}^{\left(j\right)}\right)}^{-1}{\left(G_{i,k}^{\left(j\right)}\right)}^{\mathrm{T}}$$

(38)

where

$$G_{i,k}^{\left(j\right)}=\sum _{l=0}^L{w}_{i,k}^{\left(l\right)}\left({\epsilon }_{i,k}^{\left(l\right)}-m_{i,k-1\left|k-1\right.}^{\left(j\right)}\right){\left(z_{i,k\left|k-1\right.}^{\left(l\right)}-{\eta }_{i,k\left|k-1\right.}^{\left(j\right)}\right)}^{\mathrm{T}}$$

(39)

$$S_{i,k}^{\left(j\right)}=\sum _{l=0}^L{w}_{i,k}^{\left(l\right)}\left(z_{i,k\left|k-1\right.}^{\left(l\right)}-{\eta }_{i,k\left|k-1\right.}^{\left(j\right)}\right){\left(z_{i,k\left|k-1\right.}^{\left(l\right)}-{\eta }_{i,k\left|k-1\right.}^{\left(j\right)}\right)}^{\mathrm{T}}$$

(40)

$${\eta }_{i,k\left|k-1\right.}^{\left(j\right)}=\frac{1}{1+L}\sum _{l=0}^L{z}_{i,k\left|k-1\right.}^{\left(l\right)}$$

(41)

$${\eta }_{i,k\left|k-1\right.}^{\left(j\right)}=\sum _{l=0}^L{w}_{i,k}^{\left(l\right)}{z}_{i,k\left|k-1\right.}^{\left(l\right)}$$

(42)

$$z_{i,k\left|k-1\right.}^{\left(l\right)}=h^{*}\left({\epsilon }_{i,k}^{\left(l\right)},v_{i,k}^{\left(d\right)},s_{k,d},z_k^{\left(d\right)}\right)$$

(43)

Additionally, the likelihood function ${\varphi }_{z_k^{\left(d\right)}}\left(\xi \right)$ for the measurement $z_k^{\left(d\right)}$ can be obtained by

$${\varphi }_{z_k^{\left(d\right)}}\left(\xi \right)=N\left(z_k^{\left(d\right)};{\eta }_{i,k\left|k-1\right.}^{\left(j\right)},S_{i,k}^{\left(j\right)}\right)$$

(44)

Then, the mean $m_{i,k\left|k\right.}^{\left(j\right)}$ and covariance $P_{i,k\left|k\right.}^{\left(j\right)}$ are sequentially updated by the remaining measurements in $W_k^p$. In other words, $m_{i,k\left|k\right.}^{\left(j\right)}$ and $P_{i,k\left|k\right.}^{\left(j\right)}$ substitute the predicted mean $m_{i,k-1\left|k-1\right.}^{\left(j\right)}$ and covariance $P_{i,k-1\left|k-1\right.}^{\left(j\right)}$ in Equations (37)–(39) for the subsequent measurement update. Once all the measurements in $W_k^p$ are used for the update, the mean $m_{i,k\left|k\right.}^{\left(j\right)}$ and covariance $P_{i,k\left|k\right.}^{\left(j\right)}$ are the final update results of the predicted mean $m_{i,k-1\left|k-1\right.}^{\left(j\right)}$ and covariance $P_{i,k-1\left|k-1\right.}^{\left(j\right)}$. The corresponding updated weight $w_{i,k\left|k\right.}^{\left(j\right)}$ can be derived from Equation (14), i.e.,

$$w_{i,k\left|k\right.}^{\left(j\right)}={\omega }_P\frac{e^{-\gamma \left(\xi \right)}{p}_D\left(\xi \right){\mathsf{\Phi }}_{W_k^p}\left(\xi \right)}{d_W{\displaystyle \prod _{z_k^{\left(d\right)}\in W_k^p}}{\lambda }_k{c}_k\left(z_k^{\left(d\right)}\right)}{w}_{i,k-1\left|k-1\right.}^{\left(j\right)}$$

(45)

where ${\mathsf{\Phi }}_{W_k^p}\left(\xi \right)$ denotes the likelihood function for this measurement subset $W_k^p$, which is calculated by

$${\mathsf{\Phi }}_{W_k^p}\left(\xi \right)=\prod _{z_k^{\left(d\right)}\in W_k^p}{\varphi }_{z_k^{\left(d\right)}}\left(\xi \right)$$

(46)

Reiterate Equation (37) to Equation (46) until all Gaussian mixture components of the re-initialized intensity are updated by the measurements in cell $W_k^p$. Then, the intensity $D_{i,k\left|k\right.}^D\left(\xi ,W_k^p\right)$ can be represented as follows:

$$D_{i,k\left|k\right.}^D\left(\xi ,W_k^p\right)=\sum _{j=1}^{J_{k\left|k-1\right.}}{w}_{i,k\left|k\right.}^{\left(j\right)}N\left(\xi ;m_{i,k\left|k\right.}^{\left(j\right)},P_{i,k\left|k\right.}^{\left(j\right)}\right)$$

(47)

On the basis of the measurement update equation in Equation (14), the expected form of the posterior intensity $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ can be rewritten as

$$\begin{array}{cc}\hfill D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)& =\left[1-\left(1-e^{-\gamma \left(\xi \right)}\right)p_D\left(\xi \right)\right]D_{i,k-1\left|k-1\right.}\left(\xi \right)+\hfill \\ \hfill & \sum _{P\angle Z_k}\sum _{W_k^p\in P}{D}_{i,k\left|k\right.}^D\left(\xi ,W_k^p\right)\hfill \end{array}$$

(48)

(3) Model Probability Update: The model probability for each Gaussian mixture component of the posterior intensity $D_{i,k\left|k\right.}^D\left(\xi ,W_k^p\right)$ is updated by a model likelihood function. According to Equation (27), the model probability ${\mu }_{i,k}^{\left(j\right)}$ for the $jth$ Gaussian mixture component of $D_{i,k\left|k\right.}^D\left(\xi ,W_k^p\right)$ can be rewritten as

$${\mu }_{i,k}^{\left(j\right)}=\frac{c_i^{\left(j\right)}{\mathsf{\Lambda }}_{W_{{}_k}^j}^{\left(i\right)}}{c^{\left(j\right)}}$$

(49)

$$c_i^{\left(j\right)}=\sum _{r=1}^{N_m}{\pi }_{ri}{\mu }_{i,k-1}^{\left(j\right)}$$

(50)

$$c^{\left(j\right)}=\sum _{i=1}^{N_m}\left(c_i^{\left(j\right)}{\mathsf{\Lambda }}_{W_{{}_k}^j}^{\left(i\right)}\right)$$

(51)

$${\mathsf{\Lambda }}_{W_{{}_k}^j}^{\left(i\right)}=\sum _{d=1}^{\left|W_p\right|}{\mathsf{\Lambda }}_{d,j,k}^{\left(i\right)}$$

(52)

where ${\mathsf{\Lambda }}_{d,j,k}^{\left(i\right)}$ denotes the model likelihood function of the $j^{th}$ Gaussian mixture component for model $m_k^{\left(j\right)}$, which is approximated by the likelihood function ${\varphi }_{z_k^{\left(d\right)}}\left(\xi \right)$ for the measurement $z_k^{\left(d\right)}$ in Equation (44).

(4) Estimate Fusion: Suppose that the posterior intensity $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ has the following form:

$$D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)=\sum _{j=1}^{J_{k\left|k\right.}}{w}_{i,k\left|k\right.}^{\left(j\right)}N\left(\xi ;m_{i,k\left|k\right.}^{\left(j\right)},P_{i,k\left|k\right.}^{\left(j\right)}\right)$$

(53)

Then, the mixed posterior intensity ${\widehat{D}}_{k\left|k\right.}\left(\xi \right)$ can be calculated as follows:

$${\widehat{D}}_{k\left|k\right.}\left(\xi \right)=\sum _{j=1}^{J_{k\left|k\right.}}{\widehat{w}}_{i,k\left|k\right.}^{\left(j\right)}N\left(\xi ;{\widehat{m}}_{i,k\left|k\right.}^{\left(j\right)},{\widehat{P}}_{i,k\left|k\right.}^{\left(j\right)}\right)$$

(54)

where

$${\widehat{w}}_{k\left|k\right.}^{\left(j\right)}=\sum _{i=1}^{N_m}{w}_{i,k\left|k\right.}^{\left(j\right)}{\mu }_{i,k}^{\left(j\right)}$$

(55)

$${\widehat{m}}_{k\left|k\right.}^{\left(j\right)}=\sum _{i=1}^{N_m}{m}_{i,k\left|k\right.}^{\left(j\right)}{\mu }_{i,k}^{\left(j\right)}$$

(56)

$${\widehat{P}}_{k\left|k\right.}^{\left(j\right)}=\sum _{i=1}^{N_m}\left[P_{i,k\left|k\right.}^{\left(j\right)}+\left(m_{i,k\left|k\right.}^{\left(j\right)}-{\widehat{m}}_{k\left|k\right.}^{\left(j\right)}\right){\left(m_{i,k\left|k\right.}^{\left(j\right)}-{\widehat{m}}_{k\left|k\right.}^{\left(j\right)}\right)}^{\mathrm{T}}\right]{\mu }_k^{\left(j\right)}$$

(57)

(5) Pruning and Merging: In this stage, Gaussian mixture components of the mixed posterior intensity ${\widehat{D}}_{k\left|k\right.}\left(\xi \right)$, Gaussian mixture components of posterior intensities $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ ($i\in \left[1,N_m\right]$) and model probabilities corresponding to $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ are needed to be pruned and merged, respectively. Precisely, they are described in the following three parts.

(a) Pruning and Merging of Gaussian Mixture Components of the Mixed Posterior Intensity: For the mixed posterior intensity ${\widehat{D}}_{k\left|k\right.}\left(\xi \right)$, Gaussian mixture components with weights lower than a truncation threshold $\tau$ are eliminated. Indexes of these eliminated components are denoted by a vector $\left[p_1,p_2,\cdots ,p_{N_p}\right]$, where $N_p$ is the number of eliminated components. After pruning, remaining Gaussian mixture components whose distances fall within a merging threshold U are needed to be merged. Indexes of these Gaussian mixture components are denoted by $U_m=\left[u_1,u_2,\cdots ,u_{N_u}\right]={\left\{u_g\right\}}_{g=1}^{N_u}$, where $N_u$ is the number of Gaussian mixture components to be merged. Let $w_k^{\left(f\right)}$, $m_k^{\left(f\right)}$, and $P_k^{\left(f\right)}$ denote the weight, mean, and covariance of the $f^{th}$ merged Gaussian mixture component, respectively. According to [18], $w_k^{\left(f\right)}$ and $P_k^{\left(f\right)}$ are calculated by

$$w_k^{\left(f\right)}=\sum _{u_g\in U_m}{\widehat{w}}_{k\left|k\right.}^{\left(u_g\right)}$$

(58)

$$\begin{array}{cc}\hfill P_k^{\left(f\right)}& =\frac{1}{w_k^{\left(f\right)}}\sum _{u_g\in U_m}{\widehat{w}}_{k\left|k\right.}^{\left(u_g\right)}\left({\widehat{P}}_{k\left|k\right.}^{\left(u_g\right)}+\left(m_k^{\left(f\right)}-m_{k\left|k\right.}^{\left(u_g\right)}\right){\left(m_k^{\left(f\right)}-{\widehat{m}}_{k\left|k\right.}^{\left(u_g\right)}\right)}^{\mathrm{T}}\right)\hfill \end{array}$$

(59)

The merging of means ${\widehat{m}}_{k\left|k\right.}^{\left(u_g\right)}$ includes the merging of extension states ${\widehat{m}}_{e,k\left|k\right.}^{\left(u_g\right)}$ and that of dynamic states ${\widehat{m}}_{m,k\left|k\right.}^{\left(u_g\right)}$. Suppose that $m_{e,k}^{\left(f\right)}$ and $m_{m,k}^{\left(f\right)}$ are extension states and dynamic states after merging, respectively. $m_{m,k}^{\left(f\right)}$ can be directly calculated by

$$m_{m,k}^{\left(f\right)}=\frac{1}{w_k^{\left(f\right)}}\sum _{u_g\in U_m}{\widehat{w}}_{k\left|k\right.}^{\left(u_g\right)}{\widehat{m}}_{m,k\left|k\right.}^{\left(u_g\right)}$$

(60)

Unlike the direct merging of ${\widehat{m}}_{m,k\left|k\right.}^{\left(u_g\right)}$, we merge radial functions to obtain the final $m_{e,k}^{\left(f\right)}$. Let $r\left(m_{e,k}^{\left(f\right)},{\theta }_k\right)$ denote the radial function after merging. The value of $r\left(m_{e,k}^{\left(f\right)},{\theta }_k\right)$ at each integral angle ${\theta }_k=0,1,2,3\cdots ,2\pi$ can be calculated by

$$r\left(m_{e,k}^{\left(f\right)},{\theta }_k\right)=\frac{{\displaystyle \sum _{u_g\in U_m}}R\left({\theta }_k\right){\widehat{m}}_{e,k\left|k\right.}^{\left(u_g\right)}}{N_u}$$

(61)

Then coefficients of the polar representation of $r\left(m_{e,k}^{\left(f\right)},{\theta }_k\right)$ can be obtained using the least-squares method. They are the desired extension-involved parameters $m_{e,k}^{\left(f\right)}$ after merging. The detailed least-squares method is omitted here.

(b) Pruning and Merging of Gaussian Mixture Components of Each Posterior Intensity: For each posterior intensity $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ ($i\in \left[1,N_m\right]$), Gaussian mixture components with indexes $\left[p_1,p_2,\cdots ,p_{N_p}\right]$ are eliminated, and the rest of the Gaussian mixture components with indexes $\left[u_1,u_2,\cdots ,u_{N_u}\right]$ are merged. The detailed merging procedures for the posterior intensity $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ are the same as that for the mixed posterior intensity ${\widehat{D}}_{k\left|k\right.}\left(\xi \right)$.

(c) Pruning and Merging of Model Probabilities: For each posterior intensity $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ ($i\in \left[1,N_m\right]$), its corresponding model probabilities with indexes $\left[p_1,p_2,\cdots ,p_{N_p}\right]$ are eliminated. Remaining model probabilities with indexes $\left[u_1,u_2,\cdots ,u_{N_u}\right]$ are needed to be merged.

Let ${\mu }_{i,k}^{\left(f\right)}$ denote the model probability for the $f^{th}$ Gaussian mixture component after merging and ${\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}$ represent subsets used to update Gaussian mixture components with indexes $\left[u_1,u_2,\cdots ,u_{N_u}\right]$. Following the law of total probability, ${\mu }_{i,k}^{\left(f\right)}$ can be derived as follows:

$$\begin{array}{cc}\hfill {\mu }_{i,k}^{\left(f\right)}& =P\left\{m_k^{\left(i\right)}\left|{\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\right.\right\}=\frac{P\left({\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\left|m_k^{\left(i\right)}\right.\right)P\left(m_k^{\left(i\right)}\right)}{{\displaystyle \sum _{i=1}^{N_m}}P\left({\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\left|m_k^{\left(i\right)}\right.\right)P\left(m_k^{\left(i\right)}\right)}\hfill \end{array}$$

(62)

where $P\left({\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\left|m_k^{\left(i\right)}\right.\right)$ represents the model likelihood function for these subsets. Following the knowledge of probability, $P\left({\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\left|m_k^{\left(i\right)}\right.\right)$ can be further calculated by

$$\begin{array}{cc}\hfill P\left({\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\left|m_k^{\left(i\right)}\right.\right)=& P\left(W_k^{u_1}\left|m_k^{\left(i\right)}\right.\right)+P\left(W_k^{u_2}\left|m_k^{\left(i\right)}\right.\right)+\cdots \hfill \\ \hfill & +P\left(W_k^{u_{N_u}}\left|m_k^{\left(i\right)}\right.\right)-P\left(W_k^{u_1}\cap {\left\{W_k^{u_g}\right\}}_{g=2}^{N_u}\left|m_k^{\left(i\right)}\right.\right)\hfill \\ \hfill & -P\left(W_k^{u_2}\cap {\left\{W_k^{u_g}\right\}}_{g=3}^{N_u}\left|m_k^{\left(i\right)}\right.\right)-\cdots \hfill \\ \hfill & -P\left(W_k^{u_{N_u-1}}\cap \left\{W_k^{u_{N_u}}\right\}\left|m_k^{\left(i\right)}\right.\right)\hfill \end{array}$$

(63)

where $W_k^{u_1}\cap {\left\{W_k^{u_g}\right\}}_{g=2}^{N_u}$ denotes the measurements coexisting in $W_k^{u_1}$ and ${\left\{W_k^{u_g}\right\}}_{g=2}^{N_u}$. Subsets in the same partition are regarded to be independent because no measurements in any subset will be the same. However, there may be measurements coexisting in subsets of different partitions. Considering that these measurements have little effect on the model likelihood function, $p\left({\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\left|m_k^{\left(i\right)}\right.\right)$ can be approximated by

$$P\left({\left\{W_k^{u_g}\right\}}_{g=1}^{N_u}\left|m_k^{\left(i\right)}\right.\right)\approx \sum _{g=1}^{N_u}P\left(W_k^{u_g}\left|m_k^{\left(i\right)}\right.\right)$$

(64)

Substituting Equation (64) into Equation (62) yields the following form of ${\mu }_{i,k}^{\left(f\right)}$:

$${\mu }_{i,k}^{\left(f\right)}\approx \frac{{\displaystyle \sum _{g=1}^{N_u}}P\left(W_k^{u_g}\left|m_k^{\left(i\right)}\right.\right)p\left(m_k^{\left(i\right)}\right)}{{\displaystyle \sum _{i=1}^{N_m}}\left[{\displaystyle \sum _{g=1}^{N_u}}P\left(W_k^{u_g}\left|m_k^{\left(i\right)}\right.\right)p\left(m_k^{\left(i\right)}\right)\right]}$$

(65)

where $P\left(W_k^{u_g}\left|m_k^{\left(i\right)}\right.\right)$ and $P\left(m_k^{\left(i\right)}\right)$ are calculated by Equation (25) and Equation (36), respectively. To avoid any confusion, use $c_i^{\left(f\right)}$ and ${\mathsf{\Lambda }}_{W_k^{u_g}}^{\left(i\right)}$ to replace $P\left(m_k^{\left(i\right)}\right)$ and $P\left(W_k^{u_g}\left|m_k^{\left(i\right)}\right.\right)$, respectively. Then, Equation (65) can be rewritten as

$${\mu }_{i,k}^{\left(f\right)}=\frac{c_i^{\left(f\right)}{\displaystyle \sum _{g=1}^{N_u}}{\mathsf{\Lambda }}_{W_k^{u_g}}^{\left(i\right)}}{c^{\left(f\right)}}$$

(66)

where

$$c^{\left(f\right)}=\sum _{i=1}^{N_m}\left(c_i^{\left(f\right)}\sum _{g=1}^{N_u}{\mathsf{\Lambda }}_{W_k^{u_g}}^{\left(i\right)}\right)$$

(67)

The complete pseudo codes for the pruning and merging step are summarized in

Table 1. Pseudo codes for the pruning and merging step.

Input: the mixed posterior intensity ${\widehat{D}}_{k\left|k\right.}\left(\xi \right)$,

each posterior intensity $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$ ($i\in \left[1,N_m\right]$),

and model probabilities corresponding to $D_{i,k\left|k\right.}\left(\xi \left|Z_k\right.\right)$.

$f=0,I={\left\{j\left|{\widehat{w}}_{k\left|k\right.}^{\left(j\right)}>\tau \right.\right\}}_{j=1}^{J_{k\left|k\right.}},i\in \left[1,N_m\right]$

Repeat

$f=f+1$

$p_m=\underset{j\in I}{argmax{\widehat{w}}_{i,k\left|k\right.}^{\left(j\right)}}$

$U_m=\left\{j\in I{\left|{\left({\widehat{m}}_{k\left|k\right.}^{\left(j\right)}-{\widehat{m}}_{k\left|k\right.}^{\left(p_m\right)}\right)}^{\mathrm{T}}\left({\widehat{P}}_{k\left|k\right.}^{\left(j\right)}\right)\right.}^{-1}\left({\widehat{m}}_{k\left|k\right.}^{\left(j\right)}-{\widehat{m}}_{k\left|k\right.}^{\left(p_m\right)}\right)\le U\right\}$

merging Gaussian mixture components of the mixed posterior intensity:

$w_k^{\left(f\right)}={\displaystyle \sum _{j\in U_m}}{\widehat{w}}_{k\left|k\right.}^{\left(j\right)}$

$m_{m,k}^{\left(f\right)}=\frac{1}{w_k^{\left(f\right)}}{\displaystyle \sum _{j\in U_m}}{\widehat{w}}_{k\left|k\right.}^{\left(j\right)}{\widehat{m}}_{m,k\left|k\right.}^{\left(j\right)}$

$R\left(m_{e,k}^{\left(f\right)},\theta \right)=\frac{{\displaystyle \sum _{j\in U_m}}R\left(\theta \right){\widehat{m}}_{e,k\left|k\right.}^{\left(j\right)}}{N_u}$,$\theta \in \left[0,2\pi \right]$

$m_{e,k}^{\left(f\right)}=\mathrm{leastsquare}\left(R\left(m_{e,k}^{\left(f\right)},\theta \right),\theta \right)$

$P_k^{\left(f\right)}=\frac{1}{w_k^{\left(f\right)}}{\displaystyle \sum _{j\in U_m}}{\widehat{w}}_{k\left|k\right.}^{\left(j\right)}\left({\widehat{P}}_{k\left|k\right.}^{\left(j\right)}+\left(m_k^{\left(f\right)}-m_{k\left|k\right.}^{\left(j\right)}\right){\left(m_k^{\left(f\right)}-{\widehat{m}}_{k\left|k\right.}^{\left(j\right)}\right)}^{\mathrm{T}}\right)$

merging model probabilities:

${\mu }_{i,k}^{\left(f\right)}=\frac{c_i^{\left(f\right)}{\displaystyle \sum _{j\in U_m}}{\mathsf{\Lambda }}_{W_k^j}^{\left(i\right)}}{c^{\left(f\right)}}$

${\mathsf{\Lambda }}_{W_k^j}^{\left(i\right)}={\displaystyle \sum _{d=1}^{\left|W_k^j\right|}}{\mathsf{\Lambda }}_{d,f,k}^{\left(i\right)}$

$c^{\left(f\right)}={\displaystyle \sum _{i=1}^{N_m}}{c}_i^{\left(f\right)}{\mathsf{\Lambda }}_{W_k^j}^{\left(i\right)}$

merging Gaussian mixture components of each posterior intensity:

$w_{i,k}^{\left(f\right)}={\displaystyle \sum _{j\in U_m}}{w}_{i,k\left|k\right.}^{\left(j\right)}$

$m_{i,m,k}^{\left(f\right)}=\frac{1}{w_{i,k}^{\left(f\right)}}{\displaystyle \sum _{j\in U_m}}{w}_{i,k\left|k\right.}^{\left(j\right)}{m}_{i,m,k\left|k\right.}^{\left(j\right)}$

$R\left(m_{i,e,k}^{\left(f\right)},\theta \right)=\frac{{\displaystyle \sum _{j\in U_m}}R\left(\theta \right)m_{i,e,k\left|k\right.}^{\left(j\right)}}{N_u}$,$\theta \in \left[0,2\pi \right]$

$m_{i,e,k}^{\left(f\right)}=\mathrm{leastsquare}\left(R\left(m_{i,e,k}^{\left(f\right)},\theta \right),\theta \right)$

$P_{i,k}^{\left(f\right)}=\frac{1}{w_{i,k}^{\left(f\right)}}{\displaystyle \sum _{j\in U_m}}{w}_{i,k\left|k\right.}^{\left(j\right)}\left(P_{i,k\left|k\right.}^{\left(j\right)}+\left(m_{i,k}^{\left(f\right)}-m_{i,k\left|k\right.}^{\left(j\right)}\right){\left(m_{i,k}^{\left(f\right)}-m_{i,k\left|k\right.}^{\left(j\right)}\right)}^{\mathrm{T}}\right)$

$I=I\backslash U_m$

until$I=\varnothing$

$J_k=f$

Output: the final mixed posterior intensity $D_k\left(\xi \right)={\displaystyle \sum _{j=1}^{J_k}}{w}_k^{\left(j\right)}N\left(x;m_k^{\left(j\right)},P_k^{\left(j\right)}\right)$,

each posterior intensity $D_{i,k}\left(\xi \right)={\displaystyle \sum _{j=1}^{J_k}}{w}_{i,k}^{\left(j\right)}N\left(x;m_{i,k}^{\left(j\right)},P_{i,k}^{\left(j\right)}\right)$ ($i\in \left[1,N_m\right]$),

and model probabilities corresponding to $D_{i,k}\left(\xi \right)$.

.

To summarize, the proposed IMM-GMPHD filter can effectively track multiple maneuvering extended objects with RHMs. It successfully derives the model probabilities for extended objects, which facilitates handling the uncertainty of objects’ maneuver models. This is obviously more reasonable for MMSEOT than the original IMM algorithm (of which both the maneuvering models and model probabilities are based on point target assumption). Moreover, the maneuvering models for extended objects with RHMs are applied in our filter to accurately describe the complex maneuvering behaviors of extended objects. More importantly, the geometrical significance is fully considered in the pruning and merging. This greatly facilitates the more accurate estimation of object extensions, compared with the filtering steps in the ET-GM-PHD filter based on original RHMs.

4. Simulation Results and Performance Evaluation

In this section, the proposed IMM-GMPHD filter for tracking multiple maneuvering extended objects using star-convex RHMs is evaluated by the means of moving maneuvering star-convex extended objects with different turn rates. To demonstrate the viability of this filter in handling complex maneuvering behaviors of multiple extended objects, we have designed three clutter tracking scenarios within a fixed sensor platform, which are as follows:

(1) Scenario A: a tracking scenario with two maneuvering star-convex extended objects simultaneously rotating with different turn rates in different directions (i.e., left turn and right turn);

(2) Scenario B: a tracking scenario with two maneuvering star-convex extended objects simultaneously rotating with different turn rates in the same directions (i.e., left turn);

(3) Scenario C: a tracking scenario with two maneuvering star-convex extended objects of complex crossing trajectories.

In these tracking scenarios, the number of parameters in a extension-involved state (as seen in Equation (8)) is set to be 9, i.e., $x_{k,e,i}^{\left(j\right)}=[a_k^{\left(0\right)},a_k^{\left(1\right)},b_k^{\left(1\right)},a_k^{\left(2\right)},$$b_k^{\left(2\right)},a_k^{\left(3\right)},b_k^{\left(3\right)},a_k^{\left(4\right)},b_k^{\left(4\right)}{]}^{\mathrm{T}}$. The detection probability of each object is assumed to be $p_D=0.99$. To avoid the discrimination, each object extension is initialized by a circle with the same radius (i.e., 70), and the initial centroid state of each object is $x_{k,m,i}^{\left(j\right)}={\left[0,0,0,0\right]}^{\mathrm{T}}$. Transition matrices for $x_{k,m,i}^{\left(j\right)}$ and $x_{k,e,i}^{\left(j\right)}$ can be simply calculated by Equation (16) and Equation (21), respectively. The process noises ${\mathit{w}}_{k,e,i}$ and ${\mathit{w}}_{k,m,i}$ are white Gaussian noises with covariances $Q^e=0.{01}^{2}diag\left[1,1,1,1,1,1,1,1,1\right]$ and $Q^m=0.{01}^{2}diag\left[1,1,1,1\right]$, respectively. Furthermore, it is assumed that each object randomly produces 9 measurements and there are clutters over the surveillance region at each time step. Note that measurements of each object are independent of other objects’ measurements and clutters. The sensor noise $v_k^{\left(i\right)}$ is a white Gaussian noise with the covariance $R=diag\left[5,5\right]$. To verify the feasibility of the proposed filter for MMSEOT, we use three different turn maneuver models: a constant velocity model (a special turn maneuver with the turn rate ${\omega }_A^{\left(1\right)}=0rad/s$), a constant left-turn model (with the turn rate ${\omega }_A^{\left(2\right)}=5\pi /\phantom{5\pi 180}180rad/s$), and a constant right-turn model (with the turn rate ${\omega }_A^{\left(3\right)}=-10\pi /\phantom{5\pi 180}180rad/s$). The model probabilities of these turn maneuver models were initialized as $1/\phantom{13}3$.

Generally, different object extension models correspond to different tracking algorithms. To be fair, the analogous methods for comparisons need to be related to the maneuvering object tracking using RHMs. However, it may be hard for us to find a suitable approach for comparison. Though there are several methods [32] for extended object tracking, they are not based on RHMs. Moreover, the existing algorithms based on star-convex RHMs rarely consider the maneuvering problems of objects because of the complex and unknown maneuver models of star-convex extended objects. Thus, we compare the proposed IMM-GMPHD filter with an IMM approach to demonstrate our filter’s effectiveness. The compared IMM approach applies the original non-maneuvering RHM models and filtering steps of the ET-GM-PHD filter [25] to the IMM-GMPHD recursion. By doing so, each innovative part of our filter can be compared, which facilitates the verification of our algorithms.

4.1. Tracking Performance in Scenario A

In this scenario, Object 1 and Object 2 are to be tracked, which started at time step $k=1s$ and ended at time step $k=66s$. Object 1 and Object 2 simultaneously maneuvered with different turn rates at time step $k=31s$∼$39s$. The estimation results and true trajectories of these objects are shown in Figure 3a. Gray-filled shapes are true object shapes, cross symbols denote estimations of object centers, and shapes with red and blue lines are object shape estimations by the proposed filter and the IMM algorithm based on [25], respectively.

Considering that it is hard to jointly evaluate the estimation performance of object extensions, we give detailed estimated shapes of Object 1 and Object 2 in Figure 4. From Figure 4, it is evident that object extension estimations are gradually accurate in the proposed filter. Even when these objects simultaneously maneuvered, the estimation of extensions are relatively accurate. However, the IMM algorithm based on [25] has worse estimation performance in object extensions. It illustrates that the proposed maneuver models can describe different maneuvers of the star-convex objects effectively. Furthermore, they can be easily applied in this filter to obtain good estimation performance in object extensions. Moreover, Figure 4 demonstrates that the proposed merging strategy of $x_{k,e}^{\left(j\right)}$ is reliable and capable of estimating object extensions accurately.

The estimation performance in the object number and centroid dynamics is jointly evaluated by a widely used evaluation measure known as the OSPA distance [33]. The higher the OSPA distance is, the worse the performance in the estimation of the object number and centroid positions is. From Figure 3b, the OSPA distance of the proposed filter is lower than that of the IMM algorithm based on [25]. Though the OSPA distance of the proposed filter is a little higher in the middle due to the unreasonable partitioning results of measurements, it stays low at other time steps. This implies that the proposed filter achieves more accurate estimation performance in the object number and centroid positions compared with the IMM algorithm based on [25]. This is because the proposed IMM-GMPHD filter fits well for MMSEOT to handle the complex maneuvering behaviors of star-convex extended objects.

4.2. Tracking Performance in Scenario B

In Scenario B, the star-convex extended objects to be tracked are denoted by Object 3 and Object 4, both of which started at $k=1s$ and ended at $k=90s$. In this scenario, two constant left-turn motions with different turn rates (i.e., ${\omega }_B^{\left(1\right)}=5\pi /\phantom{5\pi 180}180rad/s$ and ${\omega }_B^{\left(2\right)}=10\pi /\phantom{5\pi 180}180rad/s$) and a constant velocity model with the turn rate ${\omega }_B^{\left(3\right)}=0rad/s$ are used to further test the effectiveness of the proposed filter. The initial model probabilities of these maneuver models are the same as those of Scenario A.

The true trajectories and estimation results are illustrated in Figure 5a, and the detailed estimation results of Object 3 and Object 4 are shown in Figure 6. From these figures, it can be concluded that object shape estimations are gradually close to true shapes from initial circles. When these objects maneuvered with different turn rates simultaneously, the estimation of object shapes and orientations by the proposed filter remains fairly more accurate than that of the IMM algorithm based on [25]. This demonstrates that this proposed filter can effectively handle the complex maneuvering behaviors of multiple star-convex extended objects. It also shows that the pruning and merging step in the proposed algorithm contributes to the accurate extension estimations.

Figure 5b illustrates the OSPA distance at each time step. It shows that the estimation of the object number and dynamic states by the presented filter are more accurate than that of the IMM algorithm based on [25]. In other words, the proposed IMM-GMPHD filtering recursion and the update formulas of model probabilities are effective for tracking multiple maneuvering star-convex extended objects.

4.3. Tracking Performance in Scenario C

To further validate the feasibility of the proposed filter, a complex tracking scenario (i.e., Scenario C) with crossing trajectories was chosen. In Section C, the true trajectories and estimation results of Object 5 and Object 6 are shown in Figure 7a. It can be seen that these objects crossed twice. When Object 5 and Object 6 first crossed around time step $k=23s$∼$30s$, their turn rates were equal to ${\omega }_C^{\left(1\right)}=0rad/s$. When Object 5 and Object 6 crossed around time step $k=36s$∼$45s$, they rotated with different turn rates in different directions, i.e., Object 5 made a left-turn motion with the turn rate ${\omega }_C^{\left(2\right)}=5\pi /\phantom{5\pi 180}180rad/s$, while Object 6 made a right-turn motion with the turn rate ${\omega }_C^{\left(3\right)}=-5\pi /\phantom{5\pi 180}180rad/s$.

Figure 8 gives detailed tracking results of Object 5 and Object 6. The object extension estimations by the proposed filter are ultimately close to true object shapes. When objects made turn motions or crossed, object shapes and orientations were especially more accurate compared with the IMM algorithm based on [25]. This is due to feasible maneuver models, strict and scientific update formulas of model probabilities, and the effective merging step of extension states. The estimation performance of the object number and positions is assessed by the OSPA distance in Figure 7b. The OSPA distance of the proposed filter is relatively lower than that of the IMM algorithm based on [25]. This demonstrates that the proposed filter can achieve more accurate estimation of the object number and centroid positions than the IMM algorithm based on [25]. To summarize, our filter is capable of handling complex tracking issues in MMSEOT and obtains high tracking performances in object number, centroid dynamics, and object extensions.

5. Conclusions

To deal with the unknown and time-varying maneuvering behaviors of star-convex extended objects, this paper has proposed an IMM-GMPHD filter. By integrating linear maneuver models into this filter, various turn maneuvers of objects can be described accurately from two aspects (i.e., the dynamic evolutions of both centroids and extensions). In our filter, rigorous update and merging formulas of model probabilities have been derived for dealing with the uncertainty of maneuvering models. This is the main advantage over other extended object tracking approaches. Moreover, object states have been properly merged and pruned to obtain accurate tracking results. Specifically for extension-involved states, they turned to radial functions and merged, which further improves the estimation performance of object extensions. Simulation results have verified that the proposed filter can achieve more accurate estimation performance in object number, centroid dynamics, and extensions (e.g., orientations and shapes) compared with the IMM algorithm based on [25] with original ET-GM-PHD filtering steps and RHMs. In summary, our filter may be promising for handling a large range of problems in multiple maneuvering extended object tracking.