Joint Probabilistic Data Association Filter with Unknown Detection Probability and Clutter Rate

This paper proposes a novel joint probabilistic data association (JPDA) filter for joint target tracking and track maintenance under unknown detection probability and clutter rate. The proposed algorithm consists of two main parts: (1) the standard JPDA filter with a Poisson point process birth model for multi-object state estimation; and (2) a multi-Bernoulli filter for detection probability and clutter rate estimation. The performance of the proposed JPDA filter is evaluated through empirical tests. The results of the empirical tests show that the proposed JPDA filter has comparable performance with ideal JPDA that is assumed to have perfect knowledge of detection probability and clutter rate. Therefore, the algorithm developed is practical and could be implemented in a wide range of applications.


Introduction
There has been increasing attention on utilisation of small autonomous systems in military and civil applications. The issue is that the operations of these small autonomous systems are constrained by limited payload, as well as limited operation time and endurance. This has led to proliferation of lightweight, low-cost and energy efficient on-board sensors.
Reliable and autonomous target tracking is a fundamental aspect of situation awareness for autonomous systems [1][2][3]. Applying low-cost and lightweight on-board sensors in target tracking imposes additional challenges to the tracking problem since they are likely to contain some degree of uncertainties. Low-cost sensors are generally subject to a high clutter rate and low detection probability. When combined with the inherent uncertainties and complexity of the problem, the poor performance issue with these sensors could be significantly exacerbated in target tracking, especially in multi-target tracking (MTT) [1].
It is known that data association in MTT is a challenging issue. With a high clutter rate, data association in MTT becomes even more challenging. One of the most popular data association approaches is the multiple hypothesis tracking (MHT) proposed in [4,5]. Given high computational power, MHT is known as a powerful MTT algorithm to address the problem of measurement uncertainty. Data, in MHT, are processed recursively with a delayed logic scheme and MHT intentionally expands parent track hypotheses with received measurements to form child track hypotheses in decision-making trees. These child track hypotheses are subsequently assessed and the ones with low scores are pruned by thresholding to keep feasibility. Although MHT is known as a theoretically optimal Bayesian estimator, the exact solution of MHT is computationally intractable and thereby requires approximated implementations [6,7].
JPDA is another widely used one-scan data association method. The basic assumption of JPDA is that each measurement can originate from several candidate targets in the valid gate. Therefore, JPDA

System Models and Assumptions
This section addresses the system models that are used in the following sections. Throughout the paper, we use the symbol i to refer to a target index and j to refer to a measurement index.
The set of target states and measurements received at scan k are, respectively, defined as where N k denotes the number of targets at scan k, x i k the ith target at scan k, M k the number of measurements received at scan k, z j k (j = 0) the jth measurement received at scan k, z 0 k the dummy measurement for convenient representation of miss detection. We assume that the temporal evolution of each target is independent of others and follows a Markov transition model f (x k |x k−1 ), then, the prediction of target state of scan k is governed by As required for JPDA, we accept the assumption that each target can generate at most one measurement and each measurement can originate from at most one target. Each target-generated measurement is independent of each other and is detected with probability P D with measurement likelihood p(z |x ).
The additional assumptions made in the paper are as follows: • The clutter distribution is assumed to be unknown a priori. When there is no prior knowledge regarding the clutter measurements, the MTT problem generally assumes that the number of clutters or false alarms is locally Poisson distributed [23]. Therefore, this paper utilises a Poisson distribution as the clutter distribution. Clutters or false alarms are modelled by a local PPP with intensity λ FA (z) = N FA /V s with N FA being the average number of clutters of one scan and V s being the sensor volume. Clearly, if domain or environmental knowledge could be used to develop a more realistic clutter distribution in specific tracking tasks, the estimation performance could be improved.

•
The detection probability is assumed to be independent of the target state in filter design.
The detailed discussion will be presented in Section 4.1.
In standard JPDA filtering approach, a track is defined as a sequence of measurements that originate from the same target. The original JPDA filter assumes the number of targets is known a priori at each scan and utilise a wrapper heuristic M/N logic for track initiation and deletion. In this paper, we resort to a PPP model, similar to [24], for target birth instead of the heuristic M/N logic. In other words, we assume that target birth is modelled as a PPP with intensity λ B (x) and track is confirmed based on a random binary variable, target existence status, e i k ∈ {0, 1} with e i k = 1 being the existence of the ith target and e i k = 0 being non-existence. Then, measurements that either from new targets or false alarms can be modelled by a PPP with intensity λ E = λ B (x) p (z |x ) P D + λ FA (z). Note that, by incorporating the PPP model into JPDA filter, the posterior multi-object PDF can then be represented by the multi-Bernoulli distribution. This fact enables the application of recently proposed multi-Bernoulli RFS to JPDA for joint detection probability and clutter rate estimation. Let r i k = p e i k = 1 denote the existence probability of the ith target at scan k. Then, the time evolution of r i k can be formulated by the Markov Chain One model [9] and the Markov transition probability matrix is determined by the survival probability P S [9].

Joint Probabilistic Data Association Filter and Analysis
In this section, we first briefly review the basics of JPDA filter for the completeness and then give an intuitive analysis of the effect of the detection probability and clutter rate.

Joint Probabilistic Data Association Filter
JPDA algorithm aims to calculate the marginalized association probability based on all possible joint events for data association. A joint event is an allocation of all measurements to all tracks. In JPDA, a feasible joint event is defined as one possible mapping of the measurements to the tracks such that: (1) each measurement (except for the dummy one) is assigned to at most one target; and (2) each target is uniquely assigned to a measurement. Let Θ k = θ i k , i ∈ 1, 2, . . . , N k|k−1 , denote the joint association event. For each pre-existed target i ∈ 1, 2, . . . , N k|k−1 , θ i k ∈ {0, 1, . . . , M k } denotes the association event, where θ i k = j means the jth measurement is originated from the ith target and θ i k = 0 represents the dummy association in which the ith target is miss detected. JPDA assumes that each single association event is independent and the posterior of each target is characterised by a mixture as [5] Generally, propagation of the mixture distribution is computationally intractable due to the explosion of mixture terms. JPDA approximates this mixture term by a single Gaussian distribution through first moment matching method. More specifically, the state correction x i k|k for each pre-existed target i ∈ 1, 2, . . . , N k|k−1 and its corresponding covariance P i k|k of JPDA is obtained as , Z k is the existence conditioned marginal association probability that the jth measurement is associated with the ith target.
The hypothesis-conditioned update p x i k θ i k , e i k = 1, Z k can be calculated by standard Kalman filter algorithm. Therefore, the remaining part is how to obtain the marginal association probability p θ i k e i k = 1, Z k in a tractable way as this is the most computationally part. The posterior distribution of the joint association event p (Θ k |Z k ) consists of two parts: miss detection f m and detection f d [9]. Therefore, one can imply that Note the term 1 − r i k|k−1 in f m is for the non-existence of the ith target and r i k|k−1 (1 − P D ) is for the existence but non-detection of the ith target.
Although full enumeration is intractable in real applications due to the high demand of computational power, the marginal probability p θ i k |Z k can be approximated by m−best approximations, stochastic sampling or any other suitable approaches [10][11][12]25]. In this paper, we utilize the Gibbs sampling approach to approximate the marginal probability. This method enables fast calculation of the marginal probability with ignorable performance sacrifice. Detailed description of this algorithm can be found in our previous work [25]. For the completeness of the paper, a brief introduction to the Gibbs-JPDA implementation algorithm is provided in the Appendix.
After finding p θ i k |Z k , the joint probability p θ i k , e i k = 1 |Z k can be calculated using Bayesian rule as Based on Equations (5)- (7), the hypothesis-conditioned existence probability has the form as Finally, the posterior probability of target existence r i k and the existence-conditioned marginal probability p θ i k e i k = 1, Z k used in the estimation update can be calculated as Using r i k|k and p θ i k e i k = 1, Z k , we can now perform track update for the JPDA filter using Equations (3) and (4).
Note that, except for the updates of N k|k−1 pre-existed tracks, JPDA also creates M k new tracks for each measurement based on the birth model. After the track update, we use a similar approach, as shown in score-based approach [5], for track confirmation and deletion. A track is confirmed once its existence probability r i k exceeds a pre-defined threshold and otherwise is still tentative and needs further test to be confirmed or deleted. Meanwhile, a track is immediately deleted if its existence probability r i k is below a pre-defined threshold.

Remark 1.
Note that the utilisation of PPP birth model for target existence estimation is similar to the idea of joint integrated probabilistic data association (JIPDA) [9,26]. However, the difference lies in that the JIPDA filter assumes a stationary or constant intensity of extraneous sources, which include unknown targets and false alarms while the proposed JPDA with PPP birth model dynamically estimate this intensity. More specifically, in JIPDA derivation, Equation (7) reduces to [9] Therefore, the standard JIPDA ignores the influence of the unknown targets in marginalisation. Clearly, leveraging the PPP birth model is beneficial and more realistic in initialisation and calculating the marginal association probability.

Effect of the Detection Probability and Clutter Rate
This subsection analyses the effect of detection probability and clutter rate on the performance of JPDA filter provided that they are given in advance.
The effect of the detection probability P D on the estimation are two folds. On the one hand, it follows from Equations (6) and (7) that the the detection probability P D influences the relative ratio or weighting between the target miss detection f m and target detection f d . Therefore, artificially lower detection probability will give more penalty on the target-absence f m , leading to the delayed track maintenance. However, lower detection probability might be useful in holding targets through temporary fades. On the other hand, the fact that λ E = λ B (x) p (z |x ) P D + λ FA (z) also reveals that the detection probability also affects the track initialization process. Lower detection probability will force longer tentative tracks before a confirmation or deletion decision is made based on the score logic. This, in turn, will increase the workload of the tracking filter as more tracks remain tentative, which is especially severe in a high-clutter scenario.
From Equation (7), one can note that the clutter rate λ FA also plays an important role in the JPDA algorithm. Higher clutter rate will give more penalty on the external source term, i.e., false alarms or new targets, λ FA z j k + λ B z j k p z j k x i k P D than the target-detection term p χ i k |Z k p z j k x i k P D . Therefore, artificially higher clutter rate provides an estimation result that more tracks will be considered as false alarms by the tracking system. Similarly, lower clutter rate usually generates overestimation of the real targets.
Consequently, significant mismatch of either detection probability P D or clutter rate λ FA will result in highly biased estimation. However, both P D and λ FA are difficult or not cost-effective to obtain in practice, especially for for low-cost sensors. Therefore, it is imperative to design an improved JPDA algorithm that can accommodate this issue.

Joint Probabilistic Data Association Filter with Unknown Detection Probability and Clutter Rate
In this section, we will propose a novel JPDA filter to adjust the unknown detection probability and clutter rate by leveraging the multi-Bernoulli filter.

Detection Probability and Clutter Rate Estimator
The recently proposed multi-Bernoulli filter [18,[27][28][29] utilises the random finite set (RFS) theory to recursively propagate the posterior of multi-target state. This filter is known as a parametrised approximation of the Bayesian multi-target recursion. At each time instant k, the posterior PDF of the multi-target state is characterised as a multi-Bernoulli RFS. More specifically, each target is modelled as a single Bernoulli RFS, which is fully described by the existence probability r i k and the probability density p k (x i k ). A Bernoulli RFS is empty with probability 1 − r i k and has only one element, whose distribution is described by p(x i ), with probability r i k . The multi-Bernoulli RFS is a simple union of N k single independent Bernoulli RFSs and thus can be completely characterised by the parameter . The recursive multi-Bernoulli filter is given by two steps as [27,28]: (1) Prediction step. If the posterior multi-target density at time k − 1 is characterised as a multi-Bernoulli RFS as , then, the prediction of the multi-target density is also a multi-Bernoulli as where the first part denotes the multi-Bernoulli RFS of the birth model and the second part represents the prediction of existing targets with where α, β (2) Update step. If the predicted multi-target density at time k is characterised as a multi-Bernoulli , then, the update of the multi-target density is also a multi-Bernoulli as where the legacy (miss detection) component is given by and the measurement updated component is given by The advantage of the multi-Bernoulli filter is that it utilises the mathematically sound RFS theory and thus directly avoids the data association process in multi-target tracking, leading to computational efficiency. However, the target identity or the target label cannot be maintained during the filtering recursion. This means that the multi-Bernoulli filter cannot be used to support the target label management. Therefore, if the label information is required for practice, the application of multi-Bernoulli filter is limited. Fortunately, as our aim is to estimate the clutter rate as well as the detection probability using the multi-Bernoulli filter, there is no need to manage the clutter identity.
To apply the multi-Bernoulli filter to estimate the clutter rate, we consider each false alarm as a target with its own transition model, birth model and death model [18]. This consideration is motivated by the fact that the target and the clutter have totally different dynamics and thus can be treated separately in estimation. The true targets and clutters in recursive filtering are therefore labelled by a discrete space I = {0, 1}, where 1 is the label of real targets and 0 is for clutters. To further accommodate the detection probability estimation, we augment the state space asX = P × X × I, where P = [0, 1] denotes the space of detection probability and × represents the Cartesian product. Since the new state space is a discrete one, the integration on this space is given by In reality, the detection probability is related to the target state. This fact makes the system state propagation/prediction of the augmented system intractable. To address this problem, we further assume that the detection probability is independent of the target state, i.e., (22) Note that the target detection probability used in JPDA is typically piecewise constant and thus the independence assumption is reasonable in filter design. Furthermore, the validation presented in the following section reveals that the proposed algorithm works well in the presence of time-varying detection probability.
In general, the clutter dynamics is highly dependent on the sensors. If we have some prior information regarding the sensors, we can utilise such information to build the clutter dynamics model. Otherwise, general models can be leveraged, e.g., Gaussian distribution and random walk. Similarly, if there is no prior information, a typical choice of the detection probability model is the well-known beta distribution P D ∼ Beta (α, β) [30]. The beta distribution Beta (α, β) is a family of continuous distributions defined on the interval [0, 1] and has enough diversity to accommodate the changes of detection probability. By choosing different shaping parameters α and β, the beta distribution can accommodate different detection probability models.
By substituting the augmented state space into the original multi-Bernoulli filter, we can readily estimate the number of clutters at one scanN FA by the summation of all confirmed targets with label 0 and the clutter rate estimationλ FA is then given byλ FA =N FA /V s . As P D is part of the system states, it can directly obtained by the state estimation of multi-Bernoulli filter. Note that one can also apply the augmented system into JPDA to adapt to the unknown detection probability and clutter rate. However, the inherent data association process of JPDA will make practical implementation intractable.

Estimation Algorithm
The proposed JPDA is shown in Figure 1. At each time instant k, we utilise the original JPDA with PPP birth model as the main tracker to provide the multi-target state estimation. The required detection probability and clutter rate is obtained from an one-step multi-Bernoulli filter. As discussed earlier, by incorporating the PPP model into JPDA, the posterior of JPDA estimation can be modelled as the multi-Bernoulli RFS and thus the estimated output of the JPDA is feedback to the multi-Bernoulli filter for its initialisation at every time instant. The advantage of the proposed parallel filtering scheme is that one can fully exploit the strengths of both approaches: using the JPDA for track management and obtaining fast estimation of detection probability as well as clutter rate by the multi-Bernoulli filter.

Remark 2.
The computational complexity of the proposed estimation algorithm comes from two main parts: JPDA and multi-Bernoulli filter. Although the original marginalisation in JPDA is a #P-complete problem, leading to the fact that the computational time increases exponentially with respect to the number of targets, the approximation by Gibbs-sampling is a polynomial-time algorithm, as detailed in [25]. Additionally, since the RFS-based multi-Bernoulli filter directly avoids the typical data association procedure in multi-target tracking, it also runs in a polynomial time.

Numerical Simulations
In this section, the effectiveness of the proposed JPDA filtering algorithm is investigated through numerical simulations in a cluttered environment. Our experiments explore a scenario, involving 10 manoeuvring targets with different birth time. The ground truth of the considered scenario is depicted in Figure 2.

Simulation Setup
The state vector contains planar position and velocity. We use the well-known constant velocity (CV) model for target prediction. The CV model is defined as with where I 2×2 denotes the 2 × 2 identity matrix, T the sampling period, and w k ∼ N ·; 0, σ 2 v the Gaussian process noise.
It is assumed that we use a radar for multiple targets tracking. The target-generated nonlinear range and bearing measurements are modelled bỹ where (x T,k , y T,k ) denotes target position, (x R , y R ) radar position, and v k ∼ N (·; 0, R k ) the Gaussian measurement noise with R k = diag σ 2 r , σ 2 a . The clutter is assumed to be uniformly distributed in the surveillance region with its number being Poisson with N FA average returns at each scan. To validate the proposed algorithm in dynamic scenarios, N FA is given in a time-varying profile, as shown in Figure 3. Gating is performed with a threshold such that the gating probability is P G = 0.999. The detection probability P D peaks at the sensor origin with P D = 0.95 and exponentially decreases to 0.8 at the boundary of the sensing field-of-view. This setting is reasonable as the received signal strength decreases when the target is far away from the sensor. A tentative track is confirmed if the existence probability satisfies r i k ≥ 0.8 and a confirmed track is deleted immediately once r i k ≤ 0.1. The optimal sub-pattern assignment (OSPA) distance metric [31] is considered here for overall performance evaluation, namely, cardinality and position estimation errors. Let X and Y be the position estimation set and true target position set, respectively. The cardinality of these two sets are m and n, respectively. Then, for c > 0 and 1 ≤ p < ∞, the OSPA distance d c p (X, Y) is defined as [31], where Π n denotes the set of all permutations on {1, 2, . . . , n} for any positive integer n. d c x i , y π(i) is the cut-off Euclidean distance between two vectors as d c x i , y π(i) = min d x i , y π(i) , c with d x i , y π(i) being the Euclidean distance. The order parameter p determines the sensitivity of d c p (X, Y) in penalizing estimation outliers, while the cut-off parameter c determines the relative weighting of the penalties allocated to cardinality and localization errors. In all simulations, these two parameters are set as p = 1, c = 200.

Simulation Results
For the purpose of comparison, we also perform the ideal JPDA filter with perfect knowledge of detection probability and clutter rate, the following four non-ideal JPDA filters: (1) ideal P D ,  Table 1. The peaks of mean OSPA distance in Figure 4 are resulted from track confirmation for target birth. If the clutter rate is set lower than the real value, the JPDA filter overestimate the multi-target state and thus generating more 'confirmed' targets, leading to larger OSPA distance (e.g., non-ideal JPDA filter 1). If the clutter rate is set higher than the real value, the JPDA filter gives an underestimation the multi-target state (e.g., non-ideal JPDA filter 2). Although the results, shown in Figure 4, reveal that the JPDA filter is somewhat robust against the detection probability (e.g., non-ideal JPDA filter 3), lower detection probability gives a delayed confirmation, as we stated earlier. Moreover, when the detection probability is highly mismatched with its real value (e.g., non-ideal JPDA filter 4), the estimation is unreliable. As a comparison, the results clearly verify that the proposed JPDA filter has comparable performance as the ideal JPDA filter by utilising the multi-Bernoulli estimator in a feedback loop.

Conclusions
This paper developed an enhanced version of JPDA by incorporating with the multi-Bernoulli filter to accommodate the unknown detection probability and clutter rate. By utilising the PPP birth in JPDA filter, the multi-target estimation output can be characterised as a multi-Bernoulli RFS. This enables the application of the multi-Bernoulli filter in a feedback loop to estimate the unknown detection probability and clutter rate. Simulation results confirm that the proposed algorithm exhibits the advantages of the original JPDA filter and can correct the biased estimations induced by the unknown detection probability and clutter. Under the condition where detection probability and clutter rate are difficult to be obtained, the proposed JPDA filter can be a practical solution, providing performance comparable to the ideal JPDA filter.

Conflicts of Interest:
The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Appendix A. Approximated Marginalisation by Gibbs Sampling
According to Equation (5), the posterior distribution of the joint event can be formulated as where ϕ θ i k denotes the un-normalised PDA probability of θ i given by and which ensure that one measurement (except for the dummy measurement) can only be allocated to one target. It is clear that determining the marginal joint association probability requires full enumeration of all possible joint association events, which is a well-known #P-complete problem. To tackle the combinatorial nature of this problem, our previous work suggested a sampling based, especially Gibbs sampling-based, algorithm [25]. The key idea of utilising Gibbs sampling in JPDA is to consider each joint association event Θ k as a random variable that satisfies a distribution π (Θ k ). To ensure that the joint events with higher probability are more easily to be sampled, the sampling proposal π (Θ k ) is constructed proportional to the joint posterior probability as The issue is that direct sampling from Equation (A5) is very difficult as enumerating all possible joint events is impossible for real-time applications. Therefore, Gibbs sampling to applied for generating the samples. The transition kernel from one joint event Θ k = θ 1 k , . . . , θ N k|k−1 k to another joint eventΘ k = θ 1 k , . . . ,θ