Multi-Sensor Multi-Object Tracking With the Generalized Labeled Multi-Bernoulli Filter

This paper proposes an efficient implementation of the multi-sensor generalized labeled multi-Bernoulli (GLMB) filter. Like its single-sensor counterpart, such implementation requires truncating the GLMB sum. However the single-sensor case requires solving 2-D ranked assignment problems whereas the multi-sensor case require solving multi-dimensional ranked assignment problems, which are NP-hard. The proposed implementation exploits the GLMB joint prediction and update together with a new technique for truncating the GLMB filtering density based on Gibbs sampling. The resulting algorithm has a quadratic complexity in the number of hypothesized objects and linear in the total number of measurements from all sensors.


I. INTRODUCTION
The objective of multi-object tracking is to jointly estimate the number of objects and their trajectories from sensor data [1]- [4]. A majority of multi-object tracking techniques are developed for single sensors.
The use of multiple sensors, in principle, reduces uncertainty about the object existence as well as its states. However, this problem is computationally intractable in general, especially for more than two sensors, even though conceptually the generalization to multiple sensors can be straightforward.
The random finite set (RFS) framework developed by Mahler [3], [4] has attracted significant attention as a general systematic treatment of multi-sensor multi-object systems. This framework facilitates the development of novel filters such as the Probability Hypothesis Density (PHD) filter [5], Cardinalized PHD (CPHD) filter [6], and multi-Bernoulli filters [3], [7], [8]. While these filters were not designed to estimate the trajectories of objects, they have been successfully deployed in many applications including radar/sonar [9], [10], computer vision [11]- [13], cell biology [14], autonomous vehicle [15]- [17] automotive safety [18], [19], sensor scheduling [20]- [26], and sensor network [27]- [29]. notation f X for the product x∈X f (x), with f ∅ = 1. We denote a generalization of the Kroneker delta that takes arbitrary arguments such as sets, vectors, integers etc., by For a given set S, 1 S (·) denotes the indicator function of S, and F(S) denotes the class of finite subsets of S. Also, for notational compactness, we drop the subscript k for the current time, the next time is indicated by the subscript '+'.

C. Standard multi-object observation model
For a given multi-object state X, each (x, ℓ) ∈ X is either detected by sensor s with probability P (s) D (x, ℓ) and generates a detection z (s) ∈ Z (s) with likelihood g Assuming that, conditional on X, detections are independent of each other and clutter, the multi-object likelihood function of sensor s is given by [32], [33] where: Θ (s) is the set of positive 1-1 maps θ (s) : L → {0:|Z (s) |}, i.e. maps such that no two distinct arguments are mapped to the same positive value, Θ (s) (I) is the subset of Θ (s) with domain I; and ψ (s,j) The map θ (s) specifies which objects generated which detections from sensor s, i.e. object ℓ generates detection z θ(ℓ) ∈ Z (s) , with undetected objects assigned to 0. The positive 1-1 property means that θ (s) is 1-1 on {ℓ : θ (s) (ℓ) > 0}, the set of labels that are assigned positive values, and ensures that any detection in Z (s) is assigned to at most one object.
Assuming that the sensors are conditionally independent, the multi-sensor likelihood is given by (1) ...

D. Generalised Label Multi-Bernoulli (GLMB)
A GLMB density can written in the following form 5 where each ξ ∈ Ξ represents a history of (multi-sensor) association maps ξ = (θ 1:k ), each p (ξ) (·, ℓ) is a probability density on X, and each ω (I,ξ) is non-negative with ξ∈Ξ I⊆L ω (I,ξ) = 1. The cardinality distribution of a GLMB is given by while, the existence probability and probability density of track ℓ ∈ L are respectively Given the GLMB density (9), an intuitive multi-object estimator is the multi-Bernoulli estimator, which first determines the set of labels L ⊆ L with existence probabilities above a prescribed threshold, and second the MAP/mean estimates from the densities p(·, ℓ), ℓ ∈ L, for the states of the objects. A popular estimator is a suboptimal version of the Marginal Multi-object Estimator [3], which first determines the pair (L, ξ) with the highest weight ω (L,ξ) such that |L| coincides with the MAP cardinality estimate, and second the MAP/mean estimates from p (ξ) (·, ℓ), ℓ ∈ L, for the states of the objects.
The number of components in the δ-GLMB filtering density grows exponentially with time, and needs to be truncated at every time step, ideally, by retaining those with largest weights since this minimizes the L 1 approximation error [33].

A. Truncation by Gibbs Sampling
1) The Target Distribution: This subsection formulates the target distribution π(I + , θ + |I, ξ) for the Gibbs sampler.  For each pair (I + , θ + ) ∈ F(L + ) × Θ + (I + ), we define the array The ith row of γ is denoted as γ i .
Assuming that for all i ∈ {1:P },P whereψ (ξ,j (1) ,...,j (S) ) Z+ and j (1) , ..., j (S) are the indices of the measurements assigned to label ℓ i , with j (s) = 0 indicating that ℓ i is misdetected by sensor s, and j (1) = ... = j (S) = −1 indicating that ℓ i no longer exists (note that if a row of γ has a negative entry then the entire row consists of negative entries an hence depends on the given (I, ξ, ) and Z + , which have been omitted for compactness. The assumptions on the expected survival and detection probabilities,P (ξ) , eliminates trivial and ideal sensing scenarios, as well as ensuring η i (j (1) , ..., j (S) ) > 0.
To establish (b), assume that γ is not positive 1-1. If γn is also not positive 1-1, i.e., 1 Γ(n) (γn) = 0, then the RHS of (26) trivially equates to 0. It remains to show that even if γn is positive 1-1, the RHS of (26) still equates to 0. Since γ is not positive 1-1, there exist an s and distinct i, j such that Proof: We are interested in highlighting the functional dependence of π n (γ n |γn) on γ n , while its dependence on all other variables is aggregated into the normalizing constant: Factorizing 1 Γ (γ) using Lemma 1, gives For Hence, sampling from the conditionals π n amounts to sampling from a categorical distribution with 1+ S s=1 (M (s) +1) categories. This proceedure has an O(P S s=1 M (s) ) complexity since sampling from a categorical distribution is linear in the number of categories [55]. The Gibbs sampler is summarized in Algorithm 1, and has a complexity of O(T P 2 S s=1 M (s) ).
end Proposition 3. Starting from any initial state in Γ, the Gibbs sampler defined by the family of conditionals (27) converges to the target distribution (25) at an exponential rate. More concisely, let π j denote the jth power of the transition kernel, then where β min γ,γ ′ ∈Γ π 2 (γ ′ |γ) > 0 is the least likely 2-step transition probability.
The proposed Gibbs sampler has a short burn-in period due to its exponential convergence rate. More importantly, since we are not using the samples to approximate (25) as in an MCMC inference problem, it is not necessary to discard burn-in and wait for samples from the stationary distribution. For the purpose of approximating the GLMB filtering density, each distinct sample constitutes one term in the approximant, and reduces the L 1 approximation error by an amount proportional to its weight. Hence, regardless of their distribution, all distinct samples can be used, the larger the weights, the smaller the L 1 error between the approximant and the true GLMB. Note that this is also called a block Gibbs sampler since for each row we are sampling from the joint distribution of the elements of the row.

B. Multi-Sensor Joint Prediction and Update Implementation
A GLMB of the form (9) is completely characterized by parameters (ω (I,ξ) , p (ξ) ), (I, ξ) ∈ F(L)×Ξ, which can be enumerated as Since the GLMB (9) can now be rewritten as and implementing the GLMB filter amounts to propagating forward the parameter set Note that to be consistent with the indexing by h instead of (I, ξ), we abbreviateP p (h,t) Since we are only interested in samples that provide a good representation of the GLMB filtering density, increased efficiency (for the same H max + ) can be achieved by using annealing or tempering techniques to modify the stationary distribution so as to induce the Gibbs sampler to seek more diverse samples [56], [57] (note that the actual weights of the GLMB components are computed using the correct parameters). One example is to increase the temperature for diversity. Tempering with the birth model (e.g. by feeding the Gibbs sampler with a larger birth rate) directly induces the chain to generate more components with births. Tempering with the survival probability induces the Gibbs sampler to generate more components with object deaths and improves track termination. Tempering with parameters such as detection probabilities and clutter rate induces the Gibbs sampler to generate components that reduce the occurrence of dropped tracks.

A. 3D Linear Gaussian Scenario
A 3D linear Gaussian scenario with 3 independent sensors is considered. An unknown and time varying number objects appear (up to 10 simultaneously in total) with births, deaths and crossings. Individual object kinematics are described by a 6D state vector of position and velocity (in the x, y, z directions respectively) that follows a constant velocity model with sampling period of 1s, and process noise standard deviation σ ν = 5m/s 2 . The survival probability is P S = 0.99, and the birth model is an LMB with Sensor 1 has good resolution on the x-axis only with respective noise standard deviations σ x = 10m, σ y = 100m, σ z = 100m on each axis. Sensor 2 has good resolution on the y-axis only with noise standard deviations σ x = 100m, σ y = 10m, σ z = 100m on each axis. Sensor 3 has good resolution on the z-axis only noise standard deviations σ x = 100m, σ y = 100m, σ z = 10m on each axis. All sensors have detection probability P D = 0.66 and uniform Poisson false alarms with an average rate of λ c = 20 per scan. The multi-sensor GLMB filter is implemented via the proposed Gibbs sampling technique. The filter is run with 10000 components and the Gibbs sampler is tempered by taking the 3rd root of the cost tensor corresponding to a total of 3 sensors. Figure 1 shows the ground truths in 3D space. Figures 2, 3

B. Alternative Target Distribution
Note that (20) is not the only choice of distribution that ensures valid components with high weights are chosen more often than those with low weights. Instead of sampling from π(I + , θ + |I, ξ) ∝ ω (I,ξ,I+,θ+) Z+ , this subsection introduces an alternative target distribution for the Gibbs sampler, which can drastically reduce the complexity. Unique samples from the alternative target distribution are then reweighted according to (30).
A more expensive alternative choice is η (s) n (j (s) |j (s−1) ) =q (ξ,ℓn,s−1,s) j (s−1) , j (s) . The additional computation comes from the calculation of normalizing constants (58) rather than (59). On the other hand, this choice of η (s) n (j (s) |j (s−1) ) yields a better approximation of η (s) n (j (1) , ..., j (S) ). However, since the weights for the components will be corrected after the Gibbs sampling step, the advantage of a better approximation is not substantial.

VI. CONCLUSIONS
This paper proposed an efficient implementation of the Multi-sensor GLMB filter by integrating the prediction and update into one step along with an efficient algorithm for truncating the GLMB filtering density based on Gibbs sampling. The resulting algorithm is an on-line multi-sensor multi-object tracker with linear complexity in the number of measurements of each sensor and quadratic in the number of hypothesized tracks. This implementation is also applicable to approximations such as the labeled multi-Bernoulli (LMB) filter since this filter requires a special case of the GLMB prediction and a full GLMB update to be performed [37].