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.

1. 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 estimation, which will affect the practicability of the algorithm. In order to make the filter more suitable for practical scenarios, the robust CPHD/PHD filter [18,19], robust multi-Bernoulli filter [20] and the Cardinality-Balance MeMBer filter [21], the GLMB filter with unknown detection probability [22], and the robust PMBM filter [23] have been successively proposed. All of them model the detection probability as a Beta distribution and iteratively estimate the variable with the targets’ state. However, the performance of these filters degrades when the detection probability is low; the performance also suffers when the model initialization parameters are poor.

Target detection probability depends on the sensor, environment and the features used for detection [24,25]. Generally, the signal-to-noise ratio (SNR) and the amplitude are the most widely used features for detection [26,27]. The SNR feature is introduced to PHD/CPHD to estimate the detection probability [25], while [26] introduces the amplitude to the PHD/CPHD filter to improve filters’ performance, and [28] introduces the amplitude to multi-sensor MeMBer/CPHD to estimate the detection probability; simulation results show that the performance of the filters is improved thanks to the above.

Inspired by [25,26,28], we improve the PMBM filter with unknown detection probability using the feature used for detection. Firstly, we model the feature as an inverse gamma distribution. Then, the feature is integrated into the target kinematic state which is modeled as a Gaussian distribution. In particular, the feature is independent with the kinematic state of targets due to the feature being a relatively stable variable which will not change significantly with the small change in target location. Finally, the detection probability and the motion state can be estimated iteratively. Two outcomes of the experiment can be summarized as (a) the tracking performance of the proposed filter is similar to the standard PMBM filter, where the detection probability is exactly known, fixed and high. Whereas in the low probability scenario, the Beta Gaussian Mixture (BGM-PMBM) filter and standard PMBM filter yield slightly higher error than the proposed filter, especially the BGM-PMBM filter. (b) when the detection probability is varying with time, the proposed filter can more accurately estimate the varying detection probability and has better performance. The results presented show the robustness and effectiveness of the proposed filter.

The key contribution of this paper is the derivation of a close-form solution to the robust PMBM filter recursion which can jointly estimate the target state and detection probability. The intent is to exploit the inverse gamma (IG) component to estimate the feature used for detection and thus the detection probability. The algorithm in this paper improves the tracking performance of the target under unknown detection probability, especially when the detection probability is relatively low.

The outline of the remainder of this paper is as follows. Section 2 introduces the necessary background knowledge of the proposed algorithm, including PMBM recursion and inverse gamma (IG) distribution. Section 3 derives the proposed filter and its implementation. Section 4 gives the simulation results with a linear multi-target filtering scenario. A conclusion is provided in Section 5.

2. Background

2.1. PMBM RFS Density

Let and denote the states of all targets and all measurements at time , respectively. Then we have RFSs , which denote the multi-target state set and multi-target observation set, respectively. denotes the state space and denotes the observation space.

According to the observation, the given non-empty state set can be divided into two parts: undetected targets subsets and detected targets subsets . Hence the posterior density can be denoted as follows: where is the Poisson density, denoting the undetected targets, is the intensity function. is a multi-Bernoulli mixture denoting the potential targets which are detected at least once. is the hypothesis weight and , denotes the th Bernoulli density in the th global hypothesis, which is given by with denoting the existence probability and denoting the single target density.

2.2. PMBM Filter

For the standard PMBM filter, the recursive processes are summarized as follows [16].

Prediction step: Suppose the intensity of Poisson RFS and the MBM RFS are given at time . The prediction step can be expressed by

(a) For Poisson Component where is the intensity of the new born targets, denotes the predicted intensity. and are the transition density of state given and survival probability given , respectively.

(b) For MBM Component

Update step: Given the predicted PMBM filter with parameters and at time , the update step can be obtained based on the observations .

(a) For Poisson Component where 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 where is the likelihood function, and 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].

2.3. Gamma Distribution and Inverse Gamma Distribution

The probability density of the inverse Gamma distribution for non-negative variable can be denoted as [25] where shape parameter s > 0 and scale parameter . is Gamma function. The mode at which the probability density function is the maximum is . The mean value and variance of the IG distribution is and respectively.

The probability density of the Gamma distribution for variable is as follows where shape parameter and scale parameter . The mode and the mean value are and respectively.

3. 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.

3.1. Target State Model and Observation Model

Similar to [23,25], the feature denoted by is augmented to which denotes the kinematic state of a single target and consists of positions and velocities; let expresses the new state of a single target, i.e., . The variable denotes the SNR throughout in this paper. The detection probability and survival probability at time 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 can be denoted by the IGGM as where is the number of the IGGM components at time , and is the weight of the th IGGM component. and are the mean and covariance of Gaussian density. The Markov transition density can be expressed as

and denote the state transition matrix and process noise covariance.

Similarly, the observation state is also augmented; the new observation state is expressed as where and are the observation of position and the feature, respectively. The likelihood function of the augmented state at time can be expressed

and 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].

3.2. 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 and MBM RFS at time . The prediction step can be expressed by

(a) For Poisson Component where and can be obtained according to (25), (26). is also the IGGM form.

(b) For MBM Component where and can be obtained according to (25) and (26).

Update step: Suppose the predicted intensity of Poisson RFS and MBM RFS are given at time , the update step can be expressed by

(a) For Poisson Component

(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 is the clutter intensity.

Update for the targets detected previously where , ,,, , .

The feature can be extracted by

Then, the detection probability can be obtained by

Algorithm 1 gives the pseudo-code of the proposed algorithm.

3.3. 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 after the prediction step at time , it takes steps to calculate the updated PMBM density at time , where , and 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.

4. Simulation Setup and Results

4.1. Simulation Scenario 1

As a verification for the proposed algorithm, the simulation data in [16] was selected. The area is [110, 180] m [110, 180] m. There are four targets; Figure 1 is the targets’ real track. The targets are all born at Step 1 and alive all the time, except the blue one which dies at Step 40. Targets kinematic state include the position and velocity and each observation is a vector of position . The parameters used in (23) and (27) are given as where is the Kronecker product, = 1 is the sampling period, and . In formula (25), which denotes the prediction step of IG distribution, is set as 0.9. The parameter in likelihood function for feature d is 10. In formula (61), we set , , and , . In the simulation, we let the survival probability be 0.99, i.e., . 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 is employed to assess the performance of filters. The GOSPA is defined as

Figure 1

Multi-target scenario. There are four targets, all born at Step1. The blue target dies at Step 40 and the others are alive the entire time.

The denotes the assignment set between and . In this situation, , and denote the location error (LE), miss error (ME), and false error (FE), respectively. Besides , p.

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 and Gaussian density with mean [100;0;100;0], covariance . and denotes the parameters of IG distribution [25], thus the feature used for detection is 10. Clutter is also a Poisson process with intensity , where the clutter rate is . The parameters of IG distribution are and . The parameters are summarized in Table 1. The parameters used in the BGM distribution are the same as in [23].

Table 1

Parameters of the fixed detection probability.

parameter	k_s	ξ	p_S,k	sγi	tγi	skκ	tkκ	λ_κ
value	0.9	10	0.99	51	500	31	280	10

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. Figure 2 and Figure 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 Figure 4 and Figure 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.

Figure 2

Average of GOSPA distances in 100 Monte Carlo runs. (a) , ; (b) , .

Figure 3

The estimated average number of targets in 100 Monte Carlo runs. (a) , ; (b) , .

Figure 4

The standard deviation range value and the mean value of the estimated number of targets of IGGM-PMBM filter (a) , ; (b) , .

Figure 5

The standard deviation range value and the mean value of the estimated number of targets of BGM-PMBM filter (a) , ; (b) , .

Figure 6

The average of estimated detection probability in 100 Monte Carlo runs. (a) , ; (b) , .

Figure 7

Comparisons of (a) location error, (b) miss error and (c) false error with , .

Table 2

Root mean square GOSPA and LE, ME, and FE at all time steps with different parameters (P, λ).

(PD, λκ)	Proposed IGGM-PMBM				GM-PMBM				BGM-PMBM
	GOSPA	LE	ME	FE	GOSPA	LE	ME	FE	GOSPA	LE	ME	FE
(0.94, 10)	2.7477	1.7745	1.9277	0.8278	2.7511	1.7780	1.9277	0.8315	2.7592	1.7767	1.9181	0.8819
(0.94, 15)	2.8097	1.7897	1.9798	0.8784	2.7940	1.7876	1.9689	0.8571	2.7963	1.7895	1.9563	0.8889
(0.94, 20)	2.8417	1.7886	2.0154	0.9027	2.8534	1.7985	1.9985	0.9558	2.8538	1.7990	1.9783	0.9969
(0.94, 25)	2.8730	1.7937	2.0215	0.9750	2.9012	1.8063	2.0846	0.8992	2.9098	1.8098	2.0638	0.9655
(0.68, 10)	3.7648	2.2141	2.7499	1.3076	3.8460	2.2045	2.9187	1.1889	3.9832	2.2184	3.0123	1.3676
(0.68, 15)	4.0033	2.2626	3.0358	1.3005	3.9756	2.2422	3.0358	1.2497	4.1448	2.2610	3.1456	1.3494
(0.68, 20)	4.1351	2.2292	3.1407	1.5051	4.2651	2.1971	3.3939	1.3586	4.3428	2.2309	3.4265	1.4635
(0.68, 25)	4.2925	2.2756	3.3370	1.4530	4.3054	2.2758	3.4021	1.3240	4.4715	2.2840	3.5651	1.4380

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.

Table 3

Average computational times.

P_D	Proposed IGGM-PMBM	BGM-PMBM	GM-PMBM
0.94	3.83 s	3.71 s	3.37 s
0.68	6.07 s	5.89 s	5.16 s

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 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 are shown in Figure 8. It can be seen that the estimates of 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.

Table 4

The parameters for different detection probability.

Step	0~20	21~60	61~80
(s_γ,t_γ)	(51, 500)	(51, 385)	(51, 335)

Figure 8

Comparisons of (a) GOSPA error, (b) Location error, (c) False error, (d) Miss error, (e) the number of targets and (f) the estimated .

Figure 9

The standard deviation range value and the mean value of the estimated number of targets in the case of detection probability varying. (a) IGGM-PMBM filter, (b) BGM-PMBM filter.

4.2. 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 .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 10

Real trajectory of each target in the infrared image.

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.

Table 5

The target initial states.

Target	State	Feature	s_γ	t_γ	Survival Time (Frame)
1	[512;0;512;0]	6	51	300	[18,5400]
2	[512;0;512;0]	6	51	300	[1501,1991]
3	[512;0;512;0]	8	51	400	[1801,2139]
4	[512;0;512;0]	8.6	51	430	[2101,2329]
5	[512;0;512;0]	9.7	51	485	[2401,2584]
6	[1024;0;701;0]	9.2	41	368	[3318,5400]
7	[1024;0;960;0]	8.8	41	352	[3444,5400]
8	[523;0;1020;0]	6.7	41	268	[3940,5400]
9	[164;0;1024;0]	6	41	240	[4648,5284]

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

The average GOSPA, LE, FE, ME, and the estimate number of targets of the infrared image (a) GOSPA error, (b) Location error, (c) False error, (d) Miss error, (e) the number of targets.

5. 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.