*Article* **Evolutionary Optimization Based Set Joint Integrated Probabilistic Data Association Filter**

**Shuang Liang 1, Yun Zhu 2,\* and Hao Li 3,\***


**Abstract:** The joint integrated probabilistic data association (JIPDA) algorithm is widely used for the automatic tracking of multiple targets, but it has the well-known problem of track coalescence. By optimizing the posterior density, the accuracy of the target state estimation can be improved. Motivated by this idea, we developed a novel evolutionary optimization based joint integrated probabilistic data association (EOJIPDA) filter to overcome the coalescence problem of the JIPDA filter. The trace for the covariance matrix of the posterior density is used as the objective function for the above optimization problem. It is shown that the accuracy of the target state estimation can be improved by reducing the trace. Evolutionary optimization was employed to minimize the trace and optimize the posterior density. More specifically, we enumerated all the possible permutations of the targets and assign a unique index to each permutation. The resulting indices were randomly assigned to all possible association hypothesis events. Each assignment indicated one possible gene in the evolutionary algorithm. This process was repeated several times to arrive at the initial population. An illustrative example shows that the EOJIPDA filter can effectively improve the accuracy of state estimation. Numerical studies are presented for two challenging multi-target tracking scenarios with clutter and missed detections. The experimental results demonstrate that the EOJIPDA filter provides better tracking accuracy than traditional coalescence-avoiding methods.

**Keywords:** multi-target tracking; evolutionary optimization; random finite set; joint integrated probabilistic data association

#### **1. Introduction**

Multi-target tracking (MTT) is one of the most important low-level techniques used in radar, computer vision, Internet of things [1,2], and other surveillance systems [1–5]. Traditional MTT methods usually regard the MTT problem as the tracking of multiple single-targets. Examples of such methods are the multiple hypothesis tracking (MHT) [6,7] and joint probabilistic data association (JPDA) [8] filters. The essence of these filters is the association between the measurements and the targets being tracked. The MHT filter propagates all possible association hypothesis events over time. The target states can be estimated accurately from sufficient events. However, the main drawback of the MHT filter is its heavy computation. Compared with the MHT filter, the JPDA filter requires less computation and is effective at MTT. Nevertheless, the JPDA filter has the track coalescence problem, whereby tracks tend to coalesce when the targets are closely spaced. The JPDA filter has a well-known drawback in that it assumes that the number of targets and their initial states are known a priori. Furthermore, the number of targets remains constant during tracking.

In practice, the number and states of targets are usually unknown in advance and the number of targets may vary with time. The automatic tracking system is a natural choice

**Citation:** Liang, S.; Zhu, Y.; Li, H. Evolutionary Optimization Based Set Joint Integrated Probabilistic Data Association Filter. *Electronics* **2022**, *11*, 582. https://doi.org/10.3390/ electronics11040582

Academic Editor: Marco Mussetta

Received: 10 January 2022 Accepted: 10 February 2022 Published: 15 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to track targets in this situation. Both target and clutter measurements are used by the automatic tracking system to initiate tracks. To distinguish between false and true tracks, it is necessary to measure the quality of each track. The joint integrated probabilistic data association (JIPDA) filter uses the probability of target existence to measure the quality of the track and is effective for automatic tracking [9]. However, the JIPDA filter also suffers from the track coalescence problem.

To overcome the track coalescence problem of the JIPDA filter and further improve the tracking accuracy, we combine the advantages of the RFS theory and the standard JIPDA filter and propose a novel filter, named evolutionary optimization based joint integrated probabilistic data association (EOJIPDA). The essential reason for the track coalescence is that the overlap of the tracking gates leads to the association uncertainty problem. In this case, the posterior density becomes multimodal and the estimation of the multi-target state becomes less accurate. The trace of the covariance matrix is a measure of the multimodality of the posterior density [10,11]. To improve the accuracy of the multi-target state estimation, we used the evolutionary computation approach to minimize the trace of the covariance matrix and optimize the posterior density. First, it is necessary to define a population of possible solutions in the search space for the optimization problem. We enumerated all the possible permutations of the targets and assigned a unique index to each permutation. The resulting indices were randomly assigned to all the possible data association events. Each assignment indicated one possible gene in the evolutionary algorithm. This process was repeated several times until we obtained the initial population. Each candidate solution was assigned a fitness value based on the cost function. Fitter individuals have a higher chance of mating, producing more "fitter" individuals. The crossover and mutation are then applied to the population to generate the new offspring. This process is repeated until the stopping criterion is reached. The main contributions of the proposed EOJIPDA filter are summarized as follows.


The organization of this paper is as follows. Section 2 introduces the target and measurement models, and the necessary background on the JIPDA filter. In Section 3, the proposed EOJIPDA filter is described in detail. Section 4 evaluates the performance of the EOJIPDA filter using two different tracking scenarios. Finally, Section 5 presents concluding remarks.

#### **2. Related Work**

To date, several approaches have been proposed to overcome track coalescence, such as the exact nearest-neighbor JPDA (ENNJPDA) filter [12], the JPDA\* filter [13], the set JPDA (SJPDA) filter [10], the Kullback–Leibler SJPDA (KLSJPDA) filter [10], the multi-objective JPDA (MOJPDA) [14], and the nearest-neighbor SJPDA (NNSJPDA) filter [15]. Motivated by the JPDA\* filter [13], the JIPDA\* filter [16] has also been proposed to overcome the track coalescence problem of the traditional JIPDA filter. In the JIPDA\* filter, the target states are estimated by selecting one association hypothesis event and the other association hypothesis events are pruned. The track coalescence problem can be avoided effectively by the JIPDA\* filter, but the association hypothesis events being pruned may contain some useful information. Our recent work [11] proposes to use the information contained in all the association hypothesis events by optimizing the posterior density. The proposed algorithm strongly depends on initializations and improper initializations may result in bad local minima.

Finite set statistics (FISST)-based MTT algorithms [17], such as the probability hypothesis density (PHD) filter [18], the cardinalized PHD (CPHD) filter [19], the multi-Bernoulli (MB) filter [20], the label multi-Bernoulli (LMB) filter [21], and the generalized label multi-Bernoulli (GLMB) filter [22,23], have attracted significant attention in recent years. These filters use the random finite set (RFS) theory to model the uncertainty in systems. Among them, the PHD and CPHD filters are the moment approximations of the Bayes multi-target filter and have no analytical expression for the posterior multi-target density. Unlike the PHD and CPHD filters, the MB filter propagates parameters of the multi-Bernoulli RFS instead of the posterior multi-target density. By using the labeled RFS, the LMB and GLMB filters have also been proposed as approximations of the multi-target Bayes filter. These filters can be implemented by the Gaussian mixture (GM) method or the sequential Monte Carlo (SMC) method, depending on the characteristics of the tracking system. In Table 1, we list all the above-mentioned existing MTT methods.


**Table 1.** Existing MTT methods.


JPDA\* denotes a coalescence-avoiding version of JPDA, as shown in [13].

#### **3. Background**

*3.1. Target and Measurement Assumptions*

We consider the nearly constant velocity model in this paper. Assuming that the target state at time *k* is denoted as *xk*, the target motion is modelled as follows

$$\mathbf{x}\_{k+1} = F\_k \mathbf{x}\_k + \mathbf{u}\_k \tag{1}$$

where *F<sup>k</sup>* is the transition matrix and *u <sup>k</sup>* ∼ N (0, *Qk*) is zero-mean Gaussian distributed process noise with covariance *Qk*. If the number of all potential targets (tracks) at time *k* is *nk*, the augmented vector of labeled target states is

$$\mathbf{X}\_{k} = \left[ \left( \mathbf{x}\_{k}^{1} \right)^{\mathrm{T}}, \left( \mathbf{x}\_{k}^{2} \right)^{\mathrm{T}}, \dots, \left( \mathbf{x}\_{k}^{n\_{k}} \right)^{\mathrm{T}} \right]^{\mathrm{T}} \tag{2}$$

At time *k*, the relation between the target *x<sup>k</sup>* and its measurement *z<sup>k</sup>* is described as

$$z\_k = H\_k \mathbf{x}\_k + \mathbf{u}\_k \tag{3}$$

where *H<sup>k</sup>* is the observation matrix and *w<sup>k</sup>* ∼ N (0, *Rk*) denotes that *w<sup>k</sup>* follows the Gaussian distribution with zero-mean and covariance *Rk*. The target-originated measurements are detected with the probability of detection *PD*. Measurements may also originate from clutters. It is assumed that the number of clutter measurements follows the Poisson distribution and the average number of clutter measurements at each time step is *r* = *λ*|FoV| , where *λ* is the intensity of clutter and |FoV| denotes the sensor's field-of-view. The set of the validated measurements received at time *k* is given by

$$\mathbf{Z}\_k = \left\{ \mathbf{z}\_{k'}^1 \mathbf{z}\_{k'}^2 \cdot \cdots \cdot \mathbf{z}\_k^{m\_k} \right\} \tag{4}$$

where *mk* denotes the number of measurements. The sequence of measurement sets accumulated up to time *k* is described as *Z*1:*<sup>k</sup>* = {*Z*1,*Z*2, ··· ,*Zk*}.

#### *3.2. Joint Integrated Probabilistic Data Association Filter*

In the JIPDA filter, association hypothesis events are formed by assigning the validated measurements to the targets being tracked. From this, the target state *x<sup>t</sup> <sup>k</sup>* and covariance *Pt <sup>k</sup>* are estimated for each track *t* at time *k.* What's more, the estimation of the probability of target existence *r<sup>t</sup> <sup>k</sup>* is also taken into considered. Therefore, a track *t* can be completely described by the parameter set *rt <sup>k</sup>*, <sup>N</sup> (*x<sup>t</sup> k*,*P<sup>t</sup> k*)). , where <sup>N</sup> (*x<sup>t</sup> k*,*P<sup>t</sup> <sup>k</sup>*) is a Gaussian PDF with mean *x<sup>t</sup> <sup>k</sup>* and covariance *<sup>P</sup><sup>t</sup> <sup>k</sup>*. At time *k* + *1*, the main steps used to estimate the parameter set are as follows.

• Step 1: Predicting the target state *<sup>x</sup><sup>t</sup> <sup>k</sup>*+1|*k*, the covariance *<sup>P</sup><sup>t</sup> <sup>k</sup>*+1|*k*, and the probability of target existence *r<sup>t</sup> <sup>k</sup>*+1|*<sup>k</sup>* for each track *<sup>t</sup>*.

The predicted target state *x<sup>t</sup> <sup>k</sup>*|*k*−<sup>1</sup> and covariance *<sup>P</sup><sup>t</sup> <sup>k</sup>*|*k*−<sup>1</sup> are estimated using the Kalman filter. The propagation of the target' existence probability follows the Markov chain one,

$$r\_{k+1|k}^t = p\_{11}r\_k^t + p\_{21}(1 - r\_k^t) \tag{5}$$

where *p*<sup>11</sup> and *p*<sup>21</sup> are the Markov chain coefficients.

• Step 2: The tracking gate is generated for each target to select the validated measurements.

After the prediction, gating is performed for each track by defining an area around the predicted track. The area is named as "gate" and only the measurements falling within the gate are regarded as the validated measurements and are used to update the track. The authors of [24] present various gating techniques, among which the ellipsoidal gate [25] is widely used. The ellipsoidal gating area for track *t* at time *k* is defined as follows,

$$\mathcal{R}\_k^t = \left\{ \mathbf{z}\_k^j \in \mathbb{R}^{d\_z} \, \middle| \, d^t(\mathbf{z}\_k^j) \le G \right\} \tag{6}$$

where

$$d^t(z\_k^j) = v(k, t, j)^\mathsf{T} \left(\mathbf{S}\_k^t\right)^{-1} v(k, t, j)^\mathsf{T} \tag{7}$$

denotes the distance between predicted measurement *z<sup>t</sup> <sup>k</sup>*|*k*−<sup>1</sup> and received measurement *z j <sup>k</sup>*, *v*(*k*, *t*, *j*) = *z j <sup>k</sup>* <sup>−</sup> *<sup>z</sup><sup>t</sup> <sup>k</sup>*|*k*−<sup>1</sup> denotes the difference between *<sup>z</sup><sup>t</sup> <sup>k</sup>*|*k*−<sup>1</sup> and *<sup>z</sup> j <sup>k</sup>*, *<sup>S</sup><sup>t</sup> <sup>k</sup>* is the innovation covariance, *dz* is the measurement dimension, and *G* is the threshold that leads to a specified gating probability *PW*.

• Step 3: Association hypothesis events are formed by associating the validated measurements with the tracks.

The posterior probability of the association hypothesis event *θ<sup>h</sup>* is computed as

$$P(\theta\_h) = \mathbb{C}^{-1} \prod\_{t \in T\_0^h} \left( 1 - P\_D P\_W r\_{k-1}^t \right) \times \prod\_{t \in T\_1^h} \left( P\_D P\_W r\_{k-1}^t \frac{p\_h^t V\_k}{\hat{m}\_k} \right) \tag{8}$$

where *C* is the normalization constant; *p<sup>t</sup> <sup>h</sup>* = *<sup>f</sup> <sup>t</sup>* (*zt k Z*1:*k*−1)/*PW* with *PW* denoting the gating probability and *Z*1:*k*−<sup>1</sup> is the set of measurements accumulated up to time *k* − 1; *PD* is the probability of detection; *Vk* is the cluster volume; *m*ˆ *<sup>k</sup>* is the number of clutter measurement; *Th* <sup>0</sup> denotes the set of tracks associated with no measurement; and *<sup>T</sup><sup>h</sup>* <sup>1</sup> denotes the set of tracks associated with one measurement.

Using the posterior probability *P*(*θh*), the target existence probability of track *t* in *θ<sup>h</sup>* is computed as

$$r\_k^{t|h} = \left(1 - \sum\_{i=1}^{m\_k} \xi\_i^{t|h} \right) r\_{k,0}^{t|h} + \sum\_{i=1}^{m\_k} \xi\_i^{t|h} r\_{k,i}^{t|h} \tag{9}$$

where *r t*|*h <sup>k</sup>*,0 denotes the probability that track *t* is not associated with any measurement, *r t*|*h <sup>k</sup>*,*<sup>i</sup>* denotes the probability that track *<sup>t</sup>* is associated with measurement *<sup>z</sup><sup>i</sup> <sup>k</sup>* in *θh*, and the parameter *ξ t*|*h <sup>i</sup>* indicates whether the track *t* is associated with a measurement, i.e.,

$$\begin{array}{l} \mathcal{J}\_i^{t|h} = 1 \ \text{, if track } t \text{ is associated with measurement } \mathbf{z}\_k^i\\ \mathcal{J}\_i^{t|h} = 0 \ \text{, otherwise} \end{array} \tag{10}$$

The probabilities *r t*|*h <sup>k</sup>*,0 and *r t*|*h <sup>k</sup>*,*<sup>i</sup>* are computed as follows

$$r\_{k,0}^{t|h} = \frac{(1 - P\_D^t P\_W^t) r\_{k|k-1}^t}{1 - P\_D^t P\_W^t r\_{k|k-1}^t} P(\theta\_h), \ r\_{k,i}^{t|h} = P(\theta\_h) \tag{11}$$

Using the posterior probability *P*(*θh*) and the target existence probability *r t*|*h <sup>k</sup>* , the target state and the error covariance of track *t* in *θ<sup>h</sup>* are estimated as

$$\mathbf{x}\_{k}^{t|h} = \left(\mathbf{1} - \sum\_{i=1}^{m\_k} \boldsymbol{\xi}\_{i}^{t|h}\right) \mathbf{x}\_{k|k-1}^{t} + \sum\_{i=1}^{m\_k} \boldsymbol{\xi}\_{i}^{t|h} \mathbf{x}\_{k,i}^{t} \tag{12}$$

$$\mathbf{P}\_{k}^{t|h} = \left(1 - \sum\_{i=1}^{m\_k} \zeta\_i^{t|h}\right) \mathbf{P}\_{k|k-1}^t + \sum\_{i=1}^{m\_k} \zeta\_i^{t|h} \left[\mathbf{P}\_{k|k-1}^t - \mathbf{K}\_k^t \mathbf{S}\_k^t \left(\mathbf{K}\_k^t\right)^T\right] \tag{13}$$

where *x<sup>t</sup> <sup>k</sup>*,*<sup>i</sup>* is the updated target state using *<sup>z</sup><sup>i</sup> <sup>k</sup>* and *<sup>K</sup><sup>t</sup> <sup>k</sup>* is the filter gain.

• Step 4: The target existence probability, the target state, and the error covariance for track is approximated as [11]

$$r\_k^t = \sum\_{h=1}^{N\_{b\zeta}} r\_k^{t|h} \tag{14}$$

$$\mathbf{x}\_k^t = \frac{1}{r\_k^t} \sum\_{h=1}^{N\_{\mathcal{H}}} r\_k^{t|h} \mathbf{x}\_k^{t|h} \tag{15}$$

$$P\_k^t = \frac{1}{r\_k^t} \sum\_{h=1}^{N\_{k\ell}} r\_k^{t|h} \left[ \mathbf{P}\_k^{t|h} + \mathbf{x}\_k^{t|h} \left( \mathbf{x}\_k^{t|h} \right)^T - \mathbf{x}\_k^t \left( \mathbf{x}\_k^t \right)^T \right] \tag{16}$$

where *N*<sup>H</sup> denotes and the number of all association hypothesis events.

#### **4. Evolutionary-Optimization-Based Joint Integrated Probabilistic Data Association** *4.1. Motivations*

In automatic tracking systems, the target number may vary with time. In this paper, we assume that the targets are not discriminated from each other. Under these assumptions, the target states can be described as RFSs, in which the points are unordered and random [17]. It was shown that the posterior density can be changed into the other density within its RFS family [10]. Therefore, we propose to optimize the ordered posterior density within its RFS family to improve the tracking accuracy.

By reducing the error covariance, the problem of overlapping tracking gates can be effectively alleviated [11]. A scalar measure of the error covariance is the trace for the covariance matrix, as follows [11]

$$\begin{array}{ll} \mathsf{C}\_{k} &= \sum\_{t=1}^{n\_{k}} \mathsf{trace}(\mathbf{P}\_{k}^{t}) \\ &= \sum\_{t=1}^{n\_{k}} \sum\_{h=1}^{N\_{\mathcal{H}}} \frac{r\_{k}^{t|h}}{r\_{k}^{t}} \mathsf{trace} \left( \mathbf{P}\_{k}^{t|h} + \mathbf{x}\_{k}^{t|h} (\mathbf{x}\_{k}^{t|h})^{T} - \mathbf{x}\_{k}^{t} (\mathbf{x}\_{k}^{t})^{T} \right) \\ &= \sum\_{t=1}^{n\_{k}} \frac{1}{r\_{k}^{t}} \sum\_{h=1}^{N\_{\mathcal{H}}} r\_{k}^{t|h} \left[ \mathsf{trace} \left( \mathbf{P}\_{k}^{t|h} \right) + \left\| \mathbf{x}\_{k}^{t|h} - \mathbf{x}\_{k}^{t} \right\|^{2} \right] \end{array} \tag{17}$$

To alleviate the problem of overlapping tracking gates and improve the tracking accuracy, the trace for the covariance matrix Equation (17) is used as the cost function for the optimization of the posterior density.

#### *4.2. Evolutionary Optimization of the Posterior Density*

For the optimization of the posterior density, the target permutations under different association hypothesis events are reordered to minimize the cost function using Equation (17). We use the evolutionary algorithm [26,27] to implement the optimization. The specific steps of the proposed EOJIPDA method are as follows.

• Step 1: It is necessary to generate the initial population for the evolutionary optimization.

We enumerate all the possible permutations of the targets and assign a unique index to each permutation. For example, when the number of targets is *nk* = 3, the number of all possible permutations of the targets is *NP* = *nk*! = 3! = 6. The set of the indexes for these permutations is denoted as

$$I = \left\{ I\_1, I\_{1\prime}, \dots, I\_{n\_k!} \right\} \tag{18}$$

The indices within the set *I* are randomly assigned to each association hypothesis event. An assignment indicates a possible solution for the evolutionary algorithm. Such a solution is named as a chromosome, as follows

$$
clorumosome = [\mathbb{C}\_1, \mathbb{C}\_2, \dots, \mathbb{C}\_{N\_{\mathbb{M}}}] \tag{19}$$

where the element *Cj* ∈ *<sup>I</sup>* and *<sup>j</sup>* = 1, 2, ... , *<sup>N</sup>*H. We randomly perform the assignments and generate *Npop* chromosomes. These chromosomes are used as the initial population, i.e., the first generation, for the evolutionary optimization.

• Step 2: The value of the cost function forms the fitness value for each solution.

The fitness function is used to test how well the chromosome solves the problem. All the solutions in the population are sorted based on their fitness values.

• Step 3: The offspring solutions are generated using the selection, crossover and mutation operators.

The selection operator chooses some of the chromosomes for reproduction. For the crossover operator, we use the single-point crossover. The genes of the two-parent chromosome are interchanged after a random selected point. Along with the crossover, the mutation is performed. We random select one point and change it into a random selected point from the set using Equation (18).

The new population with all the parents and offspring is sorted again based on their fitness values and the population size is decreased to *Npop* by eliminating all the lower rank solutions.

• Step 4: Back to Step 3 and the cycle repeats.

For each generation, we record the chromosome with the highest fitness (along with the value of the fitness). The cycle is repeated until the fitness value of the "best-so-far" chromosome stabilizes.

#### *4.3. Illustrative Example*

To illustrate the procedure of the proposed EOJIPDA method, a one-dimensional example is used in this subsection. We assume that two Gaussian distributed targets generate two detections and there are two association hypothesis events *θ*<sup>1</sup> and *θ*2. The density of target *<sup>t</sup>* corresponding to *<sup>θ</sup><sup>h</sup>* is denoted as *<sup>p</sup>t*|*h*(*x*) = ) *<sup>r</sup>t*|*h*, <sup>N</sup> (*xt*|*h*, *<sup>P</sup>t*|*h*) \* . The initial posterior probability density is given as

$$\begin{aligned} \theta\_1: p^{1|1}(\mathbf{x}) &= \{0.1, \ N(1.0, 1.0)\}, \ p^{2|1}(\mathbf{x}) = \{0.4, \ N(5.0, 1.0)\} \\ \theta\_2: p^{1|2}(\mathbf{x}) &= \{0.5, \ N(4.0, 1.0)\}, \ p^{2|2}(\mathbf{x}) = \{0.5, \ N(2.0, 1.0)\} \end{aligned} \tag{20}$$

Using Equations (14)–(16), the approximated probability density for each track is described by

$$p^1(\mathbf{x}) = \{0.60, \ N(3.50, 2.25)\}, \ p^2(\mathbf{x}) = \{0.9, \ N(3.33, 3.22)\} \tag{21}$$

The posterior probability density of the joint track state and its approximated probability density of the estimated track state are shown in Figures 1a and 1b, respectively. Obviously, Figure 1b is very different with Figure 1a. Therefore, the approximation is less accurate and the trace of the covariance matrix in Equation (17) of the estimated covariance matrix is *C*<sup>0</sup> = 2.25 + 3.22 = 5.47.

To improve the accuracy of the approximation, we perform the evolutionary optimization of the posterior probability density. Five iterations are performed and the resulting posterior probability density is given as

$$\begin{aligned} \theta\_1: p^{1|1}(\mathbf{x}) &= \{0.1, \ N(1.0, 1.0)\}, \ p^{2|1}(\mathbf{x}) = \{0.4, \ N(5.0, 1.0)\} \\ \theta\_2: p^{1|2}(\mathbf{x}) &= \{0.5, \ N(2.0, 1.0)\}, \ p^{2|2}(\mathbf{x}) = \{0.5, \ N(4.0, 1.0)\} \end{aligned} \tag{22}$$

Using Equations (14)–(16), the approximated probability density for each track is described by

$$p^1(\mathbf{x}) = \{0.60, \ N(1.83, 1.14)\}, \ p^2(\mathbf{x}) = \{0.9, \ N(4.44, 1.25)\}\tag{23}$$

The resulting posterior probability density and its approximated probability density are plotted in Figure 2a,b, respectively. It is obvious that the similarity between Figure 2a,b is much higher than that between Figure 1a,b. The trace of the covariance matrix in Equation (17) of the estimated covariance matrix is reduced to *C*<sup>0</sup> = 1.14 + 1.25 = 2.39. This indicates that the accuracy of the approximation and the state estimation is greatly improved. The results show a good agreement with the theoretical analysis.

**Figure 2.** The posterior probability density and its approximated probability density of iteration 5: (**a**) posterior probability density of iteration 5 (**b**) approximated probability density of iteration 5.

#### **5. Numerical Simulation and Results**

We demonstrate the performance of the proposed EOJIPDA filter using two challenging MTT scenarios in the two-dimensional surveillance area. At time *k*, the state of the moving target is represented by the state vector

$$\mathbf{x}\_{k} = \begin{bmatrix} \mathbf{x}\_{k'} \dot{\mathbf{x}}\_{k'} \mathbf{y}\_{k'} \dot{\mathbf{y}}\_{k} \end{bmatrix}^{\mathrm{T}} \tag{24}$$

where [*xk*, *yk*] <sup>T</sup> is the target position and [ . *xk*, . *yk*] <sup>T</sup> is the target velocity. The transition matrix of the motion model in Equation (1) is

$$F = I\_2 \otimes \begin{bmatrix} 1 & T\_\sf s \\ 0 & 1 \end{bmatrix} \tag{25}$$

where ⊗ is the Kroneker product, *I*<sup>2</sup> is a 2 × 2 identity matrix, and *T*<sup>s</sup> is the sampling interval. The covariance of the Gaussian process noise *u <sup>k</sup>* is

$$\mathbf{Q} = I\_2 \otimes q \begin{bmatrix} \frac{T\_s^3}{2} & \frac{T\_s^2}{2} \\ \frac{T\_s^2}{2} & T\_s \end{bmatrix} \tag{26}$$

where *q* is a tuning parameter. In both scenarios, a sensor is used to collect the measurements in the surveillance area. For simplicity, it is assumed that the measurements collected by the sensor represent the positions of the targets. Therefore, the observation matrix of the measurement model in Equation (3) is

$$H = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \tag{27}$$

The covariance of the Gaussian measurement noise *w<sup>k</sup>* is

$$\mathbf{R} = \begin{bmatrix} \sigma\_x^2 & 0\\ 0 & \sigma\_y^2 \end{bmatrix} \tag{28}$$

where *σ<sup>x</sup>* and *σ<sup>y</sup>* are the standard deviations of the measurement noise in *x* and *y* coordinates, respectively.

The sampling interval of the sensor is fixed as *T*<sup>s</sup> = 1 s and the detection probability is *PD* = 0.92. The gating threshold is *G* = 9.21 [28], which corresponds to a twodimensional gating probability of *PW* = 0.99. The coefficients of the Markov chain one model in Equation (5) are borrowed from [9], where *p*<sup>11</sup> = 0.98 and *p*<sup>21</sup> = 0. The tracks are terminated if the target existence probability falls below the termination threshold and are confirmed if the probability exceeds the confirmation threshold. The confirmation

and termination thresholds are *PC* = 0.83 and *PT* = 0.0909 [9], respectively. The tracking performance of the filter is measured by the optimal sub-pattern assignment (OSPA) [29] multi-target miss distance, which measures the error between the estimated and true state and is widely used [17–20,30]. If *X* = {*x*1, *x*2, ··· , *xn*}, *Y* = {*y*1, *y*2, ··· , *ym*}, and *n* ≥ *m*, the OSPA multi-target miss distance is defined as [29]

$$d\_{p, \boldsymbol{\zeta}}^{\text{OSPA}}(X, Y) := \left( \frac{1}{n} \left( \min\_{\pi \in \Pi\_n} \sum\_{i=1}^{m} d^{(\varepsilon)}(\mathbf{x}\_i, \mathbf{y}\_{\pi(i)})^p + \mathbf{c}^p(n - m) \right) \right)^{1/p} \tag{29}$$

where *<sup>d</sup>*(*c*)(*x*, *<sup>y</sup>*) :<sup>=</sup> min(*c*, *x* − *y*-), *c* > 0 is the cut-off parameter, *p* ≥ 1 is an order parameter and *π* is a permutation function in the set of permutations Π*n*. If *n* < *m*, *d*OSPA *<sup>p</sup>*,*<sup>c</sup>* (*X*,*Y*) := *d*OSPA *<sup>p</sup>*,*<sup>c</sup>* (*Y*, *X*). All the experiments are tested in MATLAB R2010a and implemented on a computer with a 3.40 GHz processor.

#### *5.1. Scenario 1*

We consider the tracking of two targets in this scenario, whose trajectories are shown in Figure 3. The track coalescence of the JPDA and JIPDA filters can be clearly observed using this scenario; hence, this scenario has been widely used [10–15]. In Figure 3, the parameters are selected as: *φ* = *π*/3, *d* = 0.5 m, and *l*<sup>1</sup> = *l*<sup>2</sup> = 10 m. It is assumed that the speed of the target stays constant as *v* = 1 m/s. The parameters for the clutter measurements are *<sup>λ</sup>* <sup>=</sup> 0.36 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>m</sup>−<sup>2</sup> and FOV <sup>=</sup> 1.4 <sup>×</sup> <sup>10</sup><sup>3</sup> m2 , giving an average of five clutter returns per scan. The standard deviations of the measurement noise are *σ<sup>x</sup>* = *σ<sup>y</sup>* = 0.1 m .

**Figure 3.** Simulated scenario 1.

The estimated target positions for a single run of the JIPDA and EOJIPDA filters are shown in Figures 4a and 4b, respectively. It can be observed that the serious track coalescence problem happens in Figure 4a. When the targets are closely spaced, their position estimates coalesce. Although the targets move away from each other later, the JIPDA fails to detect their separation. Compared with the JIPDA filter, the EOJIPDA filter can accurately estimate the target positions, as shown in Figure 4b.

The average OSPA distances over 100 Monte Carlo trails for the JIPDA, exact nearest neighbor JIPDA (ENNJIPDA), JIPDA\*, and EOJIPDA filters are shown in Figure 5. Due to the track coalescence problem, the OSPA distance of the JIPDA filter is worse than in the comparative algorithms. Although the ENNJPDA filter was proposed to overcome the coalescence problem, it is sensitive to clutter and missed detections. Compared with the JIPDA\* filter, the EOJIPDA filter uses more information about the posterior density and performs better in terms of the OSPA error.

**Figure 4.** Estimated target positions for Scenario 1, in which circles represent estimated positions and dotted lines represent true trajectories: (**a**) Output of the JIPDA filter; (**b**) output of the EOJIPDA filter.

**Figure 5.** Average OSPA distances for Scenario 1. JIPDA\* denotes a coalescence-avoiding version of JIPDA, as shown in [16].

#### *5.2. Scenario 2*

The tracking of three targets is studied in this scenario. The trajectories of the targets cross each other at a crossing point, as shown in Figure 6. Each target keeps a constant speed, *v* = 1 m/s. The parameters used in the simulations are *φ*<sup>1</sup> = *φ*<sup>2</sup> = *π*/8 and *<sup>l</sup>*<sup>1</sup> <sup>=</sup> *<sup>l</sup>*<sup>2</sup> <sup>=</sup> <sup>100</sup> m, *<sup>σ</sup><sup>x</sup>* <sup>=</sup> *<sup>σ</sup><sup>y</sup>* <sup>=</sup> <sup>10</sup> m, *<sup>λ</sup>* <sup>=</sup> 2.5 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>m</sup>−<sup>2</sup> , and FOV <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup><sup>2</sup> . The two-point difference [24] is used to initiate tracks automatically and the maximum possible speed is assumed to be 200 m/s.

**Figure 6.** Simulated scenario 2.

The estimated target positions for a single run of the JIPDA and EOJIPDA filters are shown in Figures 7a and 7b, respectively. It can be observed from Figure 7 that it takes some time for the JIPDA and EOJIPDA filters to localize targets. When the targets are spaced close to each other, the JIPDA filter cannot accurately estimate the target positions and the tracks tend to coalesce. Compared with the JIPDA filter, the EOJIPDA filter can effectively overcome the track coalescence problem and provide position estimates that are close to the true values.

**Figure 7.** Estimated target positions for Scenario 2, in which circles represent estimated positions and dotted lines represent true trajectories: (**a**) Output of the JIPDA filter; (**b**) output of the EOJIPDA filter.

The average OSPA distances for the different tracking methods are shown in Figure 8. The tracking errors of all the methods are high at the initial time steps because the initial states of the targets are assumed to be unknown in this scenario. Since the track coalescence problem happens in the JIPDA filter, it can be observed that the tracking error of the JIPDA filter is high after the crossing point (at time *k* = 11). The ENNJIPDA and JIPDA\* filters can effectively reduce the tracking error of the JIPDA filter. Nevertheless, the error performance of the proposed EOJIPDA filter is better than that of other filters for almost the entire scenario.

Table 2 shows the accumulated number of confirmed tracks at the end point of the tracking by 100 MC trails for the JIPDA filter, the ENNJIPDA filter, the JIPDA\* filter and the EOJIPDA filter. The results for the different clutter rates are provided. It can be seen that the challenging tracking environment makes all the filters lose some tracks. Nevertheless, the performance of the EOJIPDA filter is superior to others under different clutter rates. In other words, the better tracking performance of the EOJIPDA filter makes the detection of target existence more reliable.

**Figure 8.** Average OSPA distances for Scenario 2.

**Table 2.** Number of confirmed targets.


#### **6. Conclusions**

The JIPDA filter an effective method for MTT, but it suffers from the serious track coalescence problem. To improve the tracking performance of the JIPDA filter, a novel MTT filter, named EOJIPDA, was proposed in this paper. We first attempted to use the evolutionary computation technique to optimize the posterior density of the JIPDA filter. When the target identity was irrelevant, we modeled the posterior density optimization problem as an evolutionary computation problem and the trace of the covariance matrix is used as the cost function. Through an indicative example, it was shown that the evolutionary computation can effectively reduce the value of the cost function, improving the accuracy of the Gaussian approximation and the state estimation.

The simulation results show that the EOJIPDA filter effectively avoids track coalescence and performs better than the traditional algorithms in terms of the OSPA error. Future work will investigate the application of more computational intelligence strategies, such as particle swarm optimization and the Cuckoo search algorithm, to improving the performances of traditional tracking methods.

**Author Contributions:** Conceptualization, S.L. and Y.Z.; methodology, S.L.; validation, Y.Z.; investigation, S.L.; writing—original draft preparation, S.L.; writing—review and editing, S.L., H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** In this work, we have used the free RFS MATLAB code at http://batuong.vo-au.com/codes.html.

**Acknowledgments:** This work was supported in part by the National Natural Science Foundation of China, under grant numbers 62007022 and 61906146; the Natural Science Foundation of Shaanxi Province, under grant number 2021JQ-209; and the Fundamental Research Funds for the Central Universities, under grant numbers GK202103082 and JB210210.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

