Optimized Sage–Husa filtering method for multitarget tracking with unknown detection profile

Conventional beta distribution

Suppose that α and β are the detection number and the missed detection number for the single target, the detection probability p D can be defined by the beta distributio B ( α , β , p D )

B ( α , β , p D ) = p D α − 1 ( 1 − p D ) β − 1 ∫ 0 1 p D α − 1 ( 1 − p D ) β − 1 d p D

(1)

Note that p D is implicitly defined by the integral calculation. Considering the discreteness of the surveillance time, we approximate p D using the mean of the beta distribution

p D = . α α + β

(2)

At time k − 1 , suppose that α k − 1 i and β k − 1 i are the cumulative number of detection and the cumulative number of missed detection for the target i, the variance of the beta distribution can be written as

( σ k − 1 i ) 2 = α k − 1 i β k − 1 i ( μ ( α k − 1 i + β k − 1 i ) + 1 ) ( α k − 1 i + β k − 1 i ) 2 ( α k − 1 i + β k − 1 i + 1 )

(3)

where μ ∈ [ 0 , 1 ) is the adjustment parameter. Note that ( σ k − 1 i ) 2 keeps unchanged when μ is 0. Otherwise, ( σ k − 1 i ) 2 has the maximum when μ approximates 1. Base on the characteristics of μ , we have the predicted equations of α k | k − 1 i

α k | k − 1 i = α k − 1 i α k − 1 i + β k − 1 i ( α k − 1 i β k − 1 i ( α k − 1 i + β k − 1 i ) 2 ( σ k − 1 i ) 2 − 1 ) = α k − 1 i α k − 1 i + β k − 1 i ( α k − 1 i β k − 1 i ( α k − 1 i + β k − 1 i ) 2 ( α k − 1 i + β k − 1 i ) 2 ( α k − 1 i + β k − 1 i + 1 ) α k − 1 i β k − 1 i ( μ ( α k − 1 i + β k − 1 i ) + 1 ) − 1 ) = α k − 1 i α k − 1 i + β k − 1 i ( α k − 1 i + β k − 1 i + 1 ) − ( μ ( α k − 1 i + β k − 1 i ) + 1 ) μ ( α k − 1 i + β k − 1 i ) + 1 = α k − 1 i α k − 1 i + β k − 1 i ( 1 − μ ) ( α k − 1 i + β k − 1 i ) μ ( α k − 1 i + β k − 1 i ) + 1 = ( 1 − μ ) α k − 1 i μ ( α k − 1 i + β k − 1 i ) + 1

(4)

Similarly, β k | k − 1 i can be derived as

β k | k − 1 i = β k − 1 i α k − 1 i + β k − 1 i ( α k − 1 i β k − 1 i ( α k − 1 i + β k − 1 i ) 2 ( σ k − 1 i ) 2 − 1 ) = ( 1 − μ ) β k − 1 i μ ( α k − 1 i + β k − 1 i ) + 1

(5)

At time k, the detection probability for the target i is updated as

p D , k i = . α k | k − 1 i α k | k − 1 i + β k | k − 1 i

(6)

Note that α k − 1 i and β k − 1 i are often more than 1 in the beta distribution. Due to the influence of prediction, α k − 1 i and β k − 1 i are less than 1 when μ approximates 1. Although the mean of the detection probability can approximate α / ( α + β ) , the short surveillance time still leads to the inaccurate estimation. Moreover, the estimated value becomes more inaccurate.

Improved estimation of detection probability

With respect to both detection and missed detection at time k, the averaged probability of detection for all targets can be updated as

p D , k = 1 n ∑ i = 1 n p D , k − 1 i ( α k | k − 1 i − β k | k − 1 i ) + β k | k − 1 i + 1 α k | k − 1 i + β k | k − 1 i + 1

(7)

Proof

At time k, the detection probability can be updated by two effective transfer paths, that is, two solid lines in Figure 1.

OPEN IN VIEWER

Figure 1.

Relationship between predicted probability and updated probability.

Note that there are two transfer probabilities from both p D , k − 1 i and ( 1 − p D , k − 1 i ) to the same p D , k i . Furthermore, ( α k | k − 1 i + 1 ) and ( β k | k − 1 i + 1 ) represent the transferring statuses of the detection and missed detection, respectively. When the cumulative number is equal to ( α k | k − 1 i + β k | k − 1 i + 1 ) , we have in hand the transfer probability:

p D , k | k − 1 i = 1 α k | k − 1 i + β k | k − 1 i + 1 × { ( α k | k − 1 i + 1 ) , detection ( β k | k − 1 i + 1 ) , missed detection

(8)

where the denominator ( α k | k − 1 i + β k | k − 1 i + 1 ) cannot be 0 when both α k − 1 i and β k − 1 i become 0. It ensures p D , k | k − 1 i can be effectively computed, which is equal to 1 simultaneously. According to the Bayesian rule, the updated probability of detection is considered as the product of the prior detection probability and its transfer probability. Subsequently, the averaged probability of detection for all targets can be written as

p D , k = 1 n ∑ i = 1 n p D , k i = 1 n ∑ i = 1 n ( p D , k − 1 i ( α k | k − 1 i + 1 ) α k | k − 1 i + β k | k − 1 i + 1 + ( 1 − p D , k − 1 i ) ( β k | k − 1 i + 1 ) α k | k − 1 i + β k | k − 1 i + 1 ) = 1 n ∑ i = 1 n p D , k − 1 i ( α k | k − 1 i + 1 ) + β k | k − 1 i + 1 α k | k − 1 i + β k | k − 1 i + 1 = 1 n ∑ i = 1 n p D , k − 1 i ( α k | k − 1 i − β k | k − 1 i ) + β k | k − 1 i + 1 α k | k − 1 i + β k | k − 1 i + 1

(9)

Obviously, the factor p D , k − 1 i ( α k | k − 1 i − β k | k − 1 i ) + β k | k − 1 i + 1 / α k | k − 1 i + β k | k − 1 i + 1 becomes 0 when β k | k − 1 i = 0 and p D , k − 1 i = 1 , so the corresponding detection probability can reach the maximal value 1. Note that p D , k approximates the unknown detection probability using the mean. Compared with the conventional beta distribution, its advantage is that the transfer probability in the space of the probability distribution can be applied. However, p D , k is determined by the transfer relationship, which is suitable for the unknown detection profile.

Model of target state and prediction

At time k − 1 , the state equation and measurement equation of the target i can be given by

x k | k − 1 i = F k | k − 1 x k − 1 i + u k − 1

(10)

z k i = H k x k | k − 1 i + w k

(11)

where F k | k − 1 is the transfer matrix of target state, u k − 1 is the matrix of process noise with the Gaussian distribution ( 0 , Q k − 1 ) , H k is the transfer matrix of measurement, and w k is the matrix of measurement noise with the Gaussian distribution ( 0 , R k ) .

Using the forgetting factor b , we have in hand the tune parameter d k

d k = 1 − b 1 − b k

(12)

where the variance Q k of u k − 1 and the variance R k of v k are given by

Q k = ( 1 − d k − 1 ) Q k − 1 + d k ( K k − 1 v k v k T ( K k − 1 ) T + P k i − F k | k − 1 P k − 1 i F k | k − 1 T )

(13)

R k = ( 1 − d k ) R k − 1 + d k − 1 ( ( I − H k K k − 1 ) v k v k T ( I − H k K k − 1 ) T + H k P k | k − 1 i H k T )

(14)

where ( · ) T is the transpose matrix, I is the unit square matrix, and the estimation covariance P k | k − 1 i is predicted as

P k | k − 1 i = F k | k − 1 P k − 1 i F k | k − 1 T + Q k

(15)

In the conventional Sage–Husa filter, the estimated R k is easy to lose the positive definiteness when the estimation accuracy of noise gradually decreases. Also, d k i becomes small and approaches to a constant value. In this event, we cannot estimate the noise with the time-varying R k when the rounding error occurs. As a result, the predicted P k | k − 1 i and updated P k i lose the positive definiteness, and then the Kalman filtering gain K k starts to diverge.

Considering the positive definiteness of P k − 1 i at time k − 1 , we have in hand the Cholesky decomposition

P k − 1 i = L k − 1 i Λ k − 1 i ( L k − 1 i ) T = ( L k − 1 i Λ k − 1 i ) ( L k − 1 i Λ k − 1 i ) T = Γ k − 1 i ( Γ k − 1 i ) T

(16)

where L k − 1 i is the lower triangular matrix, Λ k − 1 i is the diagonal matrix whose elements on the main diagonal are positive, and Γ k − 1 i = L k − 1 i Λ k − 1 i is the lower triangular matrix. As we know, the Cholesky decomposition of P k − 1 i is unique when the elements on the main diagonal are all positive. Suppose that r, s and t ( s ≤ r ) are the subscripts of the element γ , and the decomposition for Γ = ( γ st ) is given by

{ γ rr = P rr − ∑ t = 1 r − 1 ( γ rt ) 2 γ rt = 1 γ rr ( P rr − ∑ s = 1 t − 1 γ rs γ st )

(17)

Obviously, the element P rr = ∑ t = 1 r ( γ rt ) 2 meets the condition γ rr ≤ P rr . Given that the Cholesky decomposition is reversible and the element γ rr is controllable, the Cholesky decomposition keeps stable.

Suppose that the lower triangular matrix is predicted as Γ k | k − 1 i = ( F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + Q k T ) − 1 ( F k | k − 1 P k − 1 i F k | k − 1 T + Q k ) , the estimation covariance can be predicted as

P k | k − 1 i = Γ k | k − 1 i ( Γ k | k − 1 i ) T

(18)

Proof

We use the undetermined coefficients ε k − 1 i and δ k − 1 i to represent P k | k − 1 i under the Cholesky decomposition

P k | k − 1 i = F k | k − 1 P k − 1 i F k | k − 1 T + Q k = ( ε k − 1 i F k | k − 1 P k − 1 i F k | k − 1 T + δ k − 1 i Q k ) ( ε k − 1 i F k | k − 1 P k − 1 i F k | k − 1 T + δ k − 1 i Q k ) T = F k | k − 1 P k − 1 i F k | k − 1 T ( ( ε k − 1 i ) 2 F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + ε k − 1 i δ k − 1 i Q k T ) + Q k ( ε k − 1 i δ k − 1 i F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + ( δ k − 1 i ) 2 Q k T )

(19)

There are the following equations

{ ( ε k − 1 i ) 2 F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + ε k − 1 i δ k − 1 i Q k T = I ε k − 1 i δ k − 1 i F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + ( δ k − 1 i ) 2 Q k T = I

(20)

Then, ε k − 1 i and δ k − 1 i can be obtained as

ε k − 1 i = δ k − 1 i = ( F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + Q k T ) − 1

(21)

Substituting equation (21) into equation (19), we have

P k | k − 1 i = ( ( F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + Q k T ) − 1 ( F k | k − 1 P k − 1 i F k | k − 1 T + Q k ) ) ( ( F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + Q k T ) − 1 ( F k | k − 1 P k − 1 i F k | k − 1 T + Q k ) ) T = Γ k | k − 1 i ( Γ k | k − 1 i ) T

(22)

where Γ k | k − 1 i = ( F k | k − 1 ( P k − 1 i ) T F k | k − 1 T + Q k T ) − 1 ( F k | k − 1 P k − 1 i F k | k − 1 T + Q k ) meets the requirements of the positive definiteness of P k | k − 1 i that can be expressed as the product of Γ k | k − 1 i and its transposed matrix. According to the Cholesky decomposition and the element γ , P k | k − 1 i in equation (18) overcomes the divergence of the Kalman filter.

At time k, suppose that the measurement z k j ( j = 1 , 2 , … , m ) is in the set { x k | k − 1 1 , x k | k − 1 2 , … , x k | k − 1 n } ∪ { c k n + 1 , c k n + 2 , … , c k m } , p j ( c k i ) ( i = n + 1 , … , m ) is the association probability of the random clutter and p j ( x k | k − 1 i ) ( i = 1 , … , n ) is the association probability of the target. The space distribution c k i ( z k j ) of clutters obeys the Poisson distribution with the clutter rate λ

p j ( c k i ) = exp ( − λ ) × λ c k i ( z k j )

(23)

With respect to j measurements, the occurring probabilities of the state x k | k − 1 i and the clutter c k i are given by

p ( x k | k − 1 i ) = ∑ j = 1 n p j ( x k | k − 1 i ) ∑ x k | k − 1 i ∈ X ∑ j = 1 n p j ( x k | k − 1 i ) + ∑ c k i ∈ C ∑ j = n + 1 m p j ( c k i ) = 1 m ∑ j = 1 n p j ( x k | k − 1 i )

(24)

p ( c k i ) = ∑ j = n + 1 m p j ( c k i ) ∑ x k | k − 1 i ∈ X ∑ j = 1 n p j ( x k | k − 1 i ) + ∑ c k i ∈ C ∑ j = n + 1 m p j ( c k i ) = 1 m ∑ j = n + 1 m p j ( c k i )

(25)

According to the total number of all measurements, we have the sum of probabilities

∑ i = 1 n p ( x k | k − 1 i ) + ∑ i = n + 1 m p ( c k i ) = ∑ j = 1 m ( ∑ i = 1 n p j ( x k | k − 1 i ) + ∑ i = n + 1 m p j ( c k i ) ) = m

(26)

Then, the relationship between p j ( x k | k − 1 i ) and p j ( c k i ) is

∑ i = 1 n p j ( x k | k − 1 i ) = 1 − ∑ i = n + 1 m p j ( c k i )

(27)

Table 1 presents the prediction process for the multitarget tracking. From the pseudocode, the proposed steps predict the transfer probability p D , k | k − 1 i , the target state x k | k − 1 i , and the estimation covariance P k | k − 1 i .

OPEN IN VIEWER

Table 1.

Prediction process.

Algorithm 1: prediction process of the proposed filtering method Input: the estimated probability of detection p D , k − 1 i , the forgetting factor b , the clutter rate λ , the estimation covariance P k − 1 i , and the target state x k − 1 i ; Output: the beta parameters α k | k − 1 i and β k | k − 1 i , the transfer probability p D , k | k − 1 i , the target state x k | k − 1 i , and the estimation covariance P k | k − 1 i .    for i = 1 , … , n   Compute the detection probability    Compute the beta parameters α k | k − 1 i and β k | k − 1 i using equations (4) and (5);    Compute the transfer probability p D , k | k − 1 i using equation (8);    Compute the filtering parameters    Compute the tune parameter d k using equation (12);    Compute the variances Q k and R k using equations (13) and (14);    Compute the elements γ rr and γ st using equation (17);    Computing the estimation covariance P k − 1 i using equation (16);    Compute the undetermined coefficients ε k − 1 i and δ k − 1 i using equation (21);    Predict the estimation covariance P k | k − 1 i using equation (18);    Predict the target state x k | k − 1 i using equation (10);     for j = 1 , … , n     Compute the occurring probability of target state p ( x k | k − 1 i ) using equation (24);     end j     for j = n + 1 , … , m     Compute the occurring probability of clutter p ( c k i ) using equation (25);     end j    end i

Update and track gate strategy

Suppose that the innovation v k j and the covariance S k i of measurements are

v k i , j = z k j − H k x k | k − 1 i

(28)

S k i = R k | k − 1 i + H k P k | k − 1 i H k T

(29)

According to the identity of the Gaussian density function, we define the correlation degree χ k i , j based on v k j

χ k i , j = 1 2 π | S k i | exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) 2

(30)

The probability p j ( x k | k i ) is updated as

p j ( x k | k i ) = ( 1 − ∑ i = n + 1 m p j ( c k i ) ) χ j i ∑ i = 1 n χ j i = ( 1 − ∑ i = n + 1 m p j ( c k i ) ) 1 2 π | S k i | exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) 2 ∑ i = 1 n 1 2 π | S k i | exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) 2 = ( 1 − ∑ i = n + 1 m p j ( c k i ) ) exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) 2 ∑ i = 1 n exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) 2

(31)

At time k, the Kalman filtering gain K k is updated as

K k = P k | k − 1 i H k T ( S k j ) − 1

(32)

The estimation covariance P k i is updated as

P k i = ( I − K k H k ) P k | k − 1 i

(33)

According to the lower triangular matrix Γ k i = Γ k | k − 1 i ( I − ζ k i ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T ) and the undetermined coefficient ζ k i = 1 ± 1 − ( Γ k | k − 1 i H k T ) T ( S k j ) − 1 ( Γ k | k − 1 i H k T ) / ( Γ k | k − 1 i H k T ) T ( S k j ) − 1 ( Γ k | k − 1 i H k T ) , the Cholesky decomposition of P k i is given by

P k i = Γ k i ( Γ k i ) T

(34)

Proof

Substituting equation (18) into equation (33), we have in hand the following form

P k i = Γ k | k − 1 i ( I − P k | k − 1 i H k T ( S k j ) − 1 H k ) ( Γ k | k − 1 i ) T = Γ k | k − 1 i ( I − Γ k | k − 1 i H k T ( S k j ) − 1 H k ( Γ k | k − 1 i ) T ) ( Γ k | k − 1 i ) T = Γ k | k − 1 i ( I − ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T ) ( Γ k | k − 1 i ) T

(35)

By defining the undetermined coefficient ζ k i , we can rewrite the factor ( I − ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T ) in equation (35) as follows

I − ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T = ( I − ζ k i ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T ) ( I − ζ k i ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T ) T = I − 2 ζ k i ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T + ( ζ k i ) 2 ( Γ k | k − 1 i H k T ) ( Γ k | k − 1 i H k T ) T ( ( S k j ) − 1 ) 2 × ( Γ k | k − 1 i H k T ) ( Γ k | k − 1 i H k T ) T = I − ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( 2 ζ k i − ( ζ k i ) 2 ( Γ k | k − 1 i H k T ) T ( S k j ) − 1 ( Γ k | k − 1 i H k T ) ) ( Γ k | k − 1 i H k T ) T

(36)

We have in hand the following equation

( ζ k i ) 2 ( Γ k | k − 1 i H k T ) T ( S k j ) − 1 ( Γ k | k − 1 i H k T ) − 2 ζ k i + I = 0

(37)

Solving equation (37), there is

ζ k i = 1 ± 1 − ( Γ k | k − 1 i H k T ) T ( S k j ) − 1 ( Γ k | k − 1 i H k T ) ( Γ k | k − 1 i H k T ) T ( S k j ) − 1 ( Γ k | k − 1 i H k T )

(38)

Then, equation (38) can be rewritten as

P k i = Γ k | k − 1 i ( I − ζ k i ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k ) T ) ( I − ζ k i ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k ) T ) T ( Γ k | k − 1 i ) T = Γ k i ( Γ k i ) T

(39)

where Γ k i = Γ k | k − 1 i ( I − ζ k i ( Γ k | k − 1 i H k T ) ( S k j ) − 1 ( Γ k | k − 1 i H k T ) T ) .

Subsequently, we explore the track gate strategy to distinguish current measurements as follows:

1. If the number of available measurements in the predicted gate is not less than 1, according to the relationship between the predicted state x k | k − 1 i of target and the associated measurement z k j , the target-generated measurement can be computed by weighting all measurements in the gate. At time k, suppose that ρ is the radius of the circular gate, we identify z k j when the target i is in the gate based on the actual distance | | v k i , j | |

‖ v k i , j ‖ = ‖ z k j − H k x k | k − 1 i ‖

(40)

Under the condition of | | v k i , j | | ≤ ρ , z k j is related to the gate of all targets, and the total number of measurements is less than m in the set. If a measurement falls into the track gate, the measurement is from the target. The target-generated measurement is modified by weighting when j = 1 , 2 , … , n . Once the measurement j is recorded, the current measurement in the gate is updated as

z k | k j = exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) λ ( c ( z k j ) ) − 1 2 π | S k i | ( 1 − P D , k ) P D , k − 1 + ∑ i = 1 n exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) z k j

(41)

where z k | k j is the product of z k j and the corresponding target-generated probability, that is, the weighting factor of z k j . Furthermore, exp ( − 1 / 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) is the association density of the measurement j and the target i, and λ ( c ( z k j ) ) − 1 2 π | S k i | ( 1 − P D , k ) P D , k − 1 is the association density of the measurement j and the clutter i. Note that the numerator becomes large when v k i , j is small. As a result, the weight of z k j becomes large, and the component of z k | k j provided by the target i is large too. On the contrary, the component of z k | k j provided by the target i is small when v k i , j is large.

2. If z k j is outside the track gate, it is from the clutter. Similarly, the clutter-generated measurement can be updated when j = n + 1 , n + 2 , … , m

z k | k j = λ ( c ( z k j ) ) − 1 2 π | S k i | ( 1 − P D , k ) P D , k − 1 λ ( c ( z k j ) ) − 1 2 π | S k i | ( 1 − P D , k ) P D , k − 1 + ∑ i = 1 n exp ( − 1 2 ( v k i , j ) T ( S k i ) − 1 v k i , j ) z k j

(42)

where z k | k j is the product of z k j and the corresponding clutter-generated probability, that is, the weighting factor of z k j .

3. If there is no target in the track gate, the number of missed targets increases 1. According to the actual tracking strategy, we consider the target disappears when it cannot be detected within the continuous five periods. As a result, this target is deleted when the total number of missed targets during the continuous periods is more than five.

Given that each updated measurement z k | k j has impact on the target state x k | k i , the global correction for x k | k i is updated as

x k | k i = x k | k − 1 i + ∑ j = 1 m p j ( x k | k − 1 i ) K k ( z k | k j − H k x k | k − 1 i )

(43)

Furthermore, we compute the error distance to determine the target track ℓ k i , j

ℓ k i , j = ‖ z k | k j − H k x k | k i ‖

(44)

Under the condition of ℓ k i , j ≤ ρ , we delete the row and column related to ℓ k i , j , and then determine ℓ k i , j in turn. Repeating the steps, we delete the row related to the last state of the target. For example, the rest elements z k | k n + 1 , … , z k | k m express the clutter-generated measurements. Table 2 shows the update process of the proposed method. After identifying target tracks, the clutter-generated measurements are obtained.

OPEN IN VIEWER

Table 2.

Track gate strategy.

If x k | k i is remained after the mentioned operations, the current state belongs to the newborn target when a row has not been deleted using equation (44). At this time, we consider that the estimated velocity of the target keeps unchanged and the estimated position of target is the center of the track gate. Table 3 presents the update process. From the pseudocode, the steps are used to update the detection probability p D , k , the target state x k | k i , the estimation covariance P k i , and the target track ℓ k i , j .

OPEN IN VIEWER

Table 3.

Update process.

Algorithm 2: update process of the proposed filtering method Input: the beta parameters α k | k − 1 i and β k | k − 1 i , the forgetting factor b , the threshold ε , the estimation covariance P k | k − 1 i , and the target state x k | k − 1 i ; Output: the detection probability p D , k , the target state x k | k i , the estimation covariance P k i , and track label ℓ k i , j .    for i = 1 , ⋯ , n    Compute the detection probability     Update the detection probability p D , k using equation (7);    Compute the filtering parameters and state estimation     for j = 1 , … , n     Compute the innovation v k j of measurements using equation (28);     Compute the covariance S k i of measurements using equation (29);     Compute the correlation degree χ k i , j using equation (30);     Compute the probability p j ( x k | k i ) using equation (31);     Compute the Kalman filtering gain K k using equation (32);     Compute the undetermined coefficient ζ k i using equation (38);     Update the estimation covariance P k i using equation (34);     Compute the actual distance | | v k i , j | | using equation (40);       if | | v k i , j | | ≤ ε is false            then The detection number of missed targets increases 1. When the number of the missed targets is more than 5, the associated target is deleted.            else Compute the target and clutter measurements z k | k j in the associated gate using equation (41);            end         end j         for j = n + 1 , … , m          Compute the target and clutter measurements z k | k j in the associated gate using equation (42);         end j        Update the target state x k | k i using equation (43);         for j = 1 , … , m          Compute the track ℓ k i , j using equation (44) and determine it;         end j    end i