Robust PMBM Filter with Unknown Detection Probability Based on Feature Estimation

This paper provides a solution for multi-target tracking with unknown detection probability. For the standard Poisson Multi-Bernoulli Mixture (PMBM) filter, the detection probability is generally considered a priori. However, affected by sensors, the features used for detection, and other environmental factors, the detection probability is time-varying and unknown in most multi-target tracking scenarios. Therefore, the standard PMBM filter is not feasible in practical scenarios. In order to overcome these practical restrictions, we improve the PMBM filter with unknown detection probability using the feature used for detection. Specifically, the feature is modeled as an inverse gamma distribution and the target kinematic state is modeled as a Gaussian distribution; the feature is integrated into the target kinematic state to iteratively estimate the target detection probability with the motion state. Our experimental results show that the proposed method outperforms the standard PMBM filter and the robust PMBM filter based on Beta distribution in the scenarios with unknown and time-varying detection probability. Further, we apply the proposed filter to a simulated infrared image to confirm the effectiveness and robustness of the filter.


Introduction
Multi-target tracking (MTT) refers to the process of jointly estimating the state and number of targets based on the observation data obtained by sensors [1]. MMT is a central problem both in military fields and civilian fields such as defense, surveillance, biomedical research, and autonomous driving [1][2][3][4][5]. Two main types of algorithms have been developed over the past few decades: traditional MTT algorithms and MTT algorithms based on random finite set (RFS) [6,7]. The former, which simplifies the MMT to single-target tracking by solving data association problem between targets and observations, mainly includes Multiple Hypotheses Tracking (MHT) [8,9], Joint Probabilistic Data Association (JPDA) [10], etc., while the latter, which has received much attention and extensive research, due to its providing an effective solution to the data association problem, mainly consists of a Probability Hypothesis Density (PHD) [11] filter, Cardinality-PHD (CPHD) [12] filter, Multi-target multi-Bernoulli (MeMBer) [6] filter, Generalized Labeled Multi-Bernoulli (GLMB) [13][14][15] filter, and Poisson Multi-Bernoulli Mixture (PMBM) [16] filter. Among these algorithms, the PMBM filter, which is a convolution of Poisson RFS and multi-Bernoulli mixture (MBM) RFS, is considered a promising tracking algorithm thanks to its high performance and reduced computational cost [2,17].
It is worth mentioning that the standard PMBM filter requires some a priori knowledge, such as the detection probability, and it is usually considered as a fixed value. However, the detection probability is unknown and even time-varying in practical applications; a significant mismatch in detection probability can result in a significant bias or erroneous

PMBM RFS Density
Let x k,1 , . . . , x k,N(k) and z k,1 , . . . , z k,M(k) denote the states of all N(k) targets and all M(k) measurements at time k, respectively. Then we have RFSs X k = x k,1 , . . . , x k,N(k) ⊂ X , Z k = z k,1 , . . . , z k,M(k) ⊂ Z which denote the multi-target state set and multi-target observation set, respectively. X denotes the state space and Z denotes the observation space.
According to the observation, the given non-empty state set X can be divided into two parts: undetected targets subsets X u and detected targets subsets X d . Hence the posterior density can be denoted as follows: f P (X u ) = e − µ(x)dx ∏ x∈X u µ(x), (2) f mbm X d = ∑ j∈I ∑ i∈I j X i =X d ∏ i∈I j ω j,i f j,i X i , where f P (X u ) is the Poisson density, denoting the undetected targets, µ(x) is the intensity function. f mbm X d is a multi-Bernoulli mixture denoting the potential targets which are detected at least once. ω j,i is the hypothesis weight and ∑ ω j,i = 1, f j,i X i denotes the ith Bernoulli density in the jth global hypothesis, which is given by with r j,i denoting the existence probability and p j,i (x) denoting the single target density.

PMBM Filter
For the standard PMBM filter, the recursive processes are summarized as follows [16].
where γ k (x) is the intensity of the new born targets, µ k|k−1 (x) denotes the predicted intensity. f k|k−1 (x|ζ) and P S,k (ζ) are the transition density of state x given ζ and survival probability given ζ, respectively.
where P D denotes the detection probability.
(b) For MBM Component: The update step of detected targets can be divided into two types.
Update for the targets detected for the first time: where e k (z) = P D g k (z|ζ)µ k|k−1 (ζ)dζ (13)  where g k (z|x) is the likelihood function, and c(z) is the clutter intensity. Update for the targets detected previously: After the update, we get all possible new single-target hypotheses, we have to go through all possible data association hypotheses to construct the global hypotheses. In order to reduce the cost, a Gibbs sampler [29] or Murty's [14,16] algorithm can be employed to prune the number of the hypotheses to improve the computation efficiency. In this work, we use the Gibbs sampler due to the lower computation complexity [30]. The detail implementation can be found in [29,31].

Gamma Distribution and Inverse Gamma Distribution
The probability density of the inverse Gamma distribution for non-negative variable a can be denoted as [25] IG(a; s; t) = t s Γ(s) where shape parameter s > 0 and scale parameter t > 0. Γ(s) is Gamma function. The mode at which the probability density function is the maximum is t/(s + 1). The mean value and variance of the IG distribution is t/(s − 1) and t 2 / (s − 1) 2 (s − 2) respectively. The probability density of the Gamma distribution for variable a is as follows where shape parameter s > 0 and scale parameter t > 0. The mode and the mean value are (s − 1)/t and s/t respectively.

The Proposed Robust Filter with Unknown Detection Probability
In this section, the specific implementation of the robust PMBM filter based on inverse gamma Gaussian mixture (IGGM) distribution is introduced.

Target State Model and Observation Model
Similar to [23,25], the feature denoted by a is augmented to x which denotes the kinematic state of a single target and consists of positions and velocities; letx expresses the new state of a single target, i.e.,x = (x, a). The variable a denotes the SNR throughout in this paper. The detection probability and survival probability at time k can be expressed as The kinematic state and the feature are modeled as Gaussian distribution and inverse gamma distribution, respectively, so the target density at time k can be denoted by the IGGM as where J k is the number of the IGGM components at time k, and ω i k is the weight of the ith IGGM component. m and P are the mean and covariance of Gaussian density. The Markov transition density can be expressed as F k−1 and Q k−1 denote the state transition matrix and process noise covariance. Similarly, the observation state is also augmented; the new observation state is expressed asẑ = (z, h) where z and h are the observation of position and the feature, respectively. The likelihood function of the augmented state at time k can be expressed H k and R k represent the observation matrix and observation noise covariance, respectively. The likelihood of feature is gamma distribution which ensures the conjugation of the PMBM filter.
The update of the IG component can be calculated as follows [25].
IG a; s k|k−1 ; t k|k−1 g k (h|a)

The Implementtation of Proposed Algorithm
Similar to the standard PMBM filter, the proposed algorithm can be divided into Poisson components and MBM components which denote undetected targets and potentially detected targets, respectively. Besides, the iterative recursion of the Gaussian component is similar to Kalman filter (KF).
Prediction step: Given the intensity of Poisson RFS µ k−1 (x, a) and MBM RFS where s i,u k|k−1 and t i,u k|k−1 can be obtained according to (25), (26).
where m s i,u k|k−1 and t i,u k|k−1 can be obtained according to (25) and ( (b) For MBM Component: The update step of detected targets can be divided into two types Update for the targets detected for the first time where where κ k (z, h) is the clutter intensity. Update for the targets detected previously The feature can be extracted by Then, the detection probability can be obtained by Algorithm 1 gives the pseudo-code of the proposed algorithm. Algorithm 1 Description of the proposed robust filter

The Computation Complexity of the Proposed Algorithm
Before discussing the computation complexity of the proposed algorithm (IGGM-PMBM), we first analyze the computation complexity of Gaussian Mixture-PMBM (GM-PMBM) [32].
Suppose that the number of association hypotheses is A k−1|k−1 after the prediction step at time k − 1, it takes N u k|k−1 m k + A k−1|k−1 N k|k−1 m k + A k−1|k−1 (m k /N k|k−1 ) N k|k−1 steps to calculate the updated PMBM density at time k, where N u k|k−1 , N k|k−1 and m k are the number of the unknown target components, the MB's components, and the measurements after the prediction step, respectively. Thus, the complexity of the GM-PMBM update is . The detailed calculations can be found in [33]. From the results, it can be seen that the computation complexity of the PMBM filter is related to the number of the unknown target components, the MB's components, and the measurements. The framework of IGGM-PMBM filter is unchanged relative to the GM-PMBM filter. Thus, the complexity of the IGGM-PMBM filter does not increase relative to the GM-PMBM filter in theory. However, due to the need to propagate an additional function, i.e., inverse gamma distribution, the IGGM-PMBM filter has a slightly higher complexity compared to the GM-PMBM filter. The same analysis applies to BGM-PMBM filter. Besides, this IGGM-PMBM filter does not increase the complexity of the BGM-PMBM filter but improves the performance.

Simulation Scenario 1
As a verification for the proposed algorithm, the simulation data in [16]  and each observation is a vector of position z = z x ; z y . The parameters used in (23) and (27) are given as where ⊗ is the Kronecker product, T = 1 is the sampling period, and q = 0.01. In formula (25), which denotes the prediction step of IG distribution, k s is set as 0.9. The parameter ξ in likelihood function g k (h|a) for feature d is 10. In formula (61), we set δ 1 = 4, δ 2 = 2, and ε 1 = (2 − exp(−SNR th /δ 1 )) −1 , ε 2 = exp(−SNR th /δ 1 ). In the simulation, we let the survival probability be 0.99, i.e., p S,k = 0.99. The parameters used in the pruning and merging processes are the same as in [16,25]. The proposed algorithm was compared with the GM-PMBM filter and the BGM-PMBM filter which can also estimate the detection probability online. The generalized optimal sub-pattern (GOSPA) [34] assignment metric with parameters α = 2 is employed to assess the performance of filters. The GOSPA is defined as in likelihood function (ℎ| ) for feature d is 10. In formula (61), we set 1 = 4, 2 = 2,and 1 = (2 − exp (− ℎ / 1 )) −1 , 2 = exp (− ℎ / 1 ). In the simulation, we let the survival probability be 0.99, i.e., , = 0.99. The parameters used in the pruning and merging processes are the same as in [16,25]. The proposed algorithm was compared with the GM-PMBM filter and the BGM-PMBM filter which can also estimate the detection probability online. The generalized optimal sub-pattern (GOSPA) [34] assignment metric with parameters = 2 is employed to assess the performance of filters. The GOSPA is defined as  There are two cases to be considered: fixed detection probability and time-varying probability.
Case 1-Fixed detection probability: In this scenario, the target detection probability is fixed. Targets are born according to a Poisson process with intensity:  There are two cases to be considered: fixed detection probability and time-varying probability.
Case 1-Fixed detection probability: In this scenario, the target detection probability is fixed. Targets are born according to a Poisson process with intensity: where ω γ = 0.03 and Gaussian density with mean [100;0;100;0], covariance diag 150 2 , 1, 150 2 , 1 . s i γ = 51 and t i γ = 500 denotes the parameters of IG distribution [25], thus the feature used for detection is 10. Clutter is also a Poisson process with intensity κ k (z, h) = λ κ IG a; s κ k ; t κ k g(h|a), where the clutter rate is λ κ = 10. The parameters of IG distribution are s κ k = 31 and t κ k = 280. The parameters are summarized in Table 1. The parameters used in the BGM distribution are the same as in [23]. In this simulation, we set the SNR th as two different values: 9 and 5.5, thus the corresponding detection probability is 0.68 and 0.94 respectively. Figures 2 and 3 show the average results, corresponding to the performance metrics on their GOSPA error and the number of targets. In Figure 2, the result shows that the multitarget tracking performance of the proposed IGGM-PMBM filter is similar to the standard GM-PMBM filter with the similar GOSPA distance, where the detection probability is exactly known, fixed and high. Whereas in the low probability scenario, the BGM-PMBM filter and GM-PMBM filter yield slightly higher GOSPA error than IGGM-PMBM, especially the BGM-PMBM. This simulation shows that the proposed IGGM-PMBM filter outperforms the BMG-PMBM filter and GM-PMBM filter with low detection probability. The standard deviation range (StDev) values of estimated number of targets of IGGM-PMBM and BGM-PMBM are shown in Figures 4 and 5. On the other hand, it can be seen in Figure 6 that the estimation of P D of IGGM-PMBM is more accurate. Figure 7 gives the comparisons of LE, ME, and FE when the detection probability is 0.68. The performance comparisons in terms of average GOSPA, LE, ME, and FE with different parameters (detection probabilities and clutter rates) are shown in Table 2. The results show that the performance of the proposed IGGM-PMBM filter is better than the BGM-PMBM filter and the GM-PMBM filter under the same parameters.  Table 1. The parameters used in the BGM distribution are the same as in [23]. In this simulation, we set the ℎ as two different values: 9 and 5.5, thus the corresponding detection probability is 0.68 and 0.94 respectively. Figures 2 and 3 show the average results, corresponding to the performance metrics on their GOSPA error and the number of targets. In Figure 2, the result shows that the multitarget tracking performance of the proposed IGGM-PMBM filter is similar to the standard GM-PMBM filter with the similar GOSPA distance, where the detection probability is exactly known, fixed and high. Whereas in the low probability scenario, the BGM-PMBM filter and GM-PMBM filter yield slightly higher GOSPA error than IGGM-PMBM, especially the BGM-PMBM. This simulation shows that the proposed IGGM-PMBM filter outperforms the BMG-PMBM filter and GM-PMBM filter with low detection probability. The standard deviation range (StDev) values of estimated number of targets of IGGM-PMBM and BGM-PMBM are shown in Figures 4 and 5. On the other hand, it can be seen in Figure 6 that the estimation of of IGGM-PMBM is more accurate. Figure 7 gives the comparisons of LE, ME, and FE when the detection probability is 0.68. The performance comparisons in terms of average GOSPA, LE, ME, and FE with different parameters (detection probabilities and clutter rates) are shown in Table 2. The results show that the performance of the proposed IGGM-PMBM filter is better than the BGM-PMBM filter and the GM-PMBM filter under the same parameters.           The three filters are run separately on an AMD Core 3.20 GHz CPU PC with 16 GB RAM and MATLAB R2021b. The complexity is also illustrated by comparing the computational time. Based on 100 Monte Carlo runs, the average computational times of the three filters are shown in Table 3. It can be seen that the computation complexity of IGGM-PMBM is slightly higher than the GM-PMBM filter in a high detection probability scenario, and with the detection probability decreasing, the IGGM-PMBM filter costs more time to tackle the unknown target detection probability situation. Besides, the IGGM-PMBM filter has almost the same complexity as the BGM-PMBM filter. Case 2-Changing detection probability: In this scenario, the detection probability is varying with time. We set different parameters to express the different feature values for detection as Table 4. The SNR th in this case is 5.5, thus the detection probabilities are 0.94, 0.81, and 0.69 respectively. Other parameters are the same as in Case 1. The average GOSPA error, LE, FE, ME, and cardinality estimate as well as estimate of P D are shown in Figure 8. It can be seen that the estimates of P D of the IGGM-PMBM filter is more accurate than BGM-PMBM filter for each segment. Moreover, the GOSPA errors of the proposed IGGM-PMBM is smaller than that of BGM-PMBM. Figure 9 gives the StDev of the estimated number of targets of the IGGM-PMBM filter and BGM-PMBM filter under the detection probability vary.

Simulation Scenario 2
We simulate a long wave infrared image of the movement of space targets using Satellite Tool Kit (STK). The image is 1024 × 1024 pixels, and there are 5400 frames. Figure 10 is the real trajectory of each target in the infrared image. We only track the centroid of the target in this experiment. There are a total of 9 targets, of which target 1 is always in the center of the image. Figure 11 shows the birth state of each target. As shown in Table 5, the birth state of each target including the position and velocity is set according to the target initial states. The covariance is diag ([1, 1, 1, 1]).The initial value of the feature for each target is different, and the parameters s γ and t γ of IG distribution are different at the birth time. The number of frames in which the target is born and the number of frames in which it dies are also given in Table 5. In this scenario, SNR th is 5.5. Other parameters are the same as in simulation scenario1.

Simulation Scenario 2
We simulate a long wave infrared image of the movement of space targets using Satellite Tool Kit (STK). The image is 1024 × 1024 pixels, and there are 5400 frames. Figure  10 is the real trajectory of each target in the infrared image. We only track the centroid of the target in this experiment. There are a total of 9 targets, of which target 1 is always in the center of the image. Figure 11 shows the birth state of each target. As shown in Table  5, the birth state of each target including the position and velocity is set according to the target initial states. The covariance is ([1,1,1,1]).The initial value of the feature for each target is different, and the parameters and of IG distribution are different at the birth time. The number of frames in which the target is born and the number of frames in which it dies are also given in Table 5. In this scenario, ℎ is 5.5. Other parameters are the same as in simulation scenario1. Figure 12 gives the performance of the IGGM-PMBM filter used in the infrared image. Obviously, the filter is effective for the infrared image.    Figure 12 gives the performance of the IGGM-PMBM filter used in the infrared image. Obviously, the filter is effective for the infrared image. (g) (h) (i) Figure 11. The birth of each target (a) First target appears, (b) Second target appears, (c) Third target appears, (d) Fourth target appears, (e) Fifth target appears, (f) Sixth target appears, (g) Seventh target appears, (h) Eighth target appears, (i) Ninth target appears. Figure 11. The birth of each target (a) First target appears, (b) Second target appears, (c) Third target appears, (d) Fourth target appears, (e) Fifth target appears, (f) Sixth target appears, (g) Seventh target appears, (h) Eighth target appears, (i) Ninth target appears.

Conclusions
In order to solve the problem of unknown detection probability and improve the accuracy when the low detection probability is low, we propose a novel PMBM filter based on the feature estimation for multi-target tracking. We model the feature as an inverse

Conclusions
In order to solve the problem of unknown detection probability and improve the accuracy when the low detection probability is low, we propose a novel PMBM filter based on the feature estimation for multi-target tracking. We model the feature as an inverse gamma distribution and integrate it into the kinematic state of targets. The feature can be estimated along with the motion state of the target, so as to further calculate the detection probability. The feature is related to detection and can be actually measured which can improve the accuracy of the detection probability estimation. The simulation results also show that the proposed algorithm can accurately estimate the detection probability and outperforms the GM-PMBM filter and BGM-PMBM filter with slightly lower GOSPA error especially when the detection probability is low. The proposed algorithm was then used in an infrared simulation scenario, which demonstrated the robustness and the effectiveness of the proposed algorithm. With the increase in space targets and development of autonomous driving and other fields, more and more radar data and space-based infrared image data will be obtained. Thus, the verification of the proposed algorithm using real data such as radar data and infrared images data would be a worthwhile subject of future study. Extending the proposed algorithm with unknown information is also a topic worthy of study.
Author Contributions: Y.W. performed the simulation and wrote the paper; X.C. simulated the infrared image and offered some useful suggestions with regard to methodology with P.R. All authors have read and agreed to the published version of the manuscript.