Hybrid PHD-PMB Trajectory Smoothing Using Backward Simulation

Multi-object filtering concerns the estimation of the number of objects and their current states based on a sequence of noisy measurements [1, 2]. This problem is compounded by missed detections, clutter, and data association uncertainties. To address the multi-object filtering problem in a computationally efficient way, the probability hypothesis density (PHD) filter recursively propagates the first-order moment approximation of the multi-object posterior density, which is also known as the PHD [3]. By doing so, the PHD filter operates on the single-object space and, consequently, avoids the data association problem arising from multiple objects. Due to its computational efficiency, the PHD filter has been widely used in many applications [4, 5, 6, 7, 8, 9].

The estimation of multi-object states by also using measurements at later time steps is called multi-object smoothing. The PHD smoother involves first a forward PHD filtering recursion and then a backward smoothing recursion [10, 11]; the latter is performed by propagating the PHD of the multi-object smoothing density backward in time. Particle filter implementations of the PHD smoother were provided in [10, 11, 12, 13, 14] and the closed-form Gaussian implementations were presented in [15, 16]. However, it has been reported in [11, 12, 16] that, compared to the PHD filter, the PHD smoother reduces the localization error, but it does not necessarily improve the cardinality estimate. In particular, the PHD smoother suffers from premature object deaths [12] by reporting missed detections at earlier time steps when objects disappear, with a lag equal to the smoothing lag.

A common drawback of the PHD filter and smoother is that they do not maintain track continuity, and thus, cannot provide object trajectory estimates. A simple track building procedure for the PHD filter is by adding labels/tags to the object states or PHD components, and track continuity can be established by linking the object states or PHD components with the same label/tag over time [17, 18]. These methods may work well in some scenarios but can also result in track switches, and even missed and false detections, see [19].

A mathematically rigorous way to build trajectories for the PHD filter without using any labels/tags is by considering sets of trajectories [20]. The resulting filter is called the trajectory PHD filter, which directly provides object trajectory estimates by recursively propagating a Poisson multi-object density on the space of sets of trajectories [19]. Furthermore, smoothed trajectories can be obtained in forward filtering by considering smoothing-while-filtering. However, the trajectory PHD filter is mainly used for estimating trajectories of alive objects, and it cannot properly estimate the set of all trajectories, which also includes trajectories of dead objects.

The possibility of generating sets of smoothed trajectories for both alive and dead objects using a sequence of (unlabeled) multi-object filtering densities has been explored in [21]. In particular, this work showed that we can sample sets of trajectories using backward simulation [22] even when the forward multi-object filter is an unlabeled filter that does not explicitly maintain track continuity. The presented multi-trajectory smoothing solution has been applied to the Poisson multi-Bernoulli (PMB) filters [23, 24], and the resulting multi-trajectory smoother has achieved state-of-the-art multi-trajectory estimation performance. However, there are no similar methods for the PHD filters¹¹1Backward simulation was first applied to the particle-based PHD filter in [13] to generate smoothed trajectories with labels, but the algorithm developed in [13] can only track up to two objects..

In this paper, we present a hybrid PHD-PMB trajectory smoother, which involves PHD forward filtering and PMB backward smoothing. In PHD filtering, the predicted multi-object density is of the form Poisson Point Process (PPP) without approximation, but the updated multi-object density is of the form PMB [23]. The Bayesian filtering recursion is achieved by finding the Poisson approximation of the PMB density that matches its PHD. The key difference between our proposed hybrid PHD-PMB trajectory smoother and all the existing PHD smoothers [10, 11, 12, 14, 15, 16, 13] is that we perform backward smoothing on the PMB multi-object densities in the PHD filtering recursion before the Poisson approximation. Considering that the PMB representation of the multi-object posterior is an inherent intermediate result that we obtain for free in the PHD update step (as we will show in Section IV), it makes sense to perform smoothing using this more informative multi-object posterior representation.

The hybrid PHD-PMB trajectory smoother leverages the simplicity of PHD filtering [4] and the ability to generate smoothed trajectories using backward simulation [25, 21]. Importantly, the hybrid PHD-PMB trajectory smoother makes it possible to generate smoothed trajectories for all the objects that have existed in a theoretically sound manner without using any labels/tags. We compare the proposed hybrid PHD-PMB trajectory smoother to the PHD [4] and trajectory PHD [19] filters in a simulation study, and the results demonstrate that the proposed method has superior multi-object state and trajectory estimation performance. The results also show that the hybrid PHD-PMB trajectory smoother is more robust to the premature object death problem [11, 12].

The rest of the paper is organized as follows. In Section II, we introduce the background. The problem formulation is given in Section III. We present the filtering and smoothing recursions of the hybrid PHD-PMB trajectory smoother in Section IV. The results are presented in Section V, and the conclusions are drawn in Section VI.

In this section, we introduce the state variables of interest, densities on sets of objects and sets of trajectories, and the multi-object models. The set of finite subsets of a generic space $D$ is denoted by ${\cal F}(D)$, and the cardinality of a set $A\in{\cal F}(D)$ is $|A|$. The sequence of ordered positive integers $(\alpha,\alpha+1,\dots,\gamma-1,\gamma)$ is $\alpha:\gamma$. We use $\uplus$ to denote a union of sets that are mutually disjoint and $\langle f,g\rangle$ to denote the inner product $\int f(x)g(x)dx$. In addition, we use $\delta_{x}(\cdot)$ and $\delta_{x}[\cdot]$ to represent the Dirac and Kronecker delta functions, respectively, centered at $x$.

In backward smoothing for sets of trajectories, we aim to compute the multi-trajectory posterior density in a given time interval using a sequence of multi-object filtering densities and the multi-trajectory dynamic model.

We denote the multi-object state density at time step $k^{\prime}\in\{k,k+1\}$ and the multi-trajectory density in time interval $\alpha:\gamma$, both conditioned on the sequence of sets of measurements ${\bf z}_{1:k}=({\bf z}_{1},\dots,{\bf z}_{k})$ up to and including time step $k$, by $f_{k^{\prime}|k}(\cdot)$ and $\pi_{\alpha:\gamma|k}(\cdot)$, respectively. In forward filtering, the Bayes filter propagates the multi-object posterior density of $\mathbf{x}_{k}$ in time via the prediction and update steps:

where $g(\cdot|{\bf x})$ is the multi-object transition density and $\ell({\bf z}_{k}|\cdot)$ is the measurement likelihood. Given the multi-object models described in Section II-D, the predicted and posterior densities on the current set of object states are PMB mixtures (PMBMs) [23]. A PMB is a special case of PMBM with a single MB.

The multi-trajectory density $\pi_{1:K|K}(\mathbf{X}_{1:K})$ can be obtained by applying the following multi-trajectory smoothing equation recursively backwards in time [21, Theorem 2]:

where the predicted multi-trajectory density for the multi-trajectory dynamic model in Section II-E is given by [21, Theorem 1]

It can be observed that the predicted multi-trajectory density can be evaluated by multiplying the multi-trajectory density $\pi_{k:k|k}({\bf X}_{k:k})$, which takes the same value as $f_{k|k}(\tau^{k}({\bf X}_{k:k}))$, $1-p^{S}(\cdot)$ for trajectories that end at time step $k$, $g(\cdot|\cdot)p^{S}(\cdot)$ for trajectories that exist at both time step $k$ and $k+1$, and the density of newborn trajectories at time step $k+1$.

The backward smoothing equation (16), in general, cannot be computed in closed form. An effective solution to evaluate the multi-trajectory posterior is given by backward simulation, which allows us to draw samples in the multi-object trajectory space [21]. Performing backward simulation requires a backward kernel to recursively draw samples of sets of trajectories.

The backward kernel for sets of trajectories conditioned on the set ${\bf Y}$ of trajectories in the time interval $k+1:K$ and the sequence of measurement sets up to and including time step $K$, satisfies [21, Lemma 3]

That is, the backward kernel $\pi_{k:K|K}({\bf X}|{\bf Y})$ is proportional to the predicted multi-trajectory density $\pi_{k:k+1|k}({\bf X}_{k:k+1})$ if ${\bf Y}={\bf X}_{k+1:K}$, and zero otherwise. Evaluating the backward kernel $\pi_{k:K|K}({\bf X}|{\bf Y})$ considers all different ways of associating the trajectories in ${\bf X}_{k:k+1}$ to the trajectories in ${\bf Y}$.

In this paper, the forward filtering is achieved using a PHD filter [3], and the PMB multi-object filtering densities $f_{k|k}(\mathbf{x})$ before the Poisson approximation are stored at each time step. Then backward simulation is performed on these PMB filtering densities to obtain the approximate multi-trajectory posterior.

In this section, we present the forward filtering and backward smoothing recursions of the hybrid PHD-PMB trajectory smoother. An illustration of the complete recursion is shown in Fig. 1. We also discuss some practical considerations for a tractable linear Gaussian implementation.

In this section, we evaluate the performance of the proposed hybrid PHD-PMB trajectory smoother in a simulation study in comparison with the PHD filter [4] and the trajectory PHD filter [19].

We consider the same two-dimensional scenario as in [19] with 100 time steps and 4 objects. The single object state is $x=[p_{x},\dot{p_{x}},p_{y},\dot{p_{y}}]^{T}$, consisting of position and velocity. The single object dynamic model is nearly constant velocity with transition matrix $F$ and process noise covariance $Q$ given by

where $I_{2}$ is a $2\times 2$ identity matrix, $\otimes$ denotes the Kronecker product, $T_{s}=0.5\,\text{s}$ is the sampling period, and $\sigma_{q}=1.8$. Each object survives with probability $p^{S}=0.99$. The birth model is a PPP, and its Poisson intensity is a Gaussian mixture with three components. Each component has the same weight $0.1$ and covariance $\text{diag}([225,100,225,100])$; the means are $[85,0,140,0]^{T}$, $[-5,0,220,0]^{T}$, and $[7,0,50,0]^{T}$. The ground truth object trajectories are illustrated in Fig. 2.

The single object measurement model is linear and Gaussian with observation matrix $H$ and measurement noise covariance $R$ given by

where $\sigma_{r}=2$. The detection probability is $p^{D}=0.9$. The clutter is uniformly distributed in the area $[0\,\text{m},2000\,\text{m}]\times[0\,\text{m},2000\,\text{m}]$ with Poisson clutter rate $\gamma^{C}=50$.

The PHD filter [4] is implemented with pruning threshold $10^{-4}$, merging threshold $4$, and maximum number of Gaussian components $30$. Object state estimates are obtained using the estimator described in [26, Sec. 9.5.4.4]. The trajectory PHD filter [19] is implemented also with pruning threshold $10^{-4}$, absorption threshold $4$, and maximum number of Gaussian components $30$, and it performs smoothing-while-filtering without $L$-scan approximation. Its trajectory estimates are obtained using the estimator described in [19, Sec. VII-D]. The hybrid PHD-PMB trajectory smoother is implemented with $1000$ particles, and only a maximum of 100 global hypotheses are sampled. Its trajectory estimates are obtained using the estimator described in [21, App. G]. An exemplar output of the hybrid PHD-PMB trajectory smoother is shown in Fig. 3.

The multi-object state estimation performance is evaluated using the generalized optimal subpattern assignment (GOSPA) metric [30] at each time step with $\alpha=2$, $c=10$, and $p=1$. GOSPA allows for the decomposition of the estimation error into localization error, missed and false detection errors. The multi-trajectory smoothing estimation performance is evaluated for the set of all trajectories using the trajectory GOSPA (TGOSPA) metric [31] at the last time step with $c=10$, $p=1$, and $\gamma=1$. TGOSPA allows for the decomposition of the estimation error into localization error, missed and false detection errors, and track switch error. We note that the trajectory PHD filter cannot estimate the set of all trajectories accurately. To report the set of all trajectory estimates, we tag each Gaussian component and extract the trajectory estimates with unique tags via post-processing. We report the results averaged over 100 Monte Carlo runs.

Let us first study the multi-object estimation performance. The GOSPA error and its decomposition for the the PHD filter and the hybrid PHD-PMB trajectory smoother over time are shown in Fig. 4. Overall, the hybrid PHD-PMB trajectory smoother outperforms the PHD filter in terms of localization errors, missed and false detection errors by a significant margin. One exception is that the hybrid PHD-PMB trajectory smoother presents higher localization errors at the first few time steps due to the mismatch between the Gaussian means in the birth intensity and the true object birth positions. Moreover, the PHD-PMB trajectory smoother shows one apparent anomaly that it has increased missed detection error a few time steps before object disappearing. This anomaly is due to the premature object death, a problem also observed in PHD smoothers [11, 12], where missed detections happen at earlier time steps when objects disappear, with a lag equal to the smoothing lag. However, we note that in the hybrid PHD-PMB trajectory smoother this problem is less severe in the sense that it does not always happen for every disappearing object, and that missed detections due to premature object death are not propagated all the way back to earlier time steps. We also show the average GOSPA error and its decomposition per time step for different scenario parameters in Table I. We can see that the hybrid PHD-PMB trajectory smoother consistently outperforms the PHD filter by a significant margin, and it is especially robust to false detections.

We then discuss the trajectory estimation performance of the hybrid PHD-PMB trajectory smoother and the trajectory PHD filter. The TGOSPA error and its decomposition for the trajectory PHD filter and the hybrid PHD-PMB trajectory smoother over time are shown in Fig. 5. The hybrid PHD-PMB trajectory smoother surpasses the trajectory PHD filter by producing fewer false detections. However, the trajectory PHD filter exhibits slightly superior performance in localization accuracy, reduced missed detections, and fewer track switches. The average TGOSPA error and its decomposition per time step for different scenario parameters are presented in Table II. We can see that the hybrid PHD-PMB trajectory smoother generally outperforms the trajectory PHD filter, except for the scenario with low clutter rate.

In this paper, we have presented a hybrid PHD-PMB multi-object smoother, which performs backward simulation using the PMB densities obtained in the update step of forward PHD filtering to provide smoothed trajectory estimates. The hybrid PHD-PMB trajectory smoother makes it possible for the PHD filter to generate smoothed trajectory estimates for all the objects. Possible future work includes the development of a hybrid PHD-PMB trajectory smoother for non-Gaussian, non-linear scenarios using sequential Monte Carlo implementations and experiments with real data [32].