Next Article in Journal
Buried Pipeline Collapse Dynamic Evolution Processes and Their Settlement Prediction Based on PSO-LSTM
Next Article in Special Issue
Agile Attitude Maneuver Control of Micro-Satellites for Multi-Target Observation Based on Piecewise Power Reaching Law and Variable-Structure Sliding Mode Control
Previous Article in Journal
Analysis of Neural Networks Used by Artificial Intelligence in the Energy Transition with Renewable Energies
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Iterated Orthogonal Simplex Cubature Kalman Filter and Its Applications in Target Tracking

School of Space Information, Space Engineering University, Beijing 101416, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2024, 14(1), 392; https://doi.org/10.3390/app14010392
Submission received: 30 October 2023 / Revised: 23 December 2023 / Accepted: 27 December 2023 / Published: 31 December 2023
(This article belongs to the Special Issue Recent Advances and Application of Image Processing)

Abstract

:
In order to increase a nonlinear system’s state estimate precision, an iterated orthogonal simplex cubature Kalman filter (IOSCKF) is presented in this study for target tracking. The Gaussian-weighted integral is decomposed into a spherical integral and a radial integral, which are approximated using the spherical simplex-radial rule and second-order Gauss–Laguerre quadrature rule, respectively, and result in the novel simplex cubature rule. To decrease the high-order error terms, cubature points with appropriate weights are taken from the cubature rule and processed using the provided orthogonal matrix. The structure supporting the nonlinear Kalman filter incorporates the altered points and weights and the calculation steps; from this, the updated time and measurement can be inferred. The Gauss–Newton iteration is employed repeatedly to adjust the measurement update until the termination condition is met and the IOSCKF is attained. The proposed algorithms are applied in target tracking, including CV target tracking and spacecraft orbit tracking, and the simulation results reveal that the IOSCKF can achieve higher accuracy compared to the CKF, SCKF, and OSCKF. In spacecraft orbit tracking simulation, compared with the SCKF, the position tracking accuracy and velocity tracking accuracy of the OSCKF are increased by 2.21% and 1.94%, respectively, which indicates that the orthogonal transformation can improve the tracking accuracy. Furthermore, compared with the OSCKF, the position tracking accuracy and velocity tracking accuracy of the IOSCKF are increased by 2.71% and 2.97%, respectively, which indicates that the tracking accuracy can be effectively improved by introducing iterative calculation into the measurement equation, thus verifying the effectiveness of the method presented in this paper.

1. Introduction

For several decades, target tracking has had substantial applications in many fields, such as navigation guidance, military application, and sensor networks [1,2]. Target tracking is a process whereby the position and velocity of a target’s motion are estimated using to discrete measurements with random noise output by sensors. Mathematical models of target tracking are generally described in different coordinate systems, where the process model is established in the Cartesian coordinate system and the measurements are generally implemented in the spherical coordinate system [3,4,5]. Thus, target tracking can be considered to be a typical nonlinear filtering problem [6,7]. The posterior probability density function (pdf) is very important in nonlinear filtering theory. Because pdf includes every piece of information on the state vector, it provides an integral filtering solution. Therefore, nonlinear filtering should focus on the calculation of posterior pdf, as the posterior mean and covariance of the state cannot be propagated using the nonlinear function directly.
In nonlinear systems, there contain two types of approaches to estimate the posterior pdf; one is optimal estimation, while the other is suboptimal estimation [8]. The former mainly includes a particle filter [9] and makes no assumptions on the pdf, which leads to a large amount of computation and is not suitable for applications in engineering. For the latter method, the posterior pdf is assumed to have a Gaussian distribution, and the typical representation, the nonlinear Kalman filter, is primarily studied in practical applications [10].
In a nonlinear Kalman filter, the recursive estimation of the state vector is divided into two processes consisting of a time update and a measurement update, and a unified filtering calculation framework is given. The key problem in the filtering framework is in the approximation of the Gaussian-weighted integral of multi-dimensional nonlinear function [8], and two approximate approaches, including the approximation of the nonlinear function and that of the Gaussian pdf, are mainly taken. The extended Kalman filter (EKF) [11] belongs to the former, and the basic idea of EKF is to use the first-order Taylor expansion to linearize the nonlinear system equations, and then, the conventional Kalman filter is used to estimate the state [12]. It is required to calculate the Jacobian matrix of the system, which demands that the model be smooth and differentiable. The EKF filtering accuracy is reduced or even divergent for a highly nonlinear system, and an improved algorithm, named iterated EKF (IEKF) [13], is proposed. The latest measurements are fully utilized through iterative calculation in the measurement update of IEKF to improve the estimation accuracy; however, the inherent disadvantage of EKF is still present in IEKF, that is, the differential operation of the Jacobian matrix. In general, the approximation of pdf is easier than that of the nonlinear function, so the approximation is performed using the deterministic sampling points in the latter method, and results in the unscented Kalman filter (UKF) [14], Gauss–Hermite–Kalman filter (GHKF) [15], sparse-grid Kalman filter (SGKF) [16], and cubature Kalman filter (CKF) [17,18]. These filters are all derivative-free filters, that is, there is no need to calculate the Jacobian matrix to linearize the nonlinear equation [19]. Specifically, a set of sigma points is used to approximate the posterior mean and covariance in UKF, and Chang [20] proposes using the orthogonal transformation of the sigma points for improving the approximation accuracy without increasing the computational complexity, thus reducing the high-order error term. However, there exist three adjustable parameters in UKF, and the inappropriate choice of parameters may reduce the filtering accuracy and the numerical stability. The computation amount of GHKF increases exponentially with the system state, which may cause the curse of dimensionality [16]. In order to reduce the computational complexity, the sparse-grid integral rule can be utilized to derive the SGKF, and the number of integral points is increased polynomially [17]. However, the CKF is widely used in engineering [21,22,23,24]. The CKF helps the intractable Gaussian-weighted integral to decompose and is composed of a spherical integral and a radial integral, and a set of cubature points with equal weights are deduced [25]. It can achieve more filtering accuracy and numerical stability compared with UKF. Similar to the IEKF, the principle of iteration can also be expanded to these filters’ measurement update processes [26,27,28,29].
Based on the CKF, Jia [30] puts forward the high-degree cubature Kalman filter (HDCKF) by using arbitrary degree symmetric spherical rule and moment matching method, which achieve higher accuracy than the CKF. And then, Wang [31] employs a transformation group of regular simplexs to approximate the spherical integral, and proposes the spherical simplex-radial cubature Kalman filter (SSRCKF). Compared with the HDCKF, the SSRCKF has more accurate performance, whereas the moment matching method cannot guarantee optimality when it is used in the abovementioned approaches’ radial integrals. Furthermore, Singh [32] uses the high-order Gauss–Laguerre quadrature rule instead of the moment matching method to calculate the radial integral, and proposes the high-degree cubature quadrature Kalman filter (HDCQKF), which could improve the radial integral accuracy effectively, so as to enhance the filtering accuracy. The advantage of the high-degree cubature Kalman filter is that it can increase the filtering accuracy, although it will increase the computational complexity and reduce the numerical stability for the high-dimensional system, and thus limit the application of filters.
To further improve the nonlinear filtering accuracy without using the high-degree cubature rule, a new iterated orthogonal simplex cubature Kalman filter (IOSCKF) is proposed in this paper. The way to calculate the intractable integral is to combine the spherical simplex rule with the second-order Gauss–Laguerre quadrature rule. By introducing the orthogonal matrix to reduce the higher-order moments of the Taylor series expansion which appear during the approximation of the mean and covariance, the results in the novel cubature points with corresponding weights are improved. The time update is the method used to reduce the computational complexity for the linear process equation. It is equivalent to a simplification of the time update of the standard linear Kalman filter. The posterior state and covariance matrix, in the measurement update, perform Gaussian Newton iteration until the result is optimal. On target tracking problems, the IOSCKF’s performance was tested, and the simulation results verify the effectiveness of the proposed algorithm.
The remainder of this paper is organized as follows: a review of Gaussian approximation filter is described in Section 2, the improved spherical simplex-radial cubature rule is deduced in Section 3, the iterated orthogonal simplex cubature Kalman filter is derived in Section 4, the target tracking simulations and results are presented in Section 5, and the conclusions are drawn in Section 6.

2. Brief Review of Gaussian Approximation Filter

The following discrete nonlinear system model is defined:
x k = f x k 1 + w k 1 z k = h x k + v k
where k is the time step,  x k R n  denotes the state vector at time step k z k R p  represents the measurement vector at time step k, and  f  and  h  are the known nonlinear state function and measurement function, respectively. The process noise  w k  and measurement noise  v k  are uncorrelated zero-mean Gaussian white noise, with the covariance matrices being  Q k  and  R k , respectively.
The key problem in the nonlinear Kalman filter is calculating the multi-dimensional vector integral of the form “nonlinear function × Gaussian pdf” [8], that is, calculating the integral  I N = R n g x N x ; x ¯ , P x d x  in the Cartesian coordinate system, where  g x  denotes the arbitrary nonlinear function, and  N x ; x ¯ , P x  represents the Gaussian pdf, with the mean and covariance being  x ¯  and  P x , respectively.
This Gaussian-weighted integral  I N  can be approximated using various quadrature rules, and the integral with respect to  N ( x ; 0 , I )  can be approximated by the following quadrature rule:
I N = R n g x N x ; 0 , I d x i = 1 N p ω i g σ i
where I is the identity matrix,  N p  denotes the total number of points, and  σ i  and  ω i  represent the quadrature points and weights, respectively.
The integral with respect to the more general Gaussian distribution  N x ; x ¯ , P x  can be obtained using the equation below:
I N = R n g x N x ; x ¯ , P x d x = R n g x ¯ + P x x N x ; 0 , I d x i = 1 N p ω i g x ¯ + P x σ i
where  P x  denotes the square root of  P x , which can be solved by the Cholesky decomposition, and  x ¯ i = x ¯ + P x σ i  represents the integral point.
The Gaussian approximation filter, which is deduced using the numerical integration method in Equation (3), is summarized below [30].
Time Update:
The predicted state  x ^ k  and error covariance  P k  are calculated, respectively, as follows:
x ^ k = i = 1 N p ω i f x ^ k 1 i
P k = i = 1 N p ω i f x ^ k 1 i x ^ k f x ^ k 1 i x ^ k T + Q k 1
where  x ^ k 1 i = x ^ k 1 + + P k 1 + σ i  stands for the integral points generated using the numerical rule.
Measurement Update:
The updated state  x ^ k +  and corresponding error covariance  P k +  are estimated, respectively, as follows:
x ^ k + = x ^ k + K k z k z ^ k
P k + = P k + K k P z K k T
 where  K k = P x z P z 1  denotes the Kalman filtering gain,  z ^ k = i = 1 N p ω i h x ^ k i  represents the predicted measurement,  x ^ k i = x ^ k + P k σ i  denotes the integral points generated by the numerical rule,  P z = i = 1 N p ω i h x ^ k i z ^ k h x ^ k i z ^ k T + R k  denotes the innovation covariance matrix, and  P x z = i = 1 N p ω i x ^ k i x ^ k h x ^ k i z ^ k T  represents the cross-covariance matrix.

3. The Improved Spherical Simplex-Radial Cubature Rule

In this section, an improved spherical simplex-radial cubature rule is proven to be more accurate than the traditional conventional spherical-radial rule. Before introducing the improved cubature rule, the definition of the dth-degree rule is given first for the integral
Definition 1
([30]).
R n g x ω g x d x i ω i g γ i
where  x = x 1 x 2 x n T R n  and  ω g x  is a given weighing function. Equation (8) is a dth-degree rule if it is accurate for  g x , whose components are linear combinations of monomials  x 1 α 1 x 2 α 2 x n α n   with a total degree up to d.  α 1   ,  α 2   ,    ,  α n  are nonnegative integers and  α 1 + α 2 + α n d , and there is at least one monomial of degree  d + 1   for which Equation (8) is not exact.
Specifically, an intractable integral of the form  I g = R n g x e x T x d x  is considered first in the cubature rule. Let  x = r y  with  y T y = 1 , where  y  means the direction vector and the radius  r 0 . As a consequence,  x T x = r 2 ; then, the integral  I g  can be written anew as in a spherical-radial coordinate system as follows [12]:
I g = 0 U n g r y r n 1 e r 2 d σ y d r
where  U n  is the sphere surface defined by  U n = y R n : y T y = 1 , and  σ  is the area element on  U n . The integral  I g  is decomposed into the spherical integral  S r = U n g r y d σ y  and radial integral  R = 0 S r r n 1 e r 2 d r , respectively.
Because of the difficulty in achieving the analytical solutions, the approximate solution method of  I g  requires numerical integration, which should be approximated using the numerical integration method. From the third-degree spherical-radial cubature rule proposed in [17], we obtain
S r i = 1 N s ω s , i g r y i
R j = 1 N r ω r , j S r j
and  I g  is approximated using the following formula:
I g = 0 r n 1 e r 2 d r U n g r y d σ y j = 1 N r i = 1 N s ω r , j ω s , i g r j y i
where  y i  and  ω s , i  denote the points and weights to calculate the spherical integral, and the number of points is  N s r j  and  ω r , j  represent the points and weights to calculate the radial integral, and the number of points is  N r . The total number of points is  N r N s  if  r j 0 , or  N r 1 N s + 1  if one of  r j  is zero.

3.1. Spherical Simplex Rule

The proposed third-degree spherical rule and spherical simplex rule to calculate Equation (10) are listed below, respectively.
The third-degree spherical rule is
S s r i = 1 N s ω s , i g r y i = A n 2 n i = 1 n g r e i + g r e i
where  A n = 2 π n / Γ n / 2  denotes the surface area of the unit sphere, and  Γ z = 0 e t t z 1 d t  is the Gamma function. The vector  e i  represents the unit column vector with the ith element being 1. The number of cubature points needed is  N s = 2 n .
The third-degree spherical simplex rule [31] consisting of  N s = 2 n + 1  points is
S s s r = A n 2 n + 1 i = 1 n + 1 g r a i + g r a i
where  a j = a j , 1 a j , 2 a j , n T , j = 1 , 2 , , n + 1  denote the vertices of the n-dimensional simplex, and the elements are defined as follows.
a j , i = n + 1 n n i + 2 n i + 1 , i < j n + 1 n j + 1 n n j + 2 , i = j 0 , i > j
Although the aforementioned two spherical rules share the same accuracy of the third degree, through extensive simulations, it has been proven that spherical simplex rules have higher accuracy than spherical rules [31], and that is the reason the spherical simplex rule is chosen to approximate  S r  in this paper. For more details, please refer to the related references.

3.2. Radial Rule

Then, the calculation of the radial integral  R  is considered. Two methods, including the moment matching method and Gauss–Laguerre quadrature rule, may be used. To construct the (2m + 1)th-degree radial rule, the two methods lead to identical points and weights when m is odd and  m + 1 / 2  points are used. However, they may be different when m is even and  m / 2 + 1  points are used [30]. In the CKF and SSRCKF, the radial integral is calculated using the former one, and results in the following third-degree radial rule with only one point (that is, m = 1,  N r = 1 ).
R = 1 2 Γ n 2 S n 2
In this paper, the second-order Gauss–Laguerre quadrature rule is adopted to approximate the radial integral. This method improves the accuracy of radial integration.
For radial integral  R = 0 S r r n 1 e r 2 d r , let  r 2 = t  to obtain  r = t , which takes the place of the radial integral, and for  R = 1 2 0 S t t n 2 2 e t d t , let  q t = S t , and  β = n 2 1 . Then,  R  is transformed into  R = 1 2 0 q t t β e t d t , and this form can be handily approximated using the Gauss–Laguerre quadrature [32] rule as follows:
R = 1 2 0 q t t β e t d t 1 2 j = 1 m A j q t j
where  t j  denotes the quadrature points, with  A j  being the corresponding weights.
The approximation accuracy depends on the order m. The m-order Chebyshev–Laguerre polynomial’s solutions are the quadrature points, as shown in the following:
L m β t = 1 m t β e t d m d t m t β + m e t = 0
and the weights are calculated using the following equation.
A j = m ! Γ β + m + 1 t j L ˙ m β t j 2
On account of m being set at 2, the solutions are obtained are  t 1 = n / 2 + 1 + n / 2 + 1 t 2 = n / 2 + 1 n / 2 + 1 A 1 = n Γ n / 2 / 2 n + 4 + 2 2 n + 4 , and  A 2 = n Γ n / 2 / 2 n + 4 2 2 n + 4 , which are substituted into Equation (17) to deduce the following radial rule.
R = n Γ n 2 4 n + 8 + 4 2 n + 4 S n 2 + 1 + n 2 + 1 + n Γ n 2 4 n + 8 4 2 n + 4 S n 2 + 1 n 2 + 1
From Equation (20) we see that the number of points is  N r = 2 , which leads to the radial rule to achieve the fifth-degree accuracy. Compared with third-degree radial rule used in the CKF and SSRCKF Equation (16), the fifth-degree accuracy is more accurate. It is noted that the two points can be solved using the moment matching method as well, with one of the points set to zero. However, it is pointed out that the two points generated by the Gauss–Laguerre quadrature rule have higher radial accuracy than those generated using the moment matching method [30].

3.3. Spherical Simplex-Radial Cubature Rule

By substituting Equation (14) into Equation (20), we obtain the following improved spherical simplex-radial cubature rule which is composed of  N p = N r N s = 4 ( n + 1 )  points.
I N = R n g x N x ; x ¯ , P x d x = n 4 n + 1 n + 2 + 2 n + 4 i = 1 n + 1 g x ¯ + n + 2 + 2 n + 4 P x a i + g x ¯ n + 2 + 2 n + 4 P x a i + n 4 n + 1 n + 2 2 n + 4 i = 1 n + 1 g x ¯ + n + 2 2 n + 4 P x a i + g x ¯ n + 2 2 n + 4 P x a i
If we define matrix  a = a 1 a 2 a n + 1 , and the calculation methods of the cubature points and weights are extracted from cubature rule (21), we obtain the following.
x ¯ i = x ¯ + n + 2 + 2 n + 4 P x a a i , i = 1 , 2 , , 2 n + 2 x ¯ + n + 2 2 n + 4 P x a a i , i = 2 n + 3 , , 4 n + 4
ω i = n 4 n + 1 n + 2 + 2 n + 4 , i = 1 , 2 , , 2 n + 2 n 4 n + 1 n + 2 2 n + 4 , i = 2 n + 3 , , 4 n + 4
The order of the quadrature rule determines the accuracy of the cubature rule. The relationship of the number of quadrature points and the accuracy is proportional [8,32]. The radial rule proposed in Section 3.2 uses two points to achieve fifth-degree accuracy. In contrast to the third degree adopted by the CKF and SSRCKF, when only one point is used, fifth-degree is more accurate. And the spherical simplex rule performs better than the spherical rule; therefore, the improved spherical simplex-radial cubature rule, which is shown in Equation (21), can achieve higher accuracy than the conventional CKF and SSRCKF.

4. Iterated Orthogonal Simplex Cubature Kalman Filter

4.1. Orthogonal Transformation of Cubature Points

It is suggested from the invariant theory and the orthogonal transformation property that the transformed points  O × a a , with weights shown in Equation (23), remain to constitute a numerical cubature rule if  O  is an  n × n  orthogonal matrix. The orthogonally transformed points are able to reduce the higher-order terms of the Taylor series expansion that appeared while approximating the mean and covariance [20], and it is observed that more accurate estimation is obtained with the transformed points compared to the ordinary cubature points, which lead to the approximation accuracy of the cubature rule becoming further enhanced. In other words, the new set of cubature points can be modified as follows:
x ¯ i = x ¯ + n + 2 + 2 n + 4 P x γ i , i = 1 , 2 , , 2 n + 2 x ¯ + n + 2 2 n + 4 P x γ i , i = 2 n + 3 , , 4 n + 4
where  γ = γ 1 γ 2 γ 2 n + 2 = O × a a , and the orthogonal matrix  O  can be designed using different methods. As an effective design method, the  n × n  matrix  O  is constructed in this paper as  O = O 1 O 2 O n , and  O i = ( O i , 1 O i , 2 O i , n ) T , whose elements are defined as follows:
O i , 2 r 1 = 2 n cos 2 r 1 i π n O i , 2 r = 2 n sin 2 r 1 i π n , r = 1 , 2 , , max p Z p n 2
where  O i , n = 1 i / n  if  n  is odd.
The information of high-order terms can be reduced to a low level using the orthogonal matrix, which means the non-local sampling problem can be significantly reduced with a rare computational complexity increase.

4.2. Iterated Measurement Update Process

The prior state  x ^ k  and covariance matrix  P k  are used in the measurement update process to calculate the posterior state  x ^ k +  and covariance matrix  P k + , which are revised using the latest measurement  z k  via Kalman gain  K k . This means an approximation accuracy of  x ^ k +  to the real state is bound to be better than  x ^ k , and the filtering accuracy will be improved if we can use  x ^ k +  and  P k +  directly in the measurement update. Therefore, to iteratively make use of the measurement information, the Gauss–Newton iteration is introduced into the updating process of the mean and covariance of the state. Under the assumption of Gaussian variable independence, there exist  x k : N x ^ k , P k  and  z k : N h x k , R k , and the likelihood function of  x k  and  z k  is defined as follows [33].
p x k z k = C exp 1 2 x k x ^ k T P k 1 x k x ^ k + z k h x k T R k 1 z k h x k
where  C  is a constant, and the cost function is defined as follows.
J x k = x k x ^ k T P k 1 x k x ^ k + z k h x k T R k 1 z k h x k
The maximum likelihood estimate of Equation (26) is equivalent to the minimum value of the cost function (27), which is usually solved by the nonlinear iteration method. The linearized measurement equation and the Gauss–Newton iteration are used to achieve the iterative formula for solving the minimum value of  J x k  as follows:
x ^ k i + 1 = x ^ k + P k H k i T H k i P k H k i T + R k 1 z k h x ^ k i H k i x ^ k x ^ k i
where  H k i = h x k / x k x k = x ^ k i  denotes the Jacobian matrix, and i represents the iteration index.
The analogous representation that  P x z = P k H k T  and  P z = H k P k H k T + R k  is obtained using the statistical linear error approach, which is used instead of the Jacobian matrix to reduce the computational complexity and avoid the linearization error caused by the Taylor first-order expansion. Furthermore, the enhanced Gauss–Newton iteration is displayed as follows.
x ^ k i + 1 = x ^ k + P x z i P z i 1 z k h x ^ k i P x z i T P k 1 x ^ k x ^ k i

4.3. The Iterated Orthogonal Simplex Cubature Kalman Filter

The enhanced Gauss–Newton iteration and the aforementioned altered cubature points are combined to create the IOSCKF method within the nonlinear Kalman filter architecture, and the specific calculation steps are as follows.
Step 1: Filter Initialization
x ^ 0 + = E x 0 P 0 + = E x 0 x ^ 0 + x 0 x ^ 0 + T
Step 2: Time Update
x ^ k 1 +  and  P k 1 +  take the place of  x ¯  and  P x  in Equation (24) to calculate the cubature points  x ^ k 1 i , which are disseminated nonlinearly using  f , and the prior state  x ^ k  and covariance  P k  are calculated using the weights in Equation (23) on the basis of Equations (4) and (5), respectively.
Step 3: Measurement Update
Take the prior state  x ^ k 0 = x ^ k  and covariance matrix  P k 0 = P k  as the initial values of the iterated measurement update, and the estimated state and covariance matrix in the jth iteration express  x ^ k j  and  P k j  ( j = 1 , 2 , , N max ), respectively. The cubature points  x ^ k i , j 1  are calculated using  x ^ k j 1  and  P k j 1  instead of  x ¯  and  P x  in Equation (24), and propagated nonlinearly using  h . The predicted measurement  z ^ k j , the measurement covariance matrix  P z j , and the cross covariance matrix  P x z j  are calculated using the weights in Equation (23) on the basis of the measurement update process of the Gaussian approximation filter, respectively. The Kalman filtering  K k j  and covariance matrix  P k j  are estimated to be equivalent to the Gaussian approximation filter, respectively, while the posterior state  x ^ k j  is calculated using the following iterated form.
x ^ k j = x ^ k + K k j z k h x ^ k j 1 P x z j T P k 1 x ^ k x ^ k j 1
The above iteration format obtains a posterior estimate higher than the last one after each iteration step, and the system model h(·) is known and unchanged, so the whole iterative calculation process is convergent and stable.
The iteration termination condition is set as shown below.
x ^ k j x ^ k j 1 ε or   j = N max
where  ε  and  N max  are the pre-set threshold and maximum iteration times.
The number of iterations in the end is presumed to be N , whereupon the estimated state and covariance matrix at time  k  are  x ^ k + = x ^ k N  and  P k + = P k N , respectively. Note that the calculational burden of the ICKF and CKF is within the same order of magnitude if  N max  is not too large, and in general,  N max  is no more than three.
Considering that all parameters in the algorithm are definite values, the algorithm has low sensitivity and strong robustness to parameters.

5. Simulations and Results

The proposed IOSCKF algorithm is applied in CV target tracking and spacecraft orbit tracking, respectively, in this section.

5.1. CV Target Tracking Simulation

The target is assumed to be in constant velocity (CV) motion; the state equation of the CV model in a two-dimensional case is described as follows:
X k = F C V X k 1 + G w k 1
where  X k = x k x ˙ k y k y ˙ k T  represents the position and velocity of the target at time index  k w k  denotes the process noise, and  F C V  and  G  denote the state transformation matrix and the noise-driven matrix, which are defined, respectively as follows:
F C V = 1 T 0 0 0 1 0 0 0 0 1 T 0 0 0 1 ,   G = T 2 / 2 0 T 0 0 T 2 / 2 0 T
where T is the sampling interval.
In the target tracking system, the bearings-only measurement equation is written as follows.
Z k = arctan y k y r x k x r + v k
where  Z k  is the radar measurement at time  k x r y r  is the location of the radar, and  v k  is the measurement noise.
The proposed algorithm has the following two degraded forms: the cubature Kalman filter with the improved spherical simplex-radial rule (SCKF), and the SCKF with orthogonal transformation (OSCKF). In order to evaluate the performance of the proposed algorithm, it is compared to the general CKF, SSRCKF, SCKF, and OSCKF. This simulation shows the radar location  x r y r = 200   m 300   m , the simulation time is 40 s, T = 1, the initial position of the target is  x 0 y 0 = 100   m 200   m  and the initial velocity is  x ˙ 0 y ˙ 0 = 2   m / s 20   m / s . The initial state and covariance matrix in the simulation are set to  x ^ 0 + = 100   m 2   m / s 200   m 20   m / s T  and  P 0 + = d i a g 0.01 0.01 0.01 0.01 , respectively. The standard deviation of the measurement noise is 0.1 deg. The maximum iteration time  N max  is set to 3. The Monte Carlo simulation is implemented 500 times, and the position accuracy is evaluated using the root mean square error (RMSE), defined as follows:
R M S E p o s = 1 N i = 1 N x k x ^ i , k + 2 + y k y ^ i , k + 2
where N is the execution number of Monte Carlo simulations,  x k y k  is the real position of the target, and  x ^ i , k + y ^ i , k +  is the estimated position at time k in the ith Monte Carlo simulation. The velocity RMSE is defined similarly.
The simulation results are shown in Figure 1. The position RMSE and velocity RMSE of various filters are shown in Figure 1a,b, respectively, and the partial magnification is given to show the details of the curve. It can be seen from the figure that all the four filters have the ability to track the target trajectory, the tracking accuracy of the conventional CKF is significantly lower than that of the other filters, and the proposed IOSCKF achieves higher target tracking accuracy than the other filters. The average RMSE of the various filters is shown in Table 1; compared with the CKF, the estimated position accuracy of the SCKF is improved by 20.77%, indicating that the second-order Gauss–Laguerre rule used in the proposed modified spherical simplex-radial cubature rule has higher accuracy than the spherical-radial cubature rule. The accuracy of the OSCKF is 1% greater than that of the SCKF due to the orthogonal transformation of the cubature points reducing the high-order error components. While the filtering accuracy is further enhanced by the repetitive calculation of the measurement update, the IOSCKF’s accuracy is 0.5% higher than that of the OSCKF. Through the analysis of the simulation results, the validity of the proposed method is verified.

5.2. Spacecraft Target Tracking Simulation

5.2.1. Spacecraft Orbit Dynamics Model

For spacecrafts operating in low Earth orbit, Earth’s J2 term non-spherical perturbation and atmospheric drag perturbation are the most important perturbations. In the J2000 geocentric inertial coordinate system (O-XYZ), considering the above two perturbations, the orbital dynamics model of the spacecraft is as follows:
x ˙ = v x y ˙ = v y z ˙ = v z v ˙ x = μ x r 3 J 2 R e r 7.5 z 2 r 2 1.5 1 + f a x + f x v ˙ y = μ y r 3 J 2 R e r 7.5 z 2 r 2 1.5 1 + f a y + f y v ˙ z = μ z r 3 J 2 R e r 7.5 z 2 r 2 4.5 1 + f a z + f z r = x 2 + y 2 + z 2
where  X = x y z T  is the position vector of the spacecraft in O-XYZ;  X ˙ = v x v y v z T  is the velocity vector of the spacecraft in O-XYZ; J2 is the harmonic term coefficient;  μ  is the gravitational constant of the Earth;  R e  is the radius of the Earth; and  f x f y f z T  is the sum of the high-order non-spherical perturbation of the Earth, three-body gravitational perturbation, and solar luminous pressure perturbation on the three coordinate axes, which can be equivalent to Gaussian white noise with a zero mean value in this study.  f a = f a x f a y f a z T  is the atmospheric drag perturbation, and its specific expression is as follows:
f a = 1 2 c d A m ρ d v r e l v r e l
where  c d  is the atmospheric drag coefficient,  A / m  is the surface mass ratio of the spacecraft,  ρ d  is the atmospheric density, and  v r e l  is the relative speed between the spacecraft and the atmosphere, assuming that the atmosphere rotates with the Earth, as follows:
v r e l = X ˙ ω × X
where  ω = ω e × 0 0 1 T  is the angular velocity vector of the Earth rotation and  ω e  is the angular velocity of the Earth rotation. Using the fourth-order Runge–Kutta method, Equation (32) is written in the following discrete form:
X k = f X k 1 + w k 1
where  X k = x k y k z k v x , k v y , k v z , k T  is the track state at time k, and  w k  is the system noise.

5.2.2. Ground Radar Measurement Model

The radar measurement model is built in the radar horizon coordinate system ( O X h Y h Z h ), and the orbit model is built in the radar horizon coordinate system ( O X Y Z ), so it is necessary to use the WGS84 Earth fixed coordinate system ( O X e Y e Z e ) to achieve the conversion from  O X Y Z  to  O X h Y h Z h . Assuming that the orbital state of the spacecraft in  O X Y Z  is  X X ˙ T , the orbital state in  O X e Y e Z e  is  X f X ˙ f T , and the orbital state in  O X h Y h Z h  is  X h X ˙ h T , the conversion of the orbital state is completed in two steps.
Step 1: Transition from  O X Y Z  to  O X e Y e Z e
The conversion matrix is  M J W = M P w M R o M N u M P r , where  M P r  is the precession matrix,  M N u  is the nutation matrix,  M R o  is the Earth rotation matrix, and  M P w  is the pole shift matrix; then, we can obtain:
X f = M J W X
X ˙ f = M J W X ˙ + M ˙ J W X
In these formulas, the derivative of the transformation matrix  M ˙ J W  is shown as follows:
M ˙ J W = d M P w M R o M N u M P r d t M P w d M R o d t M N u M P r
where  d M R o / d t  is the derivative matrix of the Earth rotation matrix as follows:
d M R o d t = 0 ω e 0 ω e 0 0 0 0 0 M R o
Equations (36) and (37) are written in matrix form as:
X f X ˙ f = M J W 0 M ˙ J W M J W X X ˙
Step 2: Transition from  O X e Y e Z e  to  O X h Y h Z h
The conversion matrix  M W H  is:
M W H = sin λ cos λ 0 sin φ cos λ sin φ sin λ cos φ cos φ cos λ cos φ sin λ sin φ
where  λ  is the radar geocentric longitude and  φ  is the radar geocentric latitude, which can be converted to the geocentric coordinates of the radar  X c = x c y c z c T  in  O X e Y e Z e ; so, we obtain:
X h = M W H X f X c
X ˙ h = M W H X ˙ f
This is written in matrix form as below:
X h X ˙ h = M W H 0 0 M W H X f X ˙ f M W H X c 0
Substituting Equation (40) into Equation (44) yields:
X h X ˙ h = M W H 0 0 M W H M J W 0 M ˙ J W M J W X X ˙ M W H X c 0
Thus, the transformation of the orbital state from  O X Y Z  to  O X h Y h Z h  is achieved, as shown in Equation (45). In order to obtain the relationship between the measured value and the orbit state, we write  X h  and  X ˙ h  in the form of a specific vector as X h = x h y h z h T  and  X ˙ h = x ˙ h y ˙ h z ˙ h T ; then, the geometric relationship between the radar ranging value R, velocity measurement value  R ˙ , azimuth angle A, and pitch angle E and the orbit state is as follows:
R = x h 2 + y h 2 + z h 2 A = arctan y h x h E = arctan z h x h 2 + y h 2
Equation (46) can be written as the following discrete measurement equation:
Z k = h X k + v k
where  Z k = R k A k E k T  is the measured value at time k, and  v k  is the measured noise.
Thus, Equations (35) and (47) constitute the state equation and measurement equation corresponding to Equation (11) for real-time filtering orbit determination.

5.2.3. Simulation Result and Analysis

Figure 2a shows the spacecraft ground simulator on which the algorithm of High-Precision Orbit Propagation (HPOP) is run, which takes into account the perturbation effects of the Earth’s high-order non-spherical gravity, atmospheric drag, solar light pressure, three-body gravity, and tides. Among them, the 21-layer Earth model, atmospheric drag coefficient  c d = 2.2 , and spacecraft surface mass ratio  0.02   m 2 / kg  are considered; the atmospheric density model adopts the Jacchia–Roberts model; and the solar luminous pressure coefficient  c r = 1 . The orbital epoch of the spacecraft considered in the experiment is 16:00:00 (UTC) on 1 September 2016, and the six orbital roots are:
a = 7028.14   km ,   e   =   0 ,   i   =   80.7 ,   Ω   =   63.8 ,   ω   =   0 ,   f   =   107.1
As shown in Figure 2b, the lowest measured elevation angle of the radar is 10°. It is assumed that the radar ranging error is 60 m and the angle error is 0.02°.
The algorithm proposed in this paper is compared with the SCKF and OSCKF, and the root mean square error (RMSE) is used to evaluate the spacecraft tracking accuracy. The Monte-Carlo simulation is run 200 times, and the results are shown in Figure 3. The mean positioning RMSE and velocity RMSE are listed in Table 2.
It can be seen that the spacecraft tracking accuracy of the OSCKF is higher than that of the SCKF; specifically, the position tracking accuracy is increased by 2.21%, and the velocity tracking accuracy is increased by 1.94%, which indicates that the orthogonal transformation can improve the tracking accuracy. The spacecraft tracking accuracy of the IOSCKF is higher than that of the OSCKF; specifically, the position tracking accuracy is increased by 2.71%, and the velocity tracking accuracy is increased by 2.97%, which indicates that the tracking accuracy can be effectively improved by introducing iterative calculation into the measurement equation, thus verifying the effectiveness of the method presented in this paper.

6. Conclusions

A new iterated orthogonal simplex cubature Kalman filter (IOSCKF) is proposed to improve the nonlinear system filtering accuracy in target tracking. The new cubature points are determined using the orthogonal matrix after being modified using the spherical simplex-radial rule and the second-order Gauss–Laguerre quadrature rule. This lowers the high-order error terms. The Gauss–Newton iteration is used to maximize the posteriori state and covariance in the measurement update process until the iteration is finished. Target tracking simulations, including CV target tracking and spacecraft orbit tracking, are used to examine the performance of the proposed method, and the results reveal that the IOSCKF can achieve greater accuracy compared to the CKF, SCKF, and OSCKF. In spacecraft orbit tracking simulation, compared with the SCKF, the position tracking accuracy and velocity tracking accuracy of the OSCKF are increased by 2.21% and 1.94%, respectively, which indicate that the orthogonal transformation can improve the tracking accuracy. Furthermore, compared with the OSCKF, the position tracking accuracy and velocity tracking accuracy of the IOSCKF are increased by 2.71% and 2.97%, respectively, which indicate that the tracking accuracy can be effectively improved by introducing iterative calculation into the measurement equation, thus verifying the effectiveness of the method presented in this paper.

Author Contributions

Methodology, Z.L. and H.C.; Software, Z.L. and X.Y.; Investigation, X.Y. and L.L.; Resources, L.L.; Data curation, L.L.; Writing—original draft, Z.L. and H.C.; Funding acquisition, H.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China (Grant 62005320).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant 62005320).

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Wu, H.; Chen, S.; Yang, B.; Chen, K. Feedback robust cubature Kalman filter for target tracking using an angle sensor. Sensors 2016, 16, 629. [Google Scholar] [CrossRef] [PubMed]
  2. Alomair, R.A.; Hassan, S.Z.; Abdelrahman, M.A.E.; Amin, A.H.; El-Shewy, E.K. New solitary optical solutions for the NLSE with δ-potential through Brownian process. Results Phys. 2022, 40, 105814. [Google Scholar] [CrossRef]
  3. Abdelwahed, H.G.; Alsarhana, A.F.; El-Shewy, E.K.; Abdelrahman, M.A.E. Higher-Order Dispersive and Nonlinearity Modulations on the Propagating Optical Solitary Breather and Super Huge Waves. Fractal Fract. 2023, 7, 127. [Google Scholar] [CrossRef]
  4. El-Shewy, E.K.; Alharbi, Y.F.; Abdelrahman, M.A.E. On the dynamical stochastic electrostatic noise fluctuations in Zakharov model. Chaos Solitons Fractals 2023, 170, 113324. [Google Scholar] [CrossRef]
  5. Li, J.-Q.; Zhao, R.-H.; Chen, J.-L.; Zhao, C.-Y.; Zhu, Y.-P. Target tracking algorithm based on adaptive strong tracking particle filter. IET Sci. Meas. Technol. 2016, 10, 704–710. [Google Scholar]
  6. Wu, H.; Chen, S.; Yang, B.; Luo, X. Range-parameterised orthogonal simplex cubature Kalman filter for bearings-only measurements. IET Sci. Meas. Technol. 2016, 10, 370–374. [Google Scholar] [CrossRef]
  7. 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–1176. [Google Scholar] [CrossRef]
  8. Swati, S.B. Cubature quadrature Kalman filter. IET Signal Process. 2013, 7, 533–541. [Google Scholar]
  9. Ahmed, N.; Rutten, M.; Bessell, T.; Kanhere, S.S.; Gordon, N.; Jha, S. Detection and tracking using particle-filter-based wireless sensor networks. IEEE Trans. Mob. Comput. 2010, 9, 1332–1344. [Google Scholar] [CrossRef]
  10. Leong, P.H.; Arulampalam, S.; Lamahewa, T.A.; Abhayapala, T.D. Gaussian-sum cubature Kalman filter with improved robustness for bearings-only tracking. IEEE Signal Process. Lett. 2014, 21, 513–517. [Google Scholar] [CrossRef]
  11. Reif, K.; Gunther, S.; Yaz, E.; Unbehauen, R. Stochastic stability of the discrete-time extended Kalman filter. IEEE Trans. Autom. Control 1999, 44, 714–728. [Google Scholar] [CrossRef]
  12. Liu, G.; Tian, G. Square-root sigma-point information consensus filters for distributed nonlinear estimation. Sensors 2017, 17, 800. [Google Scholar] [CrossRef] [PubMed]
  13. Bell, B.M.; Cathey, F.W. The iterated Kalman filter update as a Gauss-Newton method. IEEE Trans. Autom. Control 1993, 38, 294–297. [Google Scholar] [CrossRef]
  14. Julier, S.J.; Uhlmann, J.K.; Durrant-Whyte, H.F. A new method for nonlinear transformation of means and covariances in filters and estimators. IEEE Trans. Autom. Control 2000, 45, 477–482. [Google Scholar] [CrossRef]
  15. Arasaratnam, I.; Haykin, S.; Elliott, R.J. Discrete-time nonlinear filtering algorithms using Gauss-Hermite quadrature. Proc. IEEE 2007, 95, 953–975. [Google Scholar] [CrossRef]
  16. Jia, B.; Xin, M.; Cheng, Y. Sparse-grid quadrature nonlinear filtering. Automatica 2012, 48, 327–341. [Google Scholar] [CrossRef]
  17. Arasaratnam, I.; Haykin, S. Cubature Kalman filters. IEEE Trans. Autom. Control 2009, 54, 1254–1269. [Google Scholar] [CrossRef]
  18. Arasaratnam, I.; Haykin, S.; Hurd, T.R. Cubature Kalman filtering for continuous-discrete systems: Theory and simulations. IEEE Trans. Signal Process. 2010, 58, 4977–4992. [Google Scholar] [CrossRef]
  19. Potnuru, D.; Chandra, K.P.B.; Arasaratnam, I.; Gu, D.-W.; Mary, K.A.; Ch, S.B. Derivative-free square-root cubature Kalman filter for non-linear brushless DC motors. IET Electr. Power Appl. 2016, 10, 419–429. [Google Scholar] [CrossRef]
  20. Chang, L.; Hu, B.; Li, A.; Qin, F. Transformed unscented Kalman filter. IEEE Trans. Autom. Control 2013, 58, 252–257. [Google Scholar] [CrossRef]
  21. Zhang, L.; Yang, H.; Lu, H.; Zhang, S.; Cai, H.; Qian, S. Cubature Kalman filtering for relative spacecraft attitude and position estimation. Acta Astronaut. 2014, 105, 254–264. [Google Scholar] [CrossRef]
  22. Li, W.; Jia, Y. Location of mobile station with maneuvers using an IMM-based cubature Kalman filter. IEEE Trans. Ind. Electron. 2012, 59, 4338–4348. [Google Scholar] [CrossRef]
  23. Zhang, Y.; Huang, Y.; Li, N.; Zhao, L. Embedded cubature Kalman filter with adaptive setting of free parameter. Signal Process. 2015, 114, 112–116. [Google Scholar] [CrossRef]
  24. Zhang, K.; Shan, G. Model-switched Gaussian sum cubature Kalman filter for attitude angle-aided three-dimensional target tracking. IET Radar Sonar Navig. 2015, 9, 531–539. [Google Scholar] [CrossRef]
  25. Zarei, J.; Shokri, E. Convergence analysis of non-linear filtering based on cubature Kalman filter. IET Sci. Meas. Technol. 2015, 9, 294–305. [Google Scholar] [CrossRef]
  26. Wu, P.; Li, X.; Bo, Y. Iterated square root unscented Kalman filter for maneuvering target tracking using TDOA measurements. Int. J. Control Autom. Syst. 2013, 11, 761–767. [Google Scholar] [CrossRef]
  27. Zhan, R.; Wan, J. Iterated unscented Kalman filter for passive target tracking. IEEE Trans. Aerosp. Electron. Syst. 2007, 43, 1155–1163. [Google Scholar] [CrossRef]
  28. Chang, G.; Xu, T.; Wang, Q. Alternative framework for the iterated unscented Kalman filter. IET Signal Process. 2017, 11, 258–264. [Google Scholar] [CrossRef]
  29. Cui, B.; Chen, X.; Xu, Y.; Huang, H.; Liu, X. Performance analysis of improved iterated cubature Kalman filter and its application to GNSS/INS. ISA Trans. 2017, 66, 460–468. [Google Scholar] [CrossRef]
  30. Jia, B.; Xin, M.; Cheng, Y. High-degree cubature Kalman filter. Automatica 2013, 49, 510–518. [Google Scholar] [CrossRef]
  31. Wang, S.; Feng, J.; Tse, C.K. Spherical simplex-radial cubature Kalman filter. IEEE Signal Process. Lett. 2014, 21, 43–46. [Google Scholar] [CrossRef]
  32. Singh, A.K.; Bhaumik, S. Higher degree cubature quadrature Kalman filter. Int. J. Control Autom. Syst. 2015, 13, 1097–1105. [Google Scholar] [CrossRef]
  33. Wang, C.; Zhang, J.; Mu, J. Maximum likelihood-based iterated divided difference filter for nonlinear systems from discrete noisy measurements. Sensors 2012, 12, 8912–8929. [Google Scholar] [CrossRef] [PubMed]
Figure 1. The RMSE of various filters. (a) The position RMSE. (b) The velocity RMSE.
Figure 1. The RMSE of various filters. (a) The position RMSE. (b) The velocity RMSE.
Applsci 14 00392 g001aApplsci 14 00392 g001b
Figure 2. Schematic diagram of spacecraft tracking. (a) Spacecraft simulator. (b) Radar spacecraft tracking.
Figure 2. Schematic diagram of spacecraft tracking. (a) Spacecraft simulator. (b) Radar spacecraft tracking.
Applsci 14 00392 g002
Figure 3. The spacecraft tracking RMSE of the various filters. (a) Position RMSE. (b) Velocity RMSE.
Figure 3. The spacecraft tracking RMSE of the various filters. (a) Position RMSE. (b) Velocity RMSE.
Applsci 14 00392 g003aApplsci 14 00392 g003b
Table 1. Average RMSE of various filters.
Table 1. Average RMSE of various filters.
FiltersPosition RMSE/mVelocity RMSE/(m/s)
CKF2.93790.1852
SCKF2.32770.1643
OSCKF2.30450.1622
IOSCKF2.29290.1619
Table 2. Spacecraft tracking RMSE of the various filters.
Table 2. Spacecraft tracking RMSE of the various filters.
FiltersPosition RMSE/mVelocity RMSE/(m/s)
SCKF24.860.103
OSCKF24.310.101
IOSCKF23.650.098
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Li, Z.; Yang, X.; Li, L.; Chen, H. Iterated Orthogonal Simplex Cubature Kalman Filter and Its Applications in Target Tracking. Appl. Sci. 2024, 14, 392. https://doi.org/10.3390/app14010392

AMA Style

Li Z, Yang X, Li L, Chen H. Iterated Orthogonal Simplex Cubature Kalman Filter and Its Applications in Target Tracking. Applied Sciences. 2024; 14(1):392. https://doi.org/10.3390/app14010392

Chicago/Turabian Style

Li, Zhaoming, Xinyan Yang, Lei Li, and Hang Chen. 2024. "Iterated Orthogonal Simplex Cubature Kalman Filter and Its Applications in Target Tracking" Applied Sciences 14, no. 1: 392. https://doi.org/10.3390/app14010392

APA Style

Li, Z., Yang, X., Li, L., & Chen, H. (2024). Iterated Orthogonal Simplex Cubature Kalman Filter and Its Applications in Target Tracking. Applied Sciences, 14(1), 392. https://doi.org/10.3390/app14010392

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