Robust Generalized Labeled Multi-Bernoulli Filter and Smoother for Multiple Target Tracking using Variational Bayesian

Abstract

Multiple target tracking mainly focuses on tracking unknown number of targets in the complex environment of clutter and missed detection. The generalized labeled multi-Bernoulli (GLMB) filter has been shown to be an effective approach and attracted extensive attention. However, in the scenarios where the clutter rate is high or measurement-outliers often occur, the performance of the GLMB filter will significantly decline due to the Gaussian-based likelihood function is sensitive to clutter. To solve this problem, this paper presents a robust GLMB filter and smoother to improve the tracking performance in the scenarios with high clutter rate, low detection probability, and measurement-outliers. Firstly, a Student-T distribution variational Bayesian (TDVB) filtering technology is employed to update targets' states. Then, The likelihood weight in the tracking process is deduced again. Finally, a trajectory smoothing method is proposed to improve the integrative tracking performance. The proposed method are compared with recent multiple target tracking filters, and the simulation results show that the proposed method can effectively improve tracking accuracy in the scenarios with high clutter rate, low detection rate and measurement-outliers. Code is published on GitHub.

1. Introduction

The significant difficulty of multiple target tracking (MTT) is to track an unknown number of targets in the complex scenarios with strong interference, missed detection, and clutters. In the past decade, MTT technology has been widely used in various fields, such as radar networks [31], vehicle tracking [32], robot vision [33], etc. The most important theories in MTT mainly include the random finite set (RFS) [1] and the labeled RFS (LRFS) [2]. A lot of work has been done based on these two theories, which make RFS and LRFS become hot MTT frameworks in recent years.

Based on the Finite Set Statistics (FISST), Mahler presented the RFS theory and its first-order moment, i.e., the probability hypothesis density (PHD) filter [1]. The most important implementation of the PHD framework is the Gaussian mixture PHD (GM-PHD) filter [3], which simplifies the difficulty of integral solution of a PHD by the Gaussian mixture method. The Poisson multi-Bernoulli mixture (PMBM) filter was proposed in [4], which is the closed-form solution for Poisson birth model and can provide data association result. The PMBM model has been combined with Gamma Gaussian inverse Wishart (GGIW) model and applied to extended target tracking [26]. Other important RFS based tracking methods include cardinalized PHD (CPHD) [5], particle PHD [6], etc.

However, previous studies have found that the RFS based filters cannot provide track association results, which greatly limits the application of RFS based filters. Recently, some RFS methods have been proposed to solve this problem, such as [4] and [27]. Another important solution is the LRFS method [2], which establishes a strict mathematical model for the trajectory. The important implementations of LRFS are delta GLMB (δ-GLMB) [7] and labeled multi-Bernoull (LMB) [8] filters, which can achieve more accurate tracking result than PHD filters and provide target trajectories. Therefore, the GLMB/LMB filtering technologies are state-in-art tracking frameworks in recent years. Other important discussions on GLMB/LMB can be found in [9]-[17]. Compared with RFS based filters, LRFS based filters achieve better tracking accuracy, but need more computational cost. Therefore, [18] presented an efficient GLMB (E-GLMB) filter to reduce the computational complexity using the Gibbs sampling method. The E-GLMB filter is an important work in recent years and has shown superior performance in MTT. However, in complex scenarios, such as scenarios with high clutter rate, high missed detection probability, and measurement-outliers, the performance of the E-GLMB filter will decrease significantly. Note that the measurement-outliers are generated by the sudden interference of the sensor at a certain time, i.e. the mathematical expectation of the measurement noise increases at that time.

Tracking in complex scenarios is an important problem of MTT. Particle filtering technology is a method to solve the complex or nonlinear tracking problem, and it has been discussed in GLMB [23] and PHD [24] filters. However, due to the large computational cost, it is difficult to implement a particle filter in real scenarios. Student-T distribution filtering technology is an efficient method to handle the complex linear scenarios. Compared with the Gaussian distribution, the Student-T distribution achieves better performance in dealing with outliers, thus it can improve the tracking accuracy in complex scenarios. The Student-T distribution filtering technology is a hot topic in recent years [28]-[29] and Student-T-distribution-based PHD filters have been derived in [25].

This paper proposes a Student-T distribution variational Bayesian E-GLMB (TDVB-E-GLMB) filter and a multi-target smoother. The TDVB [19] is a robust filtering technology compared with Bayesian filtering framework, thus the proposed method can improve the tracking performance in complex scenarios. For example, when the E-GLMB filter is applied to vehicle tracking, the strong light, fog and other environmental factors may lead to a large number of clutter in the radar system, resulting in the decline of tracking accuracy. The proposed method use the TDVB framework to improve the robustness of the filtering process, thus the influence of clutter on accuracy will be reduced. Moreover, a corresponding smoothing method is proposed for GLMB framework to improve the overall tracking accuracy. The contributions of this work are summarized as follows:

(1)This work applies TDVB filtering technology [19] to the E-GLMB filter. A new likelihood function is derived to replace the Gaussian likelihood. Therefore, the proposed TDVB-E-GLMB filter can significantly improve the tracking accuracy in complex scenarios.

(2)A multi-target TDVB smoothing method is proposed for GLMB framework. This method considers the missed detected targets and pseudo targets based on trajectories provided by the GLMB framework. Therefore, the proposed TDVB smoother can significantly improve the performance of a GLMB filter.

Code is published on Github [22]. Further analysis of the proposed TDVB-E-GLMB filter can be done through this code.

The paper is organized as follows. Section II introduces the E-GLMB filter, variational Bayesian inference, and TDVB filtering technology. Section III introduces the proposed TDVB-E-GLMB filter. Section IV introduces the proposed multi-target smoothing method for a GLMB filter. Section V shows the simulation results. Section VI contains our conclusions.

2. Background

2.1 Problem Statements

In general, the system function of the tracking process is described as

x_k = f(x_k-1 , v_k-1),       (1)

z_k = τ(x_k , w_k),       (2)

where f(⋅) denotes the state model and τ(⋅) denotes the observation model. v_k−1 and w_k are corresponding noises. In a standard RFS or LRFS filter, the state and observation models are usually defined as Gaussian model, due to the simple mathematical form and low computational complexity. Also, the state updating process is a Bayesian recursive process. Therefore, the state model and measurement model can be written as

x_k | x_k-1 ~ N(f(x_k-1), Q_k),       (3)

z_k | x_k ~ N(τ(x_k), R_k).       (4)

However, the Gaussian model also leads to the lack of robustness in complex scenarios such as low detection rate, high clutter rate, and measurement-outliers. The measurement-outliers is the measurement at a certain scan when expectation of measurement error increases by the strong interference. For example, the expectation of measurement noise is usually R₀ , but the expectation increases to a large R₁ at the time of strong interference. The measurement obtained by the sensor system with R₁ as the expected error calls measurement-outlier.

In the complex scenarios, using Gaussian model to describe the measurement generation will lead to the lack of robustness. Compared with Gaussian distribution, the Student-T is a distribution with lower peak and higher tail, thus it is easier to handle outliers. In [25] and [28]-[29], the Student-T distribution was employed to handle the measurement-outliers and the results show its better performance. However, the Student-T distribution can only perform well for measurement-outliers, but it is not robust to other complex scenarios. Moreover, compared with the Bayesian recursive process, the VB recursive process has more robustness. Therefore, the combination of Student-T distribution and VB process will greatly improve the tracking robustness, and this work has been discussed in [19].

In this work, we employ the TDVB model proposed in [19] to modify the E-GLMB filter. The state and measurement models are defined as the Gaussian model with the random auxiliary variable, i.e. (18)-(19). Note that this measurement model allows for measurement-outliers. The state updating model is modified by TDVB filtering technology. The proposed new likelihood is calculated by KL divergence, i.e. formula (35)-(37).

2.2 The E-GLMB Filter

The E-GLMB filter uses the Gibbs sampling and joint prediction-update approaches to reduce the high computational cost of the δ-GLMB filter. To achieve the goal of joint prediction and update, the E-GLMB filtering density is given by

\(\begin{aligned}\pi_{Z_{+}}(\mathrm{X}) \propto \Delta(\mathrm{X}) \sum_{\xi, I_{+}, \theta_{+}} w_{Z_{+}}^{\left(\xi, I_{+}, \theta_{+}\right)} \delta_{I_{+}}[L(\mathrm{X})]\left[p_{Z_{+}}^{\left(\xi, \theta_{+}\right)}\right]^{\mathrm{X}}\end{aligned}\)       (5)

where, ξ denotes the history of association. θ denotes the association map. L(Χ) denotes the label set and I denotes the set of object labels. δ_I+[L(X)] means that if the new label set I₊ = L(Χ) , the value is 1, otherwise, 0. p_Z+^{(ξ ,θ₊)} is the component function with its weight w_Z+^{(ξ,I₊,θ₊)} . Χ is the target state, and exponential notation f^Χ means that ∏_x∈Χ f(x) . ∆(Χ) means that δ_|X|[|L(X)|] , where |⋅| denotes the number of elements in a set. Each part of the formula (5) is given by

\(\begin{aligned} w_{Z_{+}}^{\left(\xi, I_{+}, \theta_{+}\right)} & =1_{\Theta_{+}\left(I_{+}\right)}\left(\theta_{+}\right)\left[1-\bar{P}_{S}^{(\xi)}\right]^{I-I_{+}}\left[\bar{P}_{S}^{(\xi)}\right]^{I \cap I_{+}} \\ & \times\left[1-r_{B,+}\right]^{\mathrm{B}_{+}-I_{+}} r_{B,+}^{\mathrm{B}_{+} \cap I_{+}}\left[\bar{\psi}_{Z_{+}}^{\left(\xi, \theta_{+}\right)}\right]^{I_{+}}\end{aligned}\)       (6)

\(\begin{aligned}\bar{P}_{S}^{(\xi)}(h)=\left\langle p^{(\xi)}(\cdot, h), P_{S}(\cdot, h)\right\rangle\end{aligned}\),       (7)

\(\begin{aligned}\bar{\psi}_{Z_{+}}^{\left(\xi, \theta_{+}\right)}\left(h_{+}\right)=\left\langle\bar{p}_{+}^{(\xi)}\left(\cdot, h_{+}\right), \psi_{Z_{+}}^{\left(\theta_{+}\left(h_{+}\right)\right)}\left(\cdot, h_{+}\right)\right\rangle\end{aligned}\),       (8)

\(\begin{aligned} \bar{p}_{+}^{(\xi)}\left(x_{+}, h_{+}\right) & =1_{C_{+}}\left(h_{+}\right) \frac{\left\langle P_{S}\left(\cdot, h_{+}\right) f_{+}\left(x_{+} \mid \cdot, h_{+}\right), p^{(\xi)}\left(\cdot, h_{+}\right)\right\rangle}{\bar{P}_{S}^{(\xi)}\left(h_{+}\right)}, \\ & +1_{\mathrm{B}_{+}}\left(h_{+}\right) p_{B,+}\left(x_{+}, h_{+}\right)\end{aligned}\)       (9)

where, I ∈ F (C), I₊ ∈ F(C₊), θ₊ ∈ Θ₊ . 1_C+(h₊) means that if h₊ ⊂ C₊ , the value is 1, otherwise, 0. Β denotes the track label space for target birth. r_B,+ is probability that a new object with label is born. ψ_Z+^{(θ₊(h₊))}(⋅,h₊) is a likelihood function which is defined in [18, section II.A]. P_S(⋅,h) denotes the target survival probability. In the specific implementation process, the GLMB component is a Gaussian mixture form

\(\begin{aligned}p_{Z_{+}}^{\left(\xi, \theta_{+}\right)}\left(x, h \mid Z_{+}\right)=\sum_{i=1}^{J^{(\xi)}(h)} \frac{w_{Z_{+}, i}^{(\xi, \theta)}(h)}{\eta_{Z_{+}}^{(\xi, \theta)}(h)} \mathrm{N}\left(x ; m_{Z_{+}, i}^{(\xi, \theta)}(h), P_{Z_{+}, i}^{(\xi, \theta)}(h)\right) \end{aligned}\)       (10)

where, N(⋅) denotes the Gaussian distribution. m_Z+,i^(ξ,θ)(h) and P_Z+,i^(ξ,θ)(h) denote the Gaussian mean and covariance matrix. w_Z+,i^(ξ,θ)(h) is the weight of the Gaussian mixture. η_Z+,i^(ξ,θ)(h) denotes the normalization.

2.3 The Variational Bayesian Inference

The key method of the VB technology is to use a function Ψ(Z) to approximate the posterior probability density g ( Z | X ) by computing the Kullback-Leibler (KL) divergence, i.e.,

ln g(X) = KL(Ψ || g) + L(Ψ).       (11)

Assuming that \(\begin{aligned}\Psi(Z)=\prod_{i} \Psi\left(z_{i}\right) \end{aligned}\) , then L(Ψ) can be written as

\(\begin{aligned} L(\Psi)= & -\mathrm{KL}\left(\Psi\left(z_{j}\right) \| \Psi^{*}\left(z_{j}\right)\right) \\ & +\prod_{i: i \neq j} H\left(\Psi\left(z_{j}\right)\right)-\text { const }\end{aligned}\).       (12)

Obviously, the function L(Ψ) will maximize when −KL(Ψ(z_j)||Ψ^*(z_j)) = 0 , thus we obtain

\(\begin{aligned}\Psi\left(z_{j}\right)=\Psi^{*}\left(z_{j}\right)=\frac{\exp \left(E_{i \neq j}[\ln g(Z, X)]\right)}{\text { normalize constant }}\end{aligned}\).       (13)

As a result, Ψ(Z) can be approximated to g( Z | X ).

2.4 The TDVB Filtering Technology

Compared with the Gaussian distribution, the Student-T distribution achieves better performance in dealing with complex scenarios with strong interference. Also, the variational Bayesian processes can deal with nonlinear cases by iteration. Therefore, the Student-T distribution variational Bayesian filtering technology is a robust approach to handle non-linear scenarios, especially in measurement-outlier scenarios.

According to [19], the system model is given by

x_k | x_k-1 ~ N(f(x_k-1), Q_k),       (14)

z_k | x_k ~ St(τ(x_k), R_k, v),       (15)

where, x_k denotes the target state. z_k denotes the measurement. St(⋅) denotes the Student-T distribution. f(⋅) and τ(⋅) are dynamic and observation models, respectively. Q_k and R_k are process and measurement noise, respectively. ν is the degree of freedom of Student-T distribution. The probability density of Student-T distribution can be written as

\(\begin{aligned}\operatorname{St}\left(\tau\left(x_{k}\right), R_{k}, v\right) \propto\left(1+\frac{1}{v} \aleph^{T} R_{k}^{-1} \aleph\right)^{-\frac{v+d}{2}}\end{aligned}\)       (16)

ℵ = (z_k - τ(x_k)),       (17)

where, d is the dimensions of physical space. The formula (16) can be approximately expressed as a hierarchical Gaussian form by employs an auxiliary random variable λ_k , thus formula (15) become

\(\begin{aligned}z_{k} \mid x_{k}, \lambda_{k} \sim \mathrm{N}\left(\tau\left(x_{k}\right), \frac{1}{\lambda_{k}} R_{k}\right) \end{aligned}\),       (18)

λ_k ~ Gamma ( v/2, v/2 ),       (19)

where, Gamma (⋅) denotes the Gamma distribution. Since (14) is a standard Gaussian form, the state prediction steps can use a standard Kalman filtering technology. According to (18)-(19), the Student-T distribution is approximated to a Gaussian distribution with the random auxiliary variable, thus the sate updating can be given by

m_k = m_k|k-1 + K_k (z_k - µ_k),       (20)

P_k = P_k|k-1 + K_kS_kK_k^T,       (21)

where,

μ_k = ∫τ(x_k)N(x_k | m_k|k-1, P_k|k-1)dx_k,       (22)

\(\begin{aligned} S_{k}=\int & \left(\tau\left(x_{k}\right)-\mu_{k}\right)\left(\tau\left(x_{k}\right)-\mu_{k}\right)^{T} \\ & \times \mathrm{N}\left(x_{k} \mid m_{k \mid k-1}, P_{k \mid k-1}\right) d x_{k}+\frac{1}{\lambda_{k}} R_{k}, \end{aligned}\)       (23)

K_k = ∫(x_k - m_k|k-1)(τ(x_k) - μ_k)^T× N(x_k | m_k|k-1, P_k|k-1)dx_k · S_k^-1,       (24)

\(\begin{aligned}\begin{array}{c}\bar{\gamma}_{k}=\operatorname{tr}\left(\int\left(z_{k}-\tau\left(x_{k}\right)\right)\left(z_{k}-\tau\left(x_{k}\right)\right)^{T},\right. \\ \left.\times \mathrm{N}\left(x_{k} \mid m_{k}, P_{k}\right) d x_{k} R_{k}^{-1}\right)\end{array}\end{aligned}\)       (25)

\(\begin{aligned}\lambda_{k}=(v+d) /\left(v+\bar{\gamma}_{k}\right),\end{aligned}\)       (26)

where, tr(⋅) denotes the trace of the matrix. d is the dimensions of physical space. Note that, formulas (20)-(26) represent an iterative process using VB filtering technology. Therefore, given the number of iterations, formulas (20)-(21) will converge through a changing λ_k .

3. The TDVB-E-GLMB Filter

The use of Gaussian component leads to the decrease of E-GLMB filter’s performance in the complex scenarios. This section proposed an implementation approach using TDVB filtering technology to enhance the robustness of the E-GLMB filer.

3.1 The Key Method

A GLMB filter uses the Kalman filtering technology to update the kinematical states of components, and uses the Gaussian distribution to calculate the likelihood. Then, this information will be employed to update the track table. The form of a track table is shown in Fig. 1. I^(h) is the set of labels, w^(h) is the weight of a track, and p^(h) denotes the densities. In the δ-GLMB filter, all the target kinematical states are stored in another memory, but the E-GLMB filter does not need to store all of them by joint prediction and update. The density p^(h) is related to the method of the filtering technology for updating target kinematical state, and its physical meaning is the spatial likelihood of the measurement.

Fig. 1. Example of track component.

In a standard GLMB filter, p^(h) is the Gaussian mixture form and the details can be found in [7, section IV.B]. However, it is found in the simulation that the robustness of Gaussian mixture form is insufficient in complex scenarios, which leads to a significant decline in tracking accuracy. Therefore, the key method of proposed approach is to employ the more robust TDVB filtering technology to update the target kinematical state. Then, deriving a corresponding p^(h)to improve the accuracy of the track table. Note that the proposed method is derived as a linear filter, and therefore can only be guaranteed to work for linear state transition and measurement models.

3.2 Implementation of the TDVB-E-GLMB Filter

Since the main contribution of this work is to implement the TDVB filtering technology and propose a new density, the trajectory density is the same as that in [18]. The system equations of the proposed method are (14)-(15). According to (14), the state prediction of TDVB approach is the Gaussian form, thus the state prediction step of the proposed TDVB-E-GLMB filter is the same as in [7, section V]. According to (15)-(17) and VB filtering technology, the density p^(h) can be written in the form of Student-T distribution, i.e.

\(\begin{aligned}\begin{array}{l}p_{k}^{(h)}(x, l \mid Z)=\sum_{i=1}^{J} \frac{w_{k \mid k, i}^{(h)}(l)}{\eta_{k \mid k}^{(h)}(l)} \operatorname{St}\left(m_{k \mid k, i}^{(h)}(l), R_{k}+P_{k \mid k, i}^{(h)}(l), v\right) \\ \approx \sum_{i=1}^{J} \frac{w_{k \mid k, i}^{(h)}(l)}{\eta_{k \mid k}^{(h)}(l)} \mathrm{N}\left(m_{k \mid k-1, i}^{(h)}(l), \frac{1}{\lambda_{k}} R_{k}+P_{k \mid k-1, i}^{(h)}(l)\right), \\ \times \operatorname{Gamma}\left(\lambda_{k} \mid \frac{v}{2}, \frac{v}{2}\right) \\\end{array}\end{aligned}\)       (27)

where,

\(\begin{aligned}m_{k \mid k, i}^{(h)}(l)=\left\{\begin{array}{ll}m_{k \mid k-1, i}^{(h)}(l)+K_{k} u_{k}^{(h)}, \text { if } \theta(l)>0 \\ m_{k \mid k-1, i}^{(h)}(l) & , \text { if } \theta(l)=0\end{array}\right.\end{aligned}\)       (28)

\(\begin{aligned}P_{k \mid k, i}^{(h)}(l)=\left\{\begin{array}{ll}P_{k \mid k-1, i}^{(h)}(l)+K_{k} S_{k}^{-1} K_{k}, & \text { if } \theta(l)>0 \\ P_{k \mid k-1, i}^{(h)}(l) & \text {, if } \theta(l)=0\end{array}\right.\end{aligned}\)       (29)

where, according to (22)-(24), the innovation parameters can be calculated by

u_k^(h) = z_k - H_km_k|k-1,i^(h)(l),       (30)

K_k = (P_k|k-1,i^(h)(l)⋅H_k^T)S_k^-1,       (31)

\(\begin{aligned}S_{k}=H_{k} P_{k \mid k-1, i}^{(h)}(l) H_{k}^{T}+\frac{1}{\bar{\lambda}_{k}} R_{k},\end{aligned}\)       (32)

where, H_k is the observation matrix. The auxiliary random variable \(\begin{aligned}\bar{\lambda}_{k}\end{aligned}\) can be obtained by

\(\begin{aligned}\bar{\lambda}_{k}=(v+d) /\left(v+\operatorname{tr}\left\{\left(u_{k}^{(h)}\left(u_{k}^{(h)}\right)^{T}\right) R_{k}^{-1}\right\}\right).\end{aligned}\)       (33)

Note that, formulas (28)-(32) need to be iterated to make \(\begin{aligned}\bar{\lambda}_{k}\end{aligned}\) converge. (P_k^(h)(x,l | Z) can be calculated by

\(\begin{aligned}w_{k \mid k, i}^{(h)}(l)=\left\{\begin{array}{ll}\frac{p_{D} w_{k \mid k-1, i}^{(h)}(l) q_{k, i}^{(h)}\left(z_{k}, \lambda_{k}, l\right)}{\kappa\left(z_{k}, l\right)} & \text { if } \theta(l)>0, \\ 1-p_{D} & \text {,if } \theta(l)=0\end{array}\right.\end{aligned}\)       (34)

where, p_D is the target missed detection probability. κ(z_k, l) , is the density of clutters. w_k|k-1,i^(h)(l) is the prediction weight and it is equal to w_k-1|k-1,i^(h)(l) (or birth weight when w_k-1|k-1,i^(h)(l) does not exist). q_k,i^(h)(z_k, λ_k, l) is the measurement likelihood. Since formulas (28)-(33) are TDVB form, q_k,i^(h)(z_k, λ_k, l) will be the same as that in [19], i.e.

q_k,i^(h)(z_k, λ_k, l) ≈ ϕ_k,i^(h)(z_k, λ_k, l)φ_k,i^(h)(λ_k,l).       (35)

where,

\(\begin{aligned}\begin{array}{l} \varphi_{k, i}^{(h)}\left(\lambda_{k}, l\right) \propto \exp \left(-\frac{1}{2} \bar{\lambda}_{k} \operatorname{tr}\left\{\left(u_{k}^{(h)}\left(u_{k}^{(h)}\right)^{T}\right) R_{k}^{-1}\right\}\right., \\ \left.+\left(\frac{v+d}{2}-1\right) \ln \bar{\lambda}_{k}-\frac{v \bar{\lambda}}{2}\right)\end{array}\\\end{aligned}\)       (36)

\(\begin{aligned}\begin{array}{l}\phi_{k, i}^{(h)}\left(z_{k}, \lambda_{k}, l\right) \propto \\ \exp \left(-\frac{1}{2} \lambda_{k}\left(z_{k}-H_{k} m_{k \mid k, i}^{(h)}(l)\right)^{T} R_{k}^{-1}\left(z_{k}-H_{k} m_{k \mid k, i}^{(h)}(l)\right). \right. \\ \left.\quad-\frac{1}{2}\left(\varepsilon_{k}^{(h)}\right)^{T}\left(P_{k \mid k-1, i}^{(h)}(l)\right)^{-1} \varepsilon_{k}^{(h)}\right)\end{array} \end{aligned}\)       (37)

The process of computing GLMB parameters

{(I₊^(h₊) , w₊^(h₊) , P₊^(h₊))}_{h+= 1}^H₊

is the same as the E-GLMB filter. Since it will take a long time to introduce the updating method of these trajectory parameters and this work has not contributed to this part, the specific formulas will not be introduced in detail in this work. The details can be found in [18, section 3 E].

In summary，this section improve the update steps of target kinematical state by formulas (28)-(33). Also, a new likelihood function is given by formulas (35)-(37). The code was published on GitHub [22] to help readers understand full details. The proof process of formula (36) and (37) can be found in the Appendix. In the next subsection, an effective multi-target TDVB smoother will be introduced.

4. The Multi-target TDVB Smoother

Smoothing the estimated state can effectively improve the performance of a tracker. The GLMB is an LRFS based tracking framework and can provide accurate trajectory information, thus it is feasible to use a smoothing method to improve the tracking performance. In this section, we proposed an MTT smoothing method based on TDVB smoothing technology.

4.1 The Key Method

The proposed smoothing method is based on the work in [19] and [20] (which can only be used in single target scenarios). The MTT scenario is different from a single target scenario, because it has missed detection, pseudo-track and data association problems. Therefore, we modify the single target smoother and apply it to the GLMB filtering framework.

The method in [19] and [20] is to smooth the estimated state after filtering. However, in MTT scenarios, the estimated state is obtained by the multiple hypotheses filtering process, which means that the estimated state is greatly affected by other targets, missed detection, and clutters. Fig. 2 shows an example in which the accuracy of multiple hypothesis filtering is reduced due to the closely-spaced targets. The gray arrow represents the movement direction. The forks and the dots represent the measurements produced by two targets at different scans, respectively. In the process of multiple hypothesis tracking, the measurement of two targets will affect the tracking accuracy. Further considering the influence of missed detection and clutters, as a result, the state estimation accuracy of filtering two targets will be lower than that of filtering a single target.

Fig. 2. Example of closely-spaced targets. Forks and the dots denote measurements of the two targets, respectively.

Therefore, in MTT scenarios, our view is not to smooth the filtered target state, but to make full use of the measurement information collected in the filtering and track association process. Our key method is: 1) storing the measurement set used to update the track during the tracking process, thus we obtain the measurement set of each target; 2) filtering the measurement set of each target, thus more accurate state estimation results can be obtained; 3) smoothing the state using TDVB smoothing technology. An example of a measurement set is all the forks in Fig. 2.

The details of the application of the proposed idea in a GLMB filter will be introduced in the next sub-section. Note that the proposed method is not only suitable for the GLMB filtering framework, but also to other filters that can provide target tracks.

4.2 The TDVB smoother for the GLMB filter

The standard GLMB filter can provide the labels and states, but not the corresponding measurement set of each target. Therefore, measurement storing steps should be included in the state update part. Formulas (23)-(27) are the steps of updating the state using a measurement z_k , thus a measurement memory \(\begin{aligned}\hat{Z}_{k}^{(h)}\end{aligned}\) of the h -th track can be included in these steps as \(\begin{aligned}\hat{Z}_{k}^{(h)} \leftarrow z_{k}\end{aligned}\) . When a GLMB component is deleted due to its small weight, the memory also needs to be deleted. As a result, at time k, all tracks output by the GLMB filter will have a corresponding measurement set \(\begin{aligned}\hat{Z}_{k}\end{aligned}\) . Note that \(\begin{aligned}\hat{Z}_{k}\end{aligned}\) is only a memory in the GLMB filtering process and does not affect other steps of the GLMB filter.

At time k, the measurement memory can be written as

\(\begin{aligned}\left\{\left\{Z_{k}^{(h)}\right\}_{k=D_{\min }^{(h)}}^{D_{\max }^{(h)}}\right\}_{h=1}^{N_{t r a j}}\end{aligned}\) ,

where, D_min^(h) and D_max^(h) are the beginning and ending time interval of the h track, respectively. N_traj is the number of all the tracks that have appeared. z_k^(h) denotes the corresponding measurement used to update the component states at time k. Fig. 3 shows an example where N_traj = 2. The memory stores the measurements corresponding to the two tracks. Note that if no measurement is used to update the component at some scan (i.e. it is a missed detection component), then z_k^(h) = ∅.

Fig. 3. Example of the measurement memory.

For each measurement set

\(\begin{aligned}\left\{Z_{k}^{(h)}\right\}_{k=D_{\min }^{(h)}}^{D_{\max }^{(h)}}\end{aligned}\) ,

the first step is to determine whether the track is generated from the clutters. Given a number N_pse , the trajectory should be deleted if

\(\begin{aligned}\left|\left\{z_{k}^{(h)}\right\}_{k=D_{\min }^{(h)}}^{D_{\max }^{(h)}}\right| \leq N_{p s e},\end{aligned}\)       (38)

otherwise, the new target states of the track can be calculated by

\(\begin{aligned}m_{k \mid k}^{\prime(h)}=\left\{\begin{array}{ll}m_{k \mid k-1}^{\prime(h)}+K_{k} u_{k}, & \text { if } z_{k}^{(h)} \neq \varnothing, \\ m_{k \mid k-1}^{\prime(h)} & , \text { if } z_{k}^{(h)}=\varnothing\end{array}\right.\end{aligned}\)       (39)

\(\begin{aligned}P_{k \mid k}^{\prime(h)}=\left\{\begin{array}{ll}P_{k \mid k-1}^{(h)}+K_{k} S_{k}^{-1} K_{k}^{T}, & \text { if } z_{k}^{(h)} \neq \varnothing, \\ P_{k \mid k-1}^{(h)} & , \text { if } z_{k}^{(h)}=\varnothing\end{array}\right.\end{aligned}\)       (40)

where, m'_k|k^(h) and P′_k|k^(h) are the new states calculated by the measurement set of the h -th track. The calculation approach of innovation parameters u_k^(h) , K_k , and S_k is the same as that in (25)-(27). The formulas (39)-(40) consider the missed detection. Therefore, formulas (38)-(40) consider the scenarios of the pseudo target and missed detection. Note that, the initialization of formulas (39)-(40) is determined by the state provided by the GLMB filter, i.e.

m'Dmin(h)(h) = Ξ({mk(h)}),       (41)

P'Dmin(h)(h) = P0(h),       (42)

where, Ξ(⋅) are defined as the first non-empty element in a set. In fact, (41) is to deal with the situation that mDmin(h)(h) provided by a GLMB filter is an empty set. P0(h) is the initial error matrix of the birth component of the GLMB filter.

The initialization of the proposed smoother is

\(\begin{aligned}m_{S, D_{\max }^{(h)}}^{(h)}=m_{D_{\max }^{(h)}}^{(h)}\end{aligned}\),       (43)

\(\begin{aligned}P_{S, D_{\max }^{(h)}}^{\prime(h)}=P_{D_{\max }^{(h)}}^{(h)}\end{aligned}\).       (44)

The smoothing process is

\(\begin{aligned}m_{S, k \mid k}^{\prime(h)}=m_{k \mid k}^{\prime(h)}+G\left(m_{S, k+1 \mid k}^{\prime(h)}-m_{k+1 \mid k}^{\prime(h)}\right)\end{aligned}\),       (45)

\(\begin{aligned}P_{S, k \mid k}^{\prime(h)}=P_{k \mid k}^{\prime(h)}-G\left(P_{k+1 \mid k}^{\prime(h)}-P_{S, k+1 \mid k}^{\prime(h)}\right) G^{T}\end{aligned}\),       (46)

\(\begin{aligned}G=P_{k \mid k}^{\prime(\ell)} F_{k}^{T}\left(P_{k+1 \mid k}^{\prime(h)}\right)^{-1}\end{aligned}\).       (47)

Given a number of iterations N_ite , the proposed TDVB smoother can be implemented to a GLMB filter. A smoother is not a real-time filter, thus a smoothing window should be set. In fact, the size of the window will affect the real-time performance and smoothing accuracy. Therefore, different smoothing windows should be selected according to different scenarios. The proposed method is based on the assumption that the window is large enough, thus the case of small window is not discussed. Research about this will be carried out in the future. Discussion on window size selection will be held in the next work. The pseudo-code is shown in Table 1 and the MATLAB code can be found in [22].

Table 1. Pseudo Code of TDVB Smoother

5. Simulation Result

In this section, the proposed methods are compared with the PHD filter [1], CPHD filter [5], and E-GLMB filter [18] with different parameters. The simulation data are generated according to [18], which simulates radar measurements. The target trajectories contain various complex cases are shown in Fig. 4. (a), and the simulated data of 100 scans (clutter rate is 50 and p_D = 0.75) is shown in Fig. 4. (b). The Scenario1 is a low interference scenario with low clutter rate, high detection probability, and no measurement-outlier. The Scenario2 is a complex scenario with a high cluter rate, low detection probability, and measurement-outlier. The section 5.3 shows the performance under different parameters.

Fig. 4. Target trajectories and simulated data.

The accuracy of each approach is judged by the OSPA metric [21], which does not consider the track association accuracy, thus it can be used to judge the performance of PHD and CPHD filters. Moreover, the OSPA(2) metric [30], which considers the track association accuracy, is used to more accurately judge the performance of GLMB filters. The simulation environment is as follows: Software: MATLAB 2019(a); Computer system: Windows 10; CPU: Intel(R) Core(TM) i7-6700HQ; RAM: 8GB.The specific tracking parameters are shown in TABLE II.

Table 2. The Parameters of Different Scenarios

In the table, u_x,0 and u_y,0 denote the coordinates of X and Y axes, respectively. The outlier noise ℜ_k with probability 10% means that: the measurement noise R_k of each target per scan has a 10% probability of being replaced by ℜ_k .

5.1 Results of Scenario1

In this subsection, the performance of each approach is judged by average result of 200 Monte Carlo runs.

Fig. 5. (a) shows the average OSPA results. The values of PHD and CPHD filters are obviously larger than GLMB filters, thus the tracking accuracy of a GLMB filter is higher than that of PHD and CPHD filters.

Fig. 5. Average performance of 200 Monte Carlo runs in Scenario1.

Fig. 5. (b) shows the average OSPA(2) results. The values of proposed TDVB-E-GLMB filter are significantly lower than that of E-GLMB filter when the data reaches the peaks. These data peaks are caused by the crossing and birth of targets, thus the proposed TDVB-E-GLMB filter is more robust in complex cases. However, the values of TDVB-E-GLMB are slightly higher than E-GLMB sometimes. It means that although the proposed method has strong robustness, the track association accuracy is slightly lower than the E-GLMB filter in ideal tracking environment. Moreover, the proposed smoother can significantly improve the accuracy of the two GLMB filters, and the smoothing results for TDVB-J-GLMB achieve the best performance.

Fig. 5. (c) shows the average estimation results of the number of targets. It can be seen that the GLMB trackers achieve high accuracy in this scenario, but the accuracy of PHD and CPHD filters is insufficient.

Fig. 5. (d) shows the average time costs. The values of the PHD and CPHD filters are significantly lower than GLMB filters. Therefore, the PHD and CPHD are low precision high speed filters. The values of proposed TDVB-E-GLMB are slightly lower than E-GLMB filter sometimes. In fact, the use of TDVB filtering technology will increase the updating computation time of each component. However, the improvement of the filtering accuracy leads to the decrease of the number of components, thus the overall time cost is slightly lower than that of the E-GLMB filter. Moreover, the proposed smoother only slightly increases the time cost of a GLMB filter (the corresponding values in the figure are obtained by filtering time + smoothing time). Therefore, the proposed smoother can greatly improve the filtering accuracy with minimal computational cost.

5.2 Results of Scenario2

Scenario2 is a complex scenario, thus it can judge the robustness of different approaches. The performances are judged by the average results of 200 Monte Carlo runs.

Fig. 6. (a) shows the average OSPA results. The GLMB filters are still better than PHD and CPHD filters in tracking accuracy.

Fig. 6. Average performance of 200 Monte Carlo runs in Scenario2.​​​​​​​

Fig. 6. (b) shows the average OSPA(2) results. The values of the TDVB-E-GLMB filter are obviously lower than that of E-GLMB, which means that the proposed method has strong robustness in complex scenarios. Moreover, the proposed smoother achieves the best performance. Therefore, the proposed smoother can handle complex scenarios as well.

Fig. 6. (c) shows the average estimation results of the number of targets. Compared with Fig. 5. (b), the complex tracking environment significantly affects the estimation accuracy of all approaches. However, the values estimated by the proposed TDVB-E-GLMB filter and smoother are closer to the real values.

Fig. 6. (d) shows the average time cost. As can be seen, the computational cost of proposed TDVB-E-GLMB filter and smoother are slightly higher than the standard E-GLMB filter.

5.3 Performance Under Different Parameters

In this scenario, we show the average OSPA(2) performance of 200 Monte Carlo runs with different parameters.

Fig. 7 shows the performance of different clutter rate and detection probability in the scenarios without measurement-outliers. It can be seen that the OSPA(2) values of the TDVB-E-GLMB filter is slightly lower than that of the E-GLMB filter. However, with the increase of the clutter rate (or the decrease of detection probability), the performance advantage of the proposed TDVB-E-GLMB filter becomes more obvious. It shows that the proposed algorithm has strong robustness under different parameters. Also, the proposed TDVB-E-GLMB filter with our smoothing approach achieves the best performance.

Fig. 7. Average OSPA(2) of 200 Monte Carlo runs with different parameters in the scenarios without measurement-outlier.​​​​​​​

Fig. 8 shows the average performance in the scenarios with measurement-outliers. Compared with Fig. 4, the performance advantage of the proposed TDVB-E-GLMB filter and the corresponding smoother is more obvious. Therefore, in the case of various parameters of different scenarios, the proposed approaches can effectively deal with measurement-outliers.

Fig. 8. Average OSPA(2) of 200 Monte Carlo runs with different parameters in the scenarios where the probability of measurement-outlier is 10%.​​​​​​​

6. Conclusion

In this paper, we propose a TDVB-E-GLMB filter and a corresponding smoother for tracking multiple targets. The first contribution is that the TDVB filtering technology is used in the E-GLMB filter to improve the accuracy and robustness. The second contribution is that we proposed a smoother for the GLMB filtering framework to improve the state estimation accuracy. The simulation results show that the performances of the proposed TDVB-E-GLMB filter and corresponding smoother are better than PHD, CPHD, and E-GLMB filters in various complex tracking scenarios. Note that, the proposed smoother uses the track information to obtain the measurement memory, and uses this memory for smoothing. Therefore, in theory, the proposed smoother can be applied to other filters which can provide tracks, rather than only to the GLMB filter.

In the future, we plan to apply the Student-T distribution to extended target tracking filters, such as [26].

Appendix