Asian Journal of Control, Vol. 17, No. 3, pp.

1–10, May 2015

Published online in Wiley Online Library (wileyonlinelibrary.com) DOI: 10.1002/asjc.926

A NEW DERIVATION OF THE CUBATURE KALMAN FILTERS

Xin-Chun Zhang and Yun-Long Teng

ABSTRACT
This paper investigates the cubature Kalman iltering (CKF) for nonlinear dynamic systems. This third-degree
rule based ilter employs a spherical-radial cubature rule to numerically compute the integrals encountered in nonlinear
iltering problems, thereby removing the requirements of explicitly computing the Jacobians. The cubature rule, however,
requires computing the intractable integrals over a high-dimensional spherical region for multidimensional applications.
Moreover, the cubature formula that has been used to construct the spherical cubature formula has some demerits,
most notably its inconvenient properties in computation and low estimation accuracy. Aimed at these issues, a general
class of CKFs that uses only cubature rules is derived in this paper. It can be shown that the conventional CKF is a
special case of the proposed algorithm. The paper also includes higher-degree CKFs, especially two representative types
of the ifth-degree CKFs. Performance of the proposed algorithms is demonstrated via two target tracking problems.
The experimental results, presented herein, illustrate the superior performance of higher-degree CKFs to conventional
nonlinear ilters.
Key Words: Nonlinear iltering, cubature Kalman ilter, high-degree cubature rule, state estimation, orbit.

I. INTRODUCTION the extended Kalman ilter (EKF) [7,8], which has been
used as the state-of-the-art technology in many engineer-
The Kalman ilter [1] is, currently, one of the most ing ields, generally linearizes all linear models via the
eficient recursive ilters for state estimation of linear Taylor series expansion. This linearization may introduce
dynamic systems. It is a widespread concept in the ields large errors if either the process or measurement model
of engineering and econometric applications from target is highly nonlinear. In addition, the computation of the
tracking and computer vision to estimation of structural Jacobian matrices is also a trivial mathematical task.
macroeconomic models [2,3]. It is also commonly used in As better alternatives to the EKF, the unscented
navigation and control of space vehicles, however, more Kalman ilter (UKF) [9,10], Gauss-Hermite ilter (GHF)
complex systems can be nonlinear. The nonlinearity can and the central difference ilter (CDF) [4], among oth-
be involved either in the process equation or the mea- ers, have been developed to cope with highly nonlinear
surement equation, which can severely degrade the ilter systems. Of all these ilters, the UKF preferred for its
performance. high accuracy and easy implementation, however, owing
Over the past several decades, considerable effort to the possible negative weight on the center point, the
has been made to develop eficient iltering algorithms for unscented transform (UT) or the UKF is believed to be
nonlinear/non-Gaussian systems. In theory, an optimal unstable [11]. As with the UKF, all the aforementioned
solution can be obtained by Bayesian probabilistic infer- ilters select a deterministic number of sampling points
ence, but the associated solution requires propagating to capture the irst two-order moments, thereby removing
the full density function, which is analytically intractable the requirements of explicitly computing the Jacobians.
[4–6]. Consequently, some assumptions and approxima- In general, they also yield higher level state estimation
tions, in Bayesian framework, have been made to design than the EKF.
practical suboptimal nonlinear ilters.
Another approach is the celebrated particle ilter
When nonlinear iltering problems are encountered,
(PF) [12] which picks a random number of parti-
the irst solution usually is to linearize. In one approach,
cles to approximate the posterior density. Similar il-
ters include Gaussian particle ilter (GPF) [13] and the
Manuscript received March 13, 2013; revised August 23, 2013; accepted Quasi-Gaussian particle ilter (QGPF) [6]. These ilters
December 21, 2013.
Xin-Chun Zhang (corresponding author, e-mail: [email protected]) is have been found to suffer from the “curse of dimen-
with the University of Electronic Science and Technology of China, Chengdu sionality” when applied to high-dimensional systems. In
611731, China.
Yun-Long Teng is with School of Energy Science and Engineering, University
addtion, some iltering algorithms, such as the H2 /H∞
of Electronic Science and Technology of China, Chengdu 611731, China. iltering [14] and the robust iltering techniques [15] are

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 Asian Journal of Control, Vol. 17, No. 3, pp. 1–10, May 2015

also often used in modern control theory and control pair of difference equations in discrete-time [11]:
systems engineering. ( )
xk = f xk−1 , uk−1 + vk−1 (1)
In 2009, Arasaratnam, Haykin et al. [11,16] pro-
posed the cubature Kalman ilter (CKF) which could ( )
solve high-dimensional nonlinear state estimation prob- zk = h xk , uk + wk (2)
lems. The third-degree cubature rule-based ilter employs
deterministic sampling strategy to capture the mean and where (1) and (2) are, respectively, the process equation
covariance. The CKF has achieved great success since its and the measurement equation, xk ∈ Rnx is the state at
inception, and has been adopted to solve many real-life time k, uk−1 ∈ Rnu is the control input, f and h are some
problems [17–19] for its more stable performance than known nonlinear functions; zk ∈ Rnz is the measure-
the conventional nonlinear ilters, such as the EKF and ment; vk−1 and wk are white noise with zero mean and
UKF [20]. covariance Qk−1 and Rk , respectively. The major task of
The CKF has been found to be not suficiently nonlinear iltering is continuous to yield the estimation
accurate in many nonlinear iltering problems [20]. To of the discrete system state at time k.
improve the CKF performance, Zhang [21] formulated a The CKF assumes that the predictive density of
square-root embedded cubature Kalman iler (SECKF), the state xk and the ilter likelihood density of the mea-
and Bhaumik and Swati [22] developed the square-root surement zk are both Gaussian, and then considers the
cubature-quadrature Kalman ilter (SR-CQKF), both of integral of the form
their algorithms deliver higher level state estimation than
∫Rn
the conventional CKF. I( f ) = f (x) exp (−xT x) dx
On the other hand, the developments of computa- ∞ (3)
∫0 ∫ Un
tional mathematics [23–27] reveal that higher estimation = f (ry) rn−1 exp(−r2 ) d� (y) dr
accuracy can be obtained using higher-degree (degree
more than three) cubature rules. Very recently, Jia and {
|
where y = ( y1 , y2 , · · · , yn )T and Un = y ∈ Rn |y21 + · · · +
Xin, et al. [20] developed a class of higher-degree CKFs } |
with more accurate state estimates. Yet their methods are y2n = 1 , �(⋅) is the spherical surface measure or the area
considered to be capable of achieving arbitrary-degree
element on Un .
accuracy in theory, only the rules of a maximum of ive
The primary contribution of [11] is to some extent
degree are taken into account. It may be noted that
the third-degree spherical-radial cubature rule, which
their ifth-degree spherical-radial cubature rule has a
makes it possible to numerically compute the integrals
complex form and is not likely to derive higher-degree
encountered in iltering problems. The rule can be sum-
CKFs in compact forms. In this paper, we present a new
marized as follows:
derivation of the CKFs, which can be easily extended to
( )
∫Rn
higher-degree CKFs. I( f ) = f (x) exp −xT x dx
An outline of the paper is as follows. The derivation
of a general class of CKFs, which use only cubature rules, √ (√ ) (4)
�n ∑
2n
n
is given in Section II. Higher-degree CKFs, especially two ≈ f [1]i
2n i=1 2
representative kinds of the ifth-degree CKFs are derived
in Section III. Section IV yields the performance com- Equation (4) consists of the spherical-radial rule
parison of the proposed algorithms with conventional
∞ ( )
∫0
nonlinear ilters. Conclusions of this paper are drawn
f (x)xn−1 exp −x2 dx ≈ w1 f (x1 ) (5)
in Section V.

and the spherical cubature formula with 2n points

∑2n

∫Un
f ( y) d� (y) ≈ w f ([u]i ) (6)
II. DERIVATION OF THE GENERAL CLASS i=1
OF CKFS where w > 0 is the weight corresponding to the genera-
tor [u].
2.1 Third-degree spherical-radial cubature rule
Since the spherical-radial cubature rule contains
Consider the iltering problem of nonlinear dyna- two basic rules (5) and (6), one should improve the
mic systems with additive process and measurement degrees of both rules simultaneously while developing
noise, whose state-space model can be expressed by the higher-degree CKFs. Moreover, the formula (6) requires

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 X.-C. Zhang and Y.-L. Teng: A New Derivation of the Cubature Kalman Filters

evaluating the integrals over an n-dimensional spherical integrals [28]

region Un , which is not always easy for high-dimensional
integrals (higher-degree spherical-radial cubature rules ̂ 0 f ([0]) + W
̂1 ∑
2n
I( f ) ≈ W f ([u1 ]i ) (8)
require higher-dimensional integrations). The third- i=1
degree rule is also found not to be accurate enough in
many real-life applications [20]. In addition, the evalua- ̂ 0 and W
where I( f ) is given by (3), W ̂ 1 are the weights
tion points of the cubature formula (6) may eventually be corresponding to the generator [0] and [u1 ]. Equation (8)
outside Rn for n suficiently large [28]. is recognized as a fully symmetric formula and it follows
the two properties [21]:

2.2 Third-degree CKFs Property 1. If the integrand f is a monomial containing

In this section, a general class of CKFs will be any odd power of coordinate variables then I(f ) = 0.
introduced. However, before proceeding further, we irst
give the deinition of an orbit, which is important in Property 2. If the integrand f is a monomial without odd
the following derivation of the CKFs that employ fully powers of coordinate variables, then I(f ) depend only on
symmetric cubature rules. the exponents.

Deinition 1. Let G a group of orthogonal transforma- By the properties above, it sufices to consider the
tions that have a ixed point at the origin, the G-orbit of even order integrals, namely, f = 1 and f = x2i for the
a point x ∈ Rn is the set {g(x)|g ∈ G}. third-degree cubature formula; the equations to be satis-
ied are
Obviously, a G-orbit of a given point is an invariant {
Ŵ 0 + 2nW ̂ 1 = I0
set with respect to G. Based on the invariant theory, we (9)
can construct fully symmetric cubature formulae that are ̂ 1 = I2
2u21 W
invariant with respect to
where
{( ) ( )}
G = x1 , · · · , xn → s1 xp1 , · · · , sn xpn (7) ( ) √
∫Rn
I0 = exp −xT x dx = � n (10)

and si = ±1 for those values of i where xpi ≠ 0, i. e.,

where {p1 , p2 , · · · , pn} is a set of the integers 1 , · · · , n,
( ) √ /
∫Rn
the G-orbit of the two-dimensional point (u, 0) is the set I2 = x2i exp −xT x dx = � n 2 (11)
{(u, 0), (0, u), (−u, 0), (0, −u)}, such point is also known as
a generator and denoted by [u] in the CKF. here (10) and (11) make use of the following relationship
Let x ∈ Rn , then there could be a great number of of the gamma function
(i)
(2n − 1)!! √
n(n − 1) · · · (n − i + 1), n(i) = n i!, and ui ≠ 0.
generators, some of them are listed / in Table I, where n =
(i) Γ(1∕2 + n) = � (12)
2n
A typical third-degree cubature formula integrating
over the n-dimensional real number ield Rn in general ̂ 0 and
It can be shown readily from (9) that u1 , W
selects the irst two orbits to evaluate the intractable ̂
W1 are three unknowns, which means the equation is

Table I. Different types of orbits.

number of number of points

generator unknowns in the orbit unknowns
[0] 1 1 weight
[u1 ] 2 2n u1 , weight
[u1 , u1 ] 2 2n(2) u1 , weight
[u1 , u2 ] 3 4n(2) u1 , u2 , weight
[u1 , u1 , u1 ] 2 8n(3) u1 , weight
[u1 , u1 , u2 ] 3 24n(3) u1 , u2 , weight
[u1 , u2 , u3 ] 4 48n(3) u1 , u2 , u3 , weight
⋮ ⋮ ⋮ ⋮
[u1 , u2 , · · · , un ] n+1 2n n! u1 , · · · , un , weight

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 Asian Journal of Control, Vol. 17, No. 3, pp. 1–10, May 2015

under-determined and there √ are ininitely many solu- simultaneously. More about constructing higher-degree
̂0
tions. Speciically, let u1 = nI2 ∕I0 , then the weight W CKFs using higher-degree spherical-radial cubature rule
vanishes, which eventually leads to the outcome W ̂1 = can be seen in [20]. Clearly, such rules are unlikely to give
I0 ∕(2n). Furthermore, we have any compact ilter. For this case, we give the derivation
of the ifth-degree CKF that uses only the ifth-degree

∫Rn
IN ( f ) = f (x)N(x; �, �) dx cubature rule.
In order to construct a fully symmetric formula with
(√ ) ( )
∫
1 fewer points, we irst introduce three consistency con-
=√ f 2x exp −xT x dx
(13) ditions that must be satisied for any fully symmetric
�n R n

formula.
≜
∑
m
�i f (� i ) Let Ni be the number of the ith orbit, the irst
i=1 consistency condition is then naturally obtained.

where � and � are the zero and identity matrix respec- Consistency Condition 1. The number of unknowns
tively,[ N(⋅; ⋅, ⋅) stands for Gaussian density. �i = 1∕(2n),
√ ] should exceed or at least no less than the number of
�i = n , i = 1, · · · , 2n. The result is in consistent with equations Numeqs
i

N0 + 2N1 + · · · + (n + 1)N1,···,n ≥ Numeqs

the conventional CKF.
Theoretically, ininitely many CKFs with 2n+1 (14)
points can be constructed by selecting different values
of u1 , but only one picking 2n points (the conventional The above condition ensures that the equations have
CKF) is available. All these demonstrate that the conven- at least one feasible solution. And the second condition
tional CKF is a special case of the proposed algorithm. that must be satisied is [29]:

Remark 1. Note that the integrals of the monomials Consistency Condition 2. For any symmetric cubature

point, namely, N0 ≤ 1.
comprised of odd power of variable(s), such as f = xi , formula, there can be at most one orbit with only one
f = x3i , f = xi x3j and f = x3i xj , are all zero according to
Property 1.
The inal consistency condition gives the lower
bound of the number of cubature points [30].
III. HIGHER-DEGREE CKFS
Consistency Condition 3. For a symmetric cubature
The importance of the third-degree cubature rules formula, the lower bound of the number of cubature
over those of higher-degree was emphasized in [11]. points is
The authors claimed that higher-degree rules “yield no
( )
Nlower ≥ 2 dim �2m+1
improvement or make the performance worse”. How-
n n+m
ever, many records, both in mathematics [23–27] and in (G) =
n
engineering [20], have shown that higher-degree cubature
( ) (15)
rules could deliver higher level estimation accuracy. Con- ∑
n−1
i+m
sider that if a integral I(f ) is exact of polynomial degree + 2i−n , if m is odd
i
i=1
at most 2m+1, where m is a positive integer, it is also exact
of degree 3. Hence, higher-degree rules can theoretically
and
yield more exact outcomes.
To gain the beneits of higher-degree cubature rules, (
)
Nlower ≥ 2 dim �2m+1
we consider to extend the CKFs to higher degrees. n n+m
(G) − 1 =
Obviously, the ifth-degree CKFs are more easier to n
( )
implement, compared with the other higher-degree ones. ∑
n−1
i+m
Hence, we give special concerns to the ifth-dgree CKFs. + (1 − 2i−n ) , if m is even
i
i=1
(16)
3.1 Type I: Fifth-degree CKF with 2n2 + 1 cubature points
n
While turning to the spherical-radial cubature rule, where �2m+1 is the n−dimensional polynomial with degree
n
one should improve the degrees of the spherical cubature no more than 2m + 1, and dim �2m+1 (G) represents the
rule and the Gauss-Laguerre formula to higher degrees dimension of the vector space.

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 X.-C. Zhang and Y.-L. Teng: A New Derivation of the Cubature Kalman Filters

It is evident that (8) satisies all the consistency According to proposition 4.2 in [11], we have the
conditions mentioned above. A third-degree cubature cubature point sets and the associated weights as follows:
formula needs only the irst two orbits. However, the
irst two orbits are not adequate for a ifth-degree ⎧ [0]i , i=1
⎪ [√ ]
cubature formula, because they don’t meet Consistency ⎪
�i = ⎨ 3 , i = 2, · · · , 2n + 1 (23)
Condition 3. i
When constructing a fully symmetric formula be ⎪ [√ √ ]
⎪ 3, 3 , i = 2(n + 1), · · · , 2n2 + 1
exact for any polynomial of degree 5, by Property 1 and ⎩ i
Property 2, it is suficient
{ } to consider all f in the set
1, xi , xi , xi xj |i ≠ j , namely, f = 1, f = x2i , f = x4i and
2 4 2 2 ⎧ (7−n)n
⎪ 1 − 18 , i = 1
f = x2i x2j . The total number of equations to be satisied is ⎪
�i = ⎨ (4 − n)∕18, i = 2, · · · , 2n + 1 (24)
four, namely
⎪
⎪ 1∕36,
N0 + 2N1 + 2N1,1 + · · · + (n + 1)N1,···,n ≥ 4
i = 2(n + 1), · · · , 2n2 + 1
(17) ⎩
By incorporating (18) and (23)-(24) in the CKF
Hence, at least one of the third or fourth orbits is framework, one can construct the ifth-degree CKF. Note
needed. Formally, the ifth or higher orbits are also possi- that the ifth-degree CKF presented herein picks the same
ble candidates, but they are far from satisfactory since the number of evaluation points but has a more concise form
cubature points required are far above the lower bound than that introduced in [20].
given by (16). A ifth-degree cubature formula with 2n2 + The other ifth-degree cubature formulae, such as
1 points that fulills the three consistency conditions is
̂ 1 ∑ f ([u1 ]i )
2n
given by I( f ) ≈ W
i=1
(25)
n≥2
∑
̂ 1 ∑ f ([u1 ]i )
2n 4n(n−1)
̂ 0 f ([0]) + W
I( f ) ≈ W ̂ 1,1
+W f ([u1 , u2 ]i ),
i=1 i=1
(18)
̂ 1,1 ∑ f ([u1 , u1 ]i ), n≥2
2n(n−1)
+W are also eligible candidates and meet all the consis-
i=1 tency conditions required, but they adopt more cubature
points. We do not consider them here.
the equations to be satisied for this case are

̂ 0 + 2nW ̂ 1 + 2n(2) W
̂ 1,1 = I0 3.2 Type II: Fifth-degree CKF with n2 + 3n + 3
⎧W cubature points
⎪ 2 ̂ 2 ̂
⎪ 2u1 W1 + 4 (n − 1) u1 W 1,1 = I2 In the previous subsection, a ifth-degree CKF with
⎨ 4 (19)
⎪ 2u1 Ŵ 1 + 4 (n − 1) u4 W ̂ = I4 2n2 + 1 cubature points is derived. Its required num-
1 1,1
⎪ 4 ber of evaluation points is signiicantly larger than the
⎩ 4u1 Ŵ 1,1 = I2,2
lower bound n2 + n + 1 when n is suficient large. More
cubature points, however, imply a heavier computational
where I0 and I2 are given by (10) and (11), load. In order to improve the computational eficiency,
we briely introduce another symmetric cubature for-
( ) √ / mula with only n2 + 3n + 3 points, which is consid-
∫Rn
I4 = x4i exp −xT x dx = 3 � n 4 (20)
ered to be one of the most eficient cubature formulae

For the case that n ≥ 4, Mysovskikh [31] derives a

presently known.
( ) √ /
∫Rn
I2,2 = x2i x2j exp −xT x dx = � n 4 (21) centrally symmetric cubature formula for the integration
region Un as follows:

The unique solution of (19) is

̂ 0 f ([0]) + W
̂a ∑
2(n+1)
I( f ) ≈ W f ([a]i )
√ √ i=1
u1 = 3∕2, W ̂ 0 = � n [1 − (7 − n) n∕18], (26)
√ √ / ∑
n(n+1)
(22) ̂b
+W f ([b]i )
̂ 1 = � n (4 − n)∕18, W
W ̂ 1,1 = � n 36
i=1

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 Asian Journal of Control, Vol. 17, No. 3, pp. 1–10, May 2015

√ √ √
̂0= 2 �n ̂ n2 (7−n) � n ̂b= 2(n−1)2 �n candidates is as follows
where W ,Wa = and W
n+2 2(n+1) (n+2)2
2
{ (n+1)2 (n+2)2
are the weights, the generator a = ± a(i) |a(i) = ∑
2n
∑
2n
̂ 0 f ([0]) + W
I( f ) ≈ W ̂1 ̂2
f ([u1 ]i ) + W f ([u2 ]i )
( )}n+1
a(r)
1
, a (r)
2
, · · · , a (r)
n , here r = 1, · · · , n + 1, and i=1 i=1
i=1
∑
2n(n−1)
∑
2n(n−1)
̂ 1,1
+W ̂ 2,2
f ([u1 , u1 ]i ) + W f ([u2 , u2 ]i )
√
⎧− (n+1)(n+2)
, i<r
i=1 i=1
⎪
n≥3
2n(n−i+2)(n−i+1) ∑
4n(n−1)(n−2)∕3
⎪√ ̂ 1,1,1
+W f ([u1 , u1 , u1 ]i ),
a(r) = ⎨ (n+1)(n+2)(n−r+1) (27)
i , i=r i=1
⎪ 2(2m+1)(n−r+2)
⎪ (30)
⎩ 0, i>r
Practically, the CKFs with degree more than seven
{ } {√ can be similarly derived by following the procedures and
( j) n(n−1)∕2 n
the generator b = ±b , where b( j) =
j=1 consistency conditions presented in Subsection 3.1. We
} 2(n−1)
can develop these ilters by adding more generators or
(a(k) + a(l) ), k < l, l = 2, · · · , n + 1 .
orbits. However, the high-degree fully symmetric cuba-
By converting the integrals to be the standard nor- ture formulae require a much larger number of sampling
mal distribution, we obtain the cubature points and the points than the lower bound given by Consistency Con-
associated weights for the ifth-degree CKF of type II dition 3, which implies expensive computational cost
for multidimensional applications. Hence, the CKFs of
⎧ [0]i , i = 1 degree more than ive are not attractive, at least for now.
⎪√
� i = ⎨ 2[a]i , i = 2, · · · , 2n + 3 (28)
⎪√ IV. SIMULATIONS
⎩ 2[b]i , i = 2 (n + 2), · · · , n2 + 3n + 3
In this section, the performance comparisons
of the proposed algorithms with conventional ilters
are demonstrated through two target tracking prob-
⎧ 2
, i=1 lems, including the bearing only tracking (BOT) model
⎪ n+2
[5,6,12,32,33], and the space vehicle tracking problem
⎪ n2 (7−n) which has been extensively investigated in [11,20,34].
�i = ⎨ 2(n+1)2 (n+2)2
, i = 2, · · · , 2n + 3
⎪
⎪ 2(n−1)2
, i = 2 (n + 2), · · · , n2 + 3n + 3
⎩ (n+1)2 (n+2)2 4.1 Bearing only tracking
(29)
The BOT problem describes tracking a moving tar-
get within a two-dimensional plane, the standard model
In a similar way to the ifth-degree CKF of type I, of which can be expressed as follows:
we can develop the ifth-degree CKF of type II.
xk = �xk−1 + 𝚪vk−1 , k = 1, · · · , N (31)
Remark 2. The less cubature points, the better compu-
where xk = (x, x, ̇ Tk , vk−1 = (vx , vy )Tk−1 , and
̇ y, y)
tational performance! Owing to the improvement in esti-
mation accuracy and the reduction of cubature points,
this ifth-degree CKF is recognized as a balance between ⎡1 1 0 0⎤ ⎡ 0.5 0 ⎤
estimation accuracy and computational load. ⎢0 1 0 0⎥ ⎢ 1 0 ⎥
�=⎢
1⎥
, and 𝚪 = ⎢
0 0 1 0 0.5 ⎥
⎢ ⎥ ⎢ ⎥
⎣0 0 0 1⎦ ⎣ 0 1 ⎦
3.3 Higher-degree CKFs
here x and y denote Cartesian coordinates of the moving
The higher-degree CKFs discussed in this section target, ẋ and ẏ are respectively the velocities in the x and y
refer to the CKFs with degree of at least seven. Actu- directions. The system noise is Gaussian distributed with
ally, by adding more generators, we can construct the vk−1 ∼ N(0, Q�� ), where �� is the 2×2 identity matrix. An
seventh-degree cubature formulae, one of the possible observer is ixed at the origin of the plane to takes noise

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 X.-C. Zhang and Y.-L. Teng: A New Derivation of the Cubature Kalman Filters

30
5th-degree CKF(Type II)
0.6
CKF
25
5th-degree CKF(Type I)
0.4
EKF
20 UKF
0.2

0 15

ln(RMSE)
y

-0.2 10

-0.4
5
True
-0.6 5th-degree CKF(Type II)
CKF 0
-0.8 5th-degree CKF(Type I)
EKF -5
-1 UKF
-0.09 -0.08 -0.07 -0.06 -0.05 -10
x x x velocity y y velocity
Axis
Fig. 1. Tracking trajectories.
Fig. 2. Averaged RMSE of all four coordinates (logarithmic
in vertical axis).
Table II. Comparison of the number of cubature
points needed for the CKFs.

Evaluations 5th-degree CKF(Type II)

0.2 CKF
RMSEpos

3rd-degree CKF 8 0.15 5th-degree CKF(Type I)

EKF
5th-degree CKF(type I) 33 0.1 UKF
5th-degree CKF(type II) 31 0.05
0
0 5 10 15 20 25 30
(a)
measurements of the target bearing 0.015
5th-degree CKF(Type II)
zk = arctan( yk ∕xk ) + wk (32) CKF
RMSEvel

0.01 5th-degree CKF(Type I)

EKF
where the measurement noise wk ∼ N(0, R). The other 0.005 UKF
reference data are generated as Q = 0.00052, R =
0.0022 and N = 30. The initial state vector is set to 0
0 5 10 15 20 25 30
be x0 = (−0.05, 0.001, 0.7, −0.055)T and the initial esti- Time,k
mate is assumed to x̂ 0|0 = x0 with covariance P0|0 = (b)
diag{0.072, 0.0052 , 0.12 , 0.012 }. For a fair comparison, we
conduct 200 independent Monte Carlo runs. The track- Fig. 3. RMSE in (a) Position; (b) Velocity. (logarithmic in
ing trajectories of all ive algorithms, including the EKF, vertical axes).
UKF, CKF and the ifth-degree CKFs (type I and type II)
are shown in Fig. 1. The evaluation points of each CKF
are listed in Table II. to quantitatively evaluate the iltering performance of
As can be seen from Fig. 1, both the ifth-degree these nonlinear ilters, where x̂ k|k is the estimate of xk at
CKFs are superior to their counterparts, but type II is the nth Monte Carlo run.
more eficiency since it employs less cubature points. And Fig. 2 presents RMSEs of all four coordinates,
it is apparent that their trajectories are much closer to the which clearly indicates that the ifth-degree CKFs have
actual one. Furthermore, we use the root-mean square smaller root-mean square errors in y axis, but slightly
error (RMSE) averaged across all time instances, which larger error in x axis than the conventional CKF. The
can be deined as RMSEs obtained by the UKF closely follow that of the
√ CKF. In contrast, the EKF has terrible iltering perfor-
N ( )2
1 ∑ mance, making it dificult to track any moving target in
RMSE = xk − x̂ k|k (33)
N k=0 a long run.

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 Asian Journal of Control, Vol. 17, No. 3, pp. 1–10, May 2015

√
The RMSE in position can be deined as 10mrad. The initial state x0 = (1000m, 300m∕s, 1000m,
√ 0m∕s, −3◦ ∕s)T , and the initial estimate is x̂ 0|0 = x0 with
N (( ) ( ))
1∑ covariance P0|0 = diag{100m2, 10m2 ∕s2 , 100m2 , 10m2 ∕s2 ,
RMSEpos (k) = xn−̂xn|n 2+ yn−̂yn|n 2 (34)
Nn=1 100mrad 2 ∕s2 }.
where (xn , yn ) and (x̂ n|n , ŷ n|n ) are the actual and estimated
position at the nth Monte Carlo run. The RMSE in Table III. Comparison of the number of
velocity can be similarly obtained. cubature points needed for the CKFs.
The RMSE in position and velocity are respectively Evaluations
shown in Fig. 3, which clearly demonstrate the over-
all higher estimation accuracy of the ifth-degree CKFs. 3rd-degree CKF 10
5th-degree CKF(type I) 51
Also, the UKF has close performance to the third-degree
5th-degree CKF(type II) 43
CKF, whereas the CKF diverges after a short period of
time.

4.2 Space vehicle tracking

24 CKF
The space vehicle tracking scenario is more complex 5th-degree CKF(Type II)
5th-degree CKF(Type I)
than the BOT model beacuse it involves measurements 23
of the range and bearing, its process model is given by
[11,34]
22
RMSEpos(m)

xk = �xk−1 + vk−1 , k = 1, · · · , N (35)

21
where xk = (x, x, ̇ Ω)Tk , vk−1 ∼ N(0, Q) is the white
̇ y, y,
noise with zero mean and a nonsigular covariance Q =
20
diag{q1M, q1 M, q2 T}, where

⎡ 1 sin ΩT 0 − 1−cos ΩT 0⎤ 19
⎢ Ω Ω
⎥
⎢ 0 cos ΩT 0 − sin ΩT 0⎥
⎢ ⎥ 50 60 70 80 90 100
� = ⎢ 0 1−cos ΩT 1 sin ΩT 0⎥, Time(sec)
⎢ Ω Ω ⎥
⎢ 0 sin ΩT 0 cos ΩT 0⎥
⎢ ⎥ Fig. 4. RMSE in Position.
⎣0 0 0 0 1⎦

⎡ T3 T2 ⎤ CKF
M=⎢ ⎥
3 2
5th-degree CKF(Type II)
⎢ T2 ⎥ 13 5th-degree CKF(Type I)
⎣ 2
T ⎦
12.5
the position and velocity in the x and y directions are
RMSEvel(m/s)

respectively denoted by x, y, ẋ and y,

̇ Ω is the turn rate, 12
T is the time interval between time k and k − 1. A ixed
radar is placed at the origin of the plane to measure the 11.5
range r and bearing �
( ) ⎛ √ 2 ⎞
11
rk xk + y2k ⎟
=⎜ + wk (36)
�k ⎜ arctan ( y ∕x ) ⎟ 10.5
⎝ k k ⎠

where the measurement noise wk ∼ 50 60 70 80 90 100

{ N(0,}R) with zero Time(sec)
mean and covariance R = diag �r2 , ��2 . The refer-
ence data is T = 1s, Ω = −3◦ s−1 , q1 = 0.1m2 ∕s3 , Fig. 5. RMSE in velocity.
�r = 10m, q2 = 1.75 × 10−4 ∕s3 , N = 100s and �� =

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 X.-C. Zhang and Y.-L. Teng: A New Derivation of the Cubature Kalman Filters

CKF
However, there are some issues that remain
1.55
5th-degree CKF(Type II) unsolved. For example, the number of cubature points
5th-degree CKF(Type I) required for the ifth-degree CKFs are much larger than
1.5
the lower bound, whether there is a more eficient algo-
1.45 rithm or not is still under investigation. Further work
RMSEomg(Deg/s)

will be carried out in this aspect.

1.4

1.35
REFERENCES
1.3 1. Kalman, R. E., “A new approach to linear ilter-
ing and prediction problems,” Trans. ASME J. Basic
1.25
Eng., Vol. 82, No. 1, pp. 35–45 (1960).
1.2 2. Castella, F. R., “An adaptive two-dimensional
50 60 70 80 90 100
Kalman tracking ilter,” IEEE Trans. Aerosp.
Time(sec) Electron. Sys, Vol. AES-16, No. 6, pp. 822–829
(1980).
Fig. 6. RMSE in turn rate. 3. Strid, I. and K. Walentin, “Block Kalman ilter-
ing for large-scale DSGE models,” Comput. Econ.,
Vol. 33, No. 3, pp. 277–304 (2009).
Table III is the cubature points used for different 4. Ito, K. and K. Xiong, “Gaussian ilters for nonlin-
CKFs, which clearly exhibits that the ifth-degree CKF of ear iltering problems,” IEEE Trans. Autom. Control.,
type II is more eficient than that of type I. To fairly eval- Vol. 45, No. 5, pp. 910–927 (2000).
uate the algorithms, we make 200 independent Monte 5. Arulampalam, M. S., S. Maskell, N. Gordon, and
Carlo runs. Figs 4–6 show the RMSEs in the time inter- T. Clapp, “A tutorial on particle ilters for online
val 50–100s. Owing to the fact that the EKF is sensitive nonlinear/non-Gaussian Bayesian tracking,” IEEE
to information available in the previous time step and the Trans. Signal Process., Vol. 50, No. 2, pp. 174–188
UKF is often found to stop running continuously, we do (2002).
not present their results in this simulation. 6. Wu, Y., D. Hu, M. Wu, and X. Hu, “Quasi-Gaussian
Figs 4–6 indicate that the ifth-degree CKFs yield particle iltering.” Proc. 6th Int. Conf. Comput. Sci.,
very close state estimations. Their RMSEs in position Reading, UK, pp. 689–696 (2006).
and velocity are dramatically smaller than the conven- 7. Sunahara, Y. and K. Yamashita, “An approximate
tional CKF as shown in Figs 4 and 5. Also, Fig. 6 method of state estimation for non-linear dynamical
illustrates that they also marginally outperform the CKF systems with state-dependent noise,” Int. J. Control,
in turn rate, which conforms to our expectation. Vol. 11, No. 6, pp. 957–972 (1970).
8. Bucy, R. S. and K. D. Senne, “Digital synthesis
of non-linear ilters,” Automatica, Vol. 7, No. 3,
V. CONCLUSIONS pp. 287–298 (1971).
9. Julier, S. J. and J. K. Uhlmann, “A new extension of
The paper presents a new derivation of the CKFs, the Kalman ilter to nonlinear systems.” Proc. SPIE,
which uses only the cubature rules. It can readily be Orlando, USA, pp. 182–193 (1997).
shown that the conventional CKF is a special case of the 10. Julier, S. J. and J. K. Uhlmann, “Unscented ilter-
proposed algorithm. Also included in the current work ing and nonlinear estimation,” Proc. IEEE, Vol. 92,
is the higher-degree extension of the third-degree CKFs. No. 3, pp. 401–422 (2004).
Three consistency conditions that must be satisied for 11. Arasaratnam, I. and S. Haykin, “Cubature Kalman
fully symmetric cubature rules are introduced for con- ilters,” IEEE Trans. Autom. Control., Vol. 54, No. 6,
structing desired CKFs, and two representative types of pp. 1254–1269 (2009).
the ifth-degree CKFs are discussed in details. 12. Gordon, N. J., D. J. Salmond, and A. F. M.
The performance comparisons of the proposed Smith, “Novel approach to nonlinear/non-Gaussian
algorithms with conventional nonlinear ilters are Bayesian state estimation,” IEE Proc. - F, Vol. 140,
demonstrated through two target tracking problems. The No. 2, pp. 107–113 (1993).
results of both experiments show the improved perfor- 13. Kotecha, J. H. and P. A. Djuric, “Gaussian parti-
mance of the higher-degree CKF over the conventional cle iltering,” IEEE Trans. Signal Process., Vol. 51,
nonlinear ilters, such as the EKF, UKF and CKF. No. 10, pp. 2592–2601 (2003).

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd
 Asian Journal of Control, Vol. 17, No. 3, pp. 1–10, May 2015

14. Palhares, R. M. and P. L. D. Peres, “LMI approach one singular point,” SIAM J. Numer. Anal., Vol. 39,
to the mixed H2 ∕H∞ iltering design for discrete-time No. 4, pp. 1132–1145 (2001).
uncertain systems,” IEEE Trans. Aerosp. Electron. 28. Mcnamee, J. and F. Stenger, “Construction of fully
Syst., Vol. 37, No. 1, pp. 292–296 (2001). symmetric numerical integration formulas,” Numer.
15. Xiong, K., C. L. Wei, and L. D. Liu, “Robust Math., Vol. 10, pp. 327–344 (1967).
unscented Kalman iltering for nonlinear uncertain 29. Mantel, F. and P. Rabinowitz, “The application
systems,” Asian J. Control, Vol. 3, No. 1, pp. 292–296 of integer programming to the computation of
(2001). fully symmetric integration formulas in two and
16. Arasaratnam, I., S. Haykin, and T. R. Hurd, three dimensions,” SIAM J. Numer. Anal., Vol. 14,
“Cubature Kalman iltering for continuous-discrete pp. 391–425 (1977).
systems: Theory and simulations,” IEEE Trans. Sig- 30. Möller, H. M., Numerische Integration, Basel,
nal Process., Vol. 58, No. 10, pp. 4977–4993 (2010). Birkhäuser (1979).
17. Li, W. and Y. Jia, “Location of mobile station with 31. Mysovskikh, I. P., Interpolatory Cubature Formulas,
maneuvers using an IMM-based cubature Kalman Moscow, Nauka (1981).
ilter,” IEEE Trans. Ind. Electron., Vol. 59, No. 11, 32. Kitagawa, G., “Non-Gaussian state-space modeling
pp. 4338–4348 (2012). of nonstationary time-series,” J. Am. Stat. Assoc.,
18. Macagnano, D. and G. T. Freitas de Abreu, “Adap- Vol. 82, No. 400, pp. 1032–1063 (1987).
tive gating for multitarget tracking with Gaussian 33. Kotecha, J. H. and P. A. Djuric, “Gaussian parti-
mixture ilters,” IEEE Trans. Signal Process., Vol. 60, cle iltering,” IEEE Trans. Signal Process., Vol. 51,
No. 3, pp. 1533–1538 (2012). No. 10, pp. 2592–2601 (2003).
19. Tang, X., J. Wei, and K. Chen, “Square-root 34. Shalom, Y. B., X. R. Li, and T. Kirubarajan, Esti-
adaptive cubature Kalman ilter with application mation with Applications to Tracking and Navigation,
to spacecraft attitude estimation.” 15th Int. Conf. Wiley, New York, USA (2001).
Inform. Fusion, Singapore, pp. 1406–1412 (2012).
20. Jia, B., M. Xin, and Y. Chen, “Higher-degree cuba-
ture Kalman ilter,” Automatica, Vol. 49, No. 2,
pp. 510–518 (2013). Xin-Chun Zhang received his B.S. degree
21. Zhang, X. C., “A novel cubature Kalman ilter in Information and Computing Science
for nonlinear state estimation.” Proc. 52nd IEEE from University of Electronic Science
Conf. Decis. Control., Florence, Italy, pp. 7797–7802 and Technology of China (UESTC),
(2013). Chengdu, China, in 2011, where he is cur-
22. Bhaumik, S. and Swati, “Square-root cubature- rently pursuing the M.S. degree in Circuits
quadrature Kalman ilter,” Asian J. Control, and Systems. He will be an R&D engineer
Vol. 16, No. 2, pp. 617–622 (2014). of Alibaba Group from 2014. Mr. Zhang has published
23. Cools, R., “Advanced in multidimensional integra- many peer-reviewed papers in Circuits, Systems and
tion,” J. Comput. Appl. Math., Vol. 149, pp. 1–12 Signal Processing, Chinese Physics B, Contol Theory &
(2002). Applications, IEEE CDC, etc. He also acts as a reviewer
24. Lyness, J. N. and U. Rude, “Cubature of integrands of many international journals, such as IET Radar,
containing derivatives,” Numer. Math., Vol. 78, Sonar & Navigation and International Journal of Systems
No. 3, pp. 439–461 (1998). Science. His research interests include nonlinear
25. Cools, R., “Extrapolation and adaptive in software iltering, radar target tracking, adaptive control and
for automatic numerical integration on a cube,” distributed computing.
Numer. Algorithms, Vol. 34, No. 2–4, pp. 259–269
(2003). Yun-Long Teng received the Ph.D. degree
26. Genz, A. and B. D. Keister, “Fully symmetric inter- in Measurement Technology and Instru-
polatory rules for multiple integrals over ininite ment from University of Electronic Sci-
regions with Gaussian weight,” J. Comput. Appl. ence and Technology of China, Chengdu,
Math., Vol. 71, No. 2, pp. 299–309 (1996). China, in 2011. His research focus is the
27. Cools, R. and E. Novak, “Spherical product algo- technology and application of satellite
rithms and the integration of smooth functions with navigation and positioning.

© 2014 Chinese Automatic Control Society and Wiley Publishing Asia Pty Ltd