Adaptive probability hypothesis density filter for multi-target tracking with unknown measurement noise statistics

Abstract

Under the Gaussian noise assumption, the probability hypothesis density (PHD) filter represents a promising tool for tracking a group of moving targets with a time-varying number. However, inaccurate prior statistics of the random noise will degrade the performance of the PHD filter in many practical applications. This paper presents an adaptive Gaussian mixture PHD (AGM-PHD) filter for the multi-target tracking (MTT) problem in the scenario where both the mean and covariance of measurement noise sequences are unknown. The conventional PHD filters are extended to jointly estimate both the multi-target state and the aforementioned measurement noise statistics. In particular, the Normal-inverse-Wishart and Gaussian distributions are first integrated to represent the joint posterior intensity by transforming the measurement model into a new formulation. Then, the updating rule for the hyperparameters of the model is derived in closed form based on variational Bayesian (VB) approximation and Bayesian conjugate prior heuristics. Finally, the dynamic system state and the noise statistics are updated sequentially in an iterative manner. Simulations results with both constant velocity and constant turn model demonstrate that the AGM-PHD filter achieves comparable performance as the ideal PHD filter with true measurement noise statistics.

Introduction

Benefit from the increasing advances in sensing techniques and computer science, multi-target tracking has found its potential applications in various disciplines, such as autonomous robotics, intelligence surveillance, remote sensing and even biomedical research.^1–3 The goal of MTT is to estimate the number of targets together with their individual states from a given sequence of measurements, which is contaminated by a number of uncertain sources. The problem of tracking a group of targets was traditionally resolved by many data association-based methods, mainly including nearest neighbor (NN), ⁴ the joint integrated probabilistic data association (JIPDA), ⁵ multiple hypothesis tracking (MHT), ⁶ and multiple-target particle filter. ⁷ However, due to its combinatorial nature, theses tracking algorithms could be computationally intensive when the target number grows.

The random finite sets (RFS) theory was first introduced into the MTT by Mahler, ⁸ where the explicit data association was avoided by representing the target states and noisy measurements as two individual RFS. Consequently, the MTT problem was solved within the Bayesian filtering framework utilizing the RFS formulation. In addition, the classical PHD filter ⁸ and its several typical extensions, such as multi-Bernoulli filters^9,10 and cardinalized PHD,^11,12 were proposed to approximate the state posterior by propagating the first-order statistical moment or the intensity of the RFS. However, the recursion of various PHD filters is still intractable in general as it involves the calculation of multiple set integrals. Two closed-form implementations of the PHD filter were proposed, including Gaussian mixture PHD (GM-PHD) filter ¹³ and particle PHD filter ¹⁴ which are based on Gaussian mixture (GM) strategy and sequential Monte Carlo (SMC), respectively. Generally, the GM-PHD filter is more efficient and reliable than particle PHD filter as it obviates the expensive clustering process in the particle-based approach. Besides, the standard Kalman filter (KF) and its nonlinear variations could be easily adopted to estimate the individual mixture components under the Gaussian model assumption.

The statistics of the measurement noise statistics, especially the noise deviation and covariance, may not always be a given premise in some practical applications. In the cases that the actual measurement noise statistics deviate from the assumed ones, those standard Gaussian filters will suffer degraded performance or even divergence. In order to realize the state estimation problem involving unknown noise statistics, several adaptive filters were introduced within the framework of Bayesian estimation. For instance, to improve the estimation accuracy for vehicle navigation with noise uncertainty, the fuzzy interacting multiple model (IMM) is employed in sensor fusion to determine the bounds of the process noise covariance. ¹⁵ Zhao et al. ¹⁶ presented a general filtering scheme based on the IMM algorithm for nonlinear systems with unknown and random nuisance parameters. Furthermore, the linear minimum mean square error (LMMSE) filter ¹⁷ was embedded into the above IMM based filtering scheme to handle ballistic target tracking with unknown ballistic coefficient efficiently. By introducing an equivalent covariance matrix into the estimator, the robustness of the standard KF is enhanced against the uncertainty in the process noise. ¹⁸ As an example for MTT, the H infinite optimality criterion is integrated with the GM-PHD filter to handle the noise parameter uncertainties. ¹⁹ However, the filter could yield a conservative result since it minimizes the possible worst-case estimation errors by treating the noise statistics as bounded quantities. In addition to the above methods, the variational Bayesian (VB) approximation is introduced into Kalman filtering to reclusively estimate the dynamic state and time-varying variance of the measurement noise for linear Gaussian models. ¹⁹ For the nonlinear system model along with time-varying measurement noise covariance, a robust unscented KF using VB approximation is proposed to improve its adaptivity. ²⁰ The VB based method has advantages over other aforementioned methods that it deals with posterior inference analytically at a low computational cost, since the posterior is assumed in the same functional form as the prior based on Bayesian conjugacy. ²¹ As a result, some adaptive multi-target tracking filters based on the VB approximation has been further addressed, and an inverse-Wishart or a product of inverse-Gamma distribution was defined as a conjugate prior for the measurement noise covariance.^22–25 However, the mean of the Gaussian measurement noise was assumed to be known as a zero vector. To deal with the cases that the mean of the noises are also unavailable, Özkan et al. ²⁶ proposed a marginalized particle filter (PF) to update the noise statistics by employing finite dimensional sufficient statistics of Normal-inverse-Wishart distributions for each particle.

When referring to the MTT scenario with unknown mean and covariance of the noise, Li et al. ²⁷ introduced a measurement preprocessing procedure prior to the multi-sensor PHD updating step to cope with the unknown statistical properties of the sensors. The measurement noise bias was handled by a Monte Carlo sampling-based debiasing approach with the multi-sensor network. ²⁸ Here, a new adaptive PHD filter is proposed to recursively estimate both the unknown mean and covariance of the measurement noise in the process of target tracking. In the proposed filter, the multi-target state is represented by a Gaussian distribution, while the statistics of the measurement noise is modeled with a Normal-inverse-Wishart distribution. The measurement equation is reformed with the residual vector, and the rule for updating the hyperparameters of the measurement noise statistics is then derived using two different methods: the VB approximation and the Bayesian conjugate prior heuristics. On this basis, the hyperparameters of the model are determined sequentially and iteratively at each measurement updating step. The effectiveness of the proposed filter is evaluated by two multi-target tracking simulations where both the mean and covariance of the measurement noise are unknown.

The organization of the paper is assigned as follows. In section 2, the model formulation of target tracking with unknown measurement noise statistics is introduced. In section 3, the updating rules for the hyperparameters of the measurement noise statistics are derived, then the details of the proposed filter are presented. Section 4 illustrates the effectiveness of the proposed adaptive PHD filter with simulations. The main conclusions are given in the last section.

Problem formulation

Consider the following discrete-time state-space model for each target

$${\mathbf{x}}_k=f({\mathbf{x}}_{k-1})+{\mathbf{v}}_k$$ (1)

$${\mathbf{z}}_k=h({\mathbf{x}}_k)+{\mathbf{w}}_k$$ (2)

where $k$ denotes the time index. ${\mathbf{x}}_k\in {\mathcal{R}}^{n_x}$ and ${\mathbf{z}}_k\in {\mathcal{R}}^{n_z}$ are the target state and the measurement, respectively. $f(•)$ is the state transition function and $h(•)$ is the measurement function. The process noise ${\mathbf{v}}_k$ is additive with the Gaussian distribution $\mathcal{N}(0,{\mathbf{Q}}_k)$ . The measurement noise ${\mathbf{w}}_k$ is also additive with the Gaussian distribution $\mathcal{N}({\mathbf{r}}_k,{\mathbf{R}}_k)$ , but with non-zero mean ${\mathbf{r}}_k$ and covariance ${\mathbf{R}}_k$ . Besides, the noises are assumed mutually independently and identically distributed.

Since the exact information about the measurement noise statistics is unavailable, not only the multi-target state but also the measurement noise statistics have to be estimated simultaneously. For this purpose, the measurement residual ${\mathbf{e}}_k$ is used to represent the non-zero mean Gaussian measurement noise as

$${\mathbf{e}}_k({\mathbf{x}}_k,{\mathbf{z}}_k)={\mathbf{z}}_k-h({\mathbf{x}}_k)$$ (3)

For multivariate Gaussian distribution of the random variable ${\mathbf{e}}_k$ with unknown mean and covariance, the Normal-inverse-Wishart distribution defines as its conjugate prior. The measurement model can thus be transformed into a hierarchical Bayesian model as (4), where the density of the mean ${\mathbf{r}}_k$ conditioned on the covariance ${\mathbf{R}}_k$ is a Gaussian distribution, while the covariance ${\mathbf{R}}_k$ is inverse-Wishart distributed, that is,

$$\begin{array}{c}p({\mathbf{e}}_k|{\mathbf{x}}_k,{\mathbf{z}}_k)=\mathcal{N}({\mathbf{e}}_k;{\mathbf{r}}_k,{\mathbf{R}}_k)\hfill \\ p({\mathbf{r}}_k|{\mathbf{R}}_k)=\mathcal{N}({\mathbf{r}}_k;{\mathbf{\mu }}_k,\frac{{\mathbf{R}}_k}{{\kappa }_k})\hfill \\ p({\mathbf{R}}_k)=\mathcal{I}{\mathcal{W}}_{{\nu }_k}({\mathbf{R}}_k;{\mathbf{\Lambda }}_k)\hfill \end{array}$$ (4)

where $\mathcal{N}({\mathbf{r}}_k;{\mathbf{\mu }}_k,\frac{{\mathbf{R}}_k}{{\kappa }_k})$ denotes a Gaussian distribution of ${\mathbf{r}}_k$ , and the hyperparameters ${\mathbf{\mu }}_k$ and ${\kappa }_k$ are the mean and scale factor, respectively. $\mathcal{I}{\mathcal{W}}_{{\nu }_k}({\mathbf{R}}_k;{\mathbf{\Lambda }}_k)$ denotes an inverse-Wishart distribution of ${\mathbf{R}}_k$ , with the scale matrix ${\Lambda }_k$ and degree of freedom ${\nu }_k$ . The probability densities of the above Gaussian distribution and inverse-Wishart distribution are as follows:

$$\begin{array}{c}p({\mathbf{r}}_k|{\mathbf{R}}_k)=\frac{{\kappa }_k}{{\left(2\pi \right)}^{n_z/2}}{\left|{\mathbf{R}}_k\right|}^{-1/2}\hfill \\ \exp \{-\frac{{\kappa }_k}{2}({\mathbf{r}}_k-{\mathbf{\mu }}_k{)}^T{\mathbf{R}}_k^{-1}({\mathbf{r}}_k-{\mathbf{\mu }}_k)\}\hfill \\ p({\mathbf{R}}_k)=\mathcal{I}{\mathcal{W}}_{{\nu }_k}({\mathbf{R}}_k;{\mathbf{\Lambda }}_k)={\left|{\mathbf{R}}_k\right|}^{-({\nu }_k+n_z+1)/2}\hfill \\ \exp \{-\frac{1}{2}\mathrm{tr}({\mathbf{\Lambda }}_k{\mathbf{R}}_k^{-1})\}\hfill \end{array}$$ (5)

To track multiple targets with a time-varying number, RFSs are used to represent both the set of the target states and the measurements. The multi-target state RFS in the state space $\mathcal{X}$ and the multi-target measurement RFS in the observation space $\mathcal{Z}$ are defined as

$$\begin{array}{c}{\mathbf{X}}_k\stackrel{·}{=}\{{\mathbf{x}}_{k,1},\cdots ,{\mathbf{x}}_{k,n_k}\}\subset \mathcal{X}\hfill \\ {\mathbf{Z}}_k\stackrel{·}{=}\{{\mathbf{z}}_{k,1},\cdots ,{\mathbf{z}}_{k,m_k}\}\subset \mathcal{Z}\hfill \end{array}$$ (6)

where $n_k$ and $m_k$ denote the number of targets and measurements, respectively.

According to the finite set statistics (FISST) theory, the estimation of multi-target state can be approximated by the PHD filter which is sub-optimal but computation tractable. Since the knowledge of the measurement noise statistics is unknown, the standard PHD filter has to be extended by augmenting the target state with the measurement noise statistics. Assume that at time $k-1$ , the joint posterior intensity of the multi-target state and the measurement noise statistics is denoted by ${\upsilon }_{k-1}({\mathbf{x}}_{k-1},{\mathbf{\theta }}_{k-1}|{\mathbf{Z}}_{1:k-1})$ , where ${\mathbf{\theta }}_{k-1}=[{\mathbf{r}}_{k-1},{\mathbf{R}}_{k-1}]$ denotes the measurement noise statistics. The predicted and updated intensities at time $k$ can be derived by the two-stage Bayesian recursion:

$$\begin{array}{c}{\upsilon }_{k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k|{\mathbf{Z}}_{1:k-1})=\int p_{s,k}p({\mathbf{x}}_k|{\mathbf{x}}_{k-1})p({\mathbf{\theta }}_k|{\mathbf{\theta }}_{k-1})\hfill \\ \times {\upsilon }_{k-1}({\mathbf{x}}_{k-1},{\mathbf{\theta }}_{k-1}|{\mathbf{Z}}_{1:k-1})d{\mathbf{x}}_{k-1}d{\mathbf{\theta }}_{k-1}+{\gamma }_k({\mathbf{x}}_k,{\mathbf{\theta }}_k)\hfill \end{array}$$ (7)

$$\begin{array}{c}{\upsilon }_k({\mathbf{x}}_k,{\mathbf{\theta }}_k|{\mathbf{Z}}_{1:k})=(1-p_{d,k}){\upsilon }_{k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k|{\mathbf{Z}}_{1:k-1})\hfill \\ +\sum _{{\mathbf{z}}_k\in {\mathbf{Z}}_k}\frac{p_{d,k}p({\mathbf{z}}_k|{\mathbf{x}}_k,{\mathbf{\theta }}_k){\upsilon }_{k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k|{\mathbf{Z}}_{1:k-1})}{{\eta }_k({\mathbf{z}}_k)+\int p_{d,k}p({\mathbf{z}}_k|{\xi }_{\mathbf{x}},{\xi }_{\mathbf{\theta }}){\upsilon }_{k|k-1}({\xi }_{\mathbf{x}},{\xi }_{\mathbf{\theta }}|{\mathbf{Z}}_{1:k-1})d{\xi }_{\mathbf{x}}d{\xi }_{\mathbf{\theta }}}\hfill \end{array}$$ (8)

where ${\gamma }_k(·)$ is the intensity of the target birth RFS. ${\eta }_k(·)$ denotes the intensity of Poisson clutter RFS. $p_{s,k}$ and $p_{d,k}$ are the target survival probability and detection probability, respectively. ${\mathbf{Z}}_{1:k}=\{{\mathbf{Z}}_1,\cdots ,{\mathbf{Z}}_k\}$ is the sequence of given measurements up to time $k$ . Besides, suppose that the evolution model of the measurement noise statistics $p({\mathbf{\theta }}_k|{\mathbf{\theta }}_{k-1})$ is independent of the given state transition model $p({\mathbf{x}}_k|{\mathbf{x}}_{k-1})$ , that is,

$$p({\mathbf{x}}_k,{\mathbf{\theta }}_k|{\mathbf{x}}_{k-1},{\mathbf{\theta }}_{k-1})=p({\mathbf{x}}_k|{\mathbf{x}}_{k-1})p({\mathbf{\theta }}_k|{\mathbf{\theta }}_{k-1})$$ (9)

It is noted that the multi-target state is coupled with the unknown noise statistics in the measurement likelihood $p({\mathbf{z}}_k|{\mathbf{x}}_k,{\mathbf{\theta }}_k)$ , which results in the problem that the integrals involved in the standard PHD recursions cannot be analytically handled.

Proposed adaptive PHD filter

Joint estimation of state and measurement noise statistics

According to the definition of the transformed measurement model, the posterior of the measurement mean and covariance is reformed as

$$\begin{array}{c}p({\mathbf{r}}_{k-1},{\mathbf{R}}_{k-1}|{\mathbf{e}}_{1:k-1})\hfill \\ =\mathcal{N}({\mathbf{r}}_{k-1};{\mathbf{\mu }}_{k-1},\frac{{\mathbf{R}}_{k-1}}{{\kappa }_{k-1}})\mathcal{I}{\mathcal{W}}_{{\nu }_{k-1}}({\mathbf{R}}_{k-1};{\mathbf{\Lambda }}_{k-1})\hfill \end{array}$$ (10)

where the quantities ${\mathbf{\mu }}_{k-1}$ , ${\kappa }_{k-1}$ , ${\nu }_{k-1}$ and ${\mathbf{\Lambda }}_{k-1}$ are the corresponding hyperparameters, and ${\mathbf{e}}_{1:k-1}=\{{\mathbf{e}}_1,{\mathbf{e}}_2,\cdots ,{\mathbf{e}}_{k-1}\}$ denotes the sequence of the measurement residuals up to the current time.

Assume that the predicted density of the measurement noise statistic is in the same form as the posterior, i.e.,

$$\begin{array}{c}p({\mathbf{r}}_{k|k-1},{\mathbf{R}}_{k|k-1}|{\mathbf{e}}_{1:k-1})=\mathcal{N}({\mathbf{r}}_{k|k-1};{\mathbf{\mu }}_{k|k-1},\hfill \\ \frac{{\mathbf{R}}_{k|k-1}}{{\kappa }_{k|k-1}})\mathcal{I}{\mathcal{W}}_{{\nu }_{k|k-1}}({\mathbf{R}}_{k|k-1};{\mathbf{\Lambda }}_{k|k-1})\hfill \end{array}$$ (11)

Here, the VB approximation ²⁹ is employed to compute the estimate of the posterior density:

$$p({\mathbf{r}}_k,{\mathbf{R}}_k|{\mathbf{e}}_{1:k})\approx Q_r({\mathbf{r}}_k)Q_R({\mathbf{R}}_k)$$ (12)

where $Q_r({\mathbf{r}}_k)$ and $Q_R({\mathbf{R}}_k)$ are the approximated densities corresponding to the unknown measurement noise mean ${\mathbf{r}}_k$ and covariance ${\mathbf{R}}_k$ , respectively. The key of the VB approximation is to determine the posterior density by calculating the minimum of Kullback-Leibler divergence (KLD), which is defined as

$$\begin{array}{c}\mathrm{KLD}(Q_r({\mathbf{r}}_k)Q_R({\mathbf{R}}_k)||p({\mathbf{r}}_k,{\mathbf{R}}_k|{\mathbf{e}}_{1:k}))\hfill \\ =\int Q_r({\mathbf{r}}_k)Q_R({\mathbf{R}}_k)\log \frac{Q_r({\mathbf{r}}_k)Q_R({\mathbf{R}}_k)}{p({\mathbf{r}}_k,{\mathbf{R}}_k|{\mathbf{e}}_{1:k})}d{\mathbf{r}}_{k}d{\mathbf{R}}_k\hfill \end{array}$$ (13)

The solution to this optimization problem is obtained by minimizing KLD with respect to the $Q_r({\mathbf{r}}_k)$ and $Q_R({\mathbf{R}}_k)$ sequentially, while retaining the other part unchanged. Finally, the hyperparameters at time step $k$ can be approximated from their corresponding values at time step $k-1$ given by

$$\begin{array}{c}{\mathbf{\mu }}_k=\frac{{\kappa }_{k|k-1}}{{\kappa }_{k|k-1}+1}{\mathbf{\mu }}_{k|k-1}+\frac{1}{{\kappa }_{k|k-1}+1}{\mathbf{e}}_k\hfill \\ {\kappa }_k={\kappa }_{k|k-1}+1\hfill \\ {\nu }_k={\nu }_{k|k-1}+1\hfill \\ {\mathbf{\Lambda }}_k={\mathbf{\Lambda }}_{k|k-1}+\frac{{\kappa }_{k|k-1}}{{\kappa }_{k|k-1}+1}({\mathbf{e}}_k-{\mathbf{\mu }}_{k|k-1})({\mathbf{e}}_k-{\mathbf{\mu }}_{k|k-1}{)}^T\hfill \end{array}$$ (14)

Proof: For the proof see appendix A.

In what follows, the measurement residuals are treated as samples to derive the updating rule for the hyperparameters of the measurement noise statistics. The Normal-inverse-Wishart posterior can be calculated from its Bayesian conjugate prior with samples of measurement data as follows ³⁰ :

$$\begin{array}{c}{\mathbf{\mu }}_k=\frac{{\kappa }_{k|k-1}}{{\kappa }_{k|k-1}+n}{\mathbf{\mu }}_{k|k-1}+\frac{n}{{\kappa }_{k|k-1}+n}{\overline{\mathbf{z}}}_k\hfill \\ {\kappa }_k={\kappa }_{k|k-1}+n\hfill \\ {\nu }_k={\nu }_{k|k-1}+n\hfill \\ {\mathbf{\Lambda }}_k={\mathbf{\Lambda }}_{k|k-1}+{\mathbf{S}}_k\hfill \\ +\frac{{\kappa }_{k|k-1}n}{{\kappa }_{k|k-1}+n}({\overline{\mathbf{z}}}_k-{\mathbf{\mu }}_{k|k-1})({\overline{\mathbf{z}}}_k-{\mathbf{\mu }}_{k|k-1}{)}^T\hfill \end{array}$$ (15)

where $n$ is the number of the data set $\{{\mathbf{z}}_i{\}}_{i=1}^n$ , ${\overline{\mathbf{z}}}_k$ is the mean of the samples, and ${\mathbf{S}}_k$ is the covariance matrix of the samples. The sample mean and covariance are defined as

$$\begin{array}{c}{\overline{\mathbf{z}}}_k=\frac{1}{n}\sum _{i=1}^n{\mathbf{z}}_i\hfill \\ {\mathbf{S}}_k=\sum _{i=1}^n({\mathbf{z}}_i-{\overline{\mathbf{z}}}_k)({\mathbf{z}}_i-{\overline{\mathbf{z}}}_k{)}^T\hfill \end{array}$$ (16)

The above general recursive equations are specialized for the update rule of the hyperparameters of the measurement noise heuristically. Considering that the measurements are processed sequentially in the update step of the filter, the sample mean is equal to the current measurement ${\mathbf{z}}_k$ and the covariance is a zero matrix. The measurement ${\mathbf{z}}_k$ is then replaced by the measurement residual ${\mathbf{e}}_k$ with a simple transformation given in equation (3). Consequently, the updated values of the hyperparameters are computed as

$$\begin{array}{c}{\mathbf{\mu }}_k={\mathbf{\mu }}_{k|k-1}+\frac{1}{{\kappa }_{k|k-1}+1}({\mathbf{e}}_k-{\mathbf{\mu }}_{k|k-1})\hfill \\ {\kappa }_k={\kappa }_{k|k-1}+1\hfill \\ {\nu }_k={\nu }_{k|k-1}+1\hfill \\ {\mathbf{\Lambda }}_k={\mathbf{\Lambda }}_{k|k-1}+\frac{{\kappa }_{k|k-1}}{{\kappa }_{k|k-1}+1}({\mathbf{e}}_k-{\mathbf{\mu }}_{k|k-1})({\mathbf{e}}_k-{\mathbf{\mu }}_{k|k-1}{)}^T\hfill \end{array}$$ (17)

Note that the above heuristic updating rules are essentially the same as the equations (17) and (19). When the hyperparameters at time step $k$ is obtained, the predicted density of the measurement residual ${\mathbf{e}}_k$ becomes a multivariate t-distribution as

$$p({\mathbf{e}}_k|{\mathbf{e}}_{1:k-1})={\mathcal{T}}_{{\nu }_k-n_z+1}({\mathbf{e}}_k;{\mathbf{\mu }}_k,\frac{{\mathbf{\Lambda }}_k({\kappa }_k+1)}{{\kappa }_k({\nu }_k-n_z+1)})$$ (18)

with center at ${\mathbf{\mu }}_k$ , precision $\frac{{\mathbf{\Lambda }}_k({\kappa }_k+1)}{{\kappa }_k({\nu }_k-n_z+1)}$ , and degree of freedom ${\nu }_k-n_z+1$ . The posterior density of the measurement noise statistics can be denoted by a single Normal-inverse-Wishart distribution as follows ³⁰ :

$$\begin{array}{c}p({\mathbf{r}}_k,{\mathbf{R}}_k|{\mathbf{e}}_{1:k})\stackrel{\mathrm{\Delta }}{=}\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_k;{\mathbf{\mu }}_k,{\kappa }_k,{\nu }_k,{\mathbf{\Lambda }}_k)\hfill \\ =\frac{1}{Z}{\left|{\mathbf{R}}_k\right|}^{-(1+\frac{{\nu }_k+n_z}{2})}\exp (-\frac{1}{2}\mathrm{tr}({\mathbf{\Lambda }}_k{\mathbf{R}}_k^{-1})\hfill \\ -\frac{{\kappa }_k}{2}({\mathbf{r}}_k-{\mathbf{\mu }}_k{)}^T{\mathbf{R}}_k^{-1}({\mathbf{r}}_k-{\mathbf{\mu }}_k))\hfill \end{array}$$ (19)

where $Z$ is a constant defined as

$$Z=\frac{2^{{\nu }_k{n}_z/2}{\Gamma }_{n_z}({\nu }_k/2){(2\pi /{\nu }_k)}^{n_z/2}}{{\left|{\mathbf{\Lambda }}_k\right|}^{{\nu }_k/2}}$$ (20)

To make full use of the Bayesian conjugate prior, both the joint prior and predicted densities of ${\mathbf{x}}_k$ and ${\mathbf{\theta }}_k$ are assumed to be a product of Normal-inverse-Wishart and Gaussian distribution. For the given transmit and measurement models, the Bayesian recursions are then employed to obtain the joint posterior density. In the prediction step, to ensure that the predicted density and the prior density are in the same formation, the predicted values of the hyperparameters are computed with a simple evolution model:

$$\begin{array}{c}{\mathbf{\mu }}_{k|k-1}={\mathbf{\mu }}_{k-1}\hfill \\ {\kappa }_{k|k-1}=\frac{1}{\rho }{\kappa }_{k-1}\hfill \\ {\nu }_{k|k-1}=\rho {\nu }_{k-1}\hfill \\ {\mathbf{\Lambda }}_{k|k-1}=\rho {\mathbf{\Lambda }}_{k-1}\hfill \end{array}$$ (21)

where $\rho \in (0,1]$ is a given forgetting factor. Meanwhile, the mean and covariance of the target state can be predicted by applying the prediction equations of Gaussian filters.

In the measurement update step, it is intractable to directly derive an exact value of the joint posterior density because the target state is coupled with the measurement noise statistics. Therefore, an iterative method is used to estimate the target state and measurement noise statistics sequentially. Specifically, the measurement noise statistics are first estimated with the updated hyperparameters, which are computed according to the updating rule presented above:

$$\begin{array}{c}{\hat{\mathbf{r}}}_k={\mathbf{\mu }}_k\hfill \\ {\hat{\mathbf{R}}}_k=\frac{{\mathbf{\Lambda }}_k}{{\nu }_k-n_z-1}\hfill \end{array}$$ (22)

The updated mean and covariance of the target state are estimated by standard update equations of Gaussian filters with the known measurement noise statistics. By substituting the updated mean of the target state into the measurement model, one can obtain the new measurement residual, which is used to update the hyperparameters of the measurement noise statistics again. These iterative steps are performed multiple times until the result is converged or the preset maximum iteration limit is reached.

Adaptive GM-PHD recursion

In this section, the standard GM-PHD recursion is extended in a sense that not only the target state is estimated, but also the measurement noise statistics are adapted at the same time. Assume that all targets follow the same state-space model defined in section 2, thus the state transition and likelihood functions are both Gaussian distributed, i.e.,

$$p({\mathbf{x}}_k^j|{\mathbf{x}}_{k-1}^j)=\mathcal{N}({\mathbf{x}}_k^j;f({\mathbf{x}}_{k-1}^j),{\mathbf{Q}}_k)$$ (23)

$$p({\mathbf{z}}_k^j|{\mathbf{x}}_k^j,{\mathbf{\theta }}_k^j)=\mathcal{N}({\mathbf{z}}_k^j;h({\mathbf{x}}_{k|k-1}^j)+{\mathbf{r}}_k,{\mathbf{R}}_k)$$ (24)

where $j$ is the target index. To ensure that the adaptive filter have closed-form solutions, the spontaneous birth are characterized by the intensities comprised of a product of Normal-inverse-Wishart and Gaussian distribution mixtures:

$$\begin{array}{c}{\gamma }_k({\mathbf{x}}_k,{\mathbf{\theta }}_k)\hfill \\ =\sum _{j=1}^{J_{\gamma ,k}}[w_{\gamma ,k}^j\mathcal{N}({\mathbf{x}}_k;{\mathbf{m}}_{\gamma ,k}^j,{\mathbf{P}}_{\gamma ,k}^j)\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_k;{\mathbf{\mu }}_{\gamma ,k}^j,{\kappa }_{\gamma ,k}^j,{\nu }_{\gamma ,k}^j,{\mathbf{\Lambda }}_{\gamma ,k}^j)]\hfill \end{array}$$ (25)

where $J_{\gamma ,k}$ is the total number of mixtures with weight $w_{\gamma ,k}^j$ , mean ${\mathbf{m}}_{\gamma ,k}^j$ and covariance ${\mathbf{P}}_{\gamma ,k}^j$ . The hyperparameters for each mixture of the birthed target are denoted by ${\mathbf{\mu }}_{\gamma ,k}^j$ , ${\kappa }_{\gamma ,k}^j$ , ${\nu }_{\gamma ,k}^j$ , and ${\mathbf{\Lambda }}_{\gamma ,k}^j$ .

By combining (7), (8), and (25) together, the proposed adaptive PHD recursions are detailed as follows.

Time Prediction: Assume that the joint prior intensity of the multi-target state and measurement noise statistics at the time $k-1$ is also represented by a product of Normal-inverse-Wishart and Gaussian mixtures:

$${\upsilon }_{k-1}({\mathbf{x}}_{k-1},{\mathbf{\theta }}_{k-1})=\sum _{j=1}^{J_{k-1}}[w_{k-1}^j\mathcal{N}({\mathbf{x}}_{k-1};{\mathbf{m}}_{k-1}^j,{\mathbf{P}}_{k-1}^j)\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_{k-1};{\mathbf{\mu }}_{\gamma ,k-1}^j,{\kappa }_{\gamma ,k-1}^j,{\nu }_{\gamma ,k-1}^j,{\mathbf{\Lambda }}_{\gamma ,k-1}^j)]$$ (26)

The joint predicted intensity is computed by summing over the two intensities with regard to existing targets and spontaneous birthed targets, that is,

$${\upsilon }_{k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k)={\upsilon }_{s,k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k)+{\gamma }_k({\mathbf{x}}_k,{\mathbf{\theta }}_k)$$ (27)

where ${\upsilon }_{s,k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k)$ denotes the predicted intensity from existing targets which is computed by

$$\begin{array}{c}{\upsilon }_{s,k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k)=p_{s,k}\hfill \\ \sum _{j=1}^{J_{k-1}}[w_{k-1}^j\mathcal{N}({\mathbf{x}}_k;{\mathbf{m}}_{s,k|k-1}^j,{\mathbf{P}}_{s,k|k-1}^j)\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_{k|k-1};{\mathbf{\mu }}_{k|k-1}^j,{\kappa }_{k|k-1}^j,{\nu }_{k|k-1}^j,{\mathbf{\Lambda }}_{k|k-1}^j)]\hfill \end{array}$$ (28)

The predicted mean and covariance matrix of the target states is estimated by

$$\begin{array}{c}{\mathbf{m}}_{s,k|k-1}^j=\int f({\mathbf{x}}_{k-1}^j)\mathcal{N}({\mathbf{x}}_{k-1}^j;{\mathbf{m}}_{k-1}^j,{\mathbf{P}}_{k-1}^j)d{\mathbf{x}}_{k-1}^j\hfill \\ {\mathbf{P}}_{s,k|k-1}^j=\int (f({\mathbf{x}}_{k-1}^j)-{\mathbf{m}}_{k-1}^j)(f({\mathbf{x}}_{k-1}^j)-{\mathbf{m}}_{k-1}^j{)}^T\mathcal{N}({\mathbf{x}}_{k-1};{\mathbf{m}}_{k-1}^j,{\mathbf{P}}_{k-1}^j)d{\mathbf{x}}_{k-1}^j+{\mathbf{Q}}_k\hfill \end{array}$$ (29)

And the predicted hyperparameters of the Normal-inverse-Wishart distribution are computed in accordance with the equation (21), that is,

$$\begin{array}{c}{\mathbf{\mu }}_{k|k-1}^j={\mathbf{\mu }}_{k-1}^j\hfill \\ {\kappa }_{k|k-1}^j=\frac{1}{\rho }{\kappa }_{k-1}^j\hfill \\ {\nu }_{k|k-1}^j=\rho {\nu }_{k-1}^j\hfill \\ {\mathbf{\Lambda }}_{k|k-1}^j=\rho {\mathbf{\Lambda }}_{k-1}^j\hfill \end{array}$$ (30)

Measurement update: Given that the joint predicted intensity at time $k$ is represented by a product of Normal-inverse-Wishart and Gaussian mixtures

$$\begin{array}{c}{\upsilon }_{k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k)\hfill \\ =\sum _{j=1}^{J_{k|k-1}}[w_{k|k-1}^j\mathcal{N}({\mathbf{x}}_k;{\mathbf{m}}_{k|k-1}^j,{\mathbf{P}}_{k|k-1}^j)\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_{k|k-1};{\mathbf{\mu }}_{k|k-1}^j,{\kappa }_{k|k-1}^j,{\nu }_{k|k-1}^j,{\mathbf{\Lambda }}_{k|k-1}^j)]\hfill \end{array}$$ (31)

Then the joint posterior intensity is computed by

$${\upsilon }_k({\mathbf{x}}_k,{\mathbf{\theta }}_k)=(1-p_{d,k}){\upsilon }_{k|k-1}({\mathbf{x}}_k,{\mathbf{\theta }}_k)+\sum _{{\mathbf{z}}_k\in {\mathbf{Z}}_k}{\upsilon }_{d,k}({\mathbf{x}}_k,{\mathbf{\theta }}_k|{\mathbf{z}}_k)$$ (32)

where the detected target intensity for each measurement ${\mathbf{z}}_k$ is given by

$$\begin{array}{c}{\upsilon }_{d,k}({\mathbf{x}}_k,{\mathbf{\theta }}_k|{\mathbf{z}}_k)=\hfill \\ \sum _{j=1}^{J_{k|k-1}}[w_k^j({\mathbf{z}}_k)\mathcal{N}({\mathbf{x}}_k;{\mathbf{m}}_k^j({\mathbf{z}}_k),{\mathbf{P}}_k^j)\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_k;{\mathbf{\mu }}_k^j,{\kappa }_k^j,{\nu }_k^j,{\mathbf{\Lambda }}_k^j)]\hfill \end{array}$$ (33)

The fixed point iteration method is employed to compute the above model parameters since it is of asymptotic optimal convergence properties. ³¹ For each new measurement ${\mathbf{z}}_k$ , an iterative PHD measurement update is performed with the initial values of the parameters given as

$$\begin{array}{c}{\mathbf{m}}_k^{j,(0)}({\mathbf{z}}_k)={\mathbf{m}}_{k|k-1}^j,{\mathbf{P}}_k^{j,(0)}={\mathbf{P}}_{k|k-1}^j,\hfill \\ {\mathbf{\mu }}_k^{j,(0)}={\mathbf{\mu }}_{k|k-1}^j,{\kappa }_k^{j,(0)}={\kappa }_{k|k-1}^j+1,\hfill \\ {\nu }_k^{j,(0)}={\nu }_{k|k-1}^j+1,{\mathbf{\Lambda }}_k^{j,(0)}={\mathbf{\Lambda }}_{k|k-1}^j\hfill \end{array}$$ (34)

where the superscript $(i)=1,2,\cdots ,N_{\max }$ is the iteration index with the preset maximum number of iterations $N_{\max }$ . The $(i+1)$ -th iteration process is depicted as follows. And the auxiliary variables, that is, the innovation vector ${\hat{\mathbf{z}}}_k^j$ , the innovation covariance matrix ${\mathbf{T}}_k^j$ and the cross-covariance matrix ${\mathbf{C}}_k^j$ , for each mixture are computed by

$$\begin{array}{c}{\hat{\mathbf{z}}}_k^j=\int h({\mathbf{x}}_k^j)\mathcal{N}({\mathbf{x}}_k^j;{\mathbf{m}}_{k|k-1}^j,{\mathbf{P}}_{k|k-1}^j)d{\mathbf{x}}_k^j\hfill \\ {\mathbf{T}}_k^j=\int (h({\mathbf{x}}_k^j)-{\hat{\mathbf{z}}}_k^j)(h({\mathbf{x}}_k^j)-{\hat{\mathbf{z}}}_k^j{)}^T\mathcal{N}({\mathbf{x}}_k^j;{\mathbf{m}}_{k|k-1}^j,{\mathbf{P}}_{k|k-1}^j)d{\mathbf{x}}_k^j\hfill \\ {\mathbf{C}}_k^j=\int ({\mathbf{x}}_k^j-{\mathbf{m}}_{k|k-1}^j)(h({\mathbf{x}}_k^j)-{\hat{\mathbf{z}}}_k^j{)}^T\mathcal{N}({\mathbf{x}}_k^j;{\mathbf{m}}_{k|k-1}^j,{\mathbf{P}}_{k|k-1}^j)d{\mathbf{x}}_k^j\hfill \end{array}$$ (35)

The estimated measurement noise statistics are updated with estimates of the hyperparameters from the previous iteration as

$$\begin{array}{c}{\hat{\mathbf{r}}}_k^{j,(i+1)}={\mathbf{\mu }}_k^{j,(i)}\hfill \\ {\hat{\mathbf{R}}}_k^{j,(i+1)}=\frac{{\mathbf{\Lambda }}_k^{j,(i)}}{{\nu }_k^{j,(i)}-n_z-1}\hfill \end{array}$$ (36)

The estimated multi-target state is updated by employing standard Gaussian filters with known measurement noise statistics as

$$\begin{array}{c}{\mathbf{S}}_k^{j,(i+1)}={\mathbf{T}}_k^j+{\hat{\mathbf{R}}}_k^{j,(i+1)}\hfill \\ {\mathbf{K}}_k^{j,(i+1)}={\mathbf{C}}_k^j[{\mathbf{S}}_k^{j,(i+1)}{]}^{-1}\hfill \\ {\mathbf{m}}_k^{j,(i+1)}({\mathbf{z}}_k)={\mathbf{m}}_{k|k-1}^j+{\mathbf{K}}_k^{j,(i+1)}({\mathbf{z}}_k-{\hat{\mathbf{z}}}_k^j-{\mathbf{r}}_k^{j,(i+1)})\hfill \\ {\mathbf{P}}_k^{j,(i+1)}={\mathbf{P}}_{k|k-1}^j-{\mathbf{K}}_k^{j,(i+1)}{\mathbf{S}}_k^{j,(i+1)}[{\mathbf{K}}_k^{j,(i+1)}{]}^T\hfill \end{array}$$ (37)

The estimated hyperparameters are updated by

$$\begin{array}{c}{\mathbf{\mu }}_k=\frac{{\kappa }_{k|k-1}^j}{{\kappa }_{k|k-1}^j+1}{\mathbf{\mu }}_{k|k-1}^j+\frac{1}{{\kappa }_{k|k-1}^j+1}{\mathbf{e}}_k^{j,(i+1)}\hfill \\ {\mathbf{\Lambda }}_k^{j,(i+1)}={\mathbf{\Lambda }}_{k|k-1}^j+\frac{{\kappa }_{k|k-1}^j}{{\kappa }_{k|k-1}^j+1}\hfill \\ ({\mathbf{e}}_k^{j,(i+1)}-{\mathbf{\mu }}_{k|k-1}^j)({\mathbf{e}}_k^{j,(i+1)}-{\mathbf{\mu }}_{k|k-1}^j{)}^T\hfill \end{array}$$ (38)

where the measurement residual is computed by

$${\mathbf{e}}_k^{j,(i+1)}={\mathbf{z}}_k-h({\mathbf{m}}_k^{j,(i+1)}({\mathbf{z}}_k))$$ (39)

Suppose that the algorithm converges after $N$ iterations, then the parameters are determined by

$$\{w_k^{j,(N)}({\mathbf{z}}_k),{\mathbf{m}}_k^{j,(N)}({\mathbf{z}}_k),{\mathbf{P}}_k^{j,(N)},{\mathbf{\mu }}_k^{j,(N)},{\kappa }_k^j,{\nu }_k^j,{\mathbf{\Lambda }}_k^{j,(N)}{\}}_{j=1}^{J_{k|k-1}}$$ (40)

where the updated weights of the posterior mixture components are calculated by

$$w_k^j({\mathbf{z}}_k)=\frac{p_{d,k}{w}_{k|k-1}^j\mathcal{N}({\mathbf{z}}_k;{\hat{\mathbf{z}}}_k^j+{\mathbf{r}}_k^j,{\mathbf{S}}_k^{j,(N)})}{{\eta }_k({\mathbf{z}}_k)+p_{d,k}\sum _{l=1}^{J_{k|k-1}}{w}_{k|k-1}^l\mathcal{N}({\mathbf{z}}_k;{\hat{\mathbf{z}}}_k^l+{\mathbf{r}}_k^j,{\mathbf{S}}_k^{l,(N)})}$$ (41)

Note that if nonlinear functions are used as dynamic and measurement models, then the general Gaussian integrals in equations (29) and (35) can be approximated by utilizing commonly used variants of non-linear Gaussian filters, such as the extended KF, the unscented KF, and the cubature KF.

Algorithm summarization

Additional steps including parameter initialization, mixture reduction, and state extraction have to be implemented when realizing the entire tracking procedure. The number of the mixture components after each PHD recursion will increase exponentially as time progresses. Consequently, similarly as in Granström and Orguner, ³² mixture reduction process is performed to ensure the tractable computational cost. This is done by first removing the components that have weights below a truncation threshold and then merging the components with reasonable Mahalanobis distance from the largest component. The multi-target state is finally extracted from the reduced mixture components that are corresponding to the peaks of the posterior cardinality. In conclusion, the pseudo-code of the proposed filter is summarized in Algorithm 1.

Algorithm 1. Pseudo-code of the proposed PHD filter

Step 1: Parameter initializationSet initial hyperparameters of the measurement noise statistics:  $$p({\mathbf{\theta }}_0)=\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_0;{\mathbf{\mu }}_0,{\kappa }_0,{\nu }_0,{\mathit{\Lambda }}_0).$$ Set initial joint PHD intensity:  $${\upsilon }_0({\mathbf{x}}_0,{\mathbf{\theta }}_0)=\sum _{j=1}^{J_0}[w_0^j\mathcal{N}({\mathbf{x}}_0;{\mathbf{m}}_0^j,{\mathbf{P}}_0^j)\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_0;{\mathbf{\mu }}_0^j,{\kappa }_0^j,{\nu }_0^j,{\mathit{\Lambda }}_0^j].$$ Step 2: PHD recursionfor $k=1:K$ doCompute the predicted PHD intensity via (26)-(30), and mixtures are represented by  $$\{w_{k|k-1}^j,{\mathbf{m}}_{k|k-1}^j,{\mathbf{P}}_{k|k-1}^j,{\mathbf{\mu }}_{k|k-1}^j,{\kappa }_{k|k-1}^j,{\nu }_{k|k-1}^j,{\mathit{\Lambda }}_{k|k-1}^j{\}}_{j=1}^{J_{k|k-1}}.$$ Set the initial conditions for the iterations via (34);repeatCompute the estimated measurement noise statistics $({\mathbf{r}}_k^j,{\mathbf{R}}_k^j)$ via (36).Compute the updated multi-target states via (37).until convergedCompute the updated mixture weights via (41), and mixtures are represented by $$\{w_k^j,{\mathbf{m}}_k^j,{\mathbf{P}}_k^j,{\mathbf{\mu }}_k^j,{\kappa }_k^j,{\nu }_k^j,{\mathit{\Lambda }}_k^j{\}}_{j=1}^{J_k}.$$ end forStep 3: Mixture reductionPrune and merge the mixtures for computational feasibility.Step 4: State extractionExtract states corresponding to the peaks of the posterior cardinality.

Simulation results

Simulation configuration

In this simulation, the performance of the proposed filter is evaluated in a 2-dimensional scenario, where the number of targets is time-varying. As shown in Figure 1, the targets are observed in a surveillance region with the scale $400\mathrm{m}\times 300\mathrm{m}$ . The target state is denoted by a vector ${\mathbf{x}}_k=[p_{x,k},{\stackrel{·}{p}}_{x,k},p_{y,k},{\stackrel{·}{p}}_{y,k}]$ , which is composed of target position $[p_{x,k},p_{y,k}]$ and velocity $[{\stackrel{·}{p}}_{x,k},{\stackrel{·}{p}}_{y,k}]$ . The true trajectories of the targets are generated by a nearly constant velocity model with the matrices

$$\begin{array}{c}{F}_k=\left[\begin{array}{cccc}1\hfill & \mathrm{\Delta }\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & \mathrm{\Delta }\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right],{\mathbf{Q}}_k={\sigma }_v^2\left[\begin{array}{cccc}\frac{{\mathrm{\Delta }}^4}{4}\hfill & \frac{{\mathrm{\Delta }}^3}{2}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{{\mathrm{\Delta }}^3}{2}}& {\mathrm{\Delta }}^2\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{{\mathrm{\Delta }}^4}{4}\hfill & \frac{{\mathrm{\Delta }}^3}{2}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{{\mathrm{\Delta }}^3}{2}\hfill & {\mathrm{\Delta }}^2\hfill \end{array}\right]\hfill \end{array}$$

where the scanning period $\mathrm{\Delta }=1s$ , and the standard deviation of the process noise ${\sigma }_v=1$  m/s. Target positions are observed by the following linear Gaussian model with unknown measurement noise statistics:

$${\mathbf{z}}_k=\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \end{array}\right]{\mathbf{x}}_k+{\mathbf{w}}_k,{\mathbf{w}}_k~\mathcal{N}({\mathbf{r}}_k,{\mathbf{R}}_k)$$

where the ground truth of the measurement noise statistics are set as ${\mathbf{r}}_{\mathrm{true}}=[9,9]$ and ${\mathbf{R}}_{\mathrm{true}}=\mathrm{diag}([1,1])$ . Note that the true value of the noise mean is beyond three standard deviations to ensure that the impact of the non-zero noise mean is not negligible. As summarized in Table 1, each target has a different survival time and initial position. There are a maximum of four moving targets during the time period [20, 80]. The spontaneous birth intensity is given by

$$\begin{array}{c}{\gamma }_k({\mathbf{x}}_k,{\mathbf{r}}_k,{\mathbf{R}}_k)=w_{\gamma ,k}\hfill \\ \sum _{j=1}^4\mathcal{N}({\mathbf{x}}_k;{\mathbf{m}}_{\gamma ,k}^j,{\mathbf{P}}_{\gamma ,k})\mathcal{N}\mathcal{I}\mathcal{W}({\mathbf{\theta }}_k;{\mathbf{\mu }}_{\gamma ,k}^j,{\kappa }_{\gamma ,k}^j,{\nu }_{\gamma ,k}^j,{\mathbf{\Lambda }}_{\gamma ,k}^j)\hfill \end{array}$$

where the weight $w_{\gamma ,k}=0.03$ , the means of Gaussian components ${\mathbf{m}}_{\gamma ,k}^1=[-140,2.5,120,-2.5]$ , ${\mathbf{m}}_{\gamma ,k}^2=[-180,4,80,-2]$ , ${\mathbf{m}}_{\gamma ,k}^3=[-180,2.5,-80,2]$ , and ${\mathbf{m}}_{\gamma ,k}^4=[-140,3,-120,1.5]$ , the identical covariance ${\mathbf{P}}_{\gamma ,k}=\mathrm{diag}([25,10,25,10])$ . The initial hyperparameters of the Normal-inverse-Wishart components are set as ${\mathbf{\mu }}_0=[4,4]$ , ${\kappa }_0=5$ , ${\nu }_0=5$ , and ${\mathbf{\Lambda }}_0=\mathrm{diag}([8,8])$ . Besides, other parameters used in the simulations are listed in Table 2.

Figure 1.

True trajectories of the targets moving with a constant velocity model.

Table 1.

Description of the targets.

Target index	Initial position	Birth time	End time
1#	(−140, 120)	1	100
2#	(−180, 80)	20	100
3#	(−180, −80)	1	80
4#	(−140, −120)	10	90

Table 2.

Other parameters used in the simulations.

Parameter description	Symbol	Values
Total time step	$K$	100
Target survival probability	$p_{s,k}$	0.99
Target detection probability	$p_{d,k}$	0.98
Average clutter number	${\lambda }_c$	10
The pruning threshold	$T_{\mathrm{th}}$	0.001
The merging threshold	$U_{\mathrm{th}}$	5
The weight threshold	$W_{\mathrm{th}}$	0.5
Maximum Gaussian components	$J_{\max }$	100
Maximum iteration number	$N$	10
The forgetting factor	$\rho$	0.97

General tracking performance

In Figure 2, the resulting target positions estimated by the AGM-PHD filter in one typical trial are presented. It can be easily seen that almost all the positions of the targets are estimated accurately in both x and y coordinates. The weighted average of the estimated measurement noise statistics from each extracted mixture is illustrated in Figures 3 and 4. The simulation results suggest that both the estimated means and standard deviations of the measurement noise converge to their true values after a few time steps. The reasons for the similar convergence speed of the means and standard deviations are that they are coupled together through one single Normal-inverse-Wishart distribution.

Figure 2.

Position estimates of AGM-PHD filter in x and y coordinates.

Figure 3.

Weighted average of estimated measurement noise means.

Figure 4.

Weighted average of estimated measurement noise standard deviations.

Comparisons with benchmark filters

To demonstrate the advantages of the AGM-PHD filter, three different variations of GM-PHD filters are also implemented and evaluated in the simulations. The highlights of these filters are described as follows.

• GMTS-PHD: the standard GM-PHD filter with the ground truth of the measurement noise statistics, that is, the means and standard deviations of the measurement noise are chosen as 9.0 and 1.0, respectively.

• GMTM-PHD: the standard GM-PHD filter with the ground truth of the measurement noise mean only, while the standard deviations of the measurement noise are chosen as 2.0 to match with the initial standard deviation of AGM-PHD.

• GM-VBPHD: the adaptive PHD filter proposed in Wu et al. ²² where only the covariance of the measurement noise is assumed unknown, while the mean of the measurement noise is chosen as 8.0 to allow a reasonable distance from the ground truth.

In all the simulations, the estimated accuracy of the target position is evaluated quantitatively by the optimal sub-pattern Assignment (OSPA). ³³ Similar to the root mean square error (RMSE) metric used to represent the differences between two vectors, the OSPA metric can be used to measure the differences between two random finite sets. In each simulation run, two hundred independent Monte Carlo trials are conducted. For the OSPA metric, the cut off factor $c=100$ , and the order is set as $p=2$ . The comparison of the averaged OSPA distance for all the four filters is presented in Figure 5. The result proves that the GMTS-PHD filter obtains the least estimation error since there is no bias in the measurement noise statistics in this case. It can also be seen from Figure 5 that the tracking performance of the AGM-PHD filter is the most comparable with that of the GMTS-PHD filter. This is due to the fact that the AGM-PHD filter can adapt to unknown measurement noise statistics with a few initial time steps. Besides, the GMTM-PHD filter achieves smaller errors than the GM-VBPHD. It implies that in the measurement update step, the mean of the measurement noise has more important impacts on state estimation than the covariance when its value exceeds three times the standard deviation. The comparison of average cardinality statistics calculated by the different filters are shown in Figure 6, where the quite similar result as the OSPA comparison can be seen. The average number of the targets determined by both the ideal GMTS-PHD filter and the AGM-PHD filter are very close to the ground truth. However, the GMTS-PHD filter and GM-VBPHD fail to be in accordance with the true number because only the mean or the covariance of the measurement noise is considered.

Figure 5.

Comparison of OSPA distances for different filters.

Figure 6.

Comparison of average cardinalities for different filters.

To compare the computational costs of these filters, the average running time of Monte Carlo trail of different filters are listed in Table 3. It is noted that the running time required for each filter are measured relative to that of the GMTS-PHD filter. As can be seen from the Table, the GMTM-PHD filter demands almost the same running time as the GMTS-PHD filter. The two VB based filters take more than 25% longer running time. This is mainly because that an additional iteration step has to be carried out to estimate the unknown model hyperparameters during the measurement update stage. The computational cost of the AGM-PHD filter is slightly higher than the GM-VBPHD filter since it involves more number of hyperparameters that need to be estimated in the AGM-PHD filter.

Table 3.

Comparison of computational costs for different filters.

Index	Filter name	Relative running time
1	GMTS-PHD	1.00
2	GMTM-PHD	1.02
3	GM-VBPHD	1.27
4	AGM-PHD	1.35

The impacts of different average clutter rate on the tracking performance are also evaluated. All the filters are carried out under six different levels of clutter rate whose value is taken from the interval [5, 30] with step 5.0. For each level of clutter rate, the average OSPA metric with regard to different filters are plotted in Figure 7. It can be observed that the estimation accuracy concerning all the filters degrades as the level of clutter rate increases. However, the differences of the OSPA distance between the AGM-PHD filter and the GMTS-PHD filter remain below 3.0 at all clutter levels. Meanwhile, the increment of the OSPA distance of the GMTM-PHD filter and GM-VBPHD filter are larger than 12.0, which is nearly 1.5 times of the other two filters. This might due to the reason that more spurious measurements are treated as targets when the nominal statistics are larger than the corresponding true values.

Figure 7.

Comparison of average OSPA metric under different average clutter rate.

Evaluation with a constant turn model

To demonstrate the performance of the AGM-PHD filter for tracking targets with variable parameters of motion model, the targets are assumed to move within the same surveillance region following a nearly constant turn (CT) model ³⁴ with known but varying turn rate. The tracking setup of the simulation keeps the same except the dynamic motion model, which is given as follows:

$${\mathbf{x}}_k=F\left({\omega }_k\right){\mathbf{x}}_{k-1}+{\mathbf{v}}_k,{\mathbf{v}}_k~\mathcal{N}(0,{\mathbf{Q}}_k)$$

where the transition matrix is given as a function of the turn rate:

$$F({\omega }_k)=\left[\begin{array}{cccc}1\hfill & \frac{\sin {\omega }_k\mathrm{\Delta }}{{\omega }_k}\hfill & 0\hfill & -\frac{1-\cos {\omega }_k\mathrm{\Delta }}{{\omega }_k}\hfill \\ \multicolumn{1}{c}{0}& \frac{1-\cos {\omega }_k\mathrm{\Delta }}{{\omega }_k}\hfill & 1\hfill & \frac{\sin {\omega }_k\mathrm{\Delta }}{{\omega }_k}\hfill \\ \multicolumn{1}{c}{0}& \cos {\omega }_k\mathrm{\Delta }\hfill & 0\hfill & -\sin {\omega }_k\mathrm{\Delta }\hfill \\ \multicolumn{1}{c}{0}& \sin {\omega }_k\mathrm{\Delta }\hfill & 0\hfill & \cos {\omega }_k\mathrm{\Delta }\hfill \end{array}\right]$$

the turn rate ${\omega }_k={\omega }_{c,k}+u_k\mathrm{\Delta }$ , $\mathrm{\Delta }=1\mathrm{s}$ and $u_k$ denotes angle rate noise with Gaussian distribution $\mathcal{N}(u_k;0,(\pi /180{)}^2)$ , ${\omega }_{c,k}$ follows a piecewise function

$${\omega }_{c,k}=\{\begin{array}{cc}0,\hfill & \mathrm{for}1\le \mathrm{k}\le 40,60<\mathrm{k}\le 100\hfill \\ \multicolumn{1}{c}{-\pi /10,}& \mathrm{for}40<\mathrm{k}\le 50\hfill \\ \multicolumn{1}{c}{\pi /10,}& \mathrm{for}50<\mathrm{k}\le 60\hfill \end{array}$$

In Figure 8, the ground truth trajectory of each target together with its start and end positions are presented. All the targets travel along nearly straight lines at the time intervals [1 s, 40 s] and [60 s, 100 s], while making two converse turns at time 40 s and 50 s. During the turning periods, tracking the targets correctly becomes more difficult as the targets’ trajectories change significantly. Same as in the first simulation, 200 independent Monte Carlo trials are conducted for all four PHD filters. Consequently, the result in Figure 9 shows that more erroneous estimates appear around the turning points. Nevertheless, the proposed adaptive PHD filter still exhibit acceptable tracking performance as a whole. This implies that the AGM-PHD filter gains certain robustness against sudden changes in the control noise of non-linear dynamic models. The performance comparison between the filters is presented in Figures 10 and 11. The OSPA distance rises during the periods when the targets make the wide turns. It is also observed in the results, the AGM-PHD filter achieves smaller OSPA distance and more accurate cardinality estimates than the other two counterparts.

Figure 8.

True trajectories of the targets moving with a constant turn model.

Figure 9.

Target’s position estimates with a constant turn model.

Figure 10.

Comparison of OSPA metric between different filters with constant turn model.

Figure 11.

Comparison of average cardinalities between different filters with constant turn model.

Conclusion

In conclusion, the classical PHD filter was devised to improve the adaptivity for MTT applications where the measurement noise statistics are unavailable. During the process of the target state estimation, the measurement noise statistics are also estimated and refined iteratively according to the derived updating rule for the hyperparameters of Normal-inverse-Wishart distribution. The updating rule can be obtained based on the VB approximation or the Bayesian conjugate prior heuristics. Simulation results show the proposed adaptive PHD filter achieves comparable tracking performance with the ideal benchmark filter under different cluttering settings. Moreover, the proposed method has a certain degree of robustness even when the motion trajectories of the targets change abruptly.