A Novel Algorithm of Joint Probability Data Association Based on Loss Function

Abstract

In this paper, a joint probabilistic data association algorithm based on loss function (LJPDA) is proposed so that the computation load and accuracy of the multi-target tracking algorithm can be guaranteed simultaneously. Firstly, data association is divided in to three cases based on the relationship among validation gates and the number of measurements in the overlapping area for validation gates. Also the contribution coefficient is employed for evaluating the contribution of a measurement to a target, and the loss function, which reflects the cost of the new proposed data association algorithm, is defined. Moreover, the equation set of optimal contribution coefficient is given by minimizing the loss function, and the optimal contribution coefficient can be attained by using the Newton-Raphson method. In this way, the weighted value of each target can be achieved, and the data association among measurements and tracks can be realized. Finally, we compare performances of LJPDA proposed and joint probabilistic data association (JPDA) algorithm via numerical simulations, and much attention is paid on real-time performance and estimation error. Theoretical analysis and experimental results reveal that the LJPDA algorithm proposed exhibits small estimation error and low computation complexity.

1. Introduction

Radar data processing is widely used in data association, smoothing, filtering, and prediction of the measured data, and estimation for the parameters of targets [1,2]. In this way, stable tracks can be attained.

There are mainly two basic methods of radar data association for multiple target tracking, the method based on the maximum likelihood principle, and the method based on the Bayes criterion [3-7]. The former shows huge computation burden and is not easy to implement [8]. The latter includes nearest neighbor data association (NNDA), probabilistic data association (PDA), multiple hypothesis tracking (MHT) and joint probabilistic data association (JPDA) algorithm [9,10]. NNDA is not appropriate for environments with severe clutters or multiple targets [11,12], and PDA is only suitable for single-target and sparse multi-target environments [13]. MHT assumes that each new measurement may come from an existing target, a new target or a clutter, and achieve multi-target tracking by hypothesis evaluation, deletion and merging [14]. MHT is seldom adopted since it depends too much on a prior information of targets and clutters. Although JPDA exhibits good performance for multi-target tracking, but it is difficult to implement due to the reason as following: on the one hand, JPDA is very sophisticated and has huge computation burden for single target case; on the other hand, computation load increases rapidly and explosively when the number of tracking targets increases [15,16]. Therefore, JPDA is not easy to implement in practical engineering [1,17]. In this paper, we consider the distribution of measurements in validation gates, and combine the Bayesian criterion with the correlation characteristics between validation gates, and introduce loss (cost) function to calculate the joint probability between the measurement and the target. In this way, a multi-target joint probabilistic data association algorithm based on loss function (LJPDA) is proposed. The effectiveness and timeliness of the proposed algorithm are verified by simulation.

2. LJPDA algorithm

In the case of multiple targets or clutters, there are likely to be multiple measurements that fall into validation gates of the target. These measurements are referred to as valid measurements [1,2]. The basic of joint probability is that each valid measurement may come from the tracked target, and different valid measurement has different probability from the tracked target. So, each valid measurement in a validation gate can make a contribution to update of the associated target. In our study, we adopt a contribution efficient to depict the contribution of a valid measurement on the associated target. Then we define Loss function based on estimation error for evaluating the cost of a data association algorithm. Further, we derive the loss function relevant to contribution efficient in different cases, and find the optimal contribution efficient by minimizing the loss function. Finally, states can be updated by using the optimal contribution efficient.

2.1 Loss function

JPDA statistically computes the impact of each valid measurement on state estimation. However, in LJPDA proposed, we define a loss function tightly related to estimation error, and consider the contribution of each valid measurement on the associated target. Then, we give the expression of loss function, and attain the optimal contribution efficient.

The filter residual (Innovation) for a tracking filter can be written as

\(V(k)=Z(k)-H X(k \mid k-1)\)       (1)

where Z(k) is the measurement vector, H is the measurement matrix, and X(k|k-1) is the prediction matrix of the state at time k given k-1 updates.

The innovation of the target t being associated with the measurement j can be expressed as

\(U_{j t}(k)=Z_{j}(k)-H X_{t}(k \mid k-1)\)

The likelihood function P_jt[1] of the measurement j being associated with the target t is given by

\(P_{j t}=\frac{1}{\left|2 \pi S_{t}(k)\right|^{1 / 2}} \exp \left\{-\frac{1}{2} U_{j t}^{\prime}(k) S_{t}^{-1}(k) U_{j t}(k)\right\}\)       (2)

where S_t(k) is often named as the innovation covariance matrix of target t.

We define β_jt as the association probability of the target t being associated with the measurement j. The two powers of Mahalanobis distance is defined as the loss function.

\(d_{j t}^{2}=U_{j t}^{\prime} S_{t}^{-1} U_{j t}\)

where \(d_{j r}^{2}\) is the two powers of the Mahalanobis distance between the measurement j and the target t.

Therefore, the association loss of target t is defined as

C_t(k) = (x% - x)'S_t(k)^-1(x% - x)

where x is the true value, x% is the estimated value, and S_t is the sample covariance matrix.

2.2 Classification

On the basis of the relationship among validation gates and the number of measurements in the overlapping area for validation gates, a data association algorithm is divided into three cases as follows.

2.2.1 Case 1

As shown in Fig. 1, there is no overlap in the validation gates or no measurement in the overlapping area.

Fig. 1. No measurement in the overlapping area or no overlapping area.

In this case, the multi-target data association problem can be considered as the single objective data association problem, and the PDA algorithm can be used for data association [2]. The combined innovation V_t(k) of the target t is given by

\(V_{t}(k)=\sum_{j=1}^{m} \frac{P_{j t} U_{j t}(k)}{\sum_{j=1}^{m} P_{j t}+B}=\sum_{j=1}^{m} \beta_{j t} U_{j t}(k) t=1,2, \cdots, N\)       (4)

where \(\beta_{j t}=\frac{P_{j t}}{\sum_{i=1}^{m} P_{i t}+B}\) is the association probability, m is the number of valid measurements in the validation gate of the target t, N is the number of targets, B = λ|2πS(k)|^1/2(1-P_DP_G)/P_G is a constant that depends on the clutter density [1], λ is the spatial density of false measurements, P_D is the target detection probability, and P_G is the gate probability which is the probability of a real target existing in the validation gate.

2.2.2 Case 2

As illustrated in Fig. 2, there is only one measurement in validation gates of the target that falls in the overlapping area.

Fig. 2. Only one measurement in validation gates and the measurement falling in the overlapping area

In Fig. 2, there is only one measurement (the measurement 1) in the validation gates for the target 1 and target 2, and the measurement 1 falls in the overlapping area. In this case, there are four kinds of possibilities as following:

I. The measurement 1 only comes from the target 1, its detection probability can be expressed as P₁=P_D1(1-P_D2).

II. The measurement 1 is only from the target 2, its detection probability can be written into P₂=P_D2(1-P_D1).

III. The measurement 1 comes from both the target 1 and target 2, its detection probability is P₃=P_D1P_D2.

IV. The measurement 1 is from false alarms, its detection probability can be expressed as P₄=(1-P_D1)(1-P_D2).

where P_Di is the detection probability of the target i, and P₁ + P₂ + P₃ + P₄ = 1.

Data association is used to establish a relationship between measurements and targets [18-20]. If the assumption in which the measurement is only from the target 1 holds, then the four situations mentioned above are real, and losses of data association (loss function of data association) can be generated in the following.

When I is the real situation, there is no association error. Loss of data association is C_І = 0. If the real situation is Ⅱ, the association of target 1 and 2 are both wrong. Loss of data association is \(C_{\mathrm{II}}=d_{11}^{2}+d_{12}^{2}\).

If Ⅲ holds, the association of target 2 is wrong. Loss of data association is \(C_{\mathrm{III}}=d_{12}^{2}\).

When the real situation is Ⅳ, the association of target 1 is wrong. Loss of data association is \(C_{\mathrm{IV}}=d_{11}^{2}\).

where \(d_{1 i}^{2}=U_{1 i}^{\prime} S_{i}^{-1} U_{1 i}\) is the two powers of the Mahalanobis distance between the measurement and the target i.

Therefore, the average loss function in the case I is given by

\(\begin{aligned} C_{10} &=P_{1} C_{I}+P_{2} C_{I I}+P_{3} C_{I I I}+P_{4} C_{I V} \\ &=P_{2}\left(d_{11}^{2}+d_{12}^{2}\right)+P_{3} d_{12}^{2}+P_{4} d_{11}^{2} \\ &=\left(1-P_{D 1}\right) d_{11}^{2}+P_{D 2} d_{12}^{2} \end{aligned}\)

Similarly, the average association loss in Ⅱ-Ⅳ are given as follows

\(\begin{gathered} C_{01}=P_{D 1} d_{11}^{2}+\left(1-P_{D 2}\right) d_{12}^{2} \\ C_{11}=\left(1-P_{D 1}\right) d_{11}^{2}+\left(1-P_{D 2}\right) d_{12}^{2} \\ C_{00}=P_{D 1} d_{11}^{2}+P_{D 2} d_{12}^{2} \end{gathered}\)

where C_ij is the average loss function, i = 1 represents the measurement 1 is associated with the target 1, and j = 1 indicates the measurement 1 is associated with the target 2.

In the case 2, the measurement should not be associated with multiple targets, otherwise it will cause trajectory overlap [21]. Consequently, we need not to consider the situation, in which a measurement is from multiple targets, and the algorithm can be simplified and extended to situations where multiple validation gates are overlapping.

We suppose that there are m overlapping gates, and only one measurement in validation gates falls in the overlapping area. We further assume events in which a measurement comes from the target i(i=1, 2, …, m) are independent, and probability of the event i is P_Di(i=1, 2, …, m). When the measurement is associated with the target n, the average association loss C_n is given by

\(\begin{aligned} C_{n} &=P_{D 1} d_{11}^{2}+P_{D 2} d_{12}^{2}+\mathrm{L}+\left(1-P_{D n}\right) d_{1 n}^{2}+\mathrm{L}+P_{D m} d_{1 m}^{2} \\ &=\left(1-2 P_{D_{n}}\right) d_{1 n}^{2}+\sum_{i=1}^{m} P_{D i} d_{1 i}^{2} \quad n=1,2, \mathrm{~L}, m \end{aligned}\)       (5)

If the measurement is regarded as clutter, the average loss function can be written into

\(C_{0}=\sum_{i=1}^{m} P_{D i} d_{1 i}^{2}\)       (6)

Then the optimal association target for the measurement is given in the following by minimized the average loss function

\(k_{0}=\underset{n}{\arg \min }\left\{C_{n}, 1 \leq n \leq m\right\}\)

According to (5), the above formula can be simplified as

\(k_{0}=\arg \min \left\{\left(1-2 P_{D n}\right) d_{1 n}^{2}, 1 \leq n \leq m\right\}\)

The combined innovation V_t(k) of the target t is given by

\(V_{t}(k)= \begin{cases}U_{1 t}(k) & t=t_{0} \\ 0 & \text { else }\end{cases}\)       (7)

2.2.3 Case 3

As shown in Fig. 3, the case 3 denotes all the possibilities except for case 1 and case 2. In this case, there are at least one measurement in the overlapping area, and there are multiple measurements in the validation gates.

Fig. 3. At least one measurement in the overlapping area and multiple measurements in the validation gates

1.Contribution coefficient

On the basis of the idea of joint probability, the measurement in the overlapping area can affect state update of all targets in the same area. We define λ_it (0 ≤ λ_it ≤ 1) as the contribution coefficient of the measurement i being associated with target t, the larger λ_it is, the greater the impact of the measurement i on the associated with target t is. Taking Fig.3 as an example, the measurement 2 affects the data association of the target 1 and the target 2, and its contribution coefficients are λ₂₁ and λ₂₂ (λ₂₁ + λ₂₂ = 1), respectively. Similarly, the measurement 3 plays impact on the association of the target 1 and the target 2, and its contribution coefficients are respectively λ₃₁ and λ₃₂. In this way, the association loss function of data association is a function of the contribution coefficient, and the optimal contribution coefficient can be attained if the average loss function is minimized.

2. Formulation

If only the measurement 2 is considered, contribution coefficients of the measurement 3 are determined. Similar to the PDA algorithm [2], the association probability of the measurement 2 being associated with the target 1 is given when the measurement 2 is really from target 1.

\(\beta_{21}^{\prime}=P_{21} /\left(P_{11}+P_{21}+\lambda_{31} P_{31}+B\right)\),

where β′_ij is called the real association probability.

If the measurement 2 is not from the target 1, the real association probability of the measurement 2 being associated with the target 1 is β′₂₁ = 0

Similarly, when the measurement 2 comes from the target 2, the real association probability of the measurement 2 being associated with the target 2 is given by

\(\beta_{22}^{\prime}=P_{22} /\left(P_{22}+\lambda_{32} P_{32}+P_{42}+\lambda_{52} P_{52}+B\right)\).

If the measurement 2 is not from the target 2, the real association probability of the measurement 2 being associated with the target 2 is β′₂₂ = 0

In engineering practice, we do not know which target a measurement is from. Therefore, the association probability for the situation in which the measurement 2 is associated with targets 1 and 2 are given as follows

\(\begin{gathered} \beta_{21}^{\prime \prime}=\lambda_{21} P_{21} /\left(P_{11}+\lambda_{21} P_{21}+\lambda_{31} P_{31}+B\right) \\ \beta_{22}^{\prime \prime}=\lambda_{22} P_{22} /\left(\lambda_{22} P_{22}+\lambda_{32} P_{32}+P_{42}+\lambda_{52} P_{52}+B\right) \end{gathered}\).

where \(\beta_{i j}^{\prime \prime}\) is called the practical association probability, and P_ij and B can refer to (4).

Considering the impact of contribution coefficient when the measurement i is associated with the target j, the association loss is given by

\(C_{i j}^{A}=\left(\beta_{i j}^{\prime \prime}-\beta_{i j}^{\prime}\right)^{2} d_{i j}^{2}\)       (8)

When the measurement 2 is only comes from the target 1, the detection probability of a target is P_A=P_D1(1-P_D2) and the real association probability is β′₂₂ = 0. In this case, the loss function can be written into

\(\begin{aligned} C_{A} &=\left(\beta_{21}^{\prime \prime}-\beta_{21}^{\prime}\right)^{2} d_{21}^{2}+\left(\beta_{22}^{\prime \prime}-\beta_{22}^{\prime}\right)^{2} d_{22}^{2} \\ &=\left(\beta_{21}^{\prime \prime}-\beta_{21}^{\prime}\right)^{2} d_{21}^{2}+\beta_{22}^{\prime \prime} d_{22}^{2} \end{aligned}\)

Similarly, if the measurement 2 only comes from the target 2, detection probability is P_B=P_D2(1-P_D1), and the association loss can be expressed as

\(\begin{aligned} C_{B} &=\left(\beta_{21}^{\prime \prime}-\beta_{21}^{\prime}\right)^{2} d_{21}^{2}+\left(\beta_{22}^{\prime \prime}-\beta_{22}^{\prime}\right)^{2} d_{22}^{2} \\ &=\beta_{21}^{\prime \prime} d_{21}^{2}+\left(\beta_{22}^{\prime \prime}-\beta_{22}^{\prime}\right)^{2} d_{22}^{2} \end{aligned}\)

If the measurement 2 is from the target 1 and the target 2, detection probability is P_c=P_D1P_D2, and the association loss is given by

\(C_{C}=\left(\beta_{21}^{\prime \prime}-\beta_{21}^{\prime}\right)^{2} d_{21}^{2}+\left(\beta_{22}^{\prime \prime}-\beta_{22}^{\prime}\right)^{2} d_{22}^{2}\)

When the measurement 2 is from false alarm, detection probability is (1-P_D1)(1-P_D2), and the association loss can be written into

\(\begin{aligned} C_{D} &=\left(\beta_{21}^{\prime \prime}-\beta_{21}^{\prime}\right)^{2} d_{21}^{2}+\left(\beta_{22}^{\prime \prime}-\beta_{22}^{\prime}\right)^{2} d_{22}^{2} \\ &=\beta_{21}^{\prime \prime 2} d_{21}^{2}+\beta_{22}^{\prime \prime} d_{22}^{2} \end{aligned}\)

Therefore, the average loss function of data association for the measurement 2 is given by

\(\begin{aligned} D_{2} &=P_{A}\left[\left(\beta_{21}^{\prime \prime}-\beta_{21}^{\prime}\right)^{2} d_{21}^{2}+\beta_{22}^{\prime \prime 2} d_{22}^{2}\right] \\ &+P_{B}\left[\beta_{21}^{\prime 2} d_{21}^{2}+\left(\beta_{22}^{\prime \prime}-\beta_{22}^{\prime}\right)^{2} d_{22}^{2}\right] \\ &+P_{C}\left[\left(\beta_{21}^{\prime \prime}-\beta_{21}^{\prime}\right)^{2} d_{21}^{2}+\left(\beta_{22}^{\prime \prime}-\beta_{22}^{\prime}\right)^{2} d_{22}^{2}\right] \\ &+\left(1-P_{D 1}\right)\left(1-P_{D 2}\right)\left(\beta_{21}^{\prime 2} d_{21}^{2}+\beta_{22}^{\prime 2} d_{22}^{2}\right) \\ =& P_{D 1} \beta_{21}^{\prime 2} d_{21}^{2}+\beta_{21}^{\prime 2} d_{21}^{2}-2 P_{D 1} \beta_{21}^{\prime} \beta_{21}^{\prime \prime} d_{21}^{2} \\ &+P_{D 2} \beta_{22}^{\prime 2} d_{22}^{2}+\beta_{22}^{\prime \prime 2} d_{22}^{2}-2 P_{D 2} \beta_{22}^{\prime} \beta_{22}^{\prime \prime} d_{22}^{2} \end{aligned}\)       (9)

where P_D1 and P_D2 are the detection probabilities of the target 1 and the target 2.

It can be seen from (9) that the average loss function of data association D₂ is a function of contribution coefficients. The optimal contribution coefficients can be obtained by minimizing D₂. Taking partial derivative of λ₂₁ on both side of D₂ and letting results be zero, we obtain

\(\begin{aligned} \frac{\partial D_{2}}{\partial \lambda_{21}}=& \frac{-2 P_{D 1} \beta_{21}^{\prime} \beta_{21}^{\prime \prime} d_{21}^{2}+2 P_{D 1} \beta_{21}^{\prime} \beta_{21}^{\prime 2} d_{21}^{2}+2 \beta_{21}^{\prime 2} d_{21}^{2}-2 \beta_{21}^{\prime 3} d_{21}^{2}}{\lambda_{21}} \\ &+\frac{2 P_{D 2} \beta_{22}^{\prime} \beta_{22}^{\prime \prime} d_{22}^{2}-2 P_{D 2} \beta_{22}^{\prime} \beta_{22}^{\prime 2} d_{22}^{2}-2 \beta_{22}^{\prime 2} d_{22}^{2}+2 \beta_{22}^{\prime 3} d_{22}^{2}}{\lambda_{22}} \end{aligned}\)

The equation set for the optimal contribution coefficients can be attained

\(\left\{\begin{array}{l} \lambda_{21}=\frac{\left(P_{D 1} \beta_{21}^{\prime} \beta_{21}^{\prime \prime} d_{21}^{2}-P_{D 1} \beta_{21}^{\prime} \beta_{21}^{\prime \prime 2} d_{21}^{2}-\beta_{21}^{\prime 2} d_{21}^{2}+\beta_{21}^{\prime 3} d_{21}^{2}\right) \lambda_{22}}{\left(P_{D 2} \beta_{22}^{\prime} \beta_{22}^{\prime \prime} d_{22}^{2}-P_{D 2} \beta_{22}^{\prime} \beta_{22}^{\prime 2} d_{22}^{2}-\beta_{22}^{\prime 2} d_{22}^{2}+\beta_{22}^{\prime 3} d_{22}^{2}\right)} \\ \lambda_{22}=1-\lambda_{21} \end{array}\right.\)       (10)

Equation (10) can be extended to situations in which multiple measurements appear in the overlapping area of multiple validation gates. We assume that the number of measurements is T, the number of validation gates is M, and the events in which the measurement i is from the target j (j=1, 2, …, M) are independent (occurrence probability of the event is P_Dj). The average loss function of data association for the situation in which the measurement i is associated with the target j is given by

\(\begin{aligned} C_{i j} &=P_{D j}\left(\beta_{i j}^{\prime \prime}-\beta_{i j}^{\prime}\right)^{2} d_{i j}^{2}+\left(1-P_{D j}\right)\left(\beta_{i j}^{\prime \prime}\right)^{2} d_{i j}^{2} \\ &=P_{D j} \beta_{i j}^{\prime 2} d_{i j}^{2}+\beta_{i j}^{\prime 2} d_{i j}^{2}-2 P_{D j} \beta_{i j}^{\prime} \beta_{i j}^{\prime \prime} d_{i j}^{2} \end{aligned}\)       (11)

The average loss function of data association for the measurement i by adopting the contribution coefficient can be given by

\(\begin{aligned} D_{i} &=\sum_{j=1}^{M} C_{i j} \\ &=\sum_{j=1}^{M}\left(P_{D j} \beta_{i j}^{\prime 2} d_{i j}^{2}+\beta_{i j}^{\prime 2} d_{i j}^{2}-2 P_{D j} \beta_{i j}^{\prime} \beta_{i j}^{\prime \prime} d_{i j}^{2}\right) \end{aligned}\)       (12)

Using Lagrange Multiplier Method, the modified function F is obtained.

\(\begin{aligned} F &=D_{i}+a\left(\sum_{j=1}^{M} \lambda_{i j}-1\right) \\ &=\sum_{j=1}^{M}\left(P_{D j} \beta_{i j}^{\prime 2} d_{i j}^{2}+\beta_{i j}^{\prime 2} d_{i j}^{2}-2 P_{D j} \beta_{i j}^{\prime} \beta_{i j}^{\prime \prime} d_{i j}^{2}\right)+a\left(\sum_{j=1}^{M} \lambda_{i j}-1\right) \end{aligned}\)       (13)

where a is the undetermined coefficient.

Taking partial derivative of λ_jt on both side of (13) and letting results be zero, we obtain

\(\frac{\partial F}{\partial \lambda_{j t}}=-2 \frac{\left(P_{D t} \beta_{j t}^{\prime}-\beta_{j t}^{\prime \prime}\right)\left(1-\beta_{j t}^{\prime \prime}\right) \beta_{j t}^{\prime \prime} d_{j t}^{2}}{\lambda_{j t}}+a\)       (14)

where P_Dt is the detection probability of the target t, \(d_{j t}^{2}\) denotes the two powers of the Mahalanobis distance between the measurement j and the target t. Additionally, \(\beta_{j t}^{\prime}\) and \( \beta_{j t}^{\prime \prime}\) are the real association probability and the practical association probability when the measurement j comes from the target t. Can be given as

\(\left\{\begin{array}{l} \beta_{j t}^{\prime}=P_{j t} /\left(\sum_{i=1}^{T} \lambda_{i t} P_{i t}-\lambda_{j t} P_{j t}+P_{j t}+B\right) \\ \beta_{j t}^{\prime \prime}=\lambda_{j t} P_{j t} /\left(\sum_{i=1}^{T} \lambda_{i t} P_{i t}+B\right) \end{array}\right.\)       (15)

It follows from (14) that

\(\left\{\begin{array}{l} a=\frac{2 D_{j 1}}{\lambda_{j 1}}=\frac{2 D_{j 2}}{\lambda_{j 2}}=\mathrm{L}=\frac{2 D_{j M}}{\lambda_{j M}} \\ 1=\sum_{i=1}^{M} \lambda_{j i} \end{array}\right.\)       (16)

where \(D_{j t}=\left(P_{D t} \beta_{j t}^{\prime}-\beta_{j t}^{\prime \prime}\right)\left(1-\beta_{j t}^{\prime \prime}\right) \beta_{j t}^{\prime \prime} d_{j t}^{2}\), and M is the number of validation gates where the measurement j is located.

If the solution of (16) does not exist in the domain of definition λ_it (0 ≤ λ_it ≤ 1), the minimum value of contribution coefficient for the modified function F is on a one-dimensional or multidimensional boundary. In this case, we let the contribution coefficient of this dimension equal to zero. Then we solve the equation set in (15). In addition, the measurement in a non overlapping area completely contributes to the corresponding target t. i.e. λ_jt = 1.

Newton-Raphson method (tangents method) [22] is adopted for solving (16), and the nonlinear equation set can be approximated as linear equation set for iterative solution. Therefore, equation (16) can be approximated to (17).

\(\left\{\begin{array}{l} \text { (1) } \beta_{j t}^{\prime}=P_{j t} /\left(\sum_{i=1}^{T} \lambda_{i t} P_{i t}-\lambda_{j t} P_{j t}+P_{j t}+B\right) \\ \text { (2) } \beta_{j t}^{\prime \prime}=\lambda_{j t} P_{j t} /\left(\sum_{i=1}^{T} \lambda_{i t} P_{i t}+B\right) \\ \text { (3) } D_{j t}=\left(P_{D t} \beta_{j t}^{\prime}-\beta_{j t}^{\prime \prime}\right)\left(1-\beta_{j t}^{\prime \prime}\right) \beta_{j t}^{\prime \prime} d_{j t}^{2} \\ \text { (4) } \frac{\partial D_{j t}}{\partial \lambda_{j t}}=\frac{\left[3 \beta_{j t}^{\prime 2}-2\left(1+P_{D t} \beta_{j t}^{\prime}\right) \beta_{j t}^{\prime \prime}+P_{D t} \beta_{j t}^{\prime}\right]\left(1-\beta_{j t}^{\prime \prime}\right) \beta_{j t}^{\prime \prime} d_{j t}^{2}}{\lambda_{j t}} \\ \text { (5) } \lambda_{j t}^{+}=\lambda_{j t}-\frac{D_{j t} \lambda_{j 1}-D_{j 1} \lambda_{j t}}{\frac{\partial D_{j t}}{\partial \lambda_{j t}} g_{j 1}-D_{j 1}} \end{array}\right.\)       (17)

where \(\frac{\partial D_{j t}}{\partial \lambda_{j t}}\) is the partial derivative of D_jt on λ_jt, λ_jt denotes the contribution of the measurement j being associated with the target t, and \(\lambda_{j t}^{+}\) is the contribution coefficient of the measurement j being associated with the target t after the iteration \(\left(\lambda_{j t}^{+} \in\{2,3, \ldots, M\}\right)\). In addition, if t is equal to 1, the contribution coefficient \(\lambda_{j t}^{+}\) after iteration is equal to \(1-\left(\sum_{i=2}^{M} \lambda_{j i}^{+}\right)\). The iterative process is summarized as follows.

Step 1 Initiation of contribution coefficients. λ_jt is set to be 1/M for the contribution coefficient in the overlapping area, and λ_jt is equal to 1 in the non-overlapping area. The fixed constant β′_jt is calculated.

Step 2 The process variables \(\beta_{j t}^{n}\), D_jt and \(\frac{\partial D_{j t}}{\partial \lambda_{j t}}\) are calculated.

Step 3 \(\lambda_{j t}^{+}\) is updated and \(\lambda_{j 1}^{+}=1-\left(\sum_{i=2}^{M} \lambda_{j i}^{+}\right)\).

Step 4 λ_jt0 and \(\lambda_{j t_{0}}^{+}\) are equal to 0, if  \(\lambda_{j t_{0}}^{+} \leq 0\).

Step 5 Making \(\lambda_{j t}=\lambda_{j t}^{+}\), and repeating from the step 2 to the step 4 until \(\left|\lambda_{j t}^{+}-\lambda_{j t}\right| \leq \sigma\) the iteration limit is reached.

It can be seen that a single iteration of the measurement j being associated with the target t, (T+17) multiplications, 3 divisions, and (T+9) additions are necessary, where T is the number of measurements in the target t (validation gate t, \(t \in\{2,3, \ldots, M\}\)). The combined innovation V_t(k) of the target t is

\(V_{t}(k)=\sum_{j=1}^{T} \frac{\lambda_{j t} P_{j t} U_{j t}(k)}{\sum_{j=1}^{T} \lambda_{j t} P_{j t}+B}=\sum_{j=1}^{T} \beta_{j t} U_{j t}(k) \quad t=1,2, \ldots, M\)       (18)

where \(\beta_{j t}=\frac{\lambda_{j t} P_{j t}}{\sum_{j=1}^{T} \lambda_{j t} P_{j t}+B}\) is association probability, M is the target number, and B is a constant depending on the clutter density.

Finally, the combined innovation V_t(k) calculated in the above three cases is used for each target to update the current state of the target (track filtering).

3. Flow diagram of the LJPDA algorithm

The multi-target tracking algorithm based on loss function is summarized as follows.

Step 1 Target state prediction and innovation Ujt(k) between measurements and targets are calculated.

Step 2 The relationship between the measurement and the target validation gate according to (2)-(3) is determined, and \(d_{j t}^{2}\) and the likelihood function P_jt are calculated.

Step 3 According to the relationship between the measurement and the target validation gate, the combined innovation V_t(k) of the target is calculated. Some specific process should be noticed.

1) If the validation gates do not overlap, or there is no measurement in the overlapping area, the combined innovation V_t(k) of the target is calculated according to (4).

2) If there is only one measurement in validation gates of the target that falls in the overlapping area, the measurement is associated with the target k₀ and \(k_{0}=\underset{n}{\arg \min }\left\{\left(1-2 P_{D n}\right) d_{1 n}^{2}, 1 \leq n \leq m\right\}\). Then combined innovation V_t(k) is calculated according to (7)

3) If there is at least one measurement in the overlapping area, also, if there are multiple measurements in the validation gates, we used iteration mothed to calculate the optimal contribution coefficient λ_jt referring to case 3. Finally, the combined innovation V_t(k) is calculated according to (18).

Step 4 Filtering is performed to update the current state of the target. LJPDA based on Kalman filter [23] is shown in Fig. 4.

Fig. 4. Flow chart of LJPDA

4. Algorithm simulation

4.1 Validity simulation

In our experiments, there are four targets moving at constant speeds simultaneously and the initial states in Cartesian coordinate system (X = [x; v_x; y; v_y]) are given as X₁ = [-6; 1; 0; 25], X₂ = [-2; 0; 0; 24], X₃ = [2; 1; 0; 28] and X₄ = [6; -1; 0; 27], respectively. Units of range and velocity are meter and meter per second, respectively. Moreover, the number of sampling points is equal to 85, standard deviation of radar ranging and angle measurements are σ_Ρ = 1.6m and σ_ѳ = 0.02rad, respectively. In addition, the number of clutter in the scanning interval follows the Poisson distribution with parameter λ = 5, the position of each clutter is uniformly distributed in the scanning interval, the detection probability P_D = 0.98, and the validation gate probability P_G = 0.989. Fig. 5 shows real and measurement tracks.

Fig. 5. Real and measurement tracks

Fig. 5 a, b, c, d exhibit tracks of four targets for real tracks, tracks prior to data association when considering the impact of noise, tracks prior to data association when considering the effect of noise and clutter, and tracks obtained after data association and track filtering considering the impact of noise and clutter, respectively. It can be seen from Fig .5a there are two intersections in real tracks and targets are very dense. It can further be seen from Fig. 5c there are too many overlapping areas in tracks when considering the impact of noise and clutter. Consequently, it is hard to perform data association successfully in our experiments with severe environment. However, Fig. 5d reveals that LJPDA can accurately associate measurements with real tracks in our experiments, where dots indicate data association is successful, and circles denote data association is unsuccessful. Extensive numerical results have verified the good performance of target tracking for LJPDA

4.2 Performance analysis

In this subsection, we analyze performance of the LJPDA algorithm proposed, and compare it with JPDA which is acknowledged as one of the ideal methods [13].

We assume there are two intersecting and uniform moving targets in our experiments, and their initial positions are X₁ = [-29500; 400; 34500; -400] and X₂ = [-26200; -200; 34500; -400], respectively.

Parameters used in our experiments are shown in Table 1, where Q and γ are noise components and validation gate parameters, respectively. Additionally, false measurements were produced uniformly in a square centered at the correct measurements with an area of A = n_c / λ = 10A_v, where λ denotes the number of false measurements per unit area, and λ = 0.004. n_c is the total number of false measurements, and     \(A_{v}=\pi \gamma|S(k)|^{1 / 2}\). Moreover, there are 70 simulation steps in every run, and totally 50 runs. The real trajectories of the two targets and measurements falling into the gates are shown in Fig. 6. The root mean square errors (RMSE) of the two algorithms are shown in Fig. 7.

Table 1. Simulation parameters

Fig. 6. The real trajectories and the measurements in validation gates

Fig. 7. RMSE of range for the target 1 and target 2

As shown in Fig. 6, many unwanted echoes (clutter) exist in the validation gate, which lead to many intersections in measurement trajectories Fig. 7 shows the root mean square error (RMSE) of ranging on the X-axis and Y-axis for the target 1 and the target 2 demonstrated in Figure 6 with σ_r =150 m (the original standard deviation of ranging is 150m). It can be observed from Fig. 7 there is only small difference for RMSE of ranging between LJPDA and JPDA, and the maximum difference is less than 1m. Extensive experiments reveal that average estimation error of JPDA and LJPDA are 73.36m and 72.74m, respectively (the average difference is less than 0.52), which are far less than the original standard deviation value of 150m. Consequently, estimation errors of JPDA and LJPDA are very close, and both JPDA and LJPDA exhibit good performance of target tracking.

Having discussed the estimation error of LJPDA, we now proceed to analyze the computation load of LJPDA proposed. In order to evaluate the computation efficiency, we run LJPDA and JPDA algorithm by using MATLAB, and compare the run time by using a general purpose CPU cores of i5-3210m operating at up to 2.5GHz with 8GB DDR3. Extensive simulations have been done by using JPDA and LJPDA, and curves of the average running time for a single-step versus the number of steps are plotted in Fig. 8. As shown in this figure, running time of LJPDA in a single step is less than 1ms in most cases, and the maximum running time is 3.2ms, which is far less than that for JPDA. Moreover, JPDA needs more than 7ms to finish data processing in most cases, and have to take more than 30ms for the 14th, 18th, 22th, 24th and 26th step. Additionally, average time consuming for LJPDA and JPDA is illustrated in the Table 2. On average, LJPDA need to take 0.997ms to perform data association whereas JPDA have to take 10.2ms on the same contexts, which have been verified by extensive numerical results. On the other hand, there is only tiny difference of average estimation error for LJPDA and JPDA. Obviously, the LJPDA algorithm proposed is much better than the JPDA algorithm.

Fig. 8. Average single-step time for data association

Table 2. The algorithm performance comparison

All in all, LJPDA proposed can track multiple targets in the clutter environment smoothly, and exhibits similar performance in tracking accuracy with JPDA. Particularly, LJPDA algorithm proposed has lower computation complexity than JPDA.

5. Conclusion

In this paper, a joint probabilistic data association algorithm (LJPDA) based on loss function is proposed. We divide data association into three cases, adopt contribution efficient to depict the contribution of a valid measurement on the associated target, and define the loss function tightly relevant to the estimation error and contribution coefficients. Then, we formulate the equation set for the optimal contribution coefficient by minimizing the loss function. Finally, the optimal contribution coefficient is attained by using Newton-Raphson method. Extensive experiments show that computation burden of the LJPDA algorithm proposed is much less than that of JPDA, and their estimation error is similar. Consequently, this new algorithm takes less run time with good performance, and is very helpful for data association of target tracking.