Next Article in Journal
An Energy Efficient Adaptive Sampling Algorithm in a Sensor Network for Automated Water Quality Monitoring
Next Article in Special Issue
Integrated Display and Simulation for Automatic Dependent Surveillance–Broadcast and Traffic Collision Avoidance System Data Fusion
Previous Article in Journal
Synergy Effect of Combining Fluorescence and Mid Infrared Fiber Spectroscopy for Kidney Tumor Diagnostics
Previous Article in Special Issue
A Novel Evidence Theory and Fuzzy Preference Approach-Based Multi-Sensor Data Fusion Technique for Fault Diagnosis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Centralized Multi-Sensor Square Root Cubature Joint Probabilistic Data Association

1
School of Electronic and Information Engineering, Beihang University, Beijing 100191, China
2
Research Institute of Information Fusion, Naval Aeronautical and Astronautical University, Yantai 264001, China
3
Department of Electronic Engineering, Tsinghua University, Beijing 100084, China
*
Authors to whom correspondence should be addressed.
Sensors 2017, 17(11), 2546; https://doi.org/10.3390/s17112546
Submission received: 8 September 2017 / Revised: 25 October 2017 / Accepted: 2 November 2017 / Published: 5 November 2017

Abstract

:
This paper focuses on the tracking problem of multiple targets with multiple sensors in a nonlinear cluttered environment. To avoid Jacobian matrix computation and scaling parameter adjustment, improve numerical stability, and acquire more accurate estimated results for centralized nonlinear tracking, a novel centralized multi-sensor square root cubature joint probabilistic data association algorithm (CMSCJPDA) is proposed. Firstly, the multi-sensor tracking problem is decomposed into several single-sensor multi-target tracking problems, which are sequentially processed during the estimation. Then, in each sensor, the assignment of its measurements to target tracks is accomplished on the basis of joint probabilistic data association (JPDA), and a weighted probability fusion method with square root version of a cubature Kalman filter (SRCKF) is utilized to estimate the targets’ state. With the measurements in all sensors processed CMSCJPDA is derived and the global estimated state is achieved. Experimental results show that CMSCJPDA is superior to the state-of-the-art algorithms in the aspects of tracking accuracy, numerical stability, and computational cost, which provides a new idea to solve multi-sensor tracking problems.

1. Introduction

The centralized state estimation method for multiple-target tracking with multiple sensors can integrate detecting information from different sensors according to certain rules, which will make the tracking results more accurate than that of a single sensor [1,2,3]. Multi-sensor fusion is an important data processing method in the field of target tracking, especially under nonlinear conditions, and the fused result with multiple sensors is often better for multiple-target tracking, which has been widely concerned by scholars both at home and abroad [3,4,5,6,7,8,9,10,11,12].
The state-of-the-art techniques for multiple-target tracking mainly contains two categories. The first category is the tracking method based on a random finite set; the other is the method based on measurement to track association. The former method can directly track multiple targets without associating the detected measurements with the tracks in the surveillance area. However, this kind of method faces complex integral operations, which are often difficult to find an exact solution. Although approximation methods have been a good choice for this problem, the computational cost is too great, and improper approximation will lead to degradation of the tracking accuracy, thus making it a long way from being extensively utilized in practical applications [4]. The latter method has been widely applied in various tracking systems, and is still the most common method for centralized tracking [5,6,7,8,9,10,11,12,13,14].
In the second category, the sequential multi-sensor joint probabilistic data association algorithm (MSJPDA) is a classical and effective method for multi-sensor tracking, which achieves better tracking performance compared with its parallel counterparts [10,11]. Unfortunately, current MSJPDA algorithms are often proposed to solve association problems under linear circumstances, and the nonlinear environment is seldom involved in the literature. As is known, the nonlinear equation can be linearized according to Taylor expansion, then a MSJPDA algorithm based on the extended Kalman filter (EKF), the so called MSJPDA-EKF is derived [12,13,14]. However, it is essential for MSJPDA-EKF to calculate the Jacobian matrix in the linearization step in each iteration, and the ignorance of higher orders may introduce large linearization errors, which may result in degraded performance, and even divergence [1,2,14]. In addition, the estimated state will affect the calculation of association probability in the next iteration, thus causing changes in the computation of weights of joint events, which may eventually makes false correlations, and even lead to divergence [15]. To conquer the shortcomings of MSJPDA-EKF, the MSJPDA algorithm based on the unscented Kalman filter (MSJPDA-UKF) is proposed. MSJPDA-UKF can acquire higher estimation accuracy than MSJPDA-EKF, and is relatively more stable [16]. However, reasonable adjustments are needed in MSJPDA-UKF to achieve the ideal performance, and a numerical instability phenomenon is easy to occur in the situation of high dimensional estimation [16,17,18].
To address the high dimensional nonlinear estimation problem, the cubature Kalman filter (CKF) based on the third-degree spherical-radial cubature rule has been recently proposed [19,20,21,22,23,24,25]. In the procedure of CKF, a certain number of typical cubature points are selected to approximate the posterior probability, then the mean and covariance are captured through a nonlinear transfer function to achieve accurate estimation of the target state [19]. CKF and UKF are both moment-matching filters which deterministically select a set of weighted sample points to approximate the posterior probability density. Compared with EKF, CKF and UKF do not have to calculate the Jacobian matrix in each iteration, thereby reducing the computational cost. At the same time, they both avoid the negative influence of the truncation error in the process of nonlinear state estimation, which is particularly effective in the case of state estimation with strong nonlinearity [19,21,22,23]. In UKF, to achieve a better estimation performance, the adjustment of the scaling parameter is an essential step, as different scaling parameters may yield completely distinct estimation results. The selection of the scaling parameter is of great significance to the performance of UKF in some sense [16,19,21,25]. However, there is no parameter to be adjusted for CKF, and the selected sampling points and weights are only related to the dimension of the target state, which can be calculated in advance to reduce the computational complexity. Furthermore, CKF is more accurate with respect to high dimensional estimation, and easy for design and implementation [19,21,22,23,24]. Especially in the case of high dimension, EKF and UKF suffer from the curse of dimensionality or divergence, or both, but CKF still performs well [19,21,25]. Unfortunately, in practice, arithmetic operations performed on finite word-length digital computers may introduce large errors; then the symmetry and positive definiteness of the covariance are often not satisfied. Moreover, perhaps the loss of positive definiteness is more perilous as it terminates the running of the CKF. Then, the square root version of CKF, the so called SRCKF is proposed to address this problem [19,21,23,24]. In SRCKF, only the square root factors are propagated during the iterations. Therefore, the square rooting operations of the matrix are avoided, thus reducing the computation cost. The symmetry and positive definiteness of the covariance are preserved, and the numerical accuracy is also improved [19,21,25].
In this paper, to deal with the target tracking problem in a nonlinear system under a cluttered environment, a novel centralized multi-sensor square root cubature joint probabilistic data association algorithm (CMSCJPDA) is proposed. In CMSCJPDA, measurements of each single sensor are sequentially processed to compute the association probabilities based on the similar rules in JPDA, then SRCKF and the validated measurements are chosen to estimate the state of targets according to weighed state fusion. Eventually, the accurate estimation of multiple targets’ states with multiple sensors is accomplished.
The rest of this paper is organized as follows: Section 2 describes the problem of tracking multiple targets with multiple sensors. Then, in Section 3, a simple overview of the numerical integral approximation method based on spherical-radial principle is introduced. The main idea of CMSCJPDA is detailed in Section 4, and an algorithm flow is given. The numerical experiments are designed in Section 5, and the comparison and analysis of CMSCJPDA against several existing methods is also presented in this section. Concluding remarks of this paper and future work are given in Section 6.

2. Problem Formulation

Consider sensors are applied to track N t targets in a cluttered environment. For arbitrary target t ( 1 t N t ) , X t ( k ) denotes the state of target t at discrete time instant k. For brevity, without the control input term, the discrete-time state equation of a nonlinear system is:
X t ( k + 1 ) = f t [ k , X t ( k ) ] + V t ( k )
where k is the discrete time instant, f t ( ) is a known nonlinear function, and V t ( k ) is independent zero-mean Gaussian process noise with covariance:
E { V t ( k ) [ V t ( l ) ] T } = Q t ( k ) δ ( k , l )
where δ ( k , j ) is the Kronecker delta function.
The measurement equation of sensor i for target t is:
Z i , t ( k ) = h i [ k , X t ( k ) ] + W i ( k )
where i = 1 , 2 , , N s represents the label of sensors. Z i , t ( k ) is the measurement vector of sensor i for target t at time instant k. h i ( ) represents a known nonlinear function. W i ( k ) is the zero-mean, independent of process noise in Equation (3) and independent from sensor to sensor, with covariance:
E { W i ( k ) [ W i ( l ) ] T } = R i ( k ) δ ( k , l )

3. Numerical Approximation of the Multi-Dimensional Weighted Integral

3.1. Third-Degree Spherical-Radial Rule

Before introducing the CMSCJPDA method, a brief description of the numerical integral approximation based on spherical-radial principle is given. Consider the following multi-dimensional Gaussian weighted integral:
I ( f ) = R n f ( x ) exp ( x T x ) d x
where f ( ) is an arbitrary function, R n denotes the integral region. As is known, the key to addressing the problem of nonlinear filtering based on Bayesian theory is to compute the first-order and second-order moments. In other words, the core of the Bayesian filter is how to compute Gaussian weighted integral which is of the form nonlinear function × Gaussian density that is illustrated in Equation (5).
Let x = r y , where y T y = 1 , So that x T x = r 2 for r [ 0 ,   ) . Then the integral Equation (5) can be rewritten as:
I ( f ) = 0 U n f ( r y ) r n 1 exp ( r 2 ) d σ ( y ) d r
where U n = { y R n | y T y = 1 } is the surface of the sphere and σ ( ) is the spherical surface measure or the area element on U n . Therefore, we can split Equation (6) into two integrals:
I ( f ) = 0 S ( r ) r n 1 exp ( r 2 ) d r
S ( r ) = U n f ( r y ) d σ ( y )
where I ( f ) is the radial weighted integral that is shown in Equation (7), S ( r ) is the spherical integral with the unit weighting function ω ( y ) = 1 that is shown in Equation (8).
Then the above integrals described in Equations (7) and (8) can be solved according to the third-degree spherical cubature rule and radical rule, respectively:
U n f ( y ) d σ ( y ) ω i = 1 2 n f [ u ] i
a b f ( x ) ω ( x ) d x i = 1 m ω i f ( x i )
where [ u ] i represents the i-th element of generator u , and u is a generator if u = ( u 1 , u 2 , , u r , 0 , 0 ) n , where u i > u i + 1 > 0 ,   i = 1 , 2 , , ( r 1 ) . For simplicity, the ( n r ) zero coordinates are ignored and notation [ u 1 , u 2 , , u r ] is used to represent a complete fully-symmetric set of points that can be obtained by permutating and changing the sign of generator u in all possible ways. For example, [ 1 ] 2 denotes the following set of points:
[ 1 ] = { ( 1 0 ) ,   ( 0 1 ) ,   ( 1 0 ) ,   ( 0 1 ) }
where ω ( x ) is a known weighting function which is non-negative on the interval [ a ,   b ] .
As is shown in Equation (9), the spherical cubature rule indicates that the spherical integral can be approximated by a weighted sum of function values at the sample points, which are located at the intersection of the unit sphere and its axes. The radical rule described in Equation (10) indicates that an m-point Gaussian quadrature is equal to the sum of (2m − 1) polynomials [21].
Especially, a standard Gaussian weighted integral can be computed based on third-degree spherical-radial rule as follows [25]:
I N ( f ) = R n f ( x ) N ( x   ;   0 , I ) d x i = 1 2 n x ω i f ( ξ i )
where:
ξ i = 2 n x 2 [ 1 ] i ,   ω i = 1 2 n x ,   i = 1 , 2 , , 2 n x
where ξ i is the cubature point, ω i is the corresponding weight. [ 1 ] i denotes the i-th row or column of the generator [ 1 ] .
Actually, the weighting function term ω ( ) of the integrand is often not subject to the standard Gaussian distribution. In other words, the Gaussian weighted integral in nonlinear filtering cannot be solved directly by Equation (12). Fortunately, the following equation makes it possible to work out the nonstandard Gaussian weighted integrals which occur in nonlinear condition:
R n f ( x ) N ( x ;   μ ,   Σ ) d x = R n f ( Σ x + μ ) N ( x ;   0 ,   I ) d x
Through Equation (14), which is proved in Appendix A, the nonstandard Gaussian weighted integrals can be transformed into the standard ones. Then, the approximation of posterior mean and error covariance can be addressed by Equation (12). Thus, an iteration of the time and the measurement updates in the Bayesian filter is completed.

3.2. Accuracy Analysis

Assume an n-dimensional vector x N ( x ¯ ,   P ) , where x ¯ is the mean of x , and P is the corresponding covariance. We extend the nonlinear function f ( x ) at x ¯ based on Taylor extension:
f ( x ) = f ( x ¯ ) + D Δ x f + D Δ x 2 f 2 ! + D Δ x 3 f 3 ! + D Δ x 4 f 4 ! +
where Δ x = x x ¯ , D Δ x f = [ Δ x ] T f ( x ) | x = x ¯ , f ( x ) is the partial derivative of f ( x ) .
(1)
We use the third-degree spherical-radial rule to approximate the multi-dimensional integral I ( f ) , then:
I CKF ( f ) = 1 2 n x i = 1 2 n x f ( x i ) = 1 2 n x i = 1 2 n x f ( x ¯ ) + D Δ x i f + D Δ x i 2 f 2 ! + D Δ x i 3 f 3 ! + D Δ x i 4 f 4 ! +
where x i is the sampled cubature point, Δ x = x i x ¯ . Through proper simplification, the final result can be written as:
I CKF ( f ) = f ( x ¯ ) + ( T P 2 ! ) f + n k 1 ( 2 k ) ! i = 1 n x ( a i 1 x 1 + + a i n x n ) 2 k f ( x ) | x = x ¯
where k = 1 , 2 , 3 , , a i j = [ P ] i j ,   i , j = 1 , 2 , , n , x i f is the partial derivative of f ( x ) at the ith component of x . The interested reader can be referred to paper [26] for details in the extension.
(2)
Like the derivation process of I CKF ( f ) , we use UKF to operate the approximation, and the final result is:
I UKF ( f ) = f ( x ¯ ) + ( T P 2 ! ) f + ( n + κ ) k 1 ( 2 k ) ! i = 1 n x ( a i 1 x 1 + + a i n x n ) 2 k f ( x ) | x = x ¯
where κ is the scaling parameter, and the remaining notations are the same as Equation (17). To capture the kurtosis of the prior density as correctly as possible, κ is suggested to be κ = 3 n x [16].
This paper is to tackle the high-dimensional estimation problem, so here we just discuss the case in which the state dimension n x > 3 . As is shown in Equations (17) and (18), the first two terms in both equations are the same, the only difference is in the last term. If n x > 3 , n k 1 > ( n + κ ) k 1 , and 2k is an even number, then I CKF ( f ) > I UKF ( f ) , which indicates that CKF is more accurate than UKF in high-dimensional estimation.
We choose the stability factor I = i | ω i | / i ω i defined as the measure of the numerical stability of the multi-dimensional integral, where ω i is the sampling weight. It is proven that the sampling rule implemented in a finite-precision arithmetic machine introduces a large amount of round-off errors when the stability factor I is larger than unity [19].
In UKF, if n x > 3 , and κ + n x = 3 , so κ = 3 n x < 0 , then ω 0 = 1 n x / 3 < 0 . The stability factor is I UKF = 2 n x 3 1 > 1 , and it scales linearly with state dimension n x , which indicates that UKF introduces significant perturbations in numerical estimates for the moment integral. However, in CKF, the stability factor I CKF always meets I CKF = 1 , which shows that CKF is more accurate than UKF for high-dimensional state estimation.

4. Centralized Multi-Sensor Square Root Cubature Joint Probabilistic Data Association

In CMSCJPDA, measurements of each sensor are processed in sequence, then the multi-sensor multi-target tracking problem can be reduced to several single-sensor multi-target tracking problems, which are easier to solve.
Assume there are m i ,   k t validated measurements for target t obtained by sensor i at time instant k. l i   ( 0 l i m i ,   k t ) is the label of validated measurements from sensor i , l i = j   ( j 0 ) represents the jth measurements in the validated region. Especially, l i = 0 indicates that there is no measurement in the validated region. Z l i i , t ( k ) represents the l i th validated measurement, Z ^ i | i 1 t ( k | k ) represents the predicted measurement, which will be discussed in detail later. β l i , i t ( k ) represents the association probability that the measurement l i originated from target t, K i t ( k ) is the filtering gain. The state estimate and the corresponding error covariance after processing the measurements of the ith sensor are denoted by X ^ i t ( k | k ) and P i t ( k | k ) , respectively. X ^ 0 t ( k | k ) and P 0 t ( k | k ) represent the initial estimation, and X ^ t ( k | k ) and P t ( k | k ) represent the final estimation at time instant k, respectively. Then, the state update is as follows:
X ^ i t ( k | k ) = X ^ i 1 t ( k | k ) + K i t ( k ) l i = 0 m i ,   k t β i ,   l i t ( k ) [ Z l i i , t ( k ) Z ^ i | i 1 t ( k | k ) ]
where:
{ X ^ 0 t ( k | k ) = X ^ t ( k | k 1 ) ,   X ^ t ( k | k ) = X ^ N s t ( k | k ) P 0 t ( k | k ) = P t ( k | k 1 ) ,   P t ( k | k ) = P N s t ( k | k )
Then, the update of the error covariance is:
P i t ( k | k ) = P i 1 t ( k | k ) [ 1 β 0 , i t ( k ) ] K i t ( k ) S i t ( k ) [ K i t ( k ) ] T + l i = 0 m i ,   k t β i ,   l i t ( k ) X ^ i ,   l i t ( k | k ) [ X ^ i ,   l i t ( k | k ) ] T X ^ i t ( k | k ) [ X ^ i t ( k | k ) ] T
Note that the key point for nonlinear single sensor multi-target tracking problem is to compute the association probability β i ,   l i t ( k ) and nonlinear state estimation X ^ i t ( k | k ) . JPDA is thought to be an effective method for multiple targets tracking with single sensor, so we choose JPDA to compute β i ,   l i t ( k ) here. The remaining problem is to obtain X ^ i t ( k | k ) . Although CKF is a good method to deal with the nonlinear state estimation problem, it is easily influenced by limited computer word-length and numerical errors, which may lead to loss of positive definiteness and symmetry of the error covariance matrix [19,21,22,23,24]. The effective way to preserve both properties and improve the numerical stability is to design a square root version of the CKF. Despite the fact that the square root cubature Kalman filter (SRCKF) is reformulated to propagate the square roots of the covariance matrices, both CKF and SRCKF are algebraically equivalent to each other, the interested reader is referred to [19]. In the proposed method, we choose SRCKF to update the state of the targets after the correlation step.

4.1. Computation of Association Probabilities

The probability β i ,   l i t ( k ) that measurement li originated from target t can be calculated by summing over all feasible events for which it is true [1]:
β l i , i t ( k ) = m = 1 n i t ( k ) ω ^ l i ,   t m ( θ m ( k ) ) Pr { θ m ( k ) | Z i t ( k ) }
where n i t ( k ) is the number of joint events, θ m ( k ) denotes the event m, and ω ^ l i ,   t m ( θ m ( k ) ) is the following binary element.
ω ^ l i ,   t m ( θ m ( k ) ) = { 1 , if   θ m ( k )   occurs 0 , else
Pr { θ m ( k ) | Z i t ( k ) } represents the posterior probability of event m, which can be computed by Equation (24):
Pr { θ m ( k ) | Z i t ( k ) } = 1 C   ϕ [ θ m ( k ) ] ! V ϕ [ θ m ( k ) ] l i = 1 m i ,   k t N l i i ,   t [ Z l i i , t ( k ) ] τ l i [ θ m ( k ) ] t = 1 N t ( P D t ) δ t [ θ m ( k ) ] ( 1 P D t ) 1 δ t [ θ m ( k ) ]
where C is a constant, ϕ [ θ m ( k ) ] is the total number of false measurements in event m. V is the volume of the entire surveillance region, and P D t is the detection probability of target t. The target detection indicator δ t [ θ m ( k ) ] is defined as:
δ t [ θ m ( k ) ] = l i = 1 m k i ω ^ l i ,   t m ( θ m ( k ) ) = { 1 ,   if   target   t   is   detected 0 ,   otherwise
which indicates whether any measurement is correlated with target t in event m. In the same way, it is also easy to define measurement association indicator:
τ l i [ θ m ( k ) ] = t = 1 N t ω ^ l i ,   t m ( θ m ( k ) ) = { 1 ,   if   measurement   l i   is   associated   with   a   target 0 ,   otherwise
which indicates whether measurement l i is correlated with a target in event m. N l i i ,   t [ Z l i i , t ( k ) ] is a conditional probability density function under Gaussian assumption, i.e.:
N l i i ,   t [ Z l i i , t ( k ) ] = 1 | 2 π S i t ( k ) | exp { 1 2 [ Z l i i , t ( k ) Z ^ i | i 1 t ( k | k ) ] T [ S i t ( k ) ] 1 [ Z l i i , t ( k ) Z ^ i | i 1 t ( k | k ) ] }

4.2. State Update

During the update of state, the data processing in first sensor is a little different from that in other sensors. Firstly, the specific procedures for the state update in sensors with label i > 1 are as follows:
Step 1: Calculate the cubature points (j = 1, 2, …, 2nx):
X ^   i | i 1 t ( k | k ) = X ^   i 1 | i 1 t ( k | k )
X j ,   i | i 1 t ( k | k ) = S i | i 1 t ( k | k ) ξ j + X ^ i | i 1 t ( k | k )
where nx is the dimension of the state. S i | i 1 t ( k | k ) is the square root factor of P i 1 t ( k | k ) , i.e.:
P i 1 t ( k | k ) = S i | i 1 t ( k | k ) [ S i | i 1 t ( k | k ) ] T
Equations (28) and (29) indicate that the estimated state X ^ i 1 | i 1 t ( k | k ) and error covariance P i 1 t ( k | k ) of the previous sensor are regarded as the predicted state X ^ i | i 1 t ( k | k ) and error covariance P i | i 1 t ( k | k ) of the next sensor, which is the greatest difference of the sensors with label i > 1 compared with the first sensor with label i = 1 in the process of state estimation. That is to say, the time update is not needed in sensors with label i > 1 , and the estimation of the previous sensor is used as the prediction of the next sensor.
Step 2: Evaluate the prediction of cubature points:
Z j , i | i 1 t ( k | k ) = h i [ k , X j ,   i | i 1 t ( k | k ) ]
Step 3: Estimate the prediction of the corresponding measurement:
Z ^ i | i 1 t ( k | k ) = 1 2 n x j = 1 2 n x Z j ,   i | i 1 t ( k | k )
Step 4: Estimate the square-root of the innovation covariance:
S z z ,   i | i 1 t ( k | k ) = Tria ( [ Z i | i 1 t ( k | k ) S R i ( k ) ] )
where S z z ,   i | i 1 t ( k | k ) is a lower triangular matrix. S R i ( k ) denotes a square root factor of R i ( k ) that R i ( k ) = S R i ( k ) [ S R i ( k ) ] T , and the weighed, centred matrix:
Z i | i 1 t ( k | k ) = 1 2 n x [ Z 1 ,   i | i 1 t ( k | k ) Z ^ i | i 1 t ( k | k ) ,   Z 2 ,   i | i 1 t ( k | k ) Z ^ i | i 1 t ( k | k ) ,   ,   Z 2 n x ,   i | i 1 t ( k | k ) Z ^ i | i 1 t ( k | k ) ]
Note that S = Tria ( A ) denotes the QR decomposition, where S is a lower triangular matrix. The relationship of matrices S and A is as follows: Let B be an upper triangular matrix obtained from the QR decomposition on A T . Then, the lower triangular matrix S is obtained as S = R T .
Step 5: Estimate the cross-covariance matrix:
P x z , i | i 1 t ( k | k ) = X i | i 1 t ( k | k ) [ Z i | i 1 t ( k | k ) ] T
where:
X i | i 1 t ( k | k ) = 1 2 n x [ X 1 ,   i | i 1 t ( k | k ) X ^ i | i 1 t ( k | k ) , X 2 ,   i | i 1 t ( k | k ) X ^ i | i 1 t ( k | k ) ,   ,   X 2 n x ,   i | i 1 t ( k | k ) X ^ i | i 1 t ( k | k ) ]
Step 6: Estimate the filter gain:
K i t ( k ) = P x z , i | i 1 t ( k | k ) / S z z ,   i | i 1 t ( k | k ) / [ S z z ,   i | i 1 t ( k | k ) ] T
Step 7: Update the nonlinear state:
X ^ l i ,   i t ( k | k ) = X ^ i 1 t ( k | k ) + K i t ( k ) [ Z l i i , t ( k ) Z ^ i | i 1 t ( k | k ) ]
Equations (28)–(38) give a solution to the estimation of the target state and error covariance at time instant k with sensor label i > 1 .
Now the estimation problem in the first sensor with label i = 1 is considered, the process of the measurement update is the same as that in sensors with label i > 1 , but the time update is quite different. In the sensor with label i = 1 , the estimation of prediction of the state and error covariance is essential. The current cubature points are:
X j t ( k 1 | k 1 ) = S t ( k 1 | k 1 ) ξ j + X ^ t ( k 1 | k 1 )
The propagated cubature points are:
X j t ( k | k 1 ) = f t [ k 1 , X j t ( k 1 | k 1 ) ]
Then the predicted state is:
X ^ t ( k | k 1 ) = 1 2 n x j = 1 2 n x X j t ( k | k 1 )
and the predicted error covariance is:
P t ( k | k 1 ) = 1 2 n x j = 1 2 n x [ X j t ( k | k 1 ) X ^ t ( k | k 1 ) ] [ X j t ( k | k 1 ) X ^ t ( k | k 1 ) ] T + Q t ( k )
Equations (39)–(42) describe the time update in the first sensor with label i = 1 . With the measurement update process composed of Equations (28)–(38) applied, the estimated state X ^ 1 t ( k | k ) and error covariance P 1 t ( k | k ) can be achieved respectively after the data of the first sensor is processed.
To state the proposed method more clearly, the algorithm flow of CMSCJPDA is illustrated in Table 1.

5. Numerical Experiments

In this section, the numerical results and the effectiveness of the proposed CMSCJPDA algorithm are discussed in two simulation scenarios: crossing target tracking and maneuvering target tracking. The performance of CMSCJPDA is compared with MSJPDA-EKF and MSJPDA-UKF in two scenarios.
For a fair comparison, we make 50 independent Monte Carlo runs. In each run, the measurement noises are generated independently, and the corresponding noisy position measurements are used to initialize the estimate, which guarantees consistency of the initialization of the filter and makes the final estimate unbiased. Except for the initialization, the process of each Monte Carlo run is the same, and the mean value of 50 runs is considered as the final result. The total number of steps per run is 100. T = 1   s is the sampling interval. The nonparametric model is used for the probability mass function (PMF) of the number of false measurements. The expected number of false measurements in the validation gate is m = 2 . The detection probabilities of all targets are assumed to be P D = 0.9 . The probability mass is assumed to be P G = 0.9997 . Due to the large initial error, we display the numerical results after the tenth step for clarity.
Targets in the surveillance area are observed by three two-dimensional sensors. The measurement equation of sensor i is:
Z i ( k ) = [ ( x ( k ) x p i ) 2 + ( y ( k ) y p i ) 2 arctan ( y ( k ) y p i x ( k ) x p i ) ] + W i ( k )
where i = 1 , 2 , 3 is the sensor label. ( x ( k ) ,   y ( k ) ) is the true state of the targets, ( x p i ,   y p i ) is the position of sensor i. Other parameters of the sensors are set as shown in Table 2.
Performance metrics: To compare the tracking performance of different algorithms, the root mean square error (RMSE) is chosen as the metric. The root mean square position error of target t at time instant k is defined as:
RMSE pos ( k ,   t ) = 1 M i = 1 M ( x i t ( k ) x ^ i t ( k | k ) ) 2
where M is the total number of Monte Carlo simulations. x i t ( k ) and x ^ i t ( k | k ) are true and estimated positions of target t in the x direction in the n-th Monte Carlo run. Similarly, we may also define formulas of the root mean square velocity error.

5.1.Crossing Targets Tracking

5.1.1. Simulation Scenario

Consider a two-crossing target tracking case. The state model is:
X ( k + 1 ) = F ( k ) X ( k ) + Γ ( k ) V ( k )
where the component of process noise q 1 = q 2 = 0.01 . The state transition matrix F ( k ) is defined as:
F ( k ) = [ 1 T 0 0 0 1 0 0 0 0 1 T 0 0 0 1 ]
and the process noise distribution matrix is:
Γ ( k ) = [ T 2 / 2 0 T 0 0 T 2 / 2 0 T ]
The initial states of two targets are X1(0) = [−29,500 m, 400 m/s, 34,500 m, 2,212,400 m/s]T, X2(0) = [−26,500 m, 296 m/s, 34,500 m, −400 m/s]T. After about 31 s, the two targets cross at the point (−17,000 m, 22,000 m).

5.1.2. Results and Analysis

The true and filtered tracks of two crossing targets are illustrated in Figure 1, which shows the tracking performance of three algorithms in a single run.
To give a more explicit picture of the estimated results, the root mean square position error of three data association algorithms in the x and y directions is shown in Figure 2 and Figure 3, while the root mean square velocity error in the x and y directions is shown in Figure 4 and Figure 5. As is shown, the three algorithms can effectively estimate the state of both targets. However, CMSCJPDA and MSJPDA-UKF, which are based on deterministic sampling, are of higher tracking accuracy than MSJPDA-EKF, and CMSCJPDA is slightly more accurate than MSJPDA-UKF, which just validates the argument in [19]. With the square root version of CKF applied in the proposed method, the square rooting operations of covariance matrix is avoided and the tracking performance is effectively enhanced.
To further validate the effectiveness of the proposed algorithm, numerical stability and computational complexity of the three algorithms are compared in Table 3, which gives the average divergence times (ADT) and average time consumption (ATC) in 50 Monte Carlo runs. As is shown, the ADT of MSJPDA-EKF is twice more than that of MSJPDA-UKF, and nearly four times more than that of CMSCJPDA, to some extent CMSCJPDA is relatively more stable in terms of convergence. As for the computational cost, the ATC of CMSCJPDA is almost half that of MSJPDA-EKF, and less than that of MSJPDA-UKF. In a word, compared with MSJPDA-EKF and MSJPDA-UKF, CMSCJPDA can achieve relatively higher tracking accuracy, and is more stable and less computationally complex under a nonlinear cluttered environment.

5.2. Maneuvering Targets Tracking

5.2.1. Simulation Scenario

Assume there are two maneuvering targets in the surveillance area. The initial states of two targets are X1(0) = [20,000 m, −600 m/s, 1800 m, 500 m/s]T, X2(0) = [4000 m, 600 m/s,1800 m, 200 m/s]T. The two targets move at a constant speed in a straight line for 30 s. Then, from 30 s to 50 s, they both move with constant accelerations of a1x = −30 m/s2, a1y = 20 m/s2, a2x = 30 m/s2, a2y = 20 m/s2 in the x and y directions, respectively. From 50 s to 70 s, they still move with constant accelerations, and the accelerations for both targets in the x and y directions are ax1 = −30 m/s2, ay1 = −20 m/s2, ax2 = 30 m/s2, ay2 = −20 m/s2, respectively. After 70 s, they both move at a constant speed again until 100 s. The true trajectories of both targets are illustrated in Figure 6.

5.2.2. Results and Analysis

The root mean square (RMS) position error of two targets is compared in Figure 7 and Figure 8. The results show that the estimation error of CMSCJPDA and MSJPDA-UKF is much smaller than that of MSJPDA-EKF in the whole simulation. When t is about 40 s, the estimation error of MSJPDA-EKF is about 200 m larger compared with CMSCJPDA and MSJPDA-UKF. Moreover, with the new association strategy for maneuvering target tracking, CMSCJPDA has relatively higher accuracy than MSJPDA-UKF. At some time instants, the RMS position error of CMSCJPDA is about 50 m smaller than that of MSJPDA-UKF, and it converges relatively faster.
Table 4 evaluates the tracking performance of the three algorithms in the aspects of mean RMS position error (MRMSE) and correct association rate (CAR). The results indicate that both CMSCJPDA and MSJPDA-UKF are more accurate than MSJPDA-UKF, and CMSCJPDA is slightly more accurate than MSJPDA-UKF. Meanwhile, the CAR of CMSCJPDA is the highest among the three algorithms, which is 11.1% higher than MSJPDA-UKF and 28.7% higher than MSJPDA-EKF, and CMSCJPDA is more time-saving.

6. Conclusions

In this paper, we propose a centralized multi-sensor square root cubature joint probabilistic data association algorithm, which is based on the spherical-radial principle and the sequential update scheme. Compared with the state-of-the-art method in two tracking scenarios, the proposed method is much more accurate than MSJPDA-EKF, and a slightly better than MSJPDA-UKF by using the square root version of cubature Kalman filter to estimate the target state. Furthermore, the experimental results also show that the proposed method is less prone to divergence than MSJPDA-EKF and MSJPDA-UKF, and is more stable with respect to its numerical characteristics. Although the proposed method is less time consuming among the compared methods, due to complicated calculation of the association probability, the joint probabilistic data association-based methods are of significant computational cost, which is still less practical in real applications. Currently, the sensors in our experiment are assumed to be synchronized sampling, and there is no time difference in measurements from different sensors. However, the actual applications can rarely meet this condition, and in the future we will pay more attention to reducing the computational cost and improving the real-time performance. At the same time, the tracking problem with asynchronous sampling and changing the number of targets also deserves consideration.

Acknowledgments

The authors give their sincere thanks to the editors and the anonymous reviewers for their constructive comments of the manuscript. This study was financially co-supported by the State Key Program for Basic Research of China (no. 613XXXXX), the National Natural Science Foundation of China (nos. 61471383, 91538201, 61531020, and 61671463), the Young Elite Scientist Sponsorship Program by CAST (YESS) (no. 2015QNRC001), and China Postdoctoral Science Foundation (2017M610032).

Author Contributions

All the listed authors contributed to the research design; Y.L., J.L., G.L. and L.Q. collected the data; Y.L., J.L. and G.L. analyzed the data; Y.L. and J.L. contributed to the formula derivation; Y.L. and G.L. designed the simulation experiments; J.L. and L.Q. designed the figures and tables; Y.L. searched the literature; Y.L., J.L. and G.L. wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Proof of Equation (14)

For simplicity, define the weighting function ω 1 ( x ) = N ( x ;   0 ,   I ) and ω 2 ( x ) = N ( x ;   μ ,   Σ ) . Consider the left side of Equation (14). Since Σ is a positive definite matrix, then Σ can be factorized to be Σ = Σ ( Σ ) T . Making a variable substitution by x = Σ y + μ , we have:
R n f ( x ) ω 2 ( x ) d x = R n f ( Σ y + μ ) 1 | 2 π Σ | exp { ( Σ y ) T Σ 1 ( Σ y ) } | Σ | d y = R n f ( Σ y + μ ) 1 | 2 π I | exp ( y T y ) d y = R n f ( Σ y + μ ) N ( y ;   0 ,   I ) d y = R n f ( Σ x + μ ) ω 1 ( x ) d x
Substituting ω 1 ( x ) and ω 2 ( x ) into Equation (A1), we have:
R n f ( x ) N ( x ;   μ ,   Σ ) d x = R n f ( Σ x + μ ) N ( x ;   0 ,   I ) d x

References

  1. Bar-Shalom, Y. Multitarget-Multisensor Tracking: Advanced Application; Artech House: Norwood, MA, USA, 1990. [Google Scholar]
  2. He, Y.; Wang, G.; Lu, D.; Peng, Y. Multisensor Information Fusion with Applications, 2nd ed.; Publishing House of Electronics Industry: Beijing, China, 2010. [Google Scholar]
  3. Zhang, K.; Shan, G.; Ji, B.; Chen, H. Maneuvering targets tracking based on pose-aided nonlinear filter algorithm. Acta Electron. Sin. 2012, 40, 1670–1675. [Google Scholar]
  4. Liu, J.; Liu, Y.; He, Y.; Sun, S. Multi-sensor fuzzy data association based on azimuth and velocity. Syst. Eng. Electron. 2016, 38, 2040–2047. [Google Scholar]
  5. Ge, Q.; Wei, Z.; Cheng, T.; Chen, C.; Wang, X. Flexible fusion structure-based performance optimization learning for multisensory target tracking. Sensors 2017, 17, 1045. [Google Scholar] [CrossRef] [PubMed]
  6. Liu, Z.; Pan, Q.; Jean, D.; Gregoire, M. Hybrid Classification System for Uncertain Data. IEEE Trans. Syst. Man Cybern. Syst. 2016, 10, 2783–2790. [Google Scholar] [CrossRef]
  7. Taghavi, E.; Tharmarasa, R.; Kirubarajan, T.; Bar-Shalom, Y.; Mcdonald, M. A practical bias estimation algorithm for multisensor-multitarget tracking. IEEE Trans. Aerosp. Electron. Syst. 2016, 52, 2–19. [Google Scholar] [CrossRef]
  8. Fantacci, C.; Papi, F. Scalable multisensor multitarget tracking using the marginalized-glmb densit. IEEE Signal Process. Lett. 2016, 23, 863–867. [Google Scholar] [CrossRef]
  9. Si, W.; Wang, L.; Qu, Z. A measurement-driven adaptive probability hypothesis density filter for multitarget tracking. Chin. J. Aeronaut. 2015, 28, 1689–1698. [Google Scholar] [CrossRef]
  10. ONeil, S.D.; Pao, L.Y. Multisensor fusion algorithm for tracking. In Proceedings of the American Control Conference, San Francisco, CA, USA, 2–4 June 1993; pp. 859–863. [Google Scholar]
  11. Pao, L.Y.; Frei, C.W. A comparison of parallel and sequential implementation of a multisensory tracking algorithm. In Proceedings of the American control Conference, Seattle, WA, USA, 21–23 June 1995; pp. 1683–1687. [Google Scholar]
  12. Bar-shalom, Y. Multitarget-Multisensor Tracking: Principles and Techniques; YBS Publishing: Stors, CT, USA, 1995. [Google Scholar]
  13. Zhang, J.; He, Y.; Xiong, W. Centralized interacted multiple model multisensory fuzzy joint probabilistic data association algorithm. Acta Electron. Sin. 2008, 36, 1655–1659. [Google Scholar]
  14. Han, C.; Zhu, Y.; Duan, Z.; Han, D.; Liu, W.; Yu, X. Multi-Source Information Fusion, 2nd ed.; Tsinghua University Press: Beijing, China, 2010. [Google Scholar]
  15. Guan, X.; Zhou, X.; Rui, G. Centralized multisensor unscented joint probabilistic data association algorithm. Syst. Eng. Electron. 2009, 31, 2602–2606. [Google Scholar]
  16. Julier, S.J.; Uhlmann, J.K. Unscented filtering and nonlinear estimation. Proc. IEEE 2004, 92, 401–422. [Google Scholar] [CrossRef]
  17. Gao, X.; Chen, J.; Tao, D.; Li, X. Multi-sensor centralized fusion without measurement noise covariance by variational Bayesian approximation. IEEE Trans. Aerosp. Electron. Syst. 2011, 47, 718–727. [Google Scholar] [CrossRef]
  18. Wang, X.; Pan, Q. Overview of deterministic sampling filtering algorithms for nonlinear system. Control Decis. 2012, 27, 801–812. [Google Scholar]
  19. Arasaratnam, I.; Haykin, S. Cubature Kalman filters. IEEE Trans. Autom. Control 2009, 54, 1254–1269. [Google Scholar] [CrossRef]
  20. Arasaratnam, I.; Haykin, S. Cubature Kalman smoothers. Automatica 2011, 47, 2245–2450. [Google Scholar] [CrossRef]
  21. Liu, Y.; He, Y.; Wang, H. Squared-root cubature information consensus filter for non-linear decentralised state estimation in sensor networks. IET Radar Sonar Navig. 2014, 8, 931–938. [Google Scholar]
  22. Leong, P.H.; Arulampalam, S.; Lamahewa, T.A.; Abhayapala, T.D. A Gaussian-sum based cubature Kalman filter for bearings-only tracking. IEEE Trans. Aerosp. Electron. Syst. 2013, 49, 1161–1174. [Google Scholar] [CrossRef]
  23. Liu, Y.; Dong, K.; Wang, H.; Liu, J.; He, Y.; Pan, L. Adaptive Gaussian sum squared-root cubature Kalman filter with split-merge scheme for state estimation. Chin. J. Aeronaut. 2014, 27, 1242–1250. [Google Scholar] [CrossRef]
  24. Liu, G.; Florentin, W.; Irene, M. Square-root sigma-point information filtering. IEEE Trans. Autom. Control 2012, 57, 2945–2950. [Google Scholar] [CrossRef]
  25. Jia, B.; Sun, T.; Xin, M. Iterative diffusion-based distributed cubature Gaussian mixture filter for multisensory estimation. Sensors 2016, 16, 1741. [Google Scholar] [CrossRef] [PubMed]
  26. Sun, F.; Tang, L. Estimation precision comparison of cubature Kalman filter and unscented Kalman filter. Control Decis. 2013, 28, 303–308. [Google Scholar]
Figure 1. True tracks and filtered tracks of two crossing targets.
Figure 1. True tracks and filtered tracks of two crossing targets.
Sensors 17 02546 g001
Figure 2. Root mean square position error of target 1.
Figure 2. Root mean square position error of target 1.
Sensors 17 02546 g002
Figure 3. Root mean square position error of target 2.
Figure 3. Root mean square position error of target 2.
Sensors 17 02546 g003
Figure 4. Root mean square velocity error of target 1.
Figure 4. Root mean square velocity error of target 1.
Sensors 17 02546 g004
Figure 5. Root mean square velocity error of target 2.
Figure 5. Root mean square velocity error of target 2.
Sensors 17 02546 g005
Figure 6. The trajectories of two cross-maneuvering targets.
Figure 6. The trajectories of two cross-maneuvering targets.
Sensors 17 02546 g006
Figure 7. Root mean square position error of target 1.
Figure 7. Root mean square position error of target 1.
Sensors 17 02546 g007
Figure 8. Root mean square position error of target 2.
Figure 8. Root mean square position error of target 2.
Sensors 17 02546 g008
Table 1. Algorithm flow of CMSCJPDA.
Table 1. Algorithm flow of CMSCJPDA.
Input: The estimated state X ^ t ( k 1 | k 1 ) and error covariance P t ( k 1 | k 1 ) at time instant k 1 ;
Output: The estimated state X ^ t ( k | k ) and error covariance P t ( k | k ) at time instant k ;
Step 1. Time update ( i = 1 )
 Compute the predicted state X ^ t ( k | k 1 ) and error covariance P t ( k | k 1 ) according to Equations (39)–(42);
Step 2. Initialization
X ^ 0 t ( k | k ) = X ^ t ( k | k 1 ) , P 0 t ( k | k ) = P t ( k | k 1 ) ;
Step 3. Measurement update
 For t = 1 : N t
  For i = 1 : N s
   Compute the estimation of state X ^ i t ( k | k ) and error covariance P i t ( k | k ) with the validated measurements of sensor i according to Equations (19)–(38);
  End
 End
Step 4. Establish the final estimation
X ^ t ( k | k ) = X ^ N s t ( k | k ) , P t ( k | k ) = P N s t ( k | k 1 )
Table 2. Sensor position and parameter setting.
Table 2. Sensor position and parameter setting.
Sensor LabeSensor Position (m)Ranging Error (m)Angle Error (rad)
1(0, 0)1000.01
2(−500, −500)2000.02
3(−500, 500)3000.03
Table 3. Performance comparison of the three algorithms.
Table 3. Performance comparison of the three algorithms.
AlgorithmsAverage Divergence TimesAverage Time Consumption (s)
MSJPDA-EKF0.870.423
MSJPDA-UKF0.420.346
CMSCJPDA0.270.257
Table 4. Performance comparison of three algorithms.
Table 4. Performance comparison of three algorithms.
AlgorithmsMRMSE (m)CAR (%)TC (s)
MSJPDA-EKF206.247.60.735
MSJPDA-UKF132.665.20.563
CMSCJPDA104.376.30.432

Share and Cite

MDPI and ACS Style

Liu, Y.; Liu, J.; Li, G.; Qi, L.; Li, Y.; He, Y. Centralized Multi-Sensor Square Root Cubature Joint Probabilistic Data Association. Sensors 2017, 17, 2546. https://doi.org/10.3390/s17112546

AMA Style

Liu Y, Liu J, Li G, Qi L, Li Y, He Y. Centralized Multi-Sensor Square Root Cubature Joint Probabilistic Data Association. Sensors. 2017; 17(11):2546. https://doi.org/10.3390/s17112546

Chicago/Turabian Style

Liu, Yu, Jun Liu, Gang Li, Lin Qi, Yaowen Li, and You He. 2017. "Centralized Multi-Sensor Square Root Cubature Joint Probabilistic Data Association" Sensors 17, no. 11: 2546. https://doi.org/10.3390/s17112546

APA Style

Liu, Y., Liu, J., Li, G., Qi, L., Li, Y., & He, Y. (2017). Centralized Multi-Sensor Square Root Cubature Joint Probabilistic Data Association. Sensors, 17(11), 2546. https://doi.org/10.3390/s17112546

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop