Innovative unscented transform–based particle cardinalized probability hypothesis density filter for multi-target tracking

Abstract

Multi-target tracking is widely applied in video surveillance systems. As we know, although the standard particle cardinalized probability hypothesis density filter can estimate state of targets, it is difficult to define the proposal distribution function in prediction stage. Since the robust particles cannot be effectively drawn, the actual tracking accuracy should be enhanced. In this paper, an innovative unscented transform–based particle cardinalized probability hypothesis density filter is derived. Considering the different state spaces, we use the auxiliary particle method and then draw robust particles from the modified distributions in order to estimate the position of targets. Simultaneously, we present the recursion of the optimized Kalman gain to improve the general unscented transform for the velocity estimates. Using the track label, we further integrate them in the framework of the jump Markov model. The simulation results show that the proposed filter has advances in the multi-target tracking scenes. Moreover, the experiments indicate that the filter can track mobile targets with satisfactory results.

Introduction

Multi-target tracking is to estimate both state and cardinality of mobile targets using available measurements contaminated by random clutters.^1–3 With respect to high nonlinearity and varying times in actual video surveillance field, multi-target tracking has become an important technology for autonomous systems.

The measurement uncertainties and filtering efficiency are always our concerns. Some scholars have researched multi-target tracking method in video surveillance system with prospective success, and many articles have been published in famous journals. First, the Kalman filter (KF), for example, predicts the future motion state with linear transformation. ⁴ As we know, the state model is generally linear, but the measurement model is frequently nonlinear. ⁵ With respect to the nonlinear system, the unscented Kalman filter (UKF) easily achieves the second-order estimation accuracy avoided the Jacobian matrix computation,^6,7 where the unscented transformation (UT) scales sigma points to incorporate more information into state estimation. Thus, the UKF becomes popular in the research community. ⁸ This kind of filter has satisfactory tracking performance only under certain assumptions. Once the assumptions are unexpected, its tracking performance severely declines. To deal with this problem, an adaptive UKF for actual state and parameter estimation was derived in Song. ⁹ Subsequently, a new adaptive UKF and its typical strategy were summarized in Pan et al. ¹⁰

With the rapid development of the random finite set theory, the cardinalized probability hypothesis density (CPHD) filter has been explored in the past decades. ¹¹ It can propagate the cardinality distribution of targets to probability hypothesis density (PHD) with the second-order moment by estimating cardinality, which is a popular solution to the multi-target tracking problem in cluttered environments. ¹² In Mahler, ¹³ a Gaussian mixture cardinalized probability hypothesis density (GMCPHD) filter was used to track mobile targets with the linear and Gaussian dynamics. An improved GMCPHD filter was proposed to track passive targets. ¹⁴ In view of the nonlinear and non-Gaussian state estimation, the particle cardinalized probability hypothesis density (PCPHD) filter is more efficient, which achieves the closed-form solution under the sequential Monte Carlo runs. In Ristic et al., ¹⁵ an extension of PCPHD filter was presented to conveniently distinguish the survival and newborn targets. In Li et al., ¹⁶ a PCPHD filter was presented based on the weight–optimization strategy. Recently, a new PCPHD filter was proposed based on the adaptive genetic method to boost estimation accuracy under the condition of the motion switching time. ¹⁷ In view of extended multi-target tracking, replacing traditional multi-measurement with a rectangular region of nonzero volume in the state space, a labeled box-PCPHD filter was presented in Zou et al., ¹⁸ which can propagate intensity function as well as cardinality distribution. In Lu et al., ¹⁹ a new PCPHD filter was proposed to complete efficient track continuity and extraction. The main idea of the mentioned filters is based on the importance sampling which limits the actual tracking performance to some extent owing to the difficulty of particle degradation. In the prediction stage, they sample the particles on the common proposal distribution. Meanwhile, it is arbitrary to draw particles from an uncertain and unknown distribution. By comparison, the auxiliary particle (AP) method allows some particles to be propagated from the effective distributions with high probability.^20,21 What is more, the single PCPHD filter cannot well complete the multi-target tracking in actual applications. Usually, one mobile target has different motion stages, or some heterogeneous targets have different motion states during the surveillance time. In order to carter the situations, the jump Markov model is more competent, which matches different dynamics under the parallel filters.^22,23 Characterizing the various maneuvers, the corresponding model is necessary to represent the multi-target motion. ²⁴ Recently, a robust and fast-tracking algorithm based on two-layer jump Markov model for curvilinear maneuvering target tracking was presented in Cai and Sun. ²⁵

Inspired by the technical specifications for multi-target tracking, this note presents an innovative UT-based PCPHD filter with its derivation. It mainly addresses the double questions: How to design an effective CPHD filter for drawing the required particles from modified proposal distribution? How to combine the enhanced CPHD filter with the enhanced UT method in the nonlinear and linear state spaces individually? Therefore, the note distinctly exhibits a novel filtering method for drawing robust particles to approximate target state based on the current measurement information. The target state is partitioned into two different spaces. The target position estimation is determined by the sampling particles drawn from the proposal distribution with an extra particle and current measurements. The target velocity estimation is just on the enhanced UKF using the optimized Kalman gain. In the framework of the jump Markov model, the heterogeneous filters are systematically integrated to track the mobile targets in an actual surveillance system.

The organization of this note is assigned as follows: In section “Background,” the basic principle and some limitations of the standard CPHD filter are briefly formulated. In section “Methodology,” we propose the improvements of the UT-based PCPHD filter and illustrate its filtering process. Subsequently, the numerical study and actual experiment are indicated with promising results to verify the overall performance of the proposed filter in section “Simulation results and discussions.” In the final section, we come to the conclusions and provide the next research plan.

Background

Suppose that ${\mathbf{x}}_{k-1}$ is the target state vector in the state space $\mathcal{X}$ at time $k-1$ , the target state vector at time k is given by^26–29

$${\mathbf{x}}_k={\mathbf{F}}_{k|k-1}{\mathbf{x}}_{k-1}+{\mathbf{\Gamma }}_k{\mathbf{u}}_k$$

(1)

where ${\mathbf{F}}_{k|k-1}$ is the state transfer matrix and ${\mathbf{\Gamma }}_k$ is the gain matrix of state noise vector ${\mathbf{u}}_k$ with the Gaussian distribution $\mathcal{N}(0,{\mathbf{Q}}_k)$ .

The measurement model in the measurement space $\mathcal{Z}$ at time k can be written as

$${\mathbf{z}}_k={\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{e}}_k$$

(2)

where ${\mathbf{z}}_k$ is the measurement matrix, ${\mathbf{H}}_k$ is the measurement transfer matrix, and ${\mathbf{e}}_k$ is the measurement noise vector with the Gaussian distribution $\mathcal{N}(0,{\mathbf{R}}_k)$ .

Remark 1

At each time, some targets may disappear, some targets may survive, and some targets may appear. Besides, the evolution of targets presents and the origin of current measurements are not often known. Considering the uncertainty and randomness characterized by equations (1) and (2), we will apply the CPHD filter to track multi-target.

Principle of standard CPHD filter

Suppose that $D_{k-1}({\mathbf{x}}_{k-1})$ is the posterior PHD at time $k-1$ , its prediction is given by^30,31

$$\begin{array}{c}\multicolumn{1}{c}{D_{k|k-1}({\mathbf{x}}_k)}\\ \multicolumn{1}{c}{=b_{k|k-1}({\mathbf{x}}_k)+\int p_{S_k|k-1}({\mathbf{x}}_{k-1})f_{k|k-1}({\mathbf{x}}_k|{\mathbf{x}}_{k-1})}\\ \multicolumn{1}{c}{D_{k-1}({\mathbf{x}}_{k-1})\mathrm{d}{\mathbf{x}}_{k-1}}\end{array}$$

(3)

where $p_{S{,}_{k|k-1}}({\mathbf{x}}_{k-1})$ is the target–survival probability, $f_{k|k-1}({\mathbf{x}}_k|{\mathbf{x}}_{k-1})$ is the Markov transfer probability, and $b_{k|k-1}({\mathbf{x}}_k)$ is the target–birth intensity.

Define the cardinality distribution ${\rho }_B(n)$ of target–birth

$${\rho }_B\left(n\right)=\int b_{k|k-1}\left({\mathbf{x}}_k\right)\mathrm{d}{\mathbf{x}}_k$$

(4)

The cardinality distribution is predicted as

$$\begin{array}{c}\multicolumn{1}{c}{{\rho }_{k|k-1}\left(n\right)=\sum _{i=0}^n{\rho }_B(n-i)}\\ \multicolumn{1}{c}{\times \sum _{n=\gamma }^{\infty }\left(\begin{array}{c}\frac{{〈p_{{S,}_{k|k-1}}\left({\mathbf{x}}_{k-1}\right),D_{k-1}\left({\mathbf{x}}_{k-1}\right)〉}^{\gamma }}{{〈1,D_{k-1}\left({\mathbf{x}}_{k-1}\right)〉}^n}{C}_n^{\gamma }{\rho }_{k-1}\left(n\right)\hfill \\ \multicolumn{1}{c}{\times {〈1-p_{{S,}_{k|k-1}}\left({\mathbf{x}}_{k-1}\right),D_{k-1}\left({\mathbf{x}}_{k-1}\right)〉}^{n-\gamma }}\end{array}\right)}\end{array}$$

(5)

where ${\rho }_{k|k-1}(n)$ is predicted cardinality distribution, $〈·,·〉$ denotes the inner product of two elements, and $C_n^{\gamma }$ is the binomial coefficient. ${\rho }_{k|k-1}(n)$ is a convolution of the cardinality distribution of the survival and newborn targets.

Subsequently, the cardinality distribution at time k is updated as

$${\rho }_k\left(n\right)=\frac{{\Upsilon }_k^0\left({\mathbf{Z}}_k\right){\rho }_{k|k-1}\left(n\right)}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}\left(n\right)〉}$$

(6)

Note that ${\rho }_k(n)$ incorporates the clutter cardinality, the measurement information, the predicted intensity and predicted cardinality distribution. Suppose that $p_D$ is the detection probability of a passive sensor and $g_k({\mathbf{z}}_k|{\mathbf{x}}_k)$ is the target likelihood function, the PHD is updated as

$$\begin{array}{c}\multicolumn{1}{c}{D_k\left({\mathbf{x}}_k\right)=\left(\begin{array}{c}(1-p_D)\frac{〈{\Upsilon }_k^1\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}\left(n\right)〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}\left(n\right)〉}\hfill \\ \multicolumn{1}{c}{+p_D\sum _{z\in {\mathbf{Z}}_k}{g}_k({\mathbf{z}}_k|{\mathbf{x}}_k)\frac{〈{\Upsilon }_k^1\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}\left(n\right)〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}\left(n\right)〉}}\end{array}\right)}\\ \multicolumn{1}{c}{D_{k|k-1}\left({\mathbf{x}}_k\right)}\end{array}$$

(7)

where the parameter ${\Upsilon }_k^u$ is

$$\begin{array}{c}\multicolumn{1}{c}{{\Upsilon }_k^u\left({\mathbf{Z}}_k\right)=\sum _{m=0}^{min(\left|{\mathbf{Z}}_k\right|,n)}}\\ \multicolumn{1}{c}{\left(\begin{array}{c}(\left|{\mathbf{Z}}_k\right|-m)!p_{c,k}(\left|{\mathbf{Z}}_k\right|-m)e_m\left({\Lambda }_k\left({\mathbf{Z}}_k\right)\right)\hfill \\ \multicolumn{1}{c}{\times \frac{P_{m+u}^n{〈1-p_D,D_{k|k-1}\left({\mathbf{x}}_k\right)〉}^{n-m-u}}{{〈1,D_{k|k-1}\left({\mathbf{x}}_k\right)〉}^n}}\end{array}\right)}\end{array}$$

(8)

where $p_{c,k}(·)$ is the cardinality distribution of random clutters, $e_m(·)$ is the symmetric function, $P_m^n$ is the permutation coefficient, and ${\Lambda }_k({\mathbf{Z}}_k)$ is partially determined by the clutter process $c_k({\mathbf{z}}_k)$

$$\begin{array}{c}\multicolumn{1}{c}{{\Lambda }_k\left({\mathbf{Z}}_k\right)}\\ \multicolumn{1}{c}{=〈D_{k|k-1}\left({\mathbf{x}}_k\right),\frac{〈1,c_k\left({\mathbf{z}}_k\right)〉}{c_k\left({\mathbf{z}}_k\right)}{p}_D{g}_k({\mathbf{z}}_k|{\mathbf{x}}_k):{\mathbf{z}}_k\in {\mathbf{Z}}_k〉}\end{array}$$

(9)

Thus, $D_k({\mathbf{x}}_k)$ is the PHD update that has one missed detection term and $\left|{\mathbf{Z}}_k\right|$ detection terms.

Using the maximum a posteriori rule, we estimate the cardinality of targets as

$${\hat{N}}_k=\arg \underset{n}{sup}{\rho }_k\left(n\right)$$

(10)

where $\arg sup(·)$ is the least upper bound. ${\hat{N}}_k$ maximal peaks $\{{\hat{\mathbf{x}}}_k^{(1)},{\hat{\mathbf{x}}}_k^{(2)},\dots ,{\hat{\mathbf{x}}}_k^{({\hat{N}}_k)}\}$ are extracted to estimate state of targets.

Problem formation

In the standard CPHD filter, the proposal distribution is often arbitrary or uncertain in prediction stage so that some particles gradually degenerate. Although they reflect the target dynamics, the actual approximation is unsatisfactory. The velocity component and position component should be individually estimated, that is, the UKF replace with the particle filter (PF) to get velocity estimates. In the common UKF, the Kalman gain needs to be optimized when estimating velocity. Besides, the number of sigma points determines tracking efficiency. As a result, the sampling of the non-local effects should be deep in thought. For these considerations, the enhanced PCPHD filter and enhanced UKF should be integrated in the jump Markov model.

Methodology

Filtering framework

Since the dimension of states is higher than the dimension of measurements, many sampling particles are required to approximate the nonlinear and linear components.^32,33 In order to avoid the mutual interference, the velocity and position components should be estimated individually. First, we distinguish them in the low-dimensional state space. A necessary truncation is performed to complete a finite-dimensional design, where the UT method is applied. As a resume, it improves the velocity estimation accuracy and then saves the computational cost compared with the single PCPHD filter.

We partition the whole state vector ${\mathbf{x}}_k$ into the nonlinear component ${\mathbf{x}}_k^{\mathrm{n}}$ and linear component ${\mathbf{x}}_k^{\mathrm{l}}$ at time k. In view of the measurement set ${\mathbf{Z}}_{0:k}$ , we have in hand the posterior probability density based on Bayesian rule^34,35

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathbf{x}}_{0:k}|{\mathbf{Z}}_{0:k})}\\ \multicolumn{1}{c}{=p({\mathbf{x}}_{0:k}^{\mathrm{n}},{\mathbf{x}}_k^{\mathrm{l}}|{\mathbf{Z}}_{0:k})}\\ \multicolumn{1}{c}{=p({\mathbf{x}}_k^{\mathrm{l}}|{\mathbf{x}}_{0:k}^{\mathrm{n}})p({\mathbf{x}}_{0:k}^{\mathrm{n}}|{\mathbf{Z}}_{0:k})}\end{array}$$

(11)

where ${\mathbf{x}}_k^{\mathrm{n}}$ contains more state information determined by ${\mathbf{Z}}_{0:k}$ . Since that the probability of ${\mathbf{x}}_k^{\mathrm{l}}$ is conditional on ${\mathbf{x}}_k^{\mathrm{n}}$ , $p({\mathbf{x}}_k^{\mathrm{l}}|{\mathbf{x}}_{0:k}^{\mathrm{n}})$ follows the Gaussian distribution

$$p\left({\mathbf{x}}_k^{\mathrm{l}}|{\mathbf{x}}_{0:k}^{\mathrm{n}}\right)=N({\mathbf{x}}_{k|k-1}^{\mathrm{l}},{\mathbf{P}}_{k|k-1})$$

(12)

where the predicted linear state ${\mathbf{x}}_{k|k-1}^{\mathrm{l}}$ and its covariance ${\mathbf{P}}_{k|k-1}$ are computed under the UT recursion.

We present a lemma of the optimized KF gain by defining a set of $\vartheta -\mathrm{dimensional}$ sigma points $\{{\mathbf{\varsigma }}_{k-1}^{(\chi )}{\}}_{\chi =0}^{\vartheta }$ . Moreover, ${\mathbf{v}}_{k-1}^{(\chi )}$ and ${\mathbf{\upsilon }}_{k-1}^{(\chi )}$ are the weights assigned to ${\mathbf{x}}_{k-1}^{\mathrm{l},(i)}$ and ${\mathbf{P}}_{k-1}$ , respectively.

Lemma 1

Let ${\mathbf{P}}_{\mathrm{XZ},k-1}$ be the covariance of ${\mathbf{x}}_{k-1}$ and ${\mathbf{z}}_{k-1}$ ; the optimized KF gain at time k is given by

$$\begin{array}{l}{\mathbf{K}}_k=({\mathbf{P}}_{k-1}{\mathbf{H}}_k^T+{\mathbf{P}}_{XZ,k-1}^T)\\ \times {({\mathbf{H}}_k{\mathbf{P}}_{k-1}{\mathbf{H}}_k^T+{\mathbf{R}}_k+{\mathbf{H}}_k{\mathbf{P}}_{XZ,k-1}+{\mathbf{P}}_{XZ,k-1}^T{\mathbf{H}}_k^T)}^{-1}\end{array}$$

(13)

Proof

We have the predicted covariance of the predicted linear state ${\mathbf{x}}_{k|k-1}^{\mathrm{l}}$ as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{P}}_{k|k-1}=\sum _{\chi =0}^{\vartheta +1}{\mathbf{\upsilon }}_{k-1}^{\left(\chi \right)}({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k|k-1}^{\mathrm{l}}){({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k|k-1}^{\mathrm{l}})}^T}\\ \multicolumn{1}{c}{=\sum _{\chi =0}^{\vartheta +1}{\mathbf{\upsilon }}_{k-1}^{\left(\chi \right)}\left(\begin{array}{c}((\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k-1}^{\mathrm{l}})-{\mathbf{K}}_k({\mathbf{z}}_k^{\left(\chi \right)}-{\mathbf{z}}_{k-1}))\hfill \\ \multicolumn{1}{c}{\times {((\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k-1}^{\mathrm{l}})-{\mathbf{K}}_k({\mathbf{z}}_k^{\left(\chi \right)}-{\mathbf{z}}_{k-1}))}^T}\end{array}\right)}\\ \multicolumn{1}{c}{=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)\sum _{\chi =0}^{\vartheta +1}{\mathbf{\upsilon }}_{k-1}^{\left(\chi \right)}({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k-1}^{\mathrm{l}}){({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k-1}^{\mathrm{l}})}^T}\\ \multicolumn{1}{c}{\times {(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)}^T+{\mathbf{K}}_k\sum _{\chi =0}^{\vartheta +1}{\mathbf{\upsilon }}_{k-1}^{\left(\chi \right)}{\mathbf{e}}_k{\mathbf{e}}_k^T{\mathbf{K}}_k^T}\\ \multicolumn{1}{c}{-(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)\sum _{\chi =0}^{\vartheta +1}{\mathbf{\upsilon }}_{k-1}^{\left(\chi \right)}({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k-1}^{\mathrm{l}}){({\mathbf{z}}_k^{\left(\chi \right)}-{\mathbf{z}}_{k-1})}^T{\mathbf{K}}_k^T}\\ \multicolumn{1}{c}{-{\mathbf{K}}_k\sum _{\chi =0}^{\vartheta +1}{\mathbf{\upsilon }}_{k-1}^{\left(\chi \right)}({\mathbf{z}}_k^{\left(\chi \right)}-{\mathbf{z}}_{k-1}){({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_{k-1}^{\mathrm{l}})}^T{(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)}^T}\\ \multicolumn{1}{c}{=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_{k-1}{(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)}^T+{\mathbf{K}}_k{\mathbf{R}}_k{\mathbf{K}}_k^T}\\ \multicolumn{1}{c}{-(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_{\mathrm{XZ},k-1}{\mathbf{K}}_k^T-{\mathbf{K}}_k{\mathbf{P}}_{\mathrm{XZ},k-1}{(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)}^T}\\ \multicolumn{1}{c}{={\mathbf{P}}_{k-1}-{\mathbf{K}}_k{\mathbf{H}}_k{\mathbf{P}}_{k-1}-{\mathbf{P}}_{k-1}{\mathbf{H}}_k^T{\mathbf{K}}_k^T+{\mathbf{K}}_k({\mathbf{H}}_k{\mathbf{P}}_{k-1}{\mathbf{H}}_k^T+{\mathbf{R}}_k){\mathbf{K}}_k^T}\\ \multicolumn{1}{c}{-{\mathbf{P}}_{\mathrm{XZ},k-1}{\mathbf{K}}_k^T+{\mathbf{K}}_k{\mathbf{H}}_k{\mathbf{P}}_{\mathrm{XZ},k-1}{\mathbf{K}}_k^T-{\mathbf{K}}_k{\mathbf{P}}_{\mathrm{XZ},k-1}+{\mathbf{K}}_k{\mathbf{P}}_{\mathrm{XZ},k-1}^T{\mathbf{H}}_k^T{\mathbf{K}}_k^T}\end{array}$$

(14)

Subsequently, we compute the following partial differential

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial {\mathbf{P}}_{k|k-1}}{\partial {\mathbf{K}}_k}=-2{\left({\mathbf{H}}_k{\mathbf{P}}_{k-1}\right)}^T+2{\mathbf{K}}_k({\mathbf{H}}_k{\mathbf{P}}_{k-1}{\mathbf{H}}_k^T+{\mathbf{R}}_k)}\\ \multicolumn{1}{c}{-2{\mathbf{P}}_{\mathrm{XZ},k-1}^T+2{\mathbf{K}}_k{\mathbf{H}}_k{\mathbf{P}}_{\mathrm{XZ},k-1}+2{\mathbf{K}}_k{\mathbf{P}}_{\mathrm{XZ},k-1}^T{\mathbf{H}}_k^T}\end{array}$$

(15)

Setting equation (15) to 0, the optimized ${\mathbf{K}}_k$ is achieved in equation (13). Recalling equation (11), we also define $p({\mathbf{x}}_{0:k}^{\mathrm{n}}|{\mathbf{Z}}_{0:k})$ as

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathbf{x}}_{0:k}^{\mathrm{n}}|{\mathbf{Z}}_{0:k})}\\ \multicolumn{1}{c}{=\frac{p({\mathbf{Z}}_k,{\mathbf{x}}_{0:k}^{\mathrm{n}}|{\mathbf{Z}}_{0:k-1})}{p({\mathbf{Z}}_k|{\mathbf{Z}}_{0:k-1})}}\\ \multicolumn{1}{c}{=\frac{p({\mathbf{Z}}_k|{\mathbf{x}}_{0:k}^{\mathrm{n}},{\mathbf{Z}}_{0:k-1})p({\mathbf{x}}_{0:k}^{\mathrm{n}}|{\mathbf{Z}}_{0:k-1})}{p({\mathbf{Z}}_k|{\mathbf{Z}}_{0:k-1})}}\\ \multicolumn{1}{c}{=\frac{p({\mathbf{Z}}_k|{\mathbf{x}}_k^{\mathrm{n}})p({\mathbf{x}}_k^{\mathrm{n}}|{\mathbf{x}}_{0:k-1}^{\mathrm{n}})p({\mathbf{x}}_{0:k-1}^{\mathrm{n}}|{\mathbf{Z}}_{0:k-1})}{p({\mathbf{Z}}_k|{\mathbf{Z}}_{0:k-1})}}\end{array}$$

(16)

where $p({\mathbf{x}}_k^{\mathrm{n}}|{\mathbf{x}}_{0:k-1}^{\mathrm{n}})$ is replaced by $p({\mathbf{x}}_k^{\mathrm{n}}|{\mathbf{x}}_{k-1}^{\mathrm{n}})$ because ${\mathbf{x}}_{k-1}^{\mathrm{n}}$ plays a leading role compared with ${\mathbf{x}}_{0:k-2}^{\mathrm{n}}$ .

Initialization stage

At time $k-1$ , we use a particle set $\{{\mathbf{w}}_{k-1}^{(i)},{\mathbf{x}}_{k-1}^{\mathrm{n},(i)},{\tau }_{k-1}^{(i)}{\}}_{i=1}^{L_{k-1}}$ to approximate posterior PHD and cardinality distribution, where ${\tau }_{k-1}^{(i)}$ is the track label. Suppose that $\delta (·)$ denotes the Dirac delta function, the posterior PHD is

$$D_{k-1}(\mathbf{x},\tau |{\mathbf{Z}}_{0:k-1})=\sum _{i=1}^{L_{k-1}}{\mathbf{w}}_{k-1}^{\left(i\right)}\delta (\mathbf{x}-{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},\tau -{\tau }_{k-1}^{\left(i\right)})$$

(17)

Let ${\varsigma }_{k-1}^{(\chi ,\vartheta )}$ be an element of the vector ${\mathbf{\varsigma }}_{k-1}^{(\chi )}$ , we have

$$\begin{array}{cc}{\mathbf{\varsigma }}_{k-1}^{\left(\chi \right)}={\mathbf{x}}_{k-1}^{\mathrm{l},\left(i\right)}+\sqrt{{\alpha }_{k-1}^2{\mathbf{P}}_{k-1}}{\varsigma }_{k-1}^{(\chi ,\vartheta )}\hfill & \chi =0,1,\dots ,\vartheta +1\hfill \end{array}$$

(18)

$${\mathbf{\upsilon }}_{k-1}^{\left(\chi \right)}=\{\begin{array}{cc}3+{\mathbf{v}}_{k-1}^{\left(0\right)}-{\alpha }_{k-1}^2\hfill & \chi =0\hfill \\ \multicolumn{1}{c}{{\mathbf{v}}_{k-1}^{\left(\chi \right)}}& \chi \ne 0\hfill \end{array}$$

(19)

$${\mathbf{v}}_{k-1}^{\left(\chi \right)}=\{\begin{array}{cc}\frac{3+{\mathbf{v}}_{k-1}^{\left(0\right)}}{{\alpha }_{k-1}^2}\hfill & \chi =0\hfill \\ \multicolumn{1}{c}{\frac{1-{\mathbf{v}}_{k-1}^{\left(0\right)}}{2^{\chi }{\alpha }_{k-1}^2}}& \chi =1,2\hfill \\ \multicolumn{1}{c}{\frac{2^{(\chi -2)}{\mathbf{v}}_{k-1}^{\left(1\right)}}{{\alpha }_{k-1}^2}}& \chi =3,4,\dots ,\vartheta +1\hfill \end{array}$$

(20)

where the UT parameter ${\alpha }_{k-1}^2$ determines the square of the spread of sigma points around ${\mathbf{x}}_{k-1}^{\mathrm{l},(i)}$ .

Prediction stage

The cardinality distribution at time $k-1$ is predicted as

$$\begin{array}{c}\multicolumn{1}{c}{{\rho }_{k|k-1}(n,{\tau }_{k-1}^{\left(i\right)})}\\ \multicolumn{1}{c}{=\sum _{m=0}^n{\rho }_{B,k}(n-m,{\tau }_{k-1}^{\left(i\right)})}\\ \multicolumn{1}{c}{\times \sum _{\ell =m}^{\infty }\left(\begin{array}{c}{C}_j^{\ell }\frac{{〈p_{S,k|k-1}\left({\mathbf{x}}_{k-1}^{n,\left(i\right)}\right),{\mathbf{w}}_{k-1}^{\left(i\right)}〉}^m}{{〈1,{\mathbf{w}}_{k-1}^{\left(i\right)}〉}^{\ell }}\hfill \\ \multicolumn{1}{c}{\times {〈1-p_{S,k|k-1}\left({\mathbf{x}}_{k-1}^{n,\left(i\right)}\right),{\mathbf{w}}_{k-1}^{\left(i\right)}〉}^{\ell -m}{\rho }_{k-1}(n,{\tau }_{k-1}^{\left(i\right)})}\end{array}\right)}\end{array}$$

(21)

In equation (21), the basic strategy is to propagate a set of weights and particles predicting the cardinality distribution recursively in time.

Given that the particle weight reflects the ability of particle to approximate the corresponding measurement, the filter should select a robust particle to propagate the cardinality from a certain distribution function. Therefore, we extend the set by inserting an extra particle, where the associated parameters are given by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{w}}_{k-1}^{(L_{k-1}+1)}=\int b_{k|k-1}({\mathbf{x}}_{k-1}^{\mathrm{n},(L_{k-1}+1)}|{\tau }_{k-1}^{(L_{k-1}+1)})\mathrm{d}}\\ \multicolumn{1}{c}{{\mathbf{x}}_{k-1}^{\mathrm{n},(L_{k-1}+1)}}\end{array}$$

(22)

$$\begin{array}{c}\multicolumn{1}{c}{f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},(L_{k-1}+1)}|{\mathbf{x}}_{k-1}^{\mathrm{n},(L_{k-1}+1)},{\tau }_{k-1}^{(L_{k-1}+1)})}\\ \multicolumn{1}{c}{=\frac{b_{k|k-1}({\mathbf{x}}_{k-1}^{\mathrm{n},(L_{k-1}+1)}|{\tau }_{k-1}^{(L_{k-1}+1)})}{{\mathbf{w}}_{k-1}^{(L_{k-1}+1)}}}\end{array}$$

(23)

$$p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},(L_{k-1}+1)}\right)=1$$

(24)

Then, the predicted PHD is

$$\begin{array}{l}{D}_{k|k-1}\left(\mathbf{x},\tau |{\mathbf{Z}}_{0:k-1}\right)=\sum _{i=1}^{L_{k-1}+1}{\mathbf{w}}_{k|k-1}^{\left(i\right)}\delta \\ (\mathbf{x}-{\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)},\tau -{\tau }_{k|k-1}^{\left(i\right)})\end{array}$$

(25)

Note that ${\mathbf{x}}_{k|k-1}^{\mathrm{n},(L_{k-1}+1)}$ is assigned to any point in an extended dimensional space by pre-selecting previous particles that reflect the current measurements.

Let the subscripts $\Gamma$ and $\mathrm{D}$ denote the undetected and detected components, the integral of a given test function $\phi (·)$ in the extended target space $\mathcal{X}\prime$ can be written as follows

$$\begin{array}{c}\multicolumn{1}{c}{\underset{\mathcal{X}\prime }{\int }\phi \left({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}\right){\mathbf{w}}_{k|k-1}^{\left(i\right)}\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}\\ \multicolumn{1}{c}{=\sum _{i=1}^{L_{k-1}+1}\underset{\mathcal{X}\prime }{\int }\phi \left({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}\right)(1-p_D)\frac{〈{\Upsilon }_k^1\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}}\\ \multicolumn{1}{c}{{\mathbf{w}}_{k|k-1,\Gamma }^{\left(i\right)}\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}\\ \multicolumn{1}{c}{+\sum _{j=1}^{\left|{\mathbf{Z}}_k\right|}\sum _{i=1}^{L_{k-1}+1}\underset{\mathcal{X}\prime }{\int }\left(\begin{array}{c}\phi \left({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}\right)p_D{g}_k({\mathbf{z}}_k^{\left(j\right)}|{\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)})\frac{〈1,c_k\left({\mathbf{z}}_k^{\left(j\right)}\right)〉}{c_k\left({\mathbf{z}}_k^{\left(j\right)}\right)}\hfill \\ \multicolumn{1}{c}{\times \frac{〈{\Upsilon }_k^1({\mathbf{Z}}_k-\left\{{\mathbf{z}}_k^{\left(j\right)}\right\}),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}{\mathbf{w}}_{k|k-1,\mathrm{D}}^{\left(i\right)}}\end{array}\right)\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}\end{array}$$

(26)

where the weights of required particles are predicted as

$${\mathbf{w}}_{k|k-1,\Gamma }^{\left(i\right)}=\frac{p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)}){\mathbf{w}}_{k-1}^{\left(i\right)}}{{\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)})}$$

(27)

$${\mathbf{w}}_{k|k-1,\mathrm{D}}^{\left(i\right)}=\frac{p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)}){\mathbf{w}}_{k-1}^{\left(i\right)}}{{\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},j)}$$

(28)

In equation (27), the proposal distribution of the undetected component ${\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)})$ is

$${\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)})={\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}){\mathrm{d}}_{\Gamma }\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)$$

(29)

where ${\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)})$ and ${\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)})$ are given by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{n,\left(i\right)})}\\ \multicolumn{1}{c}{=(1-p_D)\frac{〈{\Upsilon }_k^1\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}}\\ \multicolumn{1}{c}{\times \frac{\sum _{i=1}^{L_{k-1}+1}{\mathbf{w}}_{k-1}^{\left(i\right)}\delta (\mathbf{x}-{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)})\underset{\mathcal{X}\prime }{\int }{p}_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)})\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}{\sum _{i=1}^{L_{k-1}+1}{\mathbf{w}}_{k-1}^{\left(i\right)}\underset{\mathcal{X}\prime }{\int }{p}_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)})\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}}\end{array}$$

(30)

$$\begin{array}{l}{\mathrm{d}}_{\Gamma }\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)=\sum _{i=1}^{L_{k-1}+1}{\mathbf{w}}_{k-1}^{\left(i\right)}\underset{{\mathcal{X}}^{\prime }}{\int }{p}_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)f_{k|k-1}\\ \left({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)}\right)\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\end{array}$$

(31)

In equation (28), the proposal distribution of the detected component ${\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)},{\mathbf{x}}_{k-1}^{n,\left(i\right)},j)$ is

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},j)={\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},j)}\\ \multicolumn{1}{c}{{\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}|j){\mathrm{d}}_{\mathrm{D}}\left(j\right)}\end{array}$$

(32)

where ${\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},j)$ , ${\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}|j)$ , and ${\mathrm{d}}_{\mathrm{D}}(j)$ are given by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},j)}\\ \multicolumn{1}{c}{=\frac{p_D\sum _{i=1}^{L_{k-1}+1}{\mathbf{w}}_{k-1}^{\left(i\right)}\delta (\mathbf{x}-{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)})\underset{\mathcal{X}\prime }{\int }\left(\begin{array}{c}{g}_k({\mathbf{z}}_k^{\left(j\right)}|{\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)})p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)\hfill \\ \multicolumn{1}{c}{\times f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)})}\end{array}\right)\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}{\sum _{i=1}^{L_{k-1}+1}{\mathbf{w}}_{k-1}^{\left(i\right)}\underset{\mathcal{X}\prime }{\int }\left(\begin{array}{c}{g}_k({\mathbf{z}}_k^{\left(j\right)}|{\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)})p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)\hfill \\ \multicolumn{1}{c}{\times f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)})}\end{array}\right)\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}}\end{array}$$

(33)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}|j)=\sum _{i=1}^{L_{k-1}+1}{\mathbf{w}}_{k-1}^{\left(i\right)}\underset{\mathcal{X}\prime }{\int }}\\ \multicolumn{1}{c}{\left(\begin{array}{c}{g}_k({\mathbf{z}}_k^{\left(j\right)}|{\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)})p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}\right)\hfill \\ \multicolumn{1}{c}{\times f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)},{\tau }_{k-1}^{\left(i\right)})}\end{array}\right)\mathrm{d}{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i\right)}}\end{array}$$

(34)

$${\mathrm{d}}_{\mathrm{D}}\left(j\right)=\frac{〈1,c_k\left({\mathbf{z}}_k^{\left(j\right)}\right)〉}{c_k\left({\mathbf{z}}_k^{\left(j\right)}\right)}\frac{〈{\Upsilon }_k^1({\mathbf{Z}}_k-\left\{{\mathbf{z}}_k^{\left(j\right)}\right\}),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}$$

(35)

In the enhanced UKF, we predict the linear state and available measurement based on the sigma points

$${\mathbf{x}}_{k|k-1}^{\mathrm{l},\left(i\right)}=\sum _{\chi =0}^{\vartheta +1}{\mathbf{v}}_{k-1}^{\left(\chi \right)}{\mathbf{\varsigma }}_k^{\left(\chi \right)}$$

(36)

$${{\mathbf{z}}_k}_{|k-1}=\sum _{\chi =0}^{\vartheta +1}{\mathbf{v}}_{k-1}^{\left(\chi \right)}{{\mathbf{z}}_k}_{|k-1}^{\left(\chi \right)}$$

(37)

Recalling the optimized ${\mathbf{K}}_k$ in equation (13), we can get ${\mathbf{P}}_{k|k-1}$ defined by equation (14).

Update stage

The cardinality distribution at time k is updated as

$${\rho }_k(n,{\ell }_{k|k-1}^{\left(i\right)})=\frac{{\Upsilon }_k^0\left({\mathbf{Z}}_k\right){\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i\right)})〉}$$

(38)

We draw $L_{k-1,\Gamma }$ independent and identical distribution (iid) particles ${\mathbf{x}}_{k|k-1}^{\mathrm{n},(i_{\Gamma })}$ from ${\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},(i)},{\mathbf{x}}_{k-1}^{\mathrm{n},(i)})$ to represent the undetected target and $L_{k-1,\mathrm{D}}$ iid particles ${\mathbf{x}}_{k|k-1}^{\mathrm{n},(i_{\mathrm{D}})}$ from ${\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},(i)},{\mathbf{x}}_{k-1}^{\mathrm{n},(i)},j)$ to represent the detected target. Then, the updated weights are

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{w}}_{k|k}^{\left(i_{\Gamma }\right)}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i_{\Gamma }\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\Gamma }\right)})}\\ \multicolumn{1}{c}{=\frac{1-p_D}{L_{k-1,\Gamma }}\frac{〈{\Upsilon }_k^1\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i_{\Gamma }\right)})〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k|k-1}^{\left(i_{\Gamma }\right)})〉}}\\ \multicolumn{1}{c}{\times \frac{p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\Gamma }\right)}\right)f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i_{\Gamma }\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\Gamma }\right)}){\mathbf{w}}_{k-1}^{\left(i_{\Gamma }\right)}}{{\mathrm{d}}_{\Gamma }({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i_{\Gamma }\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\Gamma }\right)})}}\end{array}$$

(39)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{w}}_{k|k}^{\left(i_{\mathrm{D}}\right)}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)},j^{\left(i_{\mathrm{D}}\right)})}\\ \multicolumn{1}{c}{=\frac{p_D}{L_{k-1,\mathrm{D}}}{g}_k({\mathbf{z}}_k^{\left(j\right)}|{\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)})\frac{〈1,c_k\left({\mathbf{z}}_k^{\left(j\right)}\right)〉}{c_k\left({\mathbf{z}}_k^{\left(j\right)}\right)}}\\ \multicolumn{1}{c}{\times \frac{〈{\Upsilon }_k^1({\mathbf{Z}}_k-\left\{{\mathbf{z}}_k^{\left(j\right)}\right\}),{\rho }_{k|k-1}(n,{\tau }_{k-1}^{\left(i_{\mathrm{D}}\right)})〉}{〈{\Upsilon }_k^0\left({\mathbf{Z}}_k\right),{\rho }_{k|k-1}(n,{\tau }_{k-1}^{\left(i_{\mathrm{D}}\right)})〉}}\\ \multicolumn{1}{c}{\times \frac{p_{s,k|k-1}\left({\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)}\right)f_{k|k-1}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)}|{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)}){\mathbf{w}}_{k-1}^{\left(i_{\mathrm{D}}\right)}}{{\mathrm{d}}_{\mathrm{D}}({\mathbf{x}}_{k|k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)},{\mathbf{x}}_{k-1}^{\mathrm{n},\left(i_{\mathrm{D}}\right)},j)}}\end{array}$$

(40)

For limiting the number of the required particles, it can be fixed based on the effective sample size

$$L_{\mathrm{eff},k|k-1,\Gamma }={\left(\sum _{i_{\Gamma }=1}^{L_{k|k-1,\Gamma }}{\left({\mathbf{w}}_{k|k,\Gamma }^{\left(i_{\Gamma }\right)}\right)}^2\right)}^{-1}$$

(41)

$$L_{\mathrm{eff},k|k-1,\mathrm{D}}={\left(\sum _{i_{\mathrm{D}}=1}^{L_{k|k-1,\mathrm{D}}}{\left({\mathbf{w}}_{k|k,\mathrm{D}}^{\left(i_{\mathrm{D}}\right)}\right)}^2\right)}^{-1}$$

(42)

Then, the updated number of particles is

$$L_{k|k-1,\Gamma }=min(L_{\mathrm{eff},k|k-1,\Gamma },L_{k-1,\Gamma })$$

(43)

$$L_{k|k-1,\mathrm{D}}=min(L_{\mathrm{eff},k|k-1,\mathrm{D}},L_{k-1,\mathrm{D}})$$

(44)

We identify the particle index and target track label. In view of detected target, the updated particle with the maximal weight ${\mathbf{w}}_{k|k,max}^{(i)}$ represents the associated track label. We evaluate ${\mathbf{w}}_{k|k,max,\mathrm{D}}^{(i_{\mathrm{D}})}$ during five continuous periods. If they are not less than 0.5, the target is detected, and the particle index $i_{\mathrm{D}}$ and track label ${\tau }_{k|k-1}^{(i_{\mathrm{D}})}$ are confirmed. Otherwise, the parameters represent the random clutter, which should be deleted. At this time, we evaluate ${\mathbf{w}}_{k|k,max,\Gamma }^{(i_{\Gamma })}$ during five continuous periods. If they are less than 0.5, the target is undetected. Both $i_{\Gamma }$ and ${\tau }_{k|k-1}^{(i_{\Gamma })}$ are deleted. We select the effective particles based on the weights

$$\begin{array}{c}\multicolumn{1}{c}{{\{{\mathbf{w}}_{k|k}^{\left(i\right)},{\mathbf{x}}_{k|k-1}^{\left(i\right)},{\tau }_{k|k-1}^{\left(i\right)}\}}_{i=1}^{L_{k|k-1}}}\\ \multicolumn{1}{c}{={\{{\mathbf{w}}_{k|k}^{\left(i_{\Gamma }\right)},{\mathbf{x}}_{k|k-1}^{\left(i_{\Gamma }\right)},{\tau }_{k|k-1}^{\left(i_{\Gamma }\right)}\}}_{i_{\Gamma }=1}^{L_{k|k-1,\Gamma }}}\\ \multicolumn{1}{c}{\cup {\{{\mathbf{w}}_{k|k}^{\left(i_{\mathrm{D}}\right)},{\mathbf{x}}_{k|k-1}^{\left(i_{\mathrm{D}}\right)},{\tau }_{k|k-1}^{\left(i_{\mathrm{D}}\right)}\}}_{i_{\mathrm{D}}=1}^{L_{k|k-1,\mathrm{D}}}}\end{array}$$

(45)

Simultaneously, the updated linear state and its covariance are

$$\begin{array}{c}{\mathbf{x}}_k^{\mathrm{l},\left(i\right)}={\mathbf{x}}_{k|k-1}^{\mathrm{l},\left(i\right)}+{\mathbf{K}}_k({\mathbf{z}}_k^{\left(\chi \right)}-{\mathbf{z}}_{k|k-1})\hfill \end{array}$$

(46)

$${\mathbf{P}}_k={\mathbf{P}}_{k|k-1}+{\mathbf{K}}_k{\mathbf{R}}_k{\mathbf{K}}_k^T$$

(47)

Considering the unitary matrix $\mathbf{U}$ and eigenvalue matrix ${\mathbf{\Sigma }}_{\chi }$ , we have the trace of ${\mathbf{P}}_k$

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{tr}({\mathbf{P}}_k)=\mathbf{U}{\mathbf{\Sigma }}_{\chi }{\mathbf{U}}^T}\\ \multicolumn{1}{c}{={\chi }^{-1}({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_k^{\mathrm{l},\left(i\right)})({\mathbf{\varsigma }}_k^{\left(\chi \right)}-{\mathbf{x}}_k^{\mathrm{l},\left(i\right)}{)}^T}\\ \multicolumn{1}{c}{=d_k^2}\end{array}$$

(48)

where $\mathrm{tr}(·)$ is the trace of matrix. Since that $\mathrm{tr}({\mathbf{P}}_k)$ equals to the square of distance $d_k^2$ between the sigma point and its linear state, the minimal ${\alpha }_k^2$ can be written as

$$\begin{array}{c}\multicolumn{1}{c}{min\left({\alpha }_k^2\right)=\frac{d_k^2}{max{\left(d_k^{\left(\chi \right)}\right)}^2}}\\ \multicolumn{1}{c}{=\frac{\mathrm{tr}\left({\mathbf{P}}_k\right)}{max\left\{{\left(\sqrt{{\mathbf{P}}_k}{\mathbf{\varsigma }}_k^{\left(\chi \right)}\right)}^T\left(\sqrt{{\mathbf{P}}_k}{\mathbf{\varsigma }}_k^{\left(\chi \right)}\right)\right\}}}\end{array}$$

(49)

Therefore, the optimized ${\alpha }_k^2$ is given by

$${\alpha }_k^2=max\{\frac{\mathrm{tr}\left({\mathbf{P}}_k\right)}{max\left\{{\left(\sqrt{{\mathbf{P}}_k}{\mathbf{\varsigma }}_k^{\left(\chi \right)}\right)}^T\left(\sqrt{{\mathbf{P}}_k}{\mathbf{\varsigma }}_k^{\left(\chi \right)}\right)\right\}},{10}^{-8}\}$$

(50)

where ${10}^{-8}$ is a referenced value.^9,36 It is important that the optimized ${\alpha }_k^2$ is used to compute the sigma points with the minimum quantity.

Remark 2

The number of the updated particles is reduced in equations (43) and (44), and $L_{k|k-1}$ effective particles are redrawn. We can only compute the limited sigma points corresponding to the certain particle i in order to reduce computational complexity. Furthermore, we achieve all sigma points in the update stage owing to the individual operation. In other words, in the enhanced UKF, the filtering process of $(L_{k-1}-L_{k|k-1})\chi$ sigma points is saved.

Estimation stage

In this stage, a new set $\{{\mathbf{w}}_k^{(i)},{\mathbf{x}}_k^{\mathrm{n},(i)},{\tau }_k^{(i)}{\}}_{i=1}^{L_k}$ is resampled from the updated particle set $\{{\mathbf{w}}_{k|k}^{(i)},{\mathbf{x}}_{k|k-1}^{\mathrm{n},(i)},{\tau }_{k|k-1}^{(i)}{\}}_{i=1}^{L_{k|k-1}}$ . At this time, we have the following PHD

$$D_k(\mathbf{x},\tau |{\mathbf{Z}}_{1:k})=\sum _{i=1}^{L_k}{\mathbf{w}}_k^{\left(i\right)}\delta (\mathbf{x}-{\mathbf{x}}_k^{\mathrm{n},\left(i\right)},\tau -{\tau }_k^{\left(i\right)})$$

(51)

In view of the different tracks at time k, the number of the target tracks $\left|{\tau }_k^{(j)}\right|$ is estimated using the following pseudo-code

for $j=1$ to i

for $i=L_k$ to i

if ${\tau }_k^{\left(j\right)}\ne {\tau }_k^{\left(i\right)}$

then record ${\tau }_k^{\left(j\right)}$

and compute the size $\left|{\tau }_k^{\left(j\right)}\right|$

end if

end for

end for

Therefore, the cardinality of targets is estimated as

$${\hat{N}}_k=\sum _{{\tau }_k^{\left(i\right)}=1}^{\left|{\tau }_k^{\left(j\right)}\right|}\arg \underset{n}{sup}{\rho }_k(n,{\tau }_k^{\left(i\right)})$$

(52)

Here, the clustering is required to partition the particle population into ${\hat{N}}_k$ clusters, and ${\hat{\mathbf{x}}}_k^{(1)},{\hat{\mathbf{x}}}_k^{(2)},\dots ,{\hat{\mathbf{x}}}_k^{({\hat{N}}_k)}$ are selected as the multi-target state estimates

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathbf{X}}}_k=\{{\hat{\mathbf{x}}}_k^{\left(1\right)},{\hat{\mathbf{x}}}_k^{\left(2\right)},\dots ,{\hat{\mathbf{x}}}_k^{\left({\hat{N}}_k\right)}\}}\\ \multicolumn{1}{c}{=\{{\mathbf{x}}_k^{\mathrm{l},\left(1\right)},{\mathbf{x}}_k^{\mathrm{l},\left(2\right)},\dots ,{\mathbf{x}}_k^{\mathrm{l},\left({\hat{N}}_k\right)}\}}\\ \multicolumn{1}{c}{\cup \{{\mathbf{x}}_k^{n,\left(1\right)},{\mathbf{x}}_k^{\mathrm{n},\left(2\right)},\dots ,{\mathbf{x}}_k^{n,\left({\hat{N}}_k\right)}\}}\end{array}$$

(53)

Figure 1 presents the filtering framework of the proposed filter. Note that two individual filters (enhanced UKF and enhanced PCPHD filter) complete the estimation of the linear and nonlinear components in one jump Markov model. Moreover, the optimal label of the target track is predicted, and the estimated cardinality of targets is computed in the enhanced PCPHD filter. These parameters are sent to the enhanced UKF in order to get the state estimates.

Figure 1.

Filtering framework.

Simulation results and discussions

Simulation definition

In this simulation, we used a scenario to measure the performance of the proposed filter. The simulation environment was Intel^® Core™ i5, 4 GB memory and MATLAB™ R2018a.

Three targets travel in a two-dimensional surveillance region with the scale of $[-2000,2000]\times [0,2000]{\mathrm{m}}^2$ .

1. During the former 30 s, the velocity of Target 1 is (–20, –20) m/s. It executes the CV motion from the position coordinate (1000, 1500) m. During the latter 30 s, its velocity is (25, –50) m/s. It keeps the clockwise CT motion, where the turning velocity is 6°/s.

2. During 11th to 50th second, the velocity of Target 2 is (12.5, 0) m/s. It remains the clockwise CT motion from the position coordinate (768.5, 1280) m, where the turning velocity is 1°/s.

3. During the 21st to 55th second, the velocity of Target 3 is (45, 45) m/s. It keeps the CV motion from the position coordinate (–1500, 250) m.

The sensor is located on the origin coordinate (0, 0) m, the detection probability is $P_D=98\%$ , the surveillance duration is 60 s, and the sampling time is $T=1\mathrm{s}$ . Besides, the average cardinality of random clutters is ${\lambda }_k=10$ . One hundred particles are used to approximate per object in this region.

In this scenario, the birth process can be written as follows: ${\pi }_{\Gamma ,k|k-1}^{(i)}={\{(r_{\Gamma ,k|k-1}^{(i)},p_{\Gamma ,k|k-1}^{(i)})\}}_{i=1}^3$ , ${\left\{r_{\Gamma ,k|k-1}^{(i)}\right\}}_{i=1}^3=0.1$ , $p_{\Gamma ,k|k-1}^{(i)}=N(x;m_{\Gamma }^{(i)},P_{\Gamma })$ , $m_{\Gamma }^{(1)}=[1000,0,1500,0,0{]}^{\mathrm{T}}$ , $m_{\Gamma }^{(2)}={[768.5,0,1280,0,0]}^{\mathrm{T}}$ , $m_{\Gamma }^{(3)}={[-1500,0,250,0,0]}^{\mathrm{T}}$ , and $P_{\Gamma }=\mathrm{diag}(10,10,10,10,1.2\pi 180)$ , where $\mathrm{diag}(·)$ denotes the diagonal matrix. At the initial time, the probabilities of two kinds of different motion models are equivalent. The measurement vector is formulated on the sensor information of both range and azimuth

$${\mathbf{z}}_k=\left[\begin{array}{c}\arctan \left(\frac{y_k}{x_k}\right)\hfill \\ \multicolumn{1}{c}{\sqrt{x_k^2+y_k^2}}\end{array}\right]+{\mathbf{e}}_k$$

(54)

The optimal sub-pattern assignment (OSPA) distance is compared the set of estimated position and true position. Let $\mathbf{X}=\{{\mathbf{x}}_i{\}}_{i=1}^g$ and $\hat{\mathbf{X}}=\{{\hat{\mathbf{x}}}_i{\}}_{i=1}^e$ be the true and estimated track sets, the OSPA distance ${\overline{d}}_p^{\left(c\right)}(\mathbf{X},\hat{\mathbf{X}})$ is defined by

$$\begin{array}{c}{\overline{d}}_p^{\left(c\right)}(\mathbf{X},\hat{\mathbf{X}})=\hfill \\ \{\begin{array}{cc}{\left(\frac{1}{\left|\hat{\mathbf{X}}\right|}(\underset{\pi \in {\Pi }_{\left|\hat{\mathbf{X}}\right|}}{min}\sum _{i=1}^{\left|\mathbf{X}\right|}{d}^{\left(c\right)}{({\mathbf{x}}_i,{\hat{\mathbf{x}}}_{\pi \left(i\right)})}^p+c^p(\left|\hat{\mathbf{X}}\right|-\left|\mathbf{X}\right|))\right)}^{1/p}\hfill & \left|\mathbf{X}\right|\le \left|\hat{\mathbf{X}}\right|\hfill \\ \multicolumn{1}{c}{min(c,\Vert \mathbf{x}-\hat{\mathbf{x}}\Vert )}& \left|\mathbf{X}\right|>\left|\hat{\mathbf{X}}\right|\hfill \end{array}\hfill \end{array}$$

(55)

where ${\Pi }_{\left|\hat{\mathbf{X}}\right|}$ is the set of permutations in $\left|\hat{\mathbf{X}}\right|$ , c is the cut-off parameter that determines the weight of penalties assigned to cardinality error as well as localization error, and p is the order parameter that determines sensitivity to outlier. Here, $p=2$ scales the position error, and c = 50 m defines the maximal distance.

Results and discussions

Figure 2 describes three-target tracks in a cluttered environment. As seen, the surveillance region is a top half-disk, where the radius is 2000 m. The measurements are concentrated especially around the sensor. The tracks of Targets 1 and 2 are nonlinear, and the track of Target 3 is a straight line. Besides, the circle and square denote the starting and ending points, respectively.

Figure 2.

Three-target tracks and measurements.

Figure 3 shows the true tracks, measurements, and estimated positions in both x and y coordinates. Note that the estimates of the proposed filter keep close to the actual target positions. The estimated tracks of the standard PCPHD filter tend to seriously drift off the true positions. It has some errors in the estimated positions. Due to the clutter-generated measurements, the standard filter cannot distinguish the random clutters when Target 2 is tracked. Therefore, it underestimates two objects when Target 1 is strong maneuvering. By comparison, the proposed filter identifies the estimated positions in two coordinates with the particle index and track label strategy. It tracks three targets under the condition of motion state transition and strong maneuvering.

Figure 3.

Target position estimates.

In Figure 4, the statistics of the estimated cardinality are demonstrated versus time. Note that three filters can basically keep converging to the cardinality during the whole surveillance time. Obviously, the standard PCPHD filter makes an unstable estimation. On the 18th second, it overestimates an object owing to clutter around Target 2. On the 31st second, Target 1 switches its dynamics from the constant velocity (CV) motion into the constant turn (CT) motion. However, the filter cannot reflect the transformation immediately so that it misses this target. Target 1 keeps strong maneuvering during the last half of the motion period. With respect to the rapid CT motion, the filter cannot detect it on the 55th second. Although the referenced filter ²⁵ adopts the jump Markov model, it cannot well distinguish the actual target from the random clutters based on the embedded UKF-PF algorithm. Then, the over-estimated cardinality occurs. By comparison, the estimated cardinality of the proposed filter is not affected by uncertain situations because of the optimized filtering structure. Especially, the track label plays a major role during the tracking process. Since it effectively tracks all targets at each time, the inherent defects in the standard PCPHD are overcome. The enhanced UT also describes the dynamics of maneuvering targets. Obviously, the proposed filter has its confidence in the estimated cardinality with the smaller standard deviation (std).

Figure 4.

Target cardinality estimates.

Figure 5 shows the OSPA distance of three kinds of filters. Note that the OSPA distance of the proposed filter is lower, which indicates the advances in the tracking accuracy. The reason can be explained that many robust particles are drawn from the modified proposal distributions. Combined with the enhanced UKF, the estimated positions are concentrated around the true tracks. For comparison, the standard filter has estimation bias due to its inherent filtering defects. The referenced filter has compromising performance.

Figure 5.

OSPA distance.

With respect to the tracking reliability of the proposed filter in a different cluttered environment, the similar scenarios under various clutter rates ${\lambda }_k=5,10,20$ are considered. In Figure 6, the statistic cardinality information is shown. As expected, with decreasing clutter rate, the proposed filter can produce the results essentially in agreement with the true cardinality. Its estimates are uneasily influenced by the clutter and motion state.

Figure 6.

Target number estimates under different clutter rates.

Furthermore, Figure 7 gives the OSPA distance under three clutter rates. Obviously, the proposed filter has satisfactory tracking results. Recalling Figure 6, the OSPA distance is the smallest in terms of ${\lambda }_k=5$ .

Figure 7.

OSPA distance under different clutter rates.

Table 1 cites the average tracking performance of the standard PCPHD filter, referenced filter, and proposed filter, which mainly contains the average cardinality, average OSPA distance and average running time. With respect to the tracking reliability, the proposed filter keeps perfect cardinality estimates (the true cardinality is 2.250) and saves OSPA distance. Based on the two-layer model and variable structure, the referenced filter has also reduced OSPA distance. However, clutter interference is inevitable. As a result, the cardinality estimates are larger. As for the tracking efficiency, the standard PCPHD filter is of the shortest running time. Since the proposed filter has reduced the required number of particles, the running time is roughly similar to the referenced filter.

Table 1.

Comparison of average tracking performance.

	Cardinality	OSPAdistance (m)	Runningtime (s)
Proposed filter	2.2376	28.5667	2.3010
Standard PCPHDfilter	2.1199	34.7438	1.9094
Referencedfilter	2.3083	31.9356	2.2795

OSPA: optimal sub-pattern assignment; PCPHD: particle cardinalized probability hypothesis density.

Finally, we utilize the proposed filter to track some mobile targets during 60 s in the surveillance region to achieve the actual tracking results. We randomly extract five frames in Figure 8. As seen, the proposed filter applies the innovative method in the framework of jump Markov mode in order to correctly track all actual targets with various dynamics. Each target is marked by a circle and its name. Moreover, the filter can effectively identify the new target (T36) in Frame 3 because the target T15 leaves the region after Frame 2. In general, the proposed filter has satisfactory improvements due to its optimal mechanism.

Figure 8.

Experimental results of actual multi-target tracking: (a) mobile targets tracking in frame 1, (b) mobile targets tracking in frame 2, (c) mobile targets s tracking in frame 3, (d) mobile targets tracking in frame 4, and (e) mobile targets tracking in frame 5.

Conclusion

This paper has developed an innovative UT-based PCPHD filter. The primary challenges are to deal with the imprecise estimations of the standard PCPHD filter with clutter-generated measurements. We use the AP method to get certain distributions for drawing robust particles. Moreover, the UT is improved using the optimized Kalman gain. In the jump Markov model, we integrate them based on the track label strategy. From theoretical analysis and numerical study, it is revealed that the required particles can be effectively sampled from the modified proposal distributions to estimate the target position. On the other hand, the target velocity is estimated with the enhanced UT.

In addition, it is important to enhance the tracking performance of the standard PCPHD filter because of the novel filtering mechanism. The proposed filter has advanced in tracking mobile targets. However, the modified proposal distributions in extended target space bring about the additional computational load. In order to track the strong maneuvering targets, in the future research of this work, we will consider reducing computational complexity by optimizing the filter structure.