Modified prediction partitioning algorithm for Gaussian inverse Wishart probability hypothesis density filter using measurement labeling method

Abstract

The Gaussian inverse Wishart probability hypothesis density filter is a promising approach for tracking multiple extended targets. However, if targets are closely spaced and performing maneuvers, the performance of a Gaussian inverse Wishart probability hypothesis density filter will decline. The reason for this is that the measurement partitioning approaches fail to provide accurate partitions, which influences the component updating process directly. This paper presents a modified prediction partitioning algorithm for the Gaussian inverse Wishart probability hypothesis density filter in order to solve the partitioning problem of closely spaced targets. The inaccuracy of the target prediction information occurred by target maneuvers leads to the above problem, thus a modified prediction partitioning algorithm will label the components and corresponding measurements to improve the prediction accuracy. Simulation results show that the use of modified prediction partitioning can improve the performance of a Gaussian inverse Wishart probability hypothesis density filter by providing more accurate partition results when targets are closely spaced and performing maneuvers.

Introduction

Typically, multiple target tracking (MTT) is based on an assumption that each target produces at most a single measurement per scan. Much research work has been done for MTT, such as by Mahler, ¹ Vo and Ma, ² Vo et al., ³ and Chen et al. ⁴ However, along with the development of the sensor technology, high-resolution sensors, such as the laser radar, may produce more than one measurement per scan, thus typical MTT approaches are no longer applicable for these type targets, that is extended target. Recently, much research work has been done for extended target tracking (ETT) research,^5–9 and this is especially true for multiple ETT (METT).^10–15

A probability hypothesis density (PHD) filter ¹ was first presented by Mahler to solve the MTT problems. PHD is a function of the first-order moment over the target state space, and an important implementation of the PHD filter is the Gaussian mixture PHD (GM-PHD) filter. ² This filter employs the Gaussian mixture method to approximate the target state PHD, thus the target state can be estimated by the Gaussian mixture components. A particle implementation of the GM-PHD filter can be found by Vo et al., ³ and an improved GM-PHD filter to deal with state-dependent clutter can be found by Chen et al. ⁴

Tracking an extended target is quite different from tracking a single-point target. Specifically, the target extension state will be estimated by an ETT filter. A typical ETT method is the random matrix (RM) model,^5,6 which treats the target shape as an elliptical extension, thus the target extended state is defined as a symmetric positive definite matrix (i.e. covariance matrix). The RM models assume that the target extended state is Wishart distributed, therefore it can be estimated by inverse Wishart distribution. Another important ETT method is the random hypersurface model (RHM),^7–9 which assumes that the target measurements are randomly produced from a hypersurface. The target shape function is determined by this hypersurface, thus RHMs can track a target and estimate its shape by computing the parameters of the shape function.

The RM and RHM approaches cannot track multiple extended targets, thus a multiple extended target PHD framework ¹⁰ was presented by Mahler. Employing the Gaussian mixture method and the above framework, Granström et al. ¹¹ presented the extended target GM-PHD (ET-GM-PHD) filter. The ET-GM-PHD filter can track an unknown number of extended targets efficiently, but the estimation error will occur when targets are spaced closely together. The reason for this is that the distance partitioning (DP) approach fails to provide accurate partitions to update the components. Hence, a sub-partitioning (SP) approach was presented by Granström et al., ¹² which is an added partitioning process following DP to divide a measurement cell produced by two or more targets into sub-cells. However, the ET-GM-PHD filter assumes that the extended state of a target is a given value, thus it cannot estimate the target extended state. To solve this problem, a GIW-PHD filter was presented by Granström and Orguner. ¹³ This filter employs the RM approach and assumes that the last estimated PHD is approximated with an unnormalized mixture of Gaussian inverse Wishart (GIW) distributions, thus this filter can track multiple extended targets and estimate their extended states. More discussions of the GIW-PHD filter can be found by Granström et al.^14,15 However, when targets are spaced closely together and performing maneuvers, estimation errors will occur, because the partitioning approaches, such as DP-SP, expectation-maximization partitioning (EMP), and prediction partitioning (PP),^11–13 fail to provide accurate partition results. A fast partitioning approach was presented by Zhang and Ji ¹⁶ to reduce the time cost of the partitioning step, but this work does not improve the performance of the DP-SP, PP, and EMP. A shape selection partitioning (SSP) algorithm was presented by Li et al.,^17,18 which can solve the partitioning performance when targets are closely spaced and performing maneuvers. However, the SSP applied only to the scenarios where the number of targets is two, and the computation costs of this algorithm are quite larger than the DP-SP, PP, and EMP.

The measurement partitioning is an integral part of the GIW-PHD filter, therefore the improvement of the partition results can improve the performance directly. This paper modified the measurements partitioning step of the GIW-PHD filter. The contributions of this paper are as follows:

Since the GIW-PHD filter cannot provide track identities, this paper employs the component labeling technique ¹⁹ for the GIW-PHD filter to label the GIW components.

This paper presents the modified prediction partitioning (MPP) approach to solve the partitioning problems of closely spaced targets. The key method is that we not only label the GIW components, but also label their corresponding measurement cells. Then the prediction parameters used in MPP will be improved by analyzing the changes of measurement cells between time k and k + 1. As a result, the accuracy of MPP is higher than DP-SP, PP, and EMP as well as that of the SSP algorithm, but the computation costs of EMP are quite lower than SSP. Therefore, the performances of the proposed MMP are better than other approaches.

The paper is organized as follows. Section “The GIW-PHD filter” briefly introduces the GIW-PHD filter. Section “Track management for the GIW-PHD filter” presents a track labeling scheme for the GIW-PHD filter. The MPP approach is presented in the section “MPP algorithm.” Simulation results can be found in the section “Simulation results.” Conclusions are given in the section “Conclusions.”

The GIW-PHD filter

Implementation of GIW-PHD filter

The GIW-PHD filter is an important implementation of the extended target PHD framework. The states of the extended targets at time k are defined as

$$X_k=\{{\xi }_k^{(i)}{\}}_{i=1}^{N_{x,k}},{\xi }_k^{(i)}=(x_k^{(i)},X_k^{(i)})$$

(1)

where $N_{x,k}$ is the unknown number of targets, ${\xi }_k^{(i)}$ denotes a component, and $x_k^{(i)}$ is the kinematic state and $X_k^{(i)}$ is a symmetric positive definite matrix which represents an ellipsoidal extension. The dynamic model and measurement model are defined as

$$x_k^{(i)}=(f_{k|k-1}\otimes I_d)x_{k-1}^{(i)}+u_k^{(i)}$$

(2a)

$$z_k^{(j)}=(h_{k|k}\otimes I_d)x_k^{(i)}+v_k^{(i)}$$

(2b)

where $A\otimes B$ denotes the Kronecker product of matrices A and B. $f_{k|k-1}$ is a state transition matrix and $h_{k|k}=[1,0,0]$ as by Granström and Orguner. ¹³ $I_d$ denotes a d-dimensional identity matrix, $v_k^{(j)}$ is the measurement noise given by target extension $X_k^{(i)}$ . $u_k^{(i)}$ is the Gaussian distributed process noise.

The update detection PHD is

$$D_{k|k}({\xi }_k)=D_{k|k}^{ND}({\xi }_k)+\sum _{p\mathrm{\angle }{Z}_k}\sum _{W\in p}{D}_{k|k}^D({\xi }_k,W)$$

(3)

where the notation $p\mathrm{\angle }{Z}_k$ means that the summation is taken over all partitions p of the current set $Z_k$ of the measurements. $D_{k|k}^D$ denotes the missed detection PHD and the detection PHD $D_{k|k}^D$ can be approximated by a mixture of GIW distributions as follows:

$$\begin{array}{rl}& D_{k|k}^D({\xi }_k,W)=\sum _{j=1}^{J_k}{w}_{k|k}^{(j)}\mathcal{N}(x_k;m_{k|k}^{(j,W)},P_{k|k}^{(j,W)}\otimes X_k)\\ & \times IW(X_k;v_{k|k}^{(j,W)},V_{k|k}^{(j,W)})\end{array}$$

(4)

where $m_{k|k}^{(j,W)}$ and $P_{k|k}^{(j,W)}\otimes X_k$ are the mean and covariance of the jth Gaussian distribution. $V_{k|k}^{(j,W)}$ and $V_{k|k}^{(j,W)}$ are the degrees of freedom and scale matrix and of the jth GIW distribution, respectively. The updated GIW component weight is

$$w_{k|k}^{(j)}=\frac{w_p}{d_W}{e}^{-{\gamma }^{(j)}}{\left(\frac{{\gamma }^{(j)}}{{\beta }_{FA,k}}\right)}^{|W|}{p}_D{L}_k^{(j,W)}{w}_{k|k-1}^{(j)},$$

(5a)

$$d_W={\delta }_{|W|,1}+\sum _{h=1}^{J_{k|k-1}^{}}{e}^{-{\gamma }^{(h)}}{\left(\frac{{\gamma }^{(h)}}{{\beta }_{FA,k}}\right)}^{|W|}{p}_D{L}_k^{(h,W)}{w}_{k|k-1}^{(h)},$$

(5b)

$$w_p=\frac{{\prod }_{W\in p}{d}_W}{{\sum }_{p^{\mathrm{\prime }}\mathrm{\angle }{Z}^{\mathrm{\prime }}}{\prod }_{W^{\mathrm{\prime }}\in p^{\mathrm{\prime }}}{d}_{W^{\mathrm{\prime }}}},$$

(5c)

$${\delta }_{|W|,1}=\{\begin{array}{ccc}1& \mathrm{if}& |W|=1\\ 0& \mathrm{if}& |W|\ne 1\end{array}$$

(5d)

where ${\beta }_{FA,k}$ is the rate parameter of the clutter measurements per surveillance volume per scan and ${\gamma }^{(j)}$ is the mean number of the measurements produced by a target. ${\delta }_{|W|,1}$ is the Kronecker delta and $p_D$ is the detection probability. $L_k^{(j,W)}$ denotes the likelihood function of the jth component and cell W. $d_W$ denotes the probability intensity of cell W, and w_p is the weight of a partition, p.

Since we have no contribution to this part, only the basic methods and functions are introduced in this section. The details of the GIW-PHD filter can be found by Granström and Orguner. ¹³

Measurement partitioning approaches

Measurement partitioning is an integral part of a GIW-PHD filter. A division of the set of measurements, Z, into nonempty subsets or cells, W is the definition of a partition, p. ¹⁰ Measurement partitioning of closely spaced targets is a major issue of METT, and can improve the performance of a GIW-PHD filter directly.

There are three partitioning approaches used in the GIW-PHD filter, respectively DP-SP, EMP, and PP, and each method has its disadvantages. For example, the DP-SP does not consider the target extension information, thus it cannot provide accurate partition results if the target extensions are different. The EMP approach may converge to a local maximum of the likelihood function that yields a non-informative partition. The PP is sensitive to target maneuvers because of the use of target prediction information.

The PP approach uses the target prediction information to divides measurements into clusters by

$$(z_{k+1}^{(i)}-m_{k+1|k}^{(j),d}{)}^T(X_{k+1|k}^{(j)}{)}^{-1}(z_{k+1}^{(i)}-m_{k+1|k}^{(j),d})<{\mathrm{\Delta }}_d(p)$$

(6)

where $z_{k+1}^{(i)}$ is a measurement and $X_{k+1|k}^{(j)}$ is the predicted target extension matrix. Here, $m_{k+1|k}^{(j),d}$ denotes the d-dimensional target predicted position at time k + 1, and d denotes the dimension of the physical space. ${\mathrm{\Delta }}_d(p)$ is computed by using the inverse cumulative ${\chi }^2$ distribution with d degrees of freedom for probability p = 0.99. If a measurement can be divided into two clusters, chose the cluster whose component has the highest weight.

If the accuracy of the target predicted information (i.e. $m_{k+1|k}^{(j),d}$ ) can be improved, the PP can provide more accurate partitions to improve the performance of the GIW-PHD filter. Thus the MPP approach has presented in the section “MPP algorithm.”

Track management for the GIW-PHD filter

The GIW-PHD filter cannot provide trajectories of individual targets. This section proposes a track labeling approach for the GIW-PHD filter using the component labeling technique. ¹⁹ The targets can be labeled by the following steps.

Prediction

Supposing that at time $k-1(k\ge 2)$ , the labels of GIW components are denoted by

$$L_{k-1}=\{l_{k-1}^{(j)}{\}}_{j=1}^{J_{k-1}}=\{l_{k-1}^{(1)},l_{k-1}^{(2)}\cdots l_{k-1}^{(J_{k-1})}\}$$

(7)

where $J_{k-1}$ is the number of GIW components at time $k-1$ and $l_{k-1}^{(j)}$ denotes the label of the jth GIW component. The labels of predicted GIW components can be denoted by

$$L_{k|k-1}=L_{k-1}\cup L_{\beta }$$

(8a)

$$L_{\beta }=\{l_{\beta }^{(i)}{\}}_{i=1}^{J_{\beta }}=\{l_{\beta }^{(1)},l_{\beta }^{(2)}\cdots l_{\beta }^{(J_{\beta })}\}$$

(8b)

where $l_{\beta }^{(i)}$ denotes the label of birth GIW components and $J_{\beta }$ is the number of them.

Update:

Suppose ${\hat{N}}_k$ denotes the number of all cells of all partitions, that is ${\hat{N}}_k={\sum }_{i=1}^{{\hat{N}}_{p,k}}{\hat{N}}_{W,k}^{(i)}$ , where ${\hat{N}}_{p,k}$ is the number of partitions at time k, and ${\hat{N}}_{W,k}^{(i)}$ is the number of the cells of the $i$ th partition. According to (2), there would be ${\hat{N}}_k+1$ times the number of predicted GIW components, thus each predicted component gives rise to ${\hat{N}}_k+1$ corresponding updated components. The corresponding labels can be expressed as

$$L_{k|k}=L_{k|k-1}\cup \{L_{k|k-1}^{(i)}{\}}_{i=1}^{{\hat{N}}_k}$$

(9)

where $L_{k|k-1}^{(i)}=L_{k|k-1}$ .

Pruning and merging:

When the GIW components with low weight are pruned, the corresponding labels and their attributes should be also pruned. Note that if several components $\{\cdots ,{\xi }_{k|k}^{(n)},\cdots \}$ should be merged and ${\xi }_{k|k}^{(n)}$ has the largest weight, the label of the merged component is equal to the label of ${\xi }_{k|k}^{(n)}$ . Since the pruning and merging steps of a component are the same as that in the standard GIW-PHD filter, the details can be found by Granström and Orguner^{[ 13 , TABLE I]}.

MPP algorithm

Definitions and key methods of MPP

Definition 1: According to (5), the measurement cell ${\hat{W}}_{k|k}^{(j)}$ produced by partition p is the corresponding measurement cell of a component ${\xi }_{k|k}^{(j)}$ , if that $-$ the partition weight $w_p$ of the partition p is the largest; $-$ ${\hat{W}}_{k|k}^{(j)}$ get the largest likelihood $L_k^{(j,W)}$ than other cells produced by partition p.

Definition 2: ${\xi }_{k+1|k}^{(j)}$ is a closely spaced predicted component of a closely spaced component set $C_{k+1|k}^{(h)}$ at time k + 1, if there is at least one distance between ${\xi }_{k+1|k}^{(j)}$ and other components smaller than a given threshold $\vartheta$ . The set of the $h$ th group closely spaced predicted components is the closely spaced component set $C_{k+1|k}^{(h)}$ .

Definition 3: Measurement cell ${\hat{W}}_{k+1|k}^{(j)}$ is the corresponding measurement cell of a predicted component ${\xi }_{k+1|k}^{(j)}$ , if

$${\hat{W}}_{k+1|k}^{(j)}=\{z_{k+1}^{(i)}|{\eta }_{i,j}<\vartheta {\}}_{z_{k+1}^{(i)}\in Z_{k+1}}$$

(10a)

$${\eta }_{i,j}=\sqrt{{(z_{k+1}^{(i)}-m_{k+1|k}^{(j),d})}^T(z_{k+1}^{(i)}-m_{k+1|k}^{(j),d})}$$

(10b)

where ${\eta }_j$ denotes the distance between $z_{k+1}^{(i)}$ and $m_{k+1|k}^{(j),d}$ . $\vartheta$ is a given distance threshold.

According to (6), the problems of the PP are that:

Even though the consideration of ${\mathrm{\Delta }}_d(p)$ can reduce the effects of the prediction error, the PP will provide inaccurate partitions if the error is large enough, for example targets are performing maneuvers.

Moreover, if targets are spaced closely together, the method that assigning the measurement based on the weight of component will cause errors, for example all measurements of closely spaced targets might be assigned to a single cluster whose component has the highest weight.

The key method of the MPP is to reduce the error of $m_{k+1|k}^{(j),d}$ using the measurement cells ${\hat{W}}_{k|k}^{(j)}$ and ${\hat{W}}_{k+1|k}^{(j)}$ , and to modify formula (6) to improve the accuracy of the partition results.

Implementation of the MPP

For predicted GIW components with the weight $w_{k+1|k}^{(j)}>0.5$ , formula (6) will be rewritten as

$$\begin{array}{r}(z_{k+1}^{(i)}-m_{k+1|k}^{\mathrm{\prime }}(j),d{)}^T(X_{k+1|k}^{(j)}{)}^{-1}\\ (z_{k+1}^{(i)}-m_{k+1|k}^{\mathrm{\prime }}(j),d)<{\mathrm{\Delta }}_d(p)\end{array}$$

(11)

where $m_{k+1|k}^{\mathrm{\prime }}(j),d$ is the modified predicted position. Note that if a measurement can be divided into two or more clusters, chose the cluster whose component resulting in the smallest Mahalanobis distance, that is the value of the left part of (11) is the smallest. This modification can avoid the result that all measurements of the closely spaced targets are assigned to a single cluster.

According to Definition 2, if a component ${\xi }_{k+1|k}^{(j)}$ is not a closely spaced predicted component, $m_{k+1|k}^{\mathrm{\prime }}(j),d$ will be modified by

$$m_{k+1|k}^{\mathrm{\prime }}(j),d=\frac{1}{|{\hat{W}}_{k+1|k}^{(j)}|}\sum _{z_{k+1}^{(i)}\in {\hat{W}}_{k+1|k}^{(j)}}{z}_{k+1}^{(i)}$$

(12)

In fact, $m_{k+1|k}^{\mathrm{\prime }}(j),d$ is the position mean of cell ${\hat{W}}_{k+1|k}^{(j)}$ , which can avoid a significant error produced by target maneuvers. However, if ${\xi }_{k+1|k}^{(j)}\in C_{k+1|k}^{(h)}$ , formula (10) cannot ensure that ${\hat{W}}_{k+1|k}^{(j)}$ only contains the measurements produced by ${\xi }_{k+1|k}^{(j)}$ , thus we use all components $\{{\xi }_{k+1|k}^{(j)}\}$ of $C_{k+1|k}^{(h)}$ and their corresponding $\{{\hat{W}}_{k|k}^{(j)}\}$ and $\{{\hat{W}}_{k+1|k}^{(j)}\}$ to modify $\{m_{k+1|k}^{\mathrm{\prime }}(j),d\}$ collaboratively.

Suppose the corresponding measurements of all components of set $C_{k+1|k}^{(h)}$ are

$$W_{k+1|k}^{(h)}={\bigcup }_{\{j|{\xi }_{k+1|k}^{(j)}\in C_{k+1|k}^{(h)}\}}{\hat{W}}_{k+1|k}^{(j)}$$

(13a)

$$W_{k|k}^{(h)}={\bigcup }_{\{j|{\xi }_{k+1|k}^{(j)}\in C_{k+1|k}^{(h)}\}}{\hat{W}}_{k|k}^{(j)}$$

(13b)

where $W_{k+1|k}^{(h)}$ denotes the union of the corresponding measurement cells of $C_{k+1|k}^{(h)}$ at time k + 1 (gray and black dots in Figure 1). $W_{k|k}^{(h)}$ denotes the union of the corresponding measurement cells of $W_{k+1|k}^{(h)}$ at time k (gray and black squares in Figure 1). Each $m_{k+1|k}^{\mathrm{\prime }}(j),d$ can be modified by

$$m_{k+1|k}^{\mathrm{\prime }}(j),d=A_{k+1}^{(h)}(m_{k+1|k}^{(j),d}-{\overline{m}}_{k+1}^{(h)})+{\overline{z}}_{k+1}^{(h)}$$

(14a)

$${\overline{m}}_{k+1}^{(\ell )}=\frac{1}{|C_{k+1|k}^{(h)}|}{\sum }_{\{j|{\xi }_{k+1|k}^{(j)}\in C_{k+1|k}^{(h)}\}}{m}_{k+1|k}^{(j),d}$$

(14b)

$${\overline{z}}_{k+1}^{(h)}=\frac{1}{|W_{k+1|k}^{(h)}|}{\sum }_{z_{k+1}^{(i)}\in W_{k+1|k}^{(h)}}{z}_{k+1}^{(i)}$$

(14c)

where ${\overline{z}}_{k+1}^{(h)}$ is the position mean of the measurements of set $W_{k+1|k}^{(h)}$ , and ${\overline{m}}_{k+1}^{(h)}$ denotes the position mean of the predicted components of set $C_{k+1|k}^{(h)}$ . $A_{k+1}^{(h)}$ is a coordinate transformation matrix to correct the error produced by targets moving closely or separately. Formula (14) means that $-$ first the origin of coordinates will become the mean of all predicted positions $\{m_{k+1|k}^{(j),d}\}$ by displacement transformation (i.e. $(m_{k+1|k}^{(j),d}-{\overline{m}}_{k+1}^{(h)})$ ); $-$ then the relative position of each $m_{k+1|k}^{(j),d}$ will be transformed by $A_{k+1}^{(h)}$ ; $-$ finally, ${\overline{z}}_{k+1}^{(h)}$ will become the mean position of all $\{m_{k+1|k}^{(j),d}\}$ by displacement transformation.

Figure 1.

A scenario where targets are closely spaced and perform maneuvers.

The transformation matrix $A_{k+1}^{(h)}$ can be approximately calculated by the extension changes between $W_{k|k}^{(h)}$ and $W_{k+1|k}^{(h)}$ (because the positions of closely spaced components will scale or rotate along with the extension change of corresponding measurements), that is the transformation equation

$${\overline{X}}_{k+1}^{(h)}=A_{k+1}^{(h)}{\overline{X}}_k^{(h)}(A_{k+1}^{(h)}{)}^T$$

(15a)

$${\overline{X}}_k^{(h)}=\frac{{\sum }_{z_k^{(i)}\in W_{k|k}^{(h)}}{(z_k^{(i)}-{\overline{z}}_k^{(h)})}^T(z_k^{(i)}-{\overline{z}}_k^{(h)})}{|W_{k|k}^{(h)}|-d-1}$$

(15b)

$${\overline{X}}_{k+1}^{(h)}=\frac{{\sum }_{z_{k+1}^{(i)}\in W_{k+1|k}^{(h)}}{(z_{k+1}^{(i)}-{\overline{z}}_{k+1}^{(h)})}^T(z_{k+1}^{(i)}-{\overline{z}}_{k+1}^{(h)})}{|W_{k+1|k}^{(h)}|-d-1}$$

(15c)

where ${\overline{X}}_k^{(\ell )}$ and ${\overline{X}}_{k+1}^{(h)}$ are the extension matrices ¹³ of $W_{k|k}^{(h)}$ and $W_{k+1|k}^{(h)}$ . Here, ${\overline{z}}_k^{(h)}$ is the position mean of $W_{k|k}^{(h)}$ . According to Goldman, ²⁰ geometrically, ${\overline{X}}_k^{(h)}$ can changes to ${\overline{X}}_{k+1}^{(h)}$ by scaling and rotation, thus $A_{k+1}^{(h)}$ can be solved by

$$A_{k+1}^{(h)}=S_{k+1}{\mathrm{\Psi }}_{k+1}$$

(16a)

$$S_{k+1}=\mathrm{diag}\left(\left[\begin{array}{cc}{\hat{\lambda }}_{1,k+1}/{\hat{\lambda }}_{1,k}& {\hat{\lambda }}_{2,k+1}/{\hat{\lambda }}_{2,k}\end{array}\right]\right)$$

(16b)

$${\mathrm{\Psi }}_{k+1}=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]$$

(16c)

$$\theta =⟨{\alpha }_{1,k},{\alpha }_{1,k+1}⟩$$

(16d)

where $S_{k+1}$ is a scaling matrix. $\{{\hat{\lambda }}_{1,k},{\hat{\lambda }}_{2,k}\}$ and $\{{\hat{\lambda }}_{1,k+1},{\hat{\lambda }}_{2,k+1}\}$ are the characteristic roots of ${\overline{X}}_k^{(h)}$ and ${\overline{X}}_{k+1}^{(h)}$ , respectively. ${\mathrm{\Psi }}_{k+1}$ is the rotation matrix with the angle $\theta$ . Here, ${\alpha }_{1,k}$ and ${\alpha }_{1,k+1}$ denote the feature vectors of ${\overline{X}}_k^{(h)}$ and ${\overline{X}}_{k+1}^{(h)}$ , and $⟨{\alpha }_{1,k},{\alpha }_{1,k+1}⟩$ denotes the angle between them.

Figures 1 and 2 show the illumination of MPP. Figure 1 shows a scenario where two closely spaced targets perform maneuvers from time k−1 to k + 1. The target maneuvers lead to the inaccuracy of the predicted positions, thus the errors should be corrected using the measurement information. Figure 2 shows the transformation process for reducing the errors above. Figure 2(a) shows a scenario where the target predicted states are quite inaccurate. Figure 2(b) shows the result of the displacement transformation $(m_{k+1|k}^{(j),d}-{\overline{m}}_{k+1}^{(h)})$ . Then, as shown in Figure 2(c), the predicted positions will be transformed by $A_{k+1}^{(h)}$ to reduce the relative position error between each other (the distance between the triangles becomes further). Finally, as shown in Figure 2(d), all predicted positions transform to their new positions, that is $m_{k+1|k}^{\mathrm{\prime }}(j),d$ , where ${\overline{z}}_{k+1}^{(h)}$ becomes the mean of them.

Figure 2.

Illumination of modified prediction partitioning (MPP). The black squares denote the measurements at time k + 1. The gray triangles denote the predicted target positions. The solid and dotted forks denote the mean of measurements and predicted target positions, respectively.

Simulation results

The target extension is modeled by $U(B_k{\mathrm{X}}_k{B}_k^T)+\mathcal{N}(R_k)$ , ⁶ where ${\mathrm{X}}_k$ is a uniform distributed extension matrix, $B_k$ is a rotation matrix determined by the motion model, and $R_k$ is a Gaussian measurement noise. The number of measurements of each target follows the Poisson distribution.

The parameters of simulated scenarios are given as

$$T_s=1\mathrm{s},\mathcal{S}=4000\times 4000{\mathrm{m}}^2,{\beta }_{FA,k}=6.25\times {10}^{-7}$$

(17a)

$$R_k=diag([1,1]),Q_k=diag([0.5,0.5,0,0])$$

(17b)

where $T_s$ is the sensor scanning interval. $\mathcal{S}$ denotes the surveillance volume with the rate parameter, that is the Poisson mean of clutter measurements is $\mathcal{S}\times {\beta }_{FA,k}=10$ . $R_k$ and $Q_k$ are the covariances of process noise and measurement noise, respectively.

The parameters of GIW-PHD filter are

$$w_0^{}=0.1,v_0^{}=7,\tau =5,{\gamma }^{(j)}=20$$

(18a)

$$V_0^{}=diag([50,50])$$

(18b)

$$P_0^{}=diag([25,100])$$

(18c)

where $w_0$ is the weight of birth GIW component. ${\gamma }^{(j)}$ is the mean number of measurements produced by each target per time. $v_0$ , $V_0$ , and $P_0$ are inverse Wishart degrees of freedom, inverse scale matrix, and Gaussian covariance of birth GIW components, respectively.

Since the errors may occur in successive scans if only one partitioning approach was used, the DP-SP with distance thresholds $d=\{1,3,5,10,15,\dots ,30\}$ will be used as a basal partitioning approach in the simulation. Then the performance of using DP-SP with MPP will be compared with that of the original partitioning approaches, that is the DP-SP with PP and EMP.

Extended target optimal sub-pattern assignment (ET-OSPA)

The optimal sub-pattern assignment (OSPA) distance ²¹ is employed to evaluate the performance of a single-point PHD filter. For METT, an ET-OSPA was presented by Li et al. ¹⁷ with

$$\begin{array}{rl}& {\overline{d}}_p^{(c)}(X,Y)=(\frac{1}{n}(\underset{\pi \in {\prod }_n}{min}\sum _{i=1}^m{d}^{(c)}\cdots \\ & \times {(x_i,{\hat{X}}_i,y_{\pi (i)},{\hat{Y}}_{\pi (i)})}^p+c^p(n-m)){)}^{\frac{1}{p}}\end{array}$$

(19a)

$$d^{(c)}(x_i,{\hat{X}}_i,y_{\pi (i)},{\hat{Y}}_{\pi (i)}{)}^p=(x_i\text{-}{y}_{\pi (i)}{)}^p+RMS{E}_{{\hat{X}}_i,{\hat{Y}}_{\pi (i)}}^p$$

(19b)

$$RMS{E}_{{\hat{X}}_i,{\hat{Y}}_{\pi (i)}}^p=(tr[{({\hat{X}}_i-{\hat{Y}}_{\pi (i)})}^2]{)}^p$$

(19c)

where $RMS{E}_{{\tilde{X}}_i,{\tilde{Y}}_{\pi (i)}}^p$ denotes the extension error. $\{{\hat{X}}_i\}$ and $\{{\hat{Y}}_{\pi (i)}\}$ are estimated and true extension matrix respectively. The first part of (19b) is equal to the OSPA algorithm, namely the localization error. The added second part is the extension error, therefore the error function $d^{(c)}$ involves whole localization and extension errors.

Results

The target extension matrices were X⁽¹⁾ = diag([20, 4]) and X⁽²⁾ = diag([10, 2]). The targets moved closer from 0 to 20 s and moved linearly together from 21 to 40 s. They performed maneuvers from 41 to 60 s and moved linearly together again from 61 to 80 s. Finally, they separated at 81 s. The tracks are shown in Figure 3.

Figure 3.

True tracks.

Figure 4 shows the average results of 100 Monte Carlo runs. In Figure 4(a), the ET-OSPA values of using just DP-SP are significantly higher than that of other approaches, because DP-SP cannot handle the cases in which targets of various sizes are spaced closely together. The ET-OSPA values of using DP-SP with PP and EMP are lower than using just DP-SP when targets moved linearly together because EMP and PP consider the target extension information. However, when targets performed maneuvers from 41 to 60 s, the ET-OSPA values of using DP-SP with EMP and PP are roughly similar to that of using just DP-SP. The reason for this is that the target maneuvers lead to the complete failure of PP and EMP, thus the use of PP and EMP will not affect the results in these cases. The ET-OSPA values of DP-SP with MPP are significantly lower than other partitioning approaches from 41 to 60 s because MPP can handle the cases of target maneuvers by improving the prediction information. Similar results of the sum of weights can be also seen in Figure 4(b).

Figure 4.

Results of 100 Monte Carlo runs: (a) average extended target optimal sub-pattern assignment (ET-OSPA) values and (b) sum of weights.

Figure 5 shows the association results of a single trial. Figure 5(a) shows the results of using DP-SP with PP and EMP. When targets were maneuvering, a lot of incorrect and excrescent tracks occurred. The reason for this is that DP-SP with PP and EMP failed to provide accurate partitions in such a case. Then, the inaccurate partitions led to the estimation error of the filter, which resulted in the association error. Figure 5(b) shows the results of using DP-SP with MPP, which are relatively better than that in Figure 5(a). It means that even though an error may occur during an individual scan when spatially close targets performed maneuvers, MPP can generally provide better partitions than other partitioning approaches.

Figure 5.

Association results of a single trial: (a) using DP-SP with PP and EMP and (b) using DP-SP with MPP.

Figure 6 shows the average computation costs of 100 Monte Carlo runs of the GIW-PHD filter using different partitioning approaches. The computation cost of using DP-SP with MPP is larger than using just DP-SP because using DP-SP with MPP provides more partitions than using just DP-SP. Also, the computation cost of using DP-SP with MPP is smaller than using DP-SP with PP and EMP, because the use of DP-SP with PP and EMP leads to the largest number of partitions.

Figure 6.

Average computation costs of 100 Monte Carlo runs.

Conclusions

This paper presents a measurement partitioning approach, called MPP, for the GIW-PHD filter. Compared with other approaches, MPP can provide more accurate partitions with a relatively small computation cost when targets are closely spaced. In future works, we plan to apply MPP to the particle PHD filters to solve the problem of dividing particle swarms. The more accurate results of particle swarms can improve the performance of a particle PHD filter.