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:
where
k is the time step,
denotes the state vector at time step
k,
represents the measurement vector at time step
k, and
and
are the known nonlinear state function and measurement function, respectively. The process noise
and measurement noise
are uncorrelated zero-mean Gaussian white noise, with the covariance matrices being
and
, 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
in the Cartesian coordinate system, where
denotes the arbitrary nonlinear function, and
represents the Gaussian pdf, with the mean and covariance being
and
, respectively.
This Gaussian-weighted integral
can be approximated using various quadrature rules, and the integral with respect to
can be approximated by the following quadrature rule:
where
I is the identity matrix,
denotes the total number of points, and
and
represent the quadrature points and weights, respectively.
The integral with respect to the more general Gaussian distribution
can be obtained using the equation below:
where
denotes the square root of
, which can be solved by the Cholesky decomposition, and
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
and error covariance
are calculated, respectively, as follows:
where
stands for the integral points generated using the numerical rule.
Measurement Update:
The updated state
and corresponding error covariance
are estimated, respectively, as follows:
where
denotes the Kalman filtering gain,
represents the predicted measurement,
denotes the integral points generated by the numerical rule,
denotes the innovation covariance matrix, and
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])
.where and is a given weighing function. Equation (8) is a dth-degree rule if it is accurate for , whose components are linear combinations of monomials
with a total degree up to d.
,
,
, are nonnegative integers and , and there is at least one monomial of degree
for which Equation (8) is not exact. Specifically, an intractable integral of the form
is considered first in the cubature rule. Let
with
, where
means the direction vector and the radius
. As a consequence,
; then, the integral
can be written anew as in a spherical-radial coordinate system as follows [
12]:
where
is the sphere surface defined by
, and
is the area element on
. The integral
is decomposed into the spherical integral
and radial integral
, respectively.
Because of the difficulty in achieving the analytical solutions, the approximate solution method of
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
and
is approximated using the following formula:
where
and
denote the points and weights to calculate the spherical integral, and the number of points is
.
and
represent the points and weights to calculate the radial integral, and the number of points is
. The total number of points is
if
, or
if one of
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
where
denotes the surface area of the unit sphere, and
is the Gamma function. The vector
represents the unit column vector with the
ith element being 1. The number of cubature points needed is
.
The third-degree spherical simplex rule [
31] consisting of
points is
where
denote the vertices of the
n-dimensional simplex, and the elements are defined as follows.
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
in this paper. For more details, please refer to the related references.
3.2. Radial Rule
Then, the calculation of the radial integral
is considered. Two methods, including the moment matching method and Gauss–Laguerre quadrature rule, may be used. To construct the (2
m + 1)th-degree radial rule, the two methods lead to identical points and weights when
m is odd and
points are used. However, they may be different when
m is even and
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,
).
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
, let
to obtain
, which takes the place of the radial integral, and for
, let
, and
. Then,
is transformed into
, and this form can be handily approximated using the Gauss–Laguerre quadrature [
32] rule as follows:
where
denotes the quadrature points, with
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:
and the weights are calculated using the following equation.
On account of
m being set at 2, the solutions are obtained are
,
,
, and
, which are substituted into Equation (17) to deduce the following radial rule.
From Equation (20) we see that the number of points is
, 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
points.
If we define matrix
, and the calculation methods of the cubature points and weights are extracted from cubature rule (21), we obtain the following.
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
, with weights shown in Equation (23), remain to constitute a numerical cubature rule if
is an
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:
where
, and the orthogonal matrix
can be designed using different methods. As an effective design method, the
matrix
is constructed in this paper as
, and
, whose elements are defined as follows:
where
if
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
and covariance matrix
are used in the measurement update process to calculate the posterior state
and covariance matrix
, which are revised using the latest measurement
via Kalman gain
. This means an approximation accuracy of
to the real state is bound to be better than
, and the filtering accuracy will be improved if we can use
and
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
and
, and the likelihood function of
and
is defined as follows [
33].
where
is a constant, and the cost function is defined as follows.
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
as follows:
where
denotes the Jacobian matrix, and
i represents the iteration index.
The analogous representation that
and
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.
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
Step 2: Time Update
and take the place of and in Equation (24) to calculate the cubature points , which are disseminated nonlinearly using , and the prior state and covariance 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
and covariance matrix
as the initial values of the iterated measurement update, and the estimated state and covariance matrix in the
jth iteration express
and
(
), respectively. The cubature points
are calculated using
and
instead of
and
in Equation (24), and propagated nonlinearly using
. The predicted measurement
, the measurement covariance matrix
, and the cross covariance matrix
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
and covariance matrix
are estimated to be equivalent to the Gaussian approximation filter, respectively, while the posterior state
is calculated using the following iterated form.
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.
where
and
are the pre-set threshold and maximum iteration times.
The number of iterations in the end is presumed to be , whereupon the estimated state and covariance matrix at time are and , respectively. Note that the calculational burden of the ICKF and CKF is within the same order of magnitude if is not too large, and in general, 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.