Time-Matching Poisson Multi-Bernoulli Mixture Filter For Multi-Target Tracking In Sensor Scanning Mode

Abstract

In Bayesian multi-target tracking, the Poisson multi-Bernoulli mixture (PMBM) filter is a state-of-the-art filter based on the methodology of random finite set which is a conjugate prior composed of Poisson point process (PPP) and multi-Bernoulli mixture (MBM). In order to improve the random finite set-based filter utilized in multi-target tracking of sensor scanning, this paper introduces the Poisson multi-Bernoulli mixture filter into time-matching Bayesian filtering framework and derive a tractable and principled method, namely: the time-matching Poisson multi-Bernoulli mixture (TM-PMBM) filter. We also provide the Gaussian mixture implementation of the TM-PMBM filter for linear-Gaussian dynamic and measurement models. Subsequently, we compare the performance of the TM-PMBM filter with other RFS filters based on time-matching method with different birth models under directional continuous scanning and out-of-order discontinuous scanning. The results of simulation demonstrate that the proposed filter not only can effectively reduce the influence of sampling time diversity, but also improve the estimated accuracy of target state along with cardinality.

1. Introduction

Multi-target tracking (MTT) is a hot issue in sensor data fusion, which has a wide range of application domain, such as sensor network, ocean monitoring and vehicle tracking [1-3]. It is complicated to address the MTT problem because of its difficulty of correlating the target with the measurement due to factors such as false alarms and clutter in the measurement space, the dilemma to determine the new target appearance along with the old target disappearance and the limited field-of-view caused by the multi-sensor.

The global nearest neighbor (GNN) [4], the joint probabilistic data association (JPDA) [5] and the multiple hypothesis tracking (MHT) [6-7] are the earliest techniques proposed to complete MTT in the last century. Basically, these data association approaches are utilized to find the association between measurement and target state estimate, then single target filter is used to predict and update the target state estimate. However, these approaches cannot be implemented on account of their computational limitations in the strong clutter jamming circumstance [8]. In order to improve the real-time performance of the algorithm, the methods based on random finite set (RFS) is utilized to the MTT problem [9]. The multi-target state and measurement are modeled as the unified RFS and the posterior multi-target density is propagated forward to achieve the Bayesian recursive filtering of multi-target. The birth and death of the target are also modeled as Poisson process.

The MTT theory based on RFS has been continuously improved in the past twenty years. According to whether the predicted distribution and posterior distribution have the same form as the initial prior, the filter on RFS can be divided into the conjugate RFS filter and the non-conjugate RFS filter [10-11]. The probability hypothesis density (PHD) filter [12-13], the cardinalized probability hypothesis density (CPHD) filter [14], the multiple target multi-Bernoulli (MeMBer) filter [15] belong to non-conjugate RFS filter, i.e., to approximate the posterior multi-target density at low cost. The conjugate RFS filter possesses higher filtering accuracy and easier computation than the non-conjugate one due to the conjugacy property [16] and it can be divided into two types. The first conjugate prior is the Poisson multi-Bernoulli mixture (PMBM) filter which consists of the union of Poisson point process and multi-Bernoulli mixture (MBM) [16]. The addition of Poisson part can manage the number of underlying targets effectively which means that all targets that have never been detected are guaranteed. Track-oriented multiple hypothesis tracking (TOMHT) [6] can be used to complete the MBM part based on all data association hypotheses. As we can see, The PMBM filter is fascinating not only because of its conjugacy, but also nicely captures the relevant uncertainties [17]. The second conjugate prior is adding the label to the target state like delta-generalized labeled multi-Bernoulli (δ-GLMB) filter [18-19] which has been proved to be a special case of the PMBM filter [16]. Although it is capable of presenting the trajectory of target, it has a fault of long computation time.

We hope to combine theory and practice in the sensor field with the RFS-based MTT method. In the classical RFS methodology [9], the measurement is assumed to have the same sampling time because the multi-target state and measurement are modeled as a unified whole. Whereas, in the application scenarios (e.g., infrared sensor, the phased array, millimeter wave radar), it takes time to scan the whole surveillance area, for the reason that the beam width is limited and the sampling time of targets at the different positions can not be the same [20]. If we introduce method of RFS to the scanning of sensor application directly, the diversity of sampling time will lead to the mismatch between prediction step and update step in RFS theory, which seriously affects the performance of target tracking. Directing at the problem of time diversity of sensor target sampling, a time-labeled RFS is proposed, where a time-matching Bayesian filtering framework is constructed and the analytical implementations based on Gaussian mixture of PHD and GLMB are provided in refs. [21-22]. The scanning reveals the underlying nature of the RFS and thus inspires the application of time-matching.

In order to complete the time-matching framework further, we combine it with the first conjugate RFS filter, namely Poisson multi-Bernoulli mixture filter in the sensor scanning mode in this paper. We treat different kinds of targets discriminatively with combinatorial filtering performed in a parallel manner. The prediction and update formula of time-matching Poisson multi-Bernoulli mixture (TM-PMBM) filter is given and the implementation of TM-PMBM filter in Gaussian mixture model is provided. Finally, we verify the effectiveness of the proposed filter under various birth conditions and different scanning modes. Compared with other time-matching RFS filters, the TM-PMBM filter can address the mismatch problem between measurements and predicted states moreover the filtering accuracy is much higher.

The rest of this paper is organised as follows. Section 2 introduces the background of RFS and PMBM conjugate prior, the problems existing in sensor scanning is mentioned in detail. In Section 3, we provide the formula derivation and theoretical basis of TM-PMBM filter to figure out the problem of dealing with various sampling time. In Section 4, the implementation of TM-PMBM filter under linear and Gaussian model is given in detail. In Section 5, we simulate the tracking performance of TM-PMBM filter on different scenarios and compare it with other RFS filters on the same conditions. In Section 6, we draw conclusions.

2. Background

In this section, we briefly introduce the knowledge of random finite set and Poisson multi-Bernoulli mixture model. In Section 2.1, the principle of random finite set is mainly described and the problems of using random finite set theory in sensor applications are explained in Section 2.2. Section 2.3 gives the related concepts of PMBM conjugate prior.

2.1 Random finite set and Bayesian filtering

Suppose a multi-target tracking scenario on the mathematical real space ℝ containing rational and irrational numbers. Let Nx,k be the number of the targets and Nz,k be the number of measurements in the kth scan. There are target states x1k,⋯,xNx,kk ∈ 𝕏 and measurements z1k,⋯,zNz,kk ∈ ℤ where 𝕏 ⊂ ℝNx and ℤ ⊂ ℝNz denote state space and measurement space respectively. These targets and measurements are random variables with no specific order and the number is also stochastic. Therefore, the target set Xk and measurement set Zk can be regarded as the following random finite set (RFS) [9]:

Xk = {x1k,⋯,xNx,kk} ∈ F(𝕏),       (1)

Zk = {Z1k,⋯,zNz,kk} ∈ F(ℤ)       (2)

where F(𝕏) and F(ℤ) are the respective collections of all finite subsets of X_k and Z_k. Single point target state is composed by the target position, velocity, steering angle and other properties which describe the nature of the target motion completely.

The Bayesian filtering structure is applied to the scenario of random finite set to realize the tracking of multiple targets. Suppose that the dynamic state of multi-target at the end of the k-1th scan is X_k-1, the dynamic state of multi-target at the end of the kth scan is X_k, and the state transition density of multiple targets at the next scan can be considered as Markov transition density f_k|k-1(·). Besides, the likelihood function g_k(·) is used to describe the probability that the measurement set Z_k of the kth scan is generated by the multi-target X_k. Then according to the above conditions, the Chapman-Kolmogorov equation can be used to calculate the multi-target prediction [9]:

π_k|k-1(X_k|Z_1:k-1) = ∫f_k|k-1(X_k|X)π_k-1(X|Z_1:k-1)δX       (3)

According to the Bayes criterion, the multi-target posterior calculation is as follows:

\(\begin{aligned}\pi_{k}\left(X_{k} \mid Z_{1: k}\right)=\frac{g_{k}\left(Z_{k} \mid X_{k}\right) \pi_{k \mid k-1}\left(X_{k} \mid Z_{1: k-1}\right)}{\int g_{k}\left(Z_{k} \mid X_{k}\right) \pi\left(X_{k} \mid Z_{1: k-1}\right) \delta X_{k}}\end{aligned}\)       (4)

where the normalization constant is obtained by the form of set integration:

\(\begin{aligned}\int f(X) \delta X=\sum_{i=1}^{\infty} \frac{1}{i !} \int f\left(\left\{x_{1}, \cdots, x_{i}\right\}\right) d\left(x_{1}, \cdots, x_{i}\right)\end{aligned}\)        (5)

2.2 Time problem in filtering

Filtering framework is under ideal conditions in Section 2.1, nevertheless, the filter of the scenario do not follow the standard bayesian rule at the same sampling time in the real sensor application. The sensor obtains the data (e.g., location, velocity) temporarily, only when the beam scans to the targets. Take Fig. 1 as an example, the target 1 and 2 are represented as a triangle and a circle, respectively. Assume the beam is at the far right of the fan-shaped surveillance area when the scan starts. The measurements of the target 1 and 2 have been detected sequentially in process of the clockwise scanning sweeping to the left, therefore their sampling time are different. In other words, the measurements not on the same line as the sensor are bound to have different sampling times.

Fig. 1. A two-dimensional scenario is scanned from left to right at the kth period. There are two targets in the surveillance area which move in a straight line starting with states x¹_k and x²_k, ending with states x^1'_k and x^2'_k. Also measurements z¹_k and z²_k are generated within the period.

Conventional target tracking takes the end time of the scan as the same moment of prediction and update, i.e., the situation that goes through the whole scanning period to filter in (3)-(4), by the same token, the target also needs the measurement at the end of the scan to update the state estimate. However, the sampling time of the measurement is often inconsistent with the end of the scanning period due to the view within the scope, this will bring error to the update, resulting in the state estimate out of accuracy. As shown in Fig. 1, the real target states at end of the kth scanning period are provided, while targets generate the measurements within the kth scanning period only if the targets are scanned by the beam. Obviously, it is not correct to use the measurements and to update the target state throughout the whole scanning period. With the addition of scanning times, the change of sampling time in different periods will reduce the stability of state estimate and the error will increase by accumulation.

For addressing the problem of filtering above, it is necessary to consider the time information of the target to match the measured sampling time with the predicted state transition time. We hope to complete the recursive filtering of the time-matching for the standard point target with Poisson multi-Bernoulli mixture conjugate prior.

2.3 PMBM conjugate prior

PMBM conjugate prior contains two types of RFS to represent two kinds of targets [16]. On the one hand the Poisson point process (PPP) RFS is used to describe the target which has never been detected while still do exist at present, on the other hand the multi-Bernoulli mixture (MBM) RFS is established to model the potential detected target, either is previously detected, or it could be the new one currently. Two types of RFS describe the process of multi-target tracking dynamics as comprehensively as possible. Therefore, we can use concrete symbols to describe this classification method:

\(\begin{aligned}f^{p m b m}(X)=\sum_{X^{u} \uplus X^{d}=X} f^{p p p}\left(X^{u}\right) f^{m b m}\left(X^{d}\right)\end{aligned}\)       (6)

where X indicates the target set, X^u and X^d indicate the undetected target set and potential detected target set, respectively. Based on the properties of two kinds of RFS, the target set density can be expressed as:

\(\begin{aligned}f^{p p p}\left(X^{u}\right)=\exp \left(-\int \mu\left(x^{u}\right) d x^{u}\right) \prod_{x^{u} \in X^{u}} \mu\left(x^{u}\right),\end{aligned}\)       (7)

       (8)

where f^ppp(·) denotes the Poisson density, µ(·) >0 denotes the Poisson rate, f^mbm(·) denotes the MBM density, which can be represented by a series of normalized MB probability densities under the index j of MBM components. Set a total of potential detected targets, i is the index of the target, therefore w^i,j and f^i,j (X_i) repesent the weight of the target and the Bernoulli density. |·| denotes the number of variables.

Each global association history hypothesis follows a multi-Bernoulli distribution. A single global hypothesis is composed of multiple hypotheses of single potential detected target, which represents the possible association between potential detected target and measurement. It's worth noting that the relationship between the weight of Bernoulli components and the weight of MB W^j can be expressed as follows:

\(\begin{aligned}W^{j} \propto \prod_{i=1}^{\left|X^{d}\right|} w^{i, j}\end{aligned}\)       (9)

The density of the ith Bernoulli component can be expressed as follows:

\(\begin{aligned}f\left(X_{i}\right)=\left\{\begin{array}{ll}1-r^{i, j} & X_{i}=\varnothing \\ r^{i, j} p^{i, j}\left(x^{d}\right) & X_{i}=\left\{x^{d}\right\} \\ 0 & \text { otherwise }\end{array}\right.\end{aligned}\)       (10)

The potential detected targets need to reflect two aspects of uncertainty. On the one hand, the uncertainty of target existence is represented by the probability of target existence r^i,j and on the other hand, the uncertainty of target state is represented by the density of the existing state p^i,j(·). This is the reason that the Bernoulli model is the most appropriate description for the distribution of targets. In summary, PMBM propagates through the hybrid of PPP and MBM, calculating new PMBM density parameters through prediction and update of the Bayesian filtering framework, which can be referred in ref. [16].

3. Time-matching PMBM filter

In this section, we show how to realize the filtering iteration of PMBM in the sensor scanning. We describe time-matching PMBM conjugate prior in Section 3.1. The prediction and update process of the improved PMBM filter are provided in Section 3.2 and Section 3.3, respectively. Section 3.4 discusses the formation of sensor tracking under MBM filter.

3.1 Time-matching PMBM conjugate prior

In Section 2.2, we mentioned the problem that the filtering accuracy decreases due to the time mismatch between the prediction step and update step. Here, we divide the multi-target prior according to the sampling time, for the original unified transform entirety to components that can be calculated independently. This separating method can be understood as transforming the single filtering into the parallel filtering during the scanning period, where it can reduce the mutual interference between each component in the single filtering.

We want to realize time-matching in the framework of standard RFS Bayesian filters [22]. Before we start filtering, we make the following assumptions:

Assumption 1: Each target evolves independently and the target-derived measurements are independent of each other. Each target may give rise to at most one measurement theoretically in a single scanning period;

Assumption 2: For the process of scanning the whole surveillance area at single time, one region is scanned and only scanned once in the single sampling period;

Assumption 3: For the process of scanning one direction in the single sampling period, the measurement generated by the target will be scanned (detected) at most once.

Based on the above assumptions, we consider the sampling time into target tracking which means the factor of time adds to the target state, the improved target state set \(\begin{aligned}\tilde{X}_{k}\end{aligned}\) and measurement set \(\begin{aligned}\tilde{Z}_{k}\end{aligned}\) are as follow:

\(\begin{aligned}\tilde{X}_{k}=\left\{\left(X_{k}^{1}, t_{k}^{1}\right), \cdots,\left(X_{k}^{N_{k}}, t_{k}^{N_{k}}\right),\left(X_{k}^{U D}, t_{k}\right)\right\} \in F(F(\mathbb{X}) \times \mathbb{T}),\end{aligned}\)       (11)

\(\begin{aligned}\tilde{Z}_{k}=\left\{\left(Z_{k}^{1}, t_{k}^{1}\right), \cdots,\left(Z_{k}^{N_{k}}, t_{k}^{N_{k}}\right)\right\} \in F(F(\mathbb{Z}) \times \mathbb{T}),\end{aligned}\)       (12)

where tⁱ_k represents the distributed sampling time in a single scanning period, Xⁱ_k and Zⁱ_k respectively are the target state estimate and measurement sets under time tⁱ_k, X^UD_k(X^N_k+1_k) denotes missed detected target state, t_k represents the end of the time in the kth scan, 𝕏, ℤ and 𝕋 represent the target state space, measurement space and sampling time space.

We can regard the target state set with sampling time as an augmented target state set, which is similar to the label set [18]. Assumption 1 guarantees the independence between the divided components, given the posterior PMBM density \(\begin{aligned}f_{k-1}^{p m b m}\left(\tilde{X}_{k} \mid \tilde{Z}_{1: k-1}\right)\end{aligned}\) , posterior Poisson intensity f^ppp_k−1 and posterior MBM density f^mbm_k−1 at the end of the k-1th scan, we can obtain predicted PMBM density at the kth scan by substituting (6) into expression [22, (10)]:

\(\begin{aligned}f_{k \mid k-1}^{p m b m}\left(\tilde{X}_{k} \mid \tilde{Z}_{1: k-1}\right)=\sum_{X_{k}^{u} \uplus X_{k}^{d}=X_{k}}\left(\begin{array}{l}\left.\prod_{i=1}^{N_{k}^{u}}\left\langle f_{k \mid k-1}\left(\left(X_{k}^{i, u}, t_{k}^{i, u}\right) \mid(X, t)\right), f_{k-1}^{p p p}\left((X, t) \mid Z_{1: k-1}\right)\right\rangle\right) \\ \left.\prod_{j=1}^{N_{k}^{d}}\left\langle f_{k \mid k-1}\left(\left(X_{k}^{j, d}, t_{k}^{j, d}\right) \mid(X, t)\right), f_{k-1}^{m b m}\left((X, t) \mid Z_{1: k-1}\right)\right\rangle\right)\end{array}\right)\end{aligned}\)       (13)

where f_k|k-1 (·|·) is the state transition kernel. Besides, 〈·,·〉denotes an inner product operation. In addition, we can rewrite (4) to obtain the update distribution changing from the division of sampling time to the division depending on measurements:

\(\begin{aligned}\begin{array}{l}f_{k \mid k-1}^{p m b m}\left(\tilde{X}_{k} \mid \tilde{Z}_{1: k}\right)=f_{k \mid k}^{p m b m}\left(\left(X_{k}^{U D}, t_{k}\right) \mid \varnothing\right) \times \prod_{i=1}^{\left|Z_{k}\right|} f_{k \mid k}^{p m b m}\left(\left(X_{k}^{i}, t_{k}^{i}\right) \mid\left\{z_{k}^{i}\right\}\right) \\ =\prod_{i=1}^{\left|Z_{k}\right|} \frac{g_{k}\left(\left\{z_{k}^{i}\right\} \mid\left(X_{k}^{i}, t_{k}^{i}\right)\right) f_{k \mid k-1}^{p m b m}\left(\left(X_{k}^{i}, t_{k}^{i}\right) \mid Z_{1: k-1}\right)}{\left\langle g_{k}\left(\left\{z_{k}^{i}\right\} \mid\left(X_{k}^{i}, t_{k}^{i}\right)\right), f_{k \mid k-1}^{p m b m}\left(\left(X_{k}^{i}, t_{k}^{i}\right) \mid Z_{1: k-1}\right)\right\rangle},\end{array}\end{aligned}\)       (14)

where ∅ represents the empty set, i.e., the missed part and g_k(·|·) is the measurement density function whose sampling time is tⁱ_k. The proof of (14) is arranged from expressions [22, (13)-(16)] analogously. As we can see that the missed part is independent of other measurement associated part due to the lack of sampling time, therefore we should handle this kind of targets alone. After the above iteration period, it is critical that the parallel filtering makes the target achieve time-matching in the process of dynamic and measuring correction.

Different from the augmented space (e.g., [23]), the time factor is quantified by the measurement or scanning period edge, which does not need to be obtained by recursive estimation. Using the basic variable from (11)-(14), we can derive the innovative time-matching PMBM filter in the augmented space. The time-matching PMBM filter contains prediction process and update process, which are transmitted in the form of PMBM and given in Section 3.2 and Section 3.3, respectively.

3.2 Prediction Process

As we can see (13) that Poisson and Bernoulli part can be predicted separately. Given the posterior PMBM density at k-1th scanning period parameterized by Poisson intensity µ_k-1(·), MBM density {w^i,j_k-1, r^i,j_k-1, p^i,j_k-1(·)}, where w^i,j_k-1, r^i,j_k-1, p^i,j_k-1(·) are the weight, existence probability and Bernoulli density of ith target in jth global hypothesis, respectively. t' = {t¹_k-1,⋯,t^N_k-1,t^UD_k-1} is the sampling time set of k-1th scanning period. Poisson part of the prediction can be expressed using TM-PHD filter prediction equation from ref. [22]:

\(\begin{aligned}\mu_{k \mid k-1}(x, t)=\lambda_{k}^{b}(x, t)+\sum_{t^{\prime}}\left\langle p_{s, k}\left(x, t^{\prime}\right) f_{k \mid k-1}\left((x, t) \mid\left(x, t^{\prime}\right)\right), \mu_{k-1}\left(x, t^{\prime}\right)\right\rangle\end{aligned}\)       (15)

where λ^b_k is the newborn intensity of the kth scan, p_s and f_k|k-1 denote the probability of survival and state transition function, respectively.

Bernoulli part of the prediction according to expressions [11,(39)-(41)] can be attached to the sampling time expressed as follows:

_W^i,jk|k-1 = _W^i,jk-1,       (16)

\(\begin{aligned}r_{k \mid k-1}^{i, j}=r_{k-1}^{i, j} \sum_{t^{\prime}}\left\langle p_{s}\left(x, t^{\prime}\right), p_{k-1}^{i, j}\left(x, t^{\prime}\right)\right\rangle,\end{aligned}\)       (17)

\(\begin{aligned}p_{k \mid k-1}^{i, j}(x, t)=\sum_{t^{\prime}} \frac{\left\langle p_{s, k}\left(x, t^{\prime}\right) f_{k \mid k-1}\left((x, t) \mid\left(x, t^{\prime}\right)\right), p_{k-1}^{i, j}\left(x, t^{\prime}\right)\right\rangle}{\left\langle p_{s, k}\left(x, t^{\prime}\right), p_{k-1}^{i, j}\left(x, t^{\prime}\right)\right\rangle}\end{aligned}\)       (18)

The parameters of the predicted PMBM density are augmented to: μ_k|k-1(x, t) , {w^{i, j}_k|k-1, r^{i, j}_{k|k-1, p}^{i, j}_k|k-1(x, t)}. This part is mainly based on the sampling to realize the allocation of time and state in the prediction part.

3.2 Update Process

Given a predicted PMBM prior with parameter μ_k|k-1(x, t), {w^{i, j}_k|k-1, r^{i, j}_{k|k-1, p}^{i, j}_k|k-1(x, t)}, we can obtain the posterior PMBM density at the kth scan on the basis of ref. [11, Section Ⅲ. A] attached to the sampling time. Corresponding to the predicted part, , p_D,k(·) and g_k(·|·) c_k(z) denote the probability of detection, the measurement density function and the clutter intensity.

A) Undetected targets follow updated Poisson intensity:

μ_k(x, t) = (1- p_D,k(x, t))μ_k-1(x, t)       (19)

B) The Poisson components generates the potential first detected targets, i.e., a new Bernoulli component is created for each measurement:

r^p_k (z) = e_k(z) / ρ^p_k(z),       (20)

\(\begin{aligned}p_{k}^{p}(x, t \mid z)=\frac{p_{D, k}(x, t) g_{k}(z \mid x, t) \mu_{k \mid k-1}(x, t)}{e_{k}(z)},\end{aligned}\)       (21)

where

ρ^p_k(z) = e_k(z) + c_k(z),       (22)

e_k(z) = 〈p_D,k(x, t)g_k(z_k|x, t), μ_k|k-1(x, t)〉       (23)

C) Previous detected targets are subject to Bernoulli update, which is divided into two aspects: missed detection and measurement association.

For missed detected targets:

w^i,j,0_k = w^{i, j}_k|k-1ρ^i,j,0_k,       (24)

r^i,j,0_k = r^{i, j}_k|k-1ξ/ρ^i,j,0_k,       (25)

\(\begin{aligned}p_{k}^{i, j, 0}(x, t)=\frac{\left(1-p_{D, k}(x, t)\right) p_{k \mid k-1}^{i, j}(x, t)}{\rho_{k}^{i, j, 0}},\end{aligned}\)       (26)

where

ρ^i,j,0_k = 1-r^{i, j}_k|k-1 + r^{i, j}_k|k-1ξ,       (27)

ξ = 〈(1-p_D,k(x, t)), p^{i, j}_k|k-1(x,t)〉       (28)

On the other hand, for the detected targets:

r^i,j_k (z) = 1,       (29)

w^{i, j}_k = w^{i, j}_k|k-1r^{i, j}_k|k-1〈p_D,k(x, t),g_k(z|x, t)〉,       (30)

\(\begin{aligned}p_{k}^{i, j}(x, t \mid z)=\frac{p_{D, k}(x, t) g_{k}(z \mid x, t) p_{k \mid k-1}^{i, j}(x, t)}{\left\langle p_{D, k}(x, t) g_{k}(z \mid x, t), p_{k \mid k-1}^{i, j}(x, t)\right\rangle}\end{aligned}\)       (31)

After the update process, the measurement has been associated with one Bernoulli component at most in each global hypothesis and a large number of local hypotheses will be generated. The time element also becomes a part of the single target hypothesis. Consequently, the direct time correspondence relationship established by the association hypothesis can be used to deduce the complete time-matching method of PMBM filter. New global hypotheses are generated by traversing all possible combinations of the single target hypothesis by Murty algorithm [6]. The final state extraction is extracted from the global hypothesis with the highest weight, which needs pruning and selection. The specific operations of hypothesis weight are explained in Section 4.3.

3.4 Discussion

As can be seen from the above operations, the update process and the prediction process are measured with the same time, which further eliminates the problem that the sampling time is often inconsistent with the end time of the scanning period. Different types of target can also be divided independently on their sampling time to reduce the mutual influence. The first step is making a time prediction for each measurement, which overcomes the time problem, however it also brings an increase in computation. In order to solve this problem, we choose to add gating technology [24] to eliminate clutter, reducing the time of the prediction step. Since there is a one-to-one correspondence between the sampling time and the measurement, the update step of time-matching PMBM filter is equivalent to that of PMBM filter. In summary, the time-matching PMBM filter performs in parallel form and the schematic diagram of the TM-PMBM filter is shown in Fig. 2:

Fig. 2. Simple schematic diagram of the TM-PMBM filter

matching method. MBM filter turns Poisson birth model to Bernoulli or mixture Bernoulli birth model [25]. Prediction and update steps of the intensity of Poisson part set to be 0, for example the update exists without having to assign each measurement to a new potential target, because Bernoulli component is the clutter in the sense of probability. Therefore, we can transform the time component on the vanishing Poisson part into the new Bernoulli component, and the computation on the update part can be further reduced, and the matching of prediction step and update step can be achieved just like TM-PMBM. Compared with MBM birth model, Poisson birth model can cover the position of potential targets to a greater extent and is not limited by the maximum number of targets, therefore TM-PMBM also has a larger scope of application than TM-MBM.

4. Gaussian Mixture Implementation of time-matching PMBM filter

In this section, under the linear and Gaussian model, we give the implementation of the TM-PMBM filter in detail. In Section 4.1 and Section 4.2, we specifically clarify the process of the detected targets in addition to the undetected and missed targets. The final Section 4.3 describes the two-stage process of pruning global assumptions used by TM-PMBM filter to improve computational efficiency and target state extraction is also discussed.

Assumption 4: The given target follows a Gaussian linear motion model and the measurement model is also a Gaussian linear model:

f(x, t|x', t') = N(x; F(t|t')x', Q(t|t')),       (32)

h(z, t|x, t) = N(z; H(t) x, R),       (33)

where N(x; m,P) means the Gaussian density x with mean m and covariance P. F∈ℝ^n_x×n_x represents the state transition function of the single-time target, H∈ℝ^n_z×n_x represents the measurement transition function of the single-time target, Q∈ℝ^n_x×n_x represents the state transition function of the single-time target, and R∈ℝ^n_z×n_z represents the state transition function of the single-time target.

Assumption 5: The survival probability and detection probability in filter are independent of the expansion state of the target and are set as constant:

p_s,k(x', t') = p_s,k,p_D,k (x, t) = p_D,k       (34)

The posterior Poisson intensity and the posterior Bernoulli component density at the end of k-1th scanning period are set as:

\(\begin{aligned}\mu_{k-1}\left(x, t^{\prime}\right)=\sum_{q=1}^{J_{k-1}^{p}} w_{k-1}^{q}\left(t^{\prime}\right) N\left(x ; m_{k-1}^{q}\left(t^{\prime}\right), P_{k-1}^{q}\left(t^{\prime}\right)\right),\end{aligned}\)       (35)

\(\begin{aligned}p_{k-1}^{i, j}\left(x, t^{\prime}\right)=\sum_{q=1}^{J_{k-1}^{B}} w_{k-1}^{i, j, q}\left(t^{\prime}\right) N\left(x ; m_{k-1}^{i, j, q}\left(t^{\prime}\right), P_{k-1}^{i, j, q}\left(t^{\prime}\right)\right),\end{aligned}\)       (36)

where w^q and J^p denotes the weight and the number of Poisson components, w^i,j.q , r^i,j,q, J^B denotes the weight, existing probability and the number of Bernoulli components, respectively. Detected targets and undetected targets are distinguished by the superscript d and u.

In this section, we do not express the steps like traditional Bayesian prediction and update as Section 3, however do divide them according to sampling time like (13)-(14), which is more conducive for us to express the Gaussian mixture implementation of the improved filter clearly.

4.1 Detected potential targets Filtering

Prediction Step: Based on Assumption 1-5, there is a Gaussian mixture form of predicted Poisson density:

\(\begin{aligned}\mu_{k \mid k-1}^{d}\left(x, t_{m}\right)=\lambda_{k}^{b}\left(x, t_{m}\right)+\sum_{t^{\prime}} \sum_{q=1}^{J_{k-1}^{p}} w_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right) N\left(x ; m_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right), P_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right)\right),\end{aligned}\)       (37)

\(\begin{aligned}\lambda_{k}^{b}\left(x, t_{m}\right)=\sum_{q=1}^{J_{b}} w_{\gamma}^{q} N\left(x ; m_{\gamma}^{q}\left(t_{m}\right), P_{\gamma}^{q}\left(t_{m}\right)\right),\end{aligned}\)       (38)

where

w^d,q_k|k-1(t_m | t') = p_s,kw^q_k-1(t')       (39)

m^d,q_k|k-1(t_m | t') = F_k(t_m | t')m^q_k-1(t')       (40)

P^d,q_k|k-1(t_m | t') = F_k(t_m | t')P^d,q_k-1(t')F_k(t_m | t')^T + Q_k(t_m | t')       (41)

where λ^b_k(x, t), J_b and w^q_γ,k-1 denote newborn Poisson intensity, the number and the weight of newborn components, respectively.

Gaussian mixture form of predicted density with Bernoulli component:

\(\begin{aligned}r_{k \mid k-1}^{d, i, j}=r_{k-1}^{i, j} \sum_{t^{\prime}} \sum_{q=1}^{J_{k-1}^{i, j}} p_{S, k}\left(t^{\prime}\right) w_{k \mid k-1}^{d, i, j, q}\left(t^{\prime}\right),\end{aligned}\)       (42)

\(\begin{aligned}p_{k \mid k-1}^{d, i, j}\left(x, t_{m}\right)=\sum_{t^{\prime}} \sum_{q=1}^{j_{k-1}^{i, j}} w_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right) N\left(x ; m_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right), P_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right)\right),\end{aligned}\)       (43)

where

w^d,i,jq_k|k-1(t_m | t') = p_s,kw^i,j,q_k-1(t')       (44)

m^d,i,jq_k|k-1(t_m | t') = F_k(t_m | t')m^i,j,q_k-1(t')       (45)

P^d,i,j,q_k|k-1(t_m | t') = F_k(t_m | t')P^i,j,q_k-1(t')F_k(t_m | t')^T + Q_k(t_m | t')       (46)

It is important to note that, for both Poisson and Bernoulli components, the predicted time component here is assigned according to the single measurement time t_m. We add gating constraint to the prediction part of Bernoulli to reduce the number of calculated measurements, which greatly reduced the computational burden of the prediction step [24]. However, these measurements are not deleted, but they are not involved in the calculation of the prediction step. When forming a new global hypothesis, these measurements still need to be associated with the potential targets.

Update Step: The detected potential targets consist of two parts: the new Bernoulli components generated by measurements that may create new tracks and previous Bernoulli components associated with measurements.

Based on Assumption 1-5, the existence probability and the density of Bernoulli components for the first time given by (20)-(23) can be characterised by:

r^p_k(z) = e_k(z)/ρ^p_k(z),       (47)

\(\begin{aligned}p_{k}^{d}\left(x, t_{m} \mid z\right)=\sum_{q=1}^{J z} w_{k}^{d, q}\left(t_{m} \mid z\right) N\left(x ; m_{k}^{d, q}\left(t_{m} \mid z\right), P_{k}^{d, q}\left(t_{m} \mid z\right)\right),\end{aligned}\)       (48)

where

\(\begin{aligned}\begin{array}{c}\rho_{k}^{p}(z)=e_{k}(z)+c_{k}(z), \\ e_{k}(z)=\sum_{q=1}^{J_{Z}} p_{D, k} w_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right) \mu_{k \mid k-1}^{d}\left(x, t_{m}\right), \\ w_{k}^{d, q}\left(t_{m} \mid z\right)=p_{D, k} w_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right) N\left(z ; H_{k} m_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right), S_{k}^{q}\left(t_{m}\right)\right) / e_{k}(z), \\ K_{k}^{q}\left(t_{m}\right)=P_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right) H_{k}^{T}\left[S_{k}^{q}\left(t_{m}\right)\right]^{-1}, \\ \varepsilon_{k}^{q}\left(t_{m} \mid z\right)=z-H_{k} m_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right), S_{k}^{q}\left(t_{m}\right)=H_{k} P_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right) H_{k}^{T}+R_{k}, \\ m_{k}^{d, q}\left(t_{m} \mid z\right)=m_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right)+K_{k}^{q}\left(t_{m}\right) \varepsilon_{k}^{q}\left(t_{m} \mid z\right), P_{k}^{d, q}\left(t_{m} \mid z\right)=P_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right)-K_{k}^{q}\left(t_{m}\right) H_{k} P_{k \mid k-1}^{d, q}\left(t_{m} \mid t^{\prime}\right)\end{array}\end{aligned}\)       (49)

After the Bayesian update for the potential new tracks, the relevent measurements ought to be considered to connect with the previous tracks. Of course, the existence probability of hypothesis component should be 1, the weight and the density of Bernoulli components for the associated hypotheses given by (29)-(31) can be characterised by:

\(\begin{aligned}w_{k}^{d, i, j}(z)=w_{k \mid k-1}^{d, i, j} r_{k \mid k-1}^{i, j} \sum_{q=1}^{J_{k, j}^{i, j}} p_{D, k} w_{k \mid k-1}^{d, i, j, q} N\left(z ; H_{k} m_{k \mid k-1}^{d, i, j, q}\left(t_{m}\right), S_{k}^{i, j, q}\left(t_{m}\right)\right),\end{aligned}\)       (50)

\(\begin{aligned}p_{k}^{d, i, j}\left(x, t_{m} \mid z\right)=\sum_{q=1}^{J_{k \mid k-1}^{i, j}} w_{k}^{d, i, j, q}\left(t_{m} \mid z\right) N\left(x ; m_{k}^{d, i, j, q}\left(t_{m} \mid z\right), P_{k}^{d, i, j, q}\left(t_{m} \mid z\right)\right),\end{aligned}\)       (51)

where

\(\begin{aligned}\begin{array}{c}K_{k}^{i, j, q}\left(t_{m}\right)=P_{k \mid k-1}^{d, i, q}\left(t_{m} \mid t^{\prime}\right) H_{k}^{T}\left[S_{k}^{i, j, q}\left(t_{m}\right)\right]^{-1}, \\ \varepsilon_{k}^{i, j, q}\left(t_{m} \mid z\right)=z-H_{k} m_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right), S_{k}^{i, j, q}\left(t_{m}\right)=H_{k} P_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right) H_{k}^{T}+R_{k}, \\ w_{k}^{d, i, j, q}\left(t_{m} \mid z\right)=p_{D, k} w_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right) N\left(z ; H_{k} m_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right), S_{k}^{i, j, q}\left(t_{m}\right)\right) \\ \times\left[\sum_{q=1}^{J_{k, k-1}^{i, j}} p_{D, k} w_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right) N\left(z ; H_{k} m_{k \mid k-1}^{d, i, j, q}\left(t_{m}\right), S_{k}^{i, j, q}\left(t_{m}\right)\right)\right]^{-1}, \\ m_{k}^{d, i, j, q}\left(t_{m} \mid z\right)=m_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right)+K_{k}^{i, j, q}\left(t_{m}\right) \varepsilon_{k}^{i, j, q}\left(t_{m} \mid z\right), \\ P_{k}^{d, i, j, q}\left(t_{m} \mid z\right)=P_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right)-K_{k}^{i, j, q}\left(t_{m}\right) H_{k} P_{k \mid k-1}^{d, i, j, q}\left(t_{m} \mid t^{\prime}\right)\end{array}\end{aligned}\)       (52)

Similar to the Bernoulli component of the first detection, the previous detected potential targets also go through the whole Kalman update process.

4.2 Undetected and missed targets Filtering

Unlike in the case of Section 4.1, we uniformly identify the moment at the end of the kth scan period as this filtering moment t_e for both Poisson and Bernoulli components which spends the whole scanning time to filter, because it fail to find the appropriate time to match them.

Prediction Step: Gaussian mixture form of predicted Poisson density can be derived as follows:

\(\begin{aligned}\mu_{k \mid k-1}^{u}\left(x, t_{e}\right)=\lambda_{k}^{b}\left(x, t_{e}\right)+\sum_{q=1}^{J_{k-1}} w_{k \mid k-1}^{u, q}\left(t_{e} \mid t^{\prime}\right) N\left(x ; m_{k \mid k-1}^{u, q}\left(t_{e} \mid t^{\prime}\right), P_{k \mid k-1}^{u, q}\left(t_{e} \mid t^{\prime}\right)\right)\end{aligned}\)       (53)

\(\begin{aligned}\lambda_{k}^{b}\left(x, t_{e}\right)=\sum_{q=1}^{J_{b}} w_{\gamma}^{q} N\left(x ; m_{\gamma}^{q}\left(t_{e}\right), P_{\gamma}^{q}\left(t_{e}\right)\right)\end{aligned}\)       (54)

Gaussian mixture form of predicted Bernoulli density can be derived as follows:

\(\begin{aligned}r_{k \mid k-1}^{u, i, j}=r_{k-1}^{i, j} \sum_{q=1}^{\substack{J_{k-1}^{i, j}}} p_{S, k}\left(t^{\prime}\right) w_{k \mid k-1}^{u, i, j, q}\left(t^{\prime}\right),\end{aligned}\)       (55)

\(\begin{aligned}p_{k \mid k-1}^{u, i, j}\left(x, t_{e}\right)=\sum_{q=1}^{J_{k-1}^{i, j}} w_{k \mid k-1}^{u, i, j, q}\left(t_{e} \mid t^{\prime}\right) N\left(x ; m_{k \mid k-1}^{u, i, j, q}\left(t_{e} \mid t^{\prime}\right), P_{k \mid k-1}^{u, i, j, q}\left(t_{e} \mid t^{\prime}\right)\right),\end{aligned}\)       (56)

The parameters of Gaussian components have the same form of calculation just as (39)-(41) and (44)-(46). It is important that the time allocation of these targets is the end of the scan rather than the sampling time of each measurement.

Update Step: Posterior Bernoulli density with undetected targets can be stated:

\(\begin{aligned}\mu_{k}^{u}\left(x, t_{e}\right)=\sum_{q=1}^{J_{k}} w_{k}^{u, q}\left(t_{e}\right) N\left(x ; m_{k}^{u, q}\left(t_{e}\right), P_{k}^{u, q}\left(t_{e}\right)\right),\end{aligned}\)       (57)

where

w^u,q_k(t_e)=(1- p_D,k)w^u,q_k|k-1(t_e | t'),       (58)

m^u,q_k(t_e)=m^u,q_k|k-1(t_e | t'),       (59)

P^u,q_k(t_e)=P^u,q_k|k-1(t_e | t'),       (60)

where the parameters of Gaussian mixture are inherited by the predicted Poisson part so does the number J_k = J_k|k-1.

Posterior Bernoulli component density for misdetection potential targets is given by (24)-(26) which can be simplified as:

^i,j,0_k(t_e)=w^u,i,j_k|k-1(t_e | t')(1-r^u,i,j_k|k-1p_D,k),       (61)

\(\begin{aligned}r_{k}^{i, j, 0}=\frac{r_{k \mid k-1}^{u, i, j}\left(1-p_{D, k}\right)}{1-r_{k \mid k-1}^{u, i, j} p_{D, k}},\end{aligned}\)       (62)

\(\begin{aligned}p_{k}^{i, j, 0}\left(t_{e}\right)=p_{k \mid k-1}^{u, i, j, q}\left(t_{e}\right)=\sum_{q=1}^{j_{k, j}^{i, j-1}} w_{k}^{u, i, j, q}\left(t_{e}\right) N\left(x ; m_{k}^{u, i, j, q}\left(t_{e}\right), P_{k}^{u, i, j, q}\left(t_{e}\right)\right),\end{aligned}\)       (63)

where

m^u,i,j.q_k(t_e) = m^u,i,j.q_k(t_e | t'),       (64)

P^u,i,j.q_k(t_e) = P^u,i,j.q_k|k-1(t_e | t')       (65)

4.3 Pruning and State extraction

Having gone through Section 4.1 and Section 4.2, we have completed the association of the target with all the measurements, forming the available target hypothesis, and make the update in accordance with the time-matching model. In order to obtain the global hypothesis required in the kth scan period, we associate the global hypothesis j of the previous scan period through all the current data, which leads to the problem of excessive number of hypotheses. In this case, we need to determine the number of global hypotheses.

The formation and deletion of global hypotheses are mainly selected by Murty algorithm to select new global hypotheses with the highest weight. The cost matrix C_j is generated according to the associated weight of the generated potential detected targets and the measurements Z:

C_j = -ln[W_O, W_N],       (66)

where

\(\begin{aligned}\begin{array}{l}W_{O}=\left[\begin{array}{ccc}w_{O}^{1,1} & \cdots & w_{O}^{1, n} \\ \vdots & \ddots & \vdots \\ w_{O}^{m, 1} & \cdots & w_{O}^{m, n}\end{array}\right], \\ W_{N}=\operatorname{diag}\left(w_{N}^{1,1}, \cdots, w_{N}^{m, m}\right), \\ w_{O}^{i, j}=\frac{w_{k \mid k-1}^{i, j} r_{k \mid k-1}^{i, j}\left\langle p_{D, k} g_{k}\left(z \mid t_{k \mid k-1}\right), p_{k \mid k-1}^{i, j}\left(t_{k \mid k-1}\right)\right\rangle}{1-r_{k \mid k-1}^{i, j} p_{D, k}} \\\end{array}\end{aligned}\)       (67)

where W_O represents the weight matrix of W_O after removing clutter generated by updating the n old tracks in the m measurements, Then W_N is a diagonal matrix of pure measurement weights. It is clear that new k global hypotheses can be selected by minimizing tr(SC) using Murty algorithm to generate assignment matrix S of 0 or 1 entries from the diagonal matrix generated by measurements.

After that, it is still necessary to delete useless global hypotheses and Bernoulli components to save computing resources, which can be divided into four main aspects. The first is to retain the global assumption that the weight is higher than the threshold. Second, each potential track is correlated with all measurements. Obviously, a target can only generate one measurement at most and the unnecessary association hypothesis must be deleted. Third, the Bernoulli component with a lower sum of the existence probability of all associations will also increase the computational burden and these will reasonably be removed. Fourth, the global hypotheses of possible duplications generated by the deletion operation are merged into one by adding their weights and removing the duplications.

The last step of the time-matching PMBM filter is the extraction of target states. Ref. [16] puts forward three estimators. Estimator 1 is the estimation of the excess threshold extracted from the selected global hypothesis with the highest weight, Estimator 2 uses the maximum posterior estimation of the cardinality with the highest weight in GLMB. Estimator 3 uses the maximum posteriori estimation of the deterministic cardinality with the highest weight.

5. Simulation and Results

5.1 Basic simulation Settings

To prove the tracking performance of time-matching PMBM filter, we design simulation with two-dimensional challenging scenarios in this section. Suppose that there are 6 targets appear in the surveillance area and each target obeys the linear Gaussian dynamic and measurement model in accordance with (32), (33), the scenarios examined are exceptionally challenging due to constant acceleration (CA) model which is expressed as follows:

\(\begin{aligned}F_{k}\left(t \mid t^{\prime}\right)=I_{2} \otimes\left[\begin{array}{ccc}1 & t-t^{\prime} & \frac{1}{2}\left(t-t^{\prime}\right)^{2} \\ 0 & 1 & t-t^{\prime} \\ 0 & 0 & \exp \left(-\frac{t-t^{\prime}}{\tau}\right)\end{array}\right], Q_{k}\left(t \mid t^{\prime}\right)=\sigma_{w}^{2}\left(1-e^{\frac{-2\left(t-t^{\prime}\right)}{\tau}}\right) \otimes\left[\begin{array}{ccc}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 10\end{array}\right],\end{aligned}\)       (68)

\(\begin{aligned}H_{k}(t)=\left[\begin{array}{llllll}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0\end{array}\right], R_{k}=\sigma_{v}^{2}\left[\begin{array}{ll}1 & 0 \\ 0 & 1\end{array}\right],\end{aligned}\)       (69)

where I₂, t−t′ and τ = 1 represent 2×2 identity matrix, sampling time interval and maneuver correlation time, besides the standard deviation of the process and measurement noise σ_w, σ_v are 5 and 10. The state of the target in the augmented space can be written as [p_x, v_x, a_x, p_y, v_y, a_y], where p_x and p_y are the x and y axises of the target, v_x, v_y and a_x, a_y are the fractional velocities and acceleration in the x and y direction. In addition, the first and second dimensions of measurement obtained are the x and y axises of the measurement. Next are some settings for the scenario parameters. The sensor is located at (0, 0) in terms of coordinates. The size of scenario is set to [0, 2000] in range as well as [0°, 180°] in angle and the clutter generated in each scan obeys the uniform distribution and its number obeys the Poisson distribution with the parameter λ = 10. The probability of survival and detection are 0.99 and 0.85 respectively.

The initial state, birth time and death time of targets in the specific state space are shown in Table 1 and the true states of targets are reflected in Fig. 3. Specially, the PMBM filter implementation uses ellipsoidal gate, a maximum number of global hypotheses N_h =10 and Estimator 3. We use a pruning threshold of 10^-5 for the Poisson part and eliminate Bernoulli components with existence probability whose is lower than 10^-5.

Table 1. The Information of True Targets

Fig. 3. The real point tracks of the targets. Various targets are marked by discernible colors.

The key to the time-matching method is to capture the sampling time of the measurements accurately. In practical sensor applications, it is to divide the surveillance area into several fan-shaped sectors of equal interval. If the target falls into the sector, the sampling time of the target can be identified as the time in the middle of the sector whose prerequisite is that the surveillance area must be demarcated sufficiently. Consequently, we divide the monitoring area into 1000 equal parts based on the situation that the sensor is considered as the origin.

The generalized optimal subpattern assignment (GOSPA) metric is used to evaluate the performance of augmented trajectory estimation [26]. Let X, Y, γ* denote the ground truth set, its estimate and the optimal assignment. The GOSPA error d^(c,α)_p(·) is decomposed as:

d^(c,α)_p(X, Y) = [c^p_l (X, Y, γ*) + c^p_f (γ*) + c^p_m (γ*)]^1/p,       (70)

where, c^p_l⁽·) denotes localization cost, c^p_f⁽·) denotes false target cost and c^p_l⁽·) denotes missed target cost, whose parameters are set as p = 2, c = 100, α = 2. Besides, the root mean square error (RMS) of Monte Carlo runs is obtained to analyze the filter. We choose PHD filter [12], MBM filter [25], PMBM filter [16] and GLMB filter [10] as prototypes to verify the effectiveness and robustness of our fashion. GLMB filter is similar to MBM₀₁ filter structurally, which can be implemented by joining prediction and update formulation, namely joint-GLMB in ref. [19]. We arranged two scenarios with different birth model to test the RFS filters and their time-matching versions [22].

In addition of the structure of filters themselves, the scanning mode of sensor plays an important role in the tracking performance of the filters. Scanning model is mainly divided into continuous scanning and discontinuous scanning. The traditional mechanical scanning mode is that directional continuous which scans repeatedly from 180°-0° in clockwise direction or from 0°-180° in anticlockwise direction. Electrical scanning mode not only completes the mode that mechanical scanning does, but also realizes intermittent out-of-order scanning mode. In particular, the order of discontinuous scanning is shown in the Table 2 under the 20 sectors. Both modes use the velocity of 180°/s.

Table 2. The Order of non-continuous Scanning Model

5.2 Scenario 1

In the first scenario, we identify several birth places with low-weight which can be found in ref. [12]. The birth model parameterized by n_b = 3, w^b_k = r^b_k = 0.03, the mean of Gaussian components are m^b,1_k = [1000,0,0,100,0,0] , m^b,2_k = [-1500,0,0,1000,0,0] and m^b,3_k = [-1200,0,0,100,0,0] . The covariance of Gaussian components is set as P^b_k = diag([10,10,10,10,10,10])². We conduct 200 Monte Carlo runs with on Matlab 2021a 2.50 GHz Intel I5 laptop to obtain root mean square GOSPA error in Fig. 4, estimated cardinality [9] in Fig. 5 and time cost in Fig. 6. Fig. 4 can be refined as: (a)-(d) indicate the results of directional continuous scanning mode and (e)-(h) indicate the results of out-of-order mode. In this case, (a)(e) for GOSPA error, (b)(f) for localization error, (c)(g) for false target error, (d)(h) for missed target error. Fig. 5 (a) indicates the results of directional scanning mode and (b) indicates the results of out-of-order scanning mode. Fig. 6 (a) and (b) indicate the results of directional scanning mode and out-of-order scanning mode, respectively.

Fig. 4. Performance comparison among four RFS filters with the time-matching version in Scenario 1.

Fig. 5. Estimated cardinality against time of filters in Scenario 1.

Fig. 6. Time cost of filters in Scenario 1.

Under this birth model, we can see that time-matching Bayesian filtering framework can improve the tracking performance of the four RFS-based filters in MTT problems. The original filters have significantly larger localization error than the time-matching RFS filters in case of little difference among false target error, missed target error and cardinality from Fig. 4 and Fig. 5. Contrary to the increase of localization error due to the incremental number of targets in the whole simulation process over time in PMBM filter, TM-PMBM filter appears to be is more stable than PMBM filter, for the reason that the diversity of sampling time is taken into account and the additional measurement error is reduced under both scanning mode. It is obvious that PMBM filter has greater fluctuations on localization error than that of TM-PMBM filter in out-of-order discontinuous scanning. TM-Bernoulli filter is better than TM-PHD filter without error in the first moment approximation. TM-PMBM filter and TM-MBM filter have better ability to deal with global hypotheses than TM-GLMB filter for association ambiguity. Besides, TM-PMBM filter take slight advantage over TM-MBM filter in this birth model.

All the errors under out-of-order scanning are more fluctuating than one under directional continuous scanning. There are lots of high peaks in false target error and missed target error in PHD filter under out-of-order scanning compared to that under directional scanning because the former mode lead to the repeated computation due to the growing number of targets and increasing velocity. This is reflected in the Fig. 6 (b) that the number of targets filter tracks decreases dramatically after step 75. Similar to PMBM filter, the TM-PMBM filter can overcome this shortcoming maximally that the error curve is more flatter and the cardinality is more close to the real one which exhibits less impact from out-of-order scanning.

The operation of time-matching will make the TM-PMBM and TM-MBM consume more time than that of original ones due to the structure of multi-Bernoulli mixture. In spite of this, they are still faster than GLMB and TM-GLMB. Subsequently, we can address the problem by simplifying PMBM into Poisson Multi-Bernoulli (PMB) distribution to obtain the best fitting PMB filter [27] in time-matching method and handle time bottleneck in our future work.

5.3 Scenario 2

In the second scenario, we do not stare at the known areas, but take care of the wide range of area covered the region of interest which is more common in sensor application. The number of birth places decreased to one and the covariance is expanded to a certain extent. The birth model is parameterized by n_b = 1, w^b_k = r^b_k = 0.03 , the mean of Gaussian component is m^b,1_k = [0,0,0,1000,0,0], and the covariance is set as P^b_k = diag([1000,100,100,1000,100,100])².

Likewise, we conduct 200 Monte Carlo runs in both scanning mode so as to obtain root mean square GOSPA error in Fig. 7, estimated cardinality in Fig. 8 and time cost in Fig. 9. The sorts of results are the same as that in Fig. 4 to Fig. 6. The order of discontinuous scanning is the same as Scenario 1’s using the velocity of 180°/s. Similar to Scenario 1, the time-matching TM-PMBM filter exhibit significantly lower localization error than that of the PMBM filter and the TM-PHD filter behaves badly in the out-of-order scanning mode at the same way. In a nutshell, the TM-PMBM filter performs better than any other filters in this scenario under both scanning modes. For missed targets error, it is still the best.

Fig. 7. Performance comparison among four RFS filters with the time-matching version in Scenario 2.

Fig. 8. Estimated cardinality against time for the filters in Scenario 2.

Fig. 9. Time cost of the original filters and the time-matching filters in Scenario 2.

The discernible differences in the birth model can still make a great influence on results. First of all, there are high peaks in step 20, 40 on account of new birth in scenario 1 as we can see in the localization error, whereas, localization error curve in the scenario 2 is more flatter and the false target error is much lower. The total GOSPA error in scenario 2 produces more violent shaking than that in scenario 1 due to the occurrence of larger missed detection error because it lower the probability of known birth information. We can also see that all filters take time to reach the true number of the targets and TM-PMBM filter needs the least time in Fig. 8 and Fig. 9.

One thing worth mentioning i.e., the TM-PMBM filter have an edge over TM-MBM filter in term of missed target error because of its uncertainty that measurements associated with the Bernoulli components when the priori density of birth is more extensive. Creating new targets from measurements makes TM-PMBM filter exhibits lower missed target error than that of the other filters to a certain extent and TM-PMBM filter is more widely used than TM-MBM.

6. Conclusion

In this paper, aiming at the problem of dealing with sampling time diversity in sensor scanning process, we propose the Poisson multi-Bernoulli mixture filter based on time-matching method. We provide the calculation formula of different types of targets after adding sampling time and use a reasonable representation of conjugate prior for the union of PPP and MBM. Simulation based on Gaussian mixture model also demonstrate that the superiority of the structure and conjugate prior form for the time-matching Poisson multi-Bernoulli mixture filter and validity of coping with the time problem that the scanning process of prediction and update do not match. Whether in the selected birth mode or scanning mode, the time-matching Poisson multi-Bernoulli mixture filter is superior to other filters.

In the further research, we have many directions to practice. How to introduce the time-matching Poisson multi-Bernoulli mixture filter into the multi-sensor and multi-scan is worth thinking about. Another work is an extension of the work in this paper, which considers the time-matching PHD filter and PMB filter on the trajectory set, in order to achieve the extraction of complete trajectory information.