Auxiliary Truncated Unscented Kalman Filtering for Bearings-Only Maneuvering Target Tracking

Novel auxiliary truncated unscented Kalman filtering (ATUKF) is proposed for bearings-only maneuvering target tracking in this paper. In the proposed algorithm, to deal with arbitrary changes in motion models, a modified prior probability density function (PDF) is derived based on some auxiliary target characteristics and current measurements. Then, the modified prior PDF is approximated as a Gaussian density by using the statistical linear regression (SLR) to estimate the mean and covariance. In order to track bearings-only maneuvering target, the posterior PDF is jointly estimated based on the prior probability density function and the modified prior probability density function, and a practical algorithm is developed. Finally, compared with other nonlinear filtering approaches, the experimental results of the proposed algorithm show a significant improvement for both the univariate nonstationary growth model (UNGM) case and bearings-only target tracking case.


Introduction
Bearings-only maneuvering target tracking has been widely researched for decades. It is important for many applications such as maritime surveillance, navigation and aerospace, wireless sensor networks (WSN), and infrared search and track (IRST) systems [1][2][3][4][5][6]. However, while implementing this technology in unlimited situations, there remain some challenging problems, such as multiple platform tracking, uncertainty of the target model and nonlinear non-Gaussian noise. To deal with the uncertainty of the motion model, such as abrupt target maneuver, heavy clutter measurements, highly nonlinearity of dynamic models and nonlinear non-Gaussian noise, etc., the interacting multiple model (IMM) [7] based on the nonlinear filtering algorithm is a promising approach. However, to model the uncertainty of the motion model, the performance of the IMM-type algorithm is directly proportional to the number of the motion models. Generally, the more motion models we produce, the greater accuracy of the estimated state we obtain. However, the computational complexity of the algorithm becomes larger with the increase of the numbers of motion models, particularly in heavily cluttered environments. Moreover, the nonlinear filtering has been studied extensively in bearings-only maneuvering target tracking.
As is well-known, the most widely used nonlinear filtering for bearing-only tracking is to employ an extended Kalman filter (EKF) [8,9]. However, when the nonlinearity of dynamic models becomes more severe, the performance of the EKF degrades sharply. In order to solve this problem, the unscented Kalman filter (UKF) [10] and the truncated unscented Kalman filtering (TUKF) were proposed [11,12]. Compared with other conventional Kalman filter-type approaches, the TUKF can achieve better performance in the conditions of the target tracking system, and can provide very The rest of the paper is organized as follows. The proposed algorithm is given in Section 2. In Section 3, we provide the experimental results. Finally, some conclusions are given in Section 4.

Proposed Algorithm
In order to track the maneuvering target, accurate motion modeling and nonlinear filtering are two challenging problems that should not be separated. However, most research on maneuvering tracking investigates these problems separately. In this section, in the bearings-only maneuvering target tracking, novel auxiliary truncated unscented Kalman filtering is proposed. In Section 2.1, the joint prior distribution is approximately constructed. In Section 2.2, in order to track the bearings-only maneuvering target, the modified prior PDF is approximated based on statistical linear regression by introducing the target spatio-temporal information. Section 2.3 summarizes the proposed algorithm.

Joint Prior Distribution
Suppose the target dynamic system can be written as: where z n ∈ R n z denotes the observation vector at time n, x n ∈ R n x denotes the target state vector at time n, and f (·) and h(·) denote the nonlinear state transition function and observation function, respectively. Suppose that r n denotes the set of target characteristics including c independent components r n = r 1 n , r 2 n , . . . , r c n . In order to derive the proposed algorithm, there are two basic hypotheses, firstly, that the nonlinear function h n (·) in (2) is a bijective, continuous function; and secondly, that the probability density function of the measurement noise e n has bounded, connected support. p e n (v n ) = 0, e n / ∈ I e n ⊂ n z where I e n is an n z -dimensional measurement validation region. Therefore, according to the second assumption, the measurement likelihood function can be defined as follows: p( z n |x n , r n ) = p e n (z n − h n (x n ))χ I en (z n − h n (x n )) (4) p(z n x n , r n ) = p e n (z n − h n (x n ))χ I xn (z n ) (x n ) where χ I en (·) is the indicator function on the subset I e n . Therefore, the state posterior PDF can be defined as: p(x 0:n |z 1:n , r 1:n ) = p(z n |x 0:n ,z 1:n−1 ,r 1:n )p(x n |x 1:n−1 ,z 1:n−1 ,r 1:n )p(x 1:n−1 |z 1:n−1 ,r 1:n−1 ) p(z n ,r n |z 1:n−1 ,r 1:n−1 ) = p e k (z n −h(x n ))χ Ix (zn ) (x n )p(x n |xn−1,z1:n−1,r1:n)p(x1:n−1|z1:n−1,r1:n−1) p(z k ,r k |z 1:k−1 ,r 1:k−1 ) ∝ p(z n |x n )p 1 (x n |z n , x n−1 , r 1:n )p(x 1:n−1 |z 1:n−1 , r 1:n−1 ) p 1 (x n z n , x n−1 , r 1:n ) = p(x n x n−1 , z 1:n−1 , r 1:n )χ I xn (z n ) (x n )/ε 1 (8) where ε 1 is a constant. From Equation (8), we can see that the modified prior PDF is defined by incorporating the current measurement information z n . According to the conclusions in [11], if the measurement noise is informative, the modified prior p 1 (·) is not only the minimum variance of the prior p 0 (·), but also can improve the algorithm's performance. Further, to deal with the uncertainty of motion models, the joint prior distribution p(x n |x 0:n−1 , z 1:n , r 1:n ) of the proposed algorithm can be defined as follows: p(x n |x 0:n−1 , z 1:n , r 1:n ) = α n p 1 (x n |z n , x n−1 , r 1:n ) + (1 − α n )p 0 (x n |x n−1 , z 1:n−1 ) = α n p(x n x n−1 , r 1:n )χ I xn (z n ) (x n ) + (1 − α n )p 0 (x n x n−1 , z 1:n−1 ) (9) where α n ∈ [0, 1] is a proper parameter. To approximately calculate the mean and covariance of the posterior distribution, we apply a UKF update to p 0 (·) (UKF1), and another UKF update to p 1 (·) (UKF2). Finally, the posterior estimates can be approximately calculated through merging both results obtained by UKF1 and UKF2.

Approximation of p 1 (·)
In the subsection, our object is to approximate the modified prior PDF p 1 (·) as Gaussian. For this reason, we write the state vector as x n = [a n T , b n T ] T , where a n ∈ n a denotes the position components of x n , b n ∈ n b denotes the velocity components of x n , and n x = n a + n b . The derivation of the mean x p 1 ,n and covariance P p 1 ,n of p 1 (·) is the same as in the truncated unscented Kalman filter (TUKF) [11], which can be shown as follows: x p 1 ,n = µ a n ,1 µ b n ,1 , P p 1 ,n = ∑ a n ,1 ∑ a n b n ,1 (∑ a n b n ,1 ) T ∑ b n ,1 Σ a n , where µ b n ,1 , Σ b n ,1 , Σ a n b n ,1 can be found in [11]. µ a n ,1 denotes the estimated mean of state a n , R n denotes the measurement noise covariance, and H −1 n = [∇ a n h T n (a n )] T | a n =µ an,1 is the Jacobian of h T n (a n ) evaluated at µ a n ,1 . Now, how to calculate the estimated mean µ a n ,1 remains a key problem to be solved. For the passive sensor tracking system, the modeling of target dynamic system is a challenging problem when the target maneuvers, and some auxiliary target characteristics need to be used to deal with arbitrary changes in motion models. To achieve a high tracking performance, a statistical linear regression method (SLR) [25] is proposed to estimate the state mean µ a n ,1 .
Firstly, to evaluate the state mean µ a n ,1 , three approximations are used: (S1) the prior PDF p 0 (a n ) is constant over the connected region I a n (z n ); (S2) the nonlinear function h(·) can be locally linearized; and (S3) the measurement noise satisfies uniform distribution e n ∼ U I en in the connected region I e n . According to S2, the nonlinear measurement Equation (2) can be approximated as a linear estimator of z n ,ẑ n such that:ẑ n = n a n + d n (12) where n denotes a linear measurement matrix, and d n denotes a noise vector, which are derived by minimizing the objective function defined as follows: where τ n is the linearization error, τ n = z n −ẑ n . Substituting τ n into (13), and setting the partial derivative of the objective function with d n to zero, (−2)E(z n − n a n − d n ) = 0 ⇔ d n = z n − n a n (14) where a n = E(x n ) and z n = E(z n ). Substituting d n into (13) τ T τ = [(z n − z n ) − n (a n − a n )] T [(z n − z n ) − n (x n − x n )] Then, setting the gradient with respect to n to zero, (−2)E{[(z n − z n ) − n (a n − a n )][z n − z n ])} = 0 (15) Solving (15) for n , we obtain: n = P T a n z n P −1 a n a n (16) where P a n a n = E[(a n − a n )(a n − a n ) T ], P a n z n = E[(a n − a n )(z n − z n ) T ]. According to the measurement equation (12), the maximum likelihood positionâ n (z n ) of target state a n can be estimated as follows: If the maximum likelihood estimateâ n (z n ) is used to replace the estimated mean µ a n ,1 , the performance of the proposed algorithm will be improved because the current measurement information is incorporated. However, it cannot solve the uncertainty in motion models. In particular, when the target speed or the measurement sampling interval is large, the tracking performance degrades.
More recently, sophisticated techniques have been based on the target motion characteristics [6], which have been proposed to address the challenges in motion modeling. In their proposed method, three target characteristics, such as actual target speed v, time interval T of measurement and course angle θ of the target, are considered to improve the tracking performance. In (18), the relationship between three target characteristics and the state predicted error is given: where v n denotes the current estimated velocity, ∇θ denotes the estimated error of course angle, and v denotes the actual target velocity. Supposing ∇v = v − v n , we can obtain From Equation (19), we can find that the predicted error ∇σ increases monotonically with the increase of the parameters (v n , T, ∇θ). In fact, when ∇v n > 50 m/s and T > 20 s, the predicted error ∇σ will larger than 1000 m. It shows that the predicted error becomes a major reason for the performance degradation. On the other hand, when the actual target speed is relatively small, the prediction error caused by T or ∇θ is smaller than the measurement error. Therefore, to improve the performance of mean estimate µ a n ,1 , the maximum likelihood estimateâ n (z n ) is considered as the latest observation, and the modified maximum likelihood estimates that can incorporate current target characteristics is defined as follows: where λ is a constant, σ 2 e (n) denotes the variance of measurement noise, and σ 2 v (n) denotes the innovation variance. According to (20) and (16), the mean µ a n ,1 and the variance Σ a n ,1 in (10) can be approximated as: µ a n ,1 =φ(z n ), Σ a n ,1 = P a n a n (21) Finally, the modified prior PDF p 1 (·) can be approximated as a Gaussian probability density function p 1 (·) ≈ N(x n ,x p 1 ,n , P p 1 ,n ).

Summary of the Proposed Algorithm
According to the descriptions above, in order to describe clearly the proposed algorithm, the diagram of the ATUKF is shown in Figure 1. In Figure 1, it is shown that one cycle of the ATUKF algorithm consists of the following steps: (A) time update (predicted by using Kalman filtering); (B) the measurement updates based on the prior p 0 (·) and the modified prior p 1 (·); and (C) weight calculation and the joint state update. According to the derived results mentioned above and Figure 1, the detailed information of the new ATUKF algorithm can be summarized as follows. n n n n a a a n a Finally, the modified prior PDF ) ( 1  p can be approximated as a Gaussian probability density

Summary of the Proposed Algorithm
According to the descriptions above, in order to describe clearly the proposed algorithm, the diagram of the ATUKF is shown in Figure 1. In Figure 1, it is shown that one cycle of the ATUKF algorithm consists of the following steps: (A) time update (predicted by using Kalman filtering); (B) the measurement updates based on the prior    0 p and the modified prior    1 p ; and (C) weight calculation and the joint state update. According to the derived results mentioned above and Figure 1, the detailed information of the new ATUKF algorithm can be summarized as follows. Figure 1. Auxiliary truncated unscented Kalman filtering. Obtain N = 2n x + 1 sigma points χ 1 0 , χ 2 0 ..., χ N 0 and the corresponding associated weights w 1 , w 2 , ..., w N using unscented transform (UT ) based on the meanx n−1|n−1 and covariance P n−1|n−1 of the posterior PDF p 0 (x n |x n−1 , z 1:n−1 , r 1;n−1 ) , where n x denotes the dimension of state x. The predicted sigma points can be obtained by the nonlinear state function f (·):

2.
Approximate the mean and covariance of the state-predicted prior PDF p 0 (x n |x n−1 , z 1:n−1 ) 3.
Compute the predicted measurementẑ 0,n|n−1 based on the nonlinear measurement function h(·):

4.
The cross-covariance, innovation covariance and error covariance are estimated as follows:
Approximate the meanx n and covariance P n of the posterior PDF p(x n |z 1:n ) using (23) and (24).

Simulation Results
In this section, to evaluate the tracking performance of the ATUKF algorithm, two examples are employed. In Section 3.1, the univariate nonstationary growth model (UNGM) is discussed [11]. In Section 3.2, a bearings-only maneuvering tracking scenario [24], interested in defense applications, is discussed. In the first case, the EKF, UKF, the quadrature Kalman filtering (QKF) [25], the mixture truncated unscented Kalman filter (MTUKF, with three mixture components) [12] and particle filtering(PF) are utilized. In the second case, the TQKF, interacting multiple model extended Kalman filtering(IIMMEKF) and the interacting multiple model Rao-Blackwellized particle filter (IMMRBPF) are employed. In all the experiments, each simulation has been repeatedly performed 100 times.

Univariate Nonstationary Growth Model (UNGM)
In this section, due to the highly nonlinearity and non-stationarity of dynamic system, the univariate nonstationary growth model is considered. The discrete time system of this model can be written as: where the process noise m n is satisfied with Gaussian distribution with zero mean and variance one, and e n is satisfied with Gaussian distribution with zero mean and variance 0.01. α = 0.5, β = 25, γ = 8, φ 1 = 0.2 and φ 2 = 0.05 are known constants. In each Monte Carlo simulation, we assume that the initial distribution of state x 0 is uniform distribution in the interval [0 1]. The number of particles is 1000. Figure 2 shows the root-mean-square position errors (RMSE) of the EKF, UKF, QKF, PF, MTUKF and ATUKF. It is obvious from Figure 2 that the performance of the ATUKF is much better than that of the EKF, UKF and QKF. In this case, the performance of the EKF is the poorest. A reason for the poor performance of the EKF, UKF and QKF is the increase of the approximate error arising from the high Sensors 2017, 17, 972 9 of 14 nonlinearity and the non-stationarity of the dynamic system. Figure 3 shows the RMS position errors of all filters with different noise variance σ e n ∼ [0. 1 5]. From Figure 3, it is seen that whenever the measurement is informative (σ e n < 1) or the measurement is uninformative (σ e n > 1), the ATUKF is robust in all situations, its performance is similar to the MTUKF's, and it is very close to that of the PF. Moreover, among the EKF, UKF and QKF, the EKF has the poorest performance. In particular, when the measurement is very informative (σ e n < 1), the EKF yields a divergent estimate. number of particles is 1000. Figure 2 shows the root-mean-square position errors (RMSE) of the EKF, UKF, QKF, PF, MTUKF and ATUKF. It is obvious from Figure 2 that the performance of the ATUKF is much better than that of the EKF, UKF and QKF. In this case, the performance of the EKF is the poorest. A reason for the poor performance of the EKF, UKF and QKF is the increase of the approximate error arising from the high nonlinearity and the non-stationarity of the dynamic system. Figure 3 shows the RMS position errors of all filters with different noise variance   Table 1 shows the computation time statistics for all algorithms. In this case, all the experiments are performed by using MATLAB programming on Intel-Core(TM)-i2-4030U processor (1.9 GHz) based on the Windows platform. It can be seen from Table 1 that the computational load of the PF is the largest than these of other filters, such as the EKF, UKF, QKF, MTUKF and ATUKF. The ATUKF is very close to the QKF. Furthermore, the computation time for the ATUKF is much lower than the MTUKF. The main reason is that the MTUKF approximates the posterior PDF as a Gaussian mixture, and it makes the computational burden increase.

Bearings-Only Maneuvering Tracking (BOT) Scenario
In this scenario, the target makes two circular turns with rectilinear segments connecting them. Figure 4 shows the true target trajectory. The speed modulus is kept constant throughout    Table 1 shows the computation time statistics for all algorithms. In this case, all the experiments are performed by using MATLAB programming on Intel-Core(TM)-i2-4030U processor (1.9 GHz) based on the Windows platform. It can be seen from Table 1 that the computational load of the PF is the largest than these of other filters, such as the EKF, UKF, QKF, MTUKF and ATUKF. The ATUKF is very close to the QKF. Furthermore, the computation time for the ATUKF is much lower than the MTUKF. The main reason is that the MTUKF approximates the posterior PDF as a Gaussian mixture, and it makes the computational burden increase.

Bearings-Only Maneuvering Tracking (BOT) Scenario
In this scenario, the target makes two circular turns with rectilinear segments connecting them. Figure 4 shows the true target trajectory. The speed modulus is kept constant throughout (0.3 km/s). The initial position is (2 km, 8 km, 1 km), and the initial velocity is (0.15 km/s, 0.26 km/s, 0.0 km/s). The segments are defined as follows: where T denotes the sampling interval, in this paper, T is set to 1.

Model 2: Constant Turn Motion
The state transition matrix is: where w is a constant angular rate. In this paper, w is set to 0175 . 0 . The process noise covariance matrix 2 Q is the same as in Model 1.  The motion model of target is defined as follows: where x s n = (x, y, z, x , y , z ) denotes the target state vector under model x, x, y and z denote the target position coordinates, x , y and z are the target speed in x, y and z directions, respectively, F s denotes the transition matrix, and s ∈ {1, 2, . . . , M} denotes the target model index. In the maneuvering target tracking scenario, only a constant velocity model is used for the ATUKF algorithm and TQKF algorithm. In the IMMEKF and IMMRBPF, there are both clockwise-and counterclockwise-coordinated turn models that are used to simulate the target maneuvering. The details of three target motion models are defined as follows:

Model 1: Constant Velocity Motion
The state transition matrix and the process noise covariance matrix are defined by: where T denotes the sampling interval, in this paper, T is set to 1.

Model 2: Constant Turn Motion
The state transition matrix is: where w is a constant angular rate. In this paper, w is set to 0.0175. The process noise covariance matrix Q 2 is the same as in Model 1.
Model 3: w > 0 describes a clockwise turn, and Model 3 is its natural counterpart for a counterclockwise turn w < 0.
Two passive sensors are deployed in (0, −5 km, 0) and (0, 5 km, 0) respectively. Using the detection fusion architecture [24], the azimuth and elevation angles of aircraft, α i and β i respectively, measured by sensor i, are transmitted to the fusion node. The measurement function is written as: where (s i,x , s i,y , s i,z ), i = 1, 2 denote the positions of the stationary sensors. The measurement covariance can be defined as R n = 1 0 0 1 σ 2 e n .
The real initial position of the target is (2 km, 8 km, 1 km), and the initial velocity is (0.15 km/s, 0.26 km/s, 0.0 km/s). The prior PDF of state x 0 is assumed to be x 0 ∼ N x 0|0 ,P 0|0 , wherex 0|0 = [2.1 km 00. 12 kms −1 7.95 km 0. 23 kms −1 0. 95 km 0 kms −1 ] T ,P 0|0 = diag[0.144 km 2 0. 02 2 km 2 s −2 0. 144 km 2 0. 02. 2 km 2 s −2 0 0] T . The standard deviation of the process and the measurement noise are set to σ m n = 0.01 and σ e n = 0.001, and these noises are zero-mean Gaussian-distributed and independent. The number of particles is set to 200. Figure 5 shows the X RMSE, Y RMSE, Z RMSE and position RMSE of the ATUKF compared with IMMEKF, IMMRBPF and TQKF. From Figure 5, we can see that the RMSE of all algorithms abruptly increased from 30 s to 50 s, which is mainly due to the increase of the target predicted errors caused by the target maneuvering. However, in Figure 5a,c,d, it is shown that the performance of the ATUKF algorithm has outperformed the IMMEKF, IMMRBPF and TQKF. A key reason is that the proposed algorithm can incorporate the target characteristic information and current measurement information into the prior PDF, which can effectively degrade the variance of errors caused due to the maneuvering of the target. Moreover, because the flight height of target remained unchanged, from Figure 5c The real initial position of the target is (2 km, 8 km, 1 km) , and the initial velocity is (0.15 km/s, 0.26 km/s, 0.0 km/s) . The prior PDF of state 0 x is assumed to be ( ) where  Figure 5 shows the X RMSE, Y RMSE, Z RMSE and position RMSE of the ATUKF compared with IMMEKF, IMMRBPF and TQKF. From Figure 5, we can see that the RMSE of all algorithms abruptly increased from 30 s to 50 s, which is mainly due to the increase of the target predicted errors caused by the target maneuvering. However, in Figure 5a,c,d, it is shown that the performance of the Finally, the computation time statistics for all algorithms are given in Table 2. In this case, all the experiments are performed by using MATLAB programming on an Intel-Core(TM)-i2-4030U processor (1.9 GHz) based on the Windows platform. In Table 2, it is shown that the computational load of the IMMRBPF is much higher than these of the IMMEKF, TQKF and ATUKF. More importantly, the ATUKF requires much less of a computation time than the TQKF. However, the computational load of the ATUKF is nearly two times higher than that of the IMMEKF.

Conclusions
In this paper, we presented a bearings-only target tracking algorithm based on an auxiliary truncated unscented Kalman filtering (ATUKF) algorithm. Unlike the truncated unscented Kalman filtering, in the proposed algorithm, several target characteristics were introduced to construct the modified prior PDF, and the statistical linear regression was used to linearize the nonlinear non-bijective measurement function by using the sigma points. Moreover, we have developed a practical algorithm for a bearings-only target tracking system. Finally, in the simulation results, compared with the EKF, UKF, the quadrature Kalman filtering (QKF), the mixture truncated unscented Kalman filter (MTUKF) and the particle filter (PF), the ATUKF exhibits better performance. For the second case, compared with the IMMEKF algorithm, the IMMRBPF algorithm and the TQKF algorithm, the ATUKF algorithm not only improves the performance of the tracker, but significantly reduces the computation time.