Target dynamics are modeled under the nearly constant velocity assumption in this paper. The state vector of target t at discrete time step k is represented by 𝐱t,k, with following state transition equation:

x k + 1 = F ⁢ x k + u k

where F denotes the state transition matrix, and uk∼𝒩⁢(0,Q) is Gaussian process noise with zero mean and covariance Q. The measurement zk associated with target t follows a linear observation model:

z k = H ⁢ x k + w k

where H is the observation matrix, and wk∼𝒩⁢(0,R) represents zero-mean Gaussian measurement noise with covariance R.

The target-originated detections occur with probability PD, while clutter measurements follow a spatial Poisson distribution. The average number of clutter measurements per time step is r=λ⁢|FoV|, where λ is the clutter intensity and |FoV| denotes the sensor's field-of-view coverage. The validated measurement set at time k is defined as:

Z k = { z k 1 , z k 2 , ⋯ , z k m k }

where mk is the total number of measurements. The historical measurement sequence up to time k is denoted as

Z 1 : k = { Z 1 , Z 2 , ⋯ , Z k }

JIPDA generates association hypotheses by probabilistically assigning validated measurements to tracks. Each monitored target maintains a dynamic parameter tuple encapsulating its kinematic state vector xt,k, covariance matrix Pt,k, and existence probability Ps⁢(t,k) at discrete time index k. Consequently, the characterization of a track t can be achieved through the parameter set {Ps⁢(t,k),𝒩⁢(xt,k,Pt,k)}. When advancing to time step k+1, the key steps for updating this parameter set are described below.

First, in the prediction stage, the Kalman filter is used to compute the predicted target state xt,k+1|k and covariance Pt,k+1|k as follows:

𝒙 t , k + l | k = 𝑭 ⁢ 𝒙 t , k 𝑷 t , k + l | k = 𝑭 ⁢ 𝑷 t , k ⁢ 𝑭 T + 𝑸 k

where F denotes the state transition matrix, and Qk denotes process noise covariance matrix at time k.

The target existence probability propagates according to a first-order Markov chain model,

P s ⁢ ( t , k + 1 ∣ k ) = p 11 ⁢ P s ⁢ ( t , k ) + p 21 ⁢ ( 1 − P s ⁢ ( t , k ) )

where P11 and P21 are the Markov chain coefficients.

Second, analytically derived tracking gates are generated for each target to assist measurement association. As comprehensively discussed in [23], various gating methodologies exist for measurement validation, with the ellipsoidal gating approach [24] being particularly prevalent due to its optimal statistical properties under Gaussian assumptions. At time step k, for track t, the ellipsoidal gate is defined as

𝒢 t k + 1 = { 𝒛 k + 1 i ∈ ℝ d z | Δ t , j k + 1 ≤ γ }

Δ t , j k + 1 = ( 𝒛 k + 1 j − 𝑯 k ⁢ 𝒙 t , k + 1 | k ) T × ( 𝑺 t , k + 1 ) − 1 ⁢ ( 𝒛 k + 1 j − 𝑯 k ⁢ 𝒙 t , k + 1 | k )

where zk+1i denotes the ith received measurement, Si,k+1=H⁢Pi,k+1|k⁢H+R is the innovation covariance, γ is a threshold, which leads to a gating probability Pw, and dz is the measurement dimension. This gating procedure ensures that only statistically compatible measurements are considered for potential association, thereby reducing computational complexity through effective measurement pruning.

Third, association hypotheses are generated by establishing probabilistic correspondences between validated measurements and existing tracks. The posterior probability of an event θh is computed as

P ⁢ ( θ h ) = C − 1 ⁢ ∏ t ∈ T 0 h ( 1 − P D ⁢ P W ⁢ P s ⁢ ( t , k + 1 ∣ k ) ) × ∏ t ∈ T 1 h ( P D P W P s ( t , k + 1 ∣ k ) p ⁢ ( z k + 1 j ∣ θ h , Z 1 : k ) λ )

where C is the normalization constant, PD denotes probability of detection, T0h denotes the set of tracks associated without measurement, and T1h corresponds to the set of tracks associated with exactly one measurement.

The existence probability of track t is updated under the association hypothesis event θh as

P s ⁢ ( t ∣ h , k + 1 ) = ( 1 − ∑ i = 1 m k ξ i t ∣ h ) ⁢ P s , 0 ⁢ ( t ∣ h , k + 1 ) + ∑ i = 1 m k ξ i t ∣ h ⁢ P s , i ⁢ ( t ∣ h , k + 1 )

where ζit|h is a binary indicator variable, i.e.,

{ ξ i t | h = 1 if tracktis associated with measurement z k i ξ i t | h = 0 otherwise

and Ps,0⁢(t|h,k+1) denotes the probability of track t being unassociated with any measurement and Ps,i⁢(t|h,k+1) denotes the probability of track t being associated with measurement zk+1i, which are computed as

P s , 0 ⁢ ( t ∣ h , k + 1 ) = ( 1 − P D t ⁢ P W t ) ⁢ P s ⁢ ( t , k + 1 ∣ k ) 1 − P D t ⁢ P W t ⁢ P s ⁢ ( t , k + 1 ∣ k ) ⁢ P ⁢ ( θ h ) P s , i ⁢ ( t ∣ h , k + 1 ) = P ⁢ ( θ h )

The state of track t under hypothesis θh is computed as:

𝒙 t | h , k = ( 1 − ∑ i = 1 m k ξ i t | h ) ⁢ 𝒙 t , k | k − 1 + ∑ i = 1 m k ξ i t | h ⁢ 𝒙 t , k i

where xt,ki is the state updated with measurement zki. The error covariance is updated as:

𝑷 t | h , k + 1 = ( 1 − ∑ i = 1 m k ξ i t | h ) ⁢ 𝑷 t , k + 1 | k + ∑ i = 1 m k ξ i t | h ⁢ [ 𝑷 t , k + 1 | k − 𝑲 t , k + 1 ⁢ 𝑺 t , k + 1 ⁢ ( 𝑲 t , k + 1 ) T ]

where Kt,k+1=Pt,k+1|k⁢HT⁢(St,k+1)−1 denotes the Kalman filter gain.

Ultimately, the track's existence probability, state estimate and covariance are approximated by

P s ⁢ ( t , k + 1 ) = ∑ h = 1 N ℋ P s ⁢ ( t ∣ h , k + 1 )

𝒙 t , k = 1 P s ⁢ ( t , k + 1 ) ⁢ ∑ h = 1 N ℋ P s ⁢ ( t ∣ h , k + 1 ) ⁢ 𝒙 t ∣ h , k

𝑷 t , k = 1 P s ⁢ ( t , k + 1 ) ⁢ ∑ h = 1 N ℋ P s ⁢ ( t ∣ h , k + 1 ) × [ 𝑷 t ∣ h , k + 1 + 𝒙 t ∣ h , k ⁢ ( 𝒙 t ∣ h , k ) T − 𝒙 t , k ⁢ ( 𝒙 t , k ) T ]

where Nℋ denotes the number of all association hypothesis events.

Automatic target tracking systems account for time-varying number of targets. Following [16], we employ random finite sets (RFS) to characterize the inherently stochastic and disordered nature of target states. Following the work in [10], which demonstrates that posterior probability densities can be transformed within their corresponding RFS families, this paper introduces an innovative approach. Specifically, it proposes a method to improve the multi-target tracking performance by optimizing ordered posterior density functions within the RFS family framework.

Theoretical analysis and experimental studies [11] have shown that minimizing the error covariance matrix can overcome the overlapping tracking gate problem. In practical applications, we employ the covariance matrix trace [10] as the optimization criterion for posterior density refinement, formulated as

C k = ∑ t = 1 n k trace ⁢ ( P t , k )

where nk denotes the number of targets and the function trace⁢(⋅) represents the trace of the matrix.

For posterior density optimization, target permutations across distinct association hypotheses are reconfigured to minimize the cost function defined in (1). The optimization is effectuated by means of PSO [5], with the principal procedures of the PSO-JIPDA approach as follows.

Step 1: PSO setup for the optimization of the posterior density.

For PSO implementation, we adopt (1) as the global fitness function. The solution space comprises every permutation of the current target indices and is defined as:

𝒫 = { Π h ∣ h = 1 , … ⁢ N ℋ }

where Πh denotes the permutation set corresponding to the hth association hypothesis and Πh has cardinality |Πh|=nk!. We initialize a PSO swarm with Np particles. At time k, the position xpk of each particle p represents a candidate solution vector, where each element corresponds to the index (from 1 to nk) of a specific target permutation chosen for a particular hypothesis. The goal is to find the vector of permutation indices (one index per hypothesis) that optimizes a global fitness function.

Step 2: PSO Iterative Optimization.

Execute the PSO algorithm for defined generations and update the inertia weight w at time k. In each generation, for a particle, we perform the following steps.

First, compute the inertia weight. The original linear decreasing inertia weight strategy was enhanced with a nonlinear adaptive approach. The dynamic inertia weight is calculated based on the iteration progress and fitness values, as follows

w = { w min + ( w max − w min ) × ( f current − f min ) ( f mean − f min ) if f current < f mean w max otherwise

where wmax and wmin represent the maximum and minimum inertia weights, respectively; f⁢(⋅) denotes fitness function; fcurrent, fmean, and fmin denote the current fitness function, the mean fitness function and the minimum fitness function, respectively.

Figure 1.

Flowchart of PSO-JIPDA algorithm.

To illustrate the PSO-JIPDA algorithm's operational process, we consider a simplified one-dimensional case with two Gaussian-distributed targets generating two measurements. Assume two association hypothesis events ψ1 and ψ2 are established, with density of target t under hypothesis ψh is modeled as pt|h⁢(x)={rt|h,𝒩⁢(xt|h,pt|h)}. The initial posterior density is specified by

ψ 1 : { p 1 | 1 ⁢ ( x ) = { 0.1 , 𝒩 ⁢ ( 1 , 1 ) } p 2 | 1 ⁢ ( x ) = { 0.4 , 𝒩 ⁢ ( 4 , 1 ) } ψ 2 : { p 1 | 2 ⁢ ( x ) = { 0.5 , 𝒩 ⁢ ( 5 , 1 ) } p 2 | 2 ⁢ ( x ) = { 0.5 , 𝒩 ⁢ ( 2.3 , 1 ) }

Based on (15)-(17), the Gaussian-approximated probability density for each target is characterized as

p 1 ⁢ ( x ) = { 0.6 , N ⁢ ( 3.5 , 2.25 ) } p 2 ⁢ ( x ) = { 0.9 , N ⁢ ( 3.5 , 2.80 ) }

Figure 2(a) and 2(b) present the posterior probability density of the joint track state and its Gaussian approximation, respectively. It can be observed that a considerable difference exists between these two densities, demonstrating the limited precision of the Gaussian approximation. In this case, the trace of the covariance matrix is

C 0 = 2.25 + 2.80 = 5.05

which can be computed using (1).

(a) Posterior and (b) approximated probability densities at iteration 0.

The posterior probability density is then optimized using the particle swarm optimization (PSO) algorithm. Following five iterations, the final optimized result is derived as follows:

ψ 1 : { p 11 ⁢ ( x ) = { 0.4 , N ⁢ ( 1 , 1 ) } p 2 | 1 ⁢ ( x ) = { 0.1 , N ⁢ ( 5 , 1 ) } derived as follows ψ 2 : { p 1 | 2 ⁢ ( x ) = { 0.5 , N ⁢ ( 2.3 , 1 ) } p 2 | 2 ⁢ ( x ) = { 0.5 , N ⁢ ( 4 , 1 ) }

Using (15)-(17), the approximated probability density for each track is given as

p 1 ⁢ ( x ) = { 0.60 , N ⁢ ( 4.44 , 1.25 ) } p 2 ⁢ ( x ) = { 0.90 , N ⁢ ( 2.08 , 1.23 ) }

Figure 3(a) displays the initial posterior probability density, while its approximated single Gaussian density is presented in Figure 3(b). Obviously, Figures 3(a) and 3(b) demonstrate a strong correlation, whereas Figures 2(a) and 2(b) exhibit a notable divergence. In this case, the trace of the covariance matrix is

C 0 = 1.25 + 1.23 = 2.48

After PSO optimization, the trace decreases by 50.89% relative to the initial trace (9), demonstrating a significant improvement in Gaussian approximation and state estimation precision. These statistical analyses exhibit strong concordance with theoretical expectations.

(a) Posterior and (b) approximated probability densities at iteration 5.