Journal of Global Positioning Systems (2007)
Vol.6, No.1: 65-73
Modified Gaussian Sum Filtering Methods for INS/GPS Integration
Yukihiro Kubo, Takuya Sato and Sueo Sugimoto
Department of El ec t rical and Electronic Engineering, Ritsumeikan University, Shiga, Japan, 525-8577
Abstract. In INS (Inertial Navigation System) /GPS
(Global Positioning System) integration, nonlinear
models should be properly handled. The most popular
and commonly used method is the Extended Kalman
Filter (EKF) which approximates the nonlinear state and
measurement equations using the first order Taylor series
expansion. On the other hand, recently, some nonlinear
filtering methods such as Gaussian Sum filter, particle
filter and unscented Kalman filter have been applied to
the integrated systems. In this paper, we propose a
modified Gaussian Sum filtering method and apply it to
land-vehicle INS/GPS integrated navigation as well as the
in-motion alignment systems. The modification of
Gaussian Sum filter is based on a combination of
Gaussian Sum filter and so-called unscented
transformation which is utilized in the unscented Kalman
filter in order to improve the treatment of the non linearity
in Gaussian Sum filter. In this paper, the performance of
modified Gaussian Sum filter based integrated systems is
compared with other filters in numerical simulations.
From simulation results, it was found that the proposed
filter can improve transient responses of the filter under
large initial estimation errors.
Key words: INS , GPS, integration, nonlinear filter
1 Introduction
In the INS/GPS integrated system, the complementary
characteristics of INS and GPS are exploited. INS
provides position, velocity and attitude information at a
high update rate with the continuous availability, and the
long term accuracy of position and velocity information
of GPS prevents th e growing navigation er rors of INS. In
other words, the navigation errors of INS are estimated
and corrected by using GPS measurements (Siouris,
1993; Grewal, 2001).
For many years, the extended Kalman filter (EKF) has
been widely utilized as the estimator in the integrated
navigation systems (Maybeck, 1979; Gelb, 1974).
Additionally, in the case of conventional navigation
systems, the initialization of INS navigation states is
completed prior to vehicle motion and then the ordinary
integrated navigation is implemented. Usually, this
initialization method needs 5 to 10 minutes and the
vehicle must be stopped. It is, however, inconvenient and
impractical when there is not enough time to stop at a
start point. Thus it is motivated to develop in-motion
alignment and navigation algorithms which can provide
the accurate attitude information while moving. Because
the initial attitude of the land-vehicle is unknown, the
attitude is usually assumed to b e 0. Thus, when the in itial
heading error is large, the nonlinear character of the INS
error equations is emphasized for in motion alignment
(Rogers, 2001). Therefore several nonlinear filtering
methods such as Monte Carlo filter (Kitagawa, 1996;
Doucet, 2000), Quasi-linear optimal filter (Sunahara,
1970), Gaussian Sum filter (Alspach, 1972) and
unscented Kalman filter (Julier, 2000), have been applied
to the integrated navigation systems. The performance
comparisons of the nonlinear filters in the integrated
navigation systems also have been reported by the
authors (Tanikawara, 2004; Fujioka, 2005; Nishiyama,
2006).
According to (Nishiyama, 2006), although Gaussian Sum
filter (GSF) works well with large uncertainties in the
initial attitude information, the linearization technique is
employed similarly to the extended Kalman filter. On the
other hand, the unscented Kalman filter (UKF) has been
recently paid much attention in the area of the integrated
navigation (Yi, 2005; An, 2005; Shin, 2007). The
unscented Kalman filter calculates the predicting mean
and covariance of the state vector from a set of samples
that are called sigma points by means of so-called
unscented transformation. In this paper, we try to
combine the GSF and the unscented transformation in
order to improve the treatment of the nonlinearity in the
GSF. With this combination, it is expected that the
transient response of the filter can be improved under
large initial estimation errors.
In this paper, firstly we briefly review the algorithms of
the nonlinear filters that are applied in this paper. Then,
the modified Gaussian Sum filtering algorithm is derived
66 Journal of Global Positioning Systems
by utilizing the unscented transformation. Finally, the
performance of EKF, GSF, UKF based and modified
Gaussian Sum filter based integrated systems is compared
in numerical simulations.
2 Description of the system
In this work, closed-loop, tightly coupled mechanization
is adopted for the INS/GPS integration. Fig. 1 shows the
architecture of the integration with major data paths
between the system components. The components of the
system are strapdown INS and GPS receiver. The INS
contains IMU (Inertial Measurement Unit: accelerometer
and gyro). Based on the measured acceleration and
angular rate, the INS computes the position, velocity and
attitude of the vehicle relative to th eir initial value at h igh
frequency. But there exist unbounded position errors that
grow slowly with time. The concept of the integrated
navigation system of Fig. 1 is to reduce the INS errors by
using some external measurement from a GPS receiver.
In this research, GPS double differenced carrier phase
and undifferenced Doppler measurements are employed
as external measurements to remove the INS errors. The
nonlinear filter estimates the errors in the navigation and
attitude information using the raw GPS data.
Fig. 1. Description of the system
2.1 Coordinate systems
To integrate the navigation systems, definitions of
coordinate systems that the navigation systems or
included sensors refer to are important. This section
defines the coordinate frames used in this paper and
represents the angular relationship between them. The
coordinate frames are defined as follows:
1) The E frame (,,)
EE E
X
YZ is the right-handed earth
fixed coordinate frame. It has the origin at the center
of the earth; the E
Z
- axis is directed toward the
North Pole; the E
X
- and E
Y- axes are located in the
equatorial plane, whereby the E
X
- axis is directed
toward the Greenwich Meridian. It is used for the
definition of position location such as latitude and
longitude.
2) The L frame (,,)
LL L
X
YZ is the right-handed
locally level coordinate frame. The L
X
- and L
Y-
axes are directed toward local north and east
respectively; L
Z
- axis is downward vertical at the
local earth surface referenced position location. It is
used for defining the angular orientation of the local
vertical in the E frame.
3) The C frame (,,)
CC C
X
YZ is the right-handed
computer frame that is defined by rotating the L
frame around negative L
Z
- axis by the “wander
angle”
α
; toward the negative L
Y- axis and the C
Z
-
axis is directed toward the negative L
Z
- axis
(upward vertical). It is used for integrating
acceleration into velocity, and used as the reference
for describing the strapdown sensor coordinate frame
orientation.
4) The B frame (,,)
BBB
YZ is the strapdown inertial
sensor coordinate frame (body frame). The B
X
- axis
is directed toward the head of the vehicle; the B
Y-
axis is the right-hand of the vehicle; the B
Z
- axis is
downward vertical to the B
X
-B
Y plane. The frame
is fixed on the vehicle and rotates with th e motion of
the vehicle.
Fig. 2. Coordinate frames
Fig. 2 shows the spatial image of the E, L and C frames,
where
λ
and
ϕ
represent the longitude and the latitude
respectively. In the inertial computations, the acceleration
sensed with respect to th e B frame have to be transformed
onto the C frame. The velocity and po sition of the ve hicle
are then computed with respect to the C frame. Such a
transformation is known as the Euler angle
transformation. We define the product of direction cosine
matrix for this transformation as C
B
T. Then the
coordinates (,,)
BBB
x
yz in the B frame are transformed
into (,,)
CCC
x
yz in the C frame as follows:
Kubo et al.: Modified Gaussian Sum Filtering Methods for INS/GPS Integration 67
CB
C
CBB
CB
x
x
yTy
zz
⎡⎤ ⎡⎤
⎢⎥ ⎢⎥
=
⎢⎥ ⎢⎥
⎢⎥ ⎢⎥
⎣⎦ ⎣⎦
(1)
where C
B
T is the direction cosine matrix.
3 INS error model
The direction cosine matrix C
E
T is represents the
transformation from the E frame to the C frame, and it
can be decomposed as follows.
CCL
ELE
TTT= (2)
On the other hand, the computed matrices C
L
T and L
E
T
contain errors C
L
T
δ
and L
E
T
δ
respectively. The error
E
C
T
δ
can be formulated as
{[( )]( )}
CCC
EEE
CL CL
LE LE
CCL
LLLLE
L
E
TTT
TT TT
TIrTr T
T
δ
δδ δ
≡−
=−
=−×−×
=R
(3)
where L
r
δ
[,Lx
r
δ
,,Ly
r
δ
,0]T is horizontal angular
position error, and the relation of L
E
T=[()]
L
LE
I
rT
δ
−× is
used in the calculation of equation (3) with the
assumption that ,Lx
r
δ
and ,Ly
r
δ
are small. Also, ()a
×
for 31× vector a= [
x
a,y
a,z
a]Tis the skew- symmetric
matrix defined by
0
()0 0
zy
zx
yx
aa
aa a
aa
⎡⎤
⎢⎥
×≡ −
⎢⎥
⎢⎥
⎣⎦
(4)
And R is the position error matrix defined as follows
(Rogers 2003, 2001).
,,
,,
,,
cos sincossin
sin coscossin
0
Cy Cx
Cx Cy
Ly Lx
rr
rr
rr
δ
αδαδδαδ δα
δ
αδαδδαδδα
δδ
⎡⎤
−−
⎢⎥
≡− −−−
⎢⎥
⎢⎥
⎢⎥
⎣⎦
R (5)
where
sinsin() sin
δ
ααδαα
≡+− (6)
coscos() cos
δ
ααδαα
≡+− (7)
According to (Rogers, 2001, 2003), we have
///
()()( )()
LCCCC
EL LL
EC EC
TT
ωω ω
=×−×++×
RRR (8)
where the dot above a letter denotes differentiation with
respect to time, the vector /
L
EL
ω
is the rotation rate of the
L frame with respect to the E frame in the L frame
coordinate system, and the vector /
C
EC
ω
is similarly
defined. From equation (8), the position error (,Cx
r
δ
,
,Cy
r
δ
) as well as azimuth error (
δ
α
) equations can be
derived.
3.1 Velocity error model
The computed velocity C
v also contains the velocity
error C
δ
v such that
CC C
vv v
δ
=+ (9)
and the velocity equation is given by
(2)
CC CCCC
vf vg
ρ
=
−+Ω×+ (10)
where C
f
is non-gravitational specific force vector,
ρ
is
relative rate vector, and C
Ω
is earth rate vector. The
specific force is proportional to the inertial acceleration
of the system due to all forces except gravity measured
by the accelerometer. C
g
is the gravity vector, positive
toward the centre of the earth in the C frame. From
equations (9) and (10), the velocity error is modelled by
(2)
(2)
(2)
CCC CCCC
CCC
CCCC
vbf v
v
vg
δδθδρδ
ρδ
δρ δδ δ
=
+× +×+Ω
−+Ω×
−+Ω×+
(11)
where C
δ
θ
[,Cx
δ
θ
,,Cy
δ
θ
,,Cz
δ
θ
]T is the attitude error.
3.2 Attitude error model
The attitude error C
δ
θ
causes the error of the
transformation matrix C
B
T. The computed matrix C
B
T
which contains the attitude error is formulated by
[( )]
CC
BCB
TI T
δθ
=−× (12)
Therefore, we have following attitud e error model.
// /
CC C
CCC
ECIEIC d
δθδωδωδθ ω
=
++×+
(13)
where C
d denotes gyro drift.
68 Journal of Global Positioning Systems
3.3 Sensor error model
In this paper, the accelerometer bias B
b and gyro bias
B
d are modelled as the first order Markov processes
respectively as follows:
1
()() ()
1
()()()
BBb
b
BBd
d
btbt ut
dtdt ut
τ
τ
=− +
=− +
(14)
where b
τ
and d
τ
are the correlation time constants and
()
b
ut, ()
d
ut are zero mean Gaussian white noise
processes.
3.4 State equation
In order to implement the nonlinear filtering for
integrated navigation, here, we define the state vector.
Because the double differenced carrier phases are used as
the measurements in this paper, the unknown integer
ambiguities should be simultaneously estimated.
Therefore the state vector is defined such that it includes
the INS errors as well as the integer ambiguities as
follows:
,, , , , ,
,, ,,,,,
1,
1,2 1,3
,, ,
,,,,,,,
,,,, ,, ,,,
,,,,
CxCyCx Cy Cx CyC
Cz BxByBzBxByBz
ns
ku kuku
rrvv h
vbbbd dd
ctN NN
δδδδδθ δθδ
δ
γβ
δ
=
x
T
(15)
where 1,2
,ku
N denotes the double differenced integer
ambiguity of the satellites 1, 2 and the receivers k, u, and
s
n is the number of visible satellites.
β
and
γ
are
defined as follows:
cos 1
sin
β
δα
γδα
≡−
The descriptions of the state vector components are listed
in Table 1. Then, from equations (8), (11), (13) and (14),
the state equation can be formulated by
()( (),)()
x
tfxttt
η
=+
(16)
where (•, )
f
t is the time-varying nonlinear function, and
the process noise ()t
η
is assumed to be mutually
independent zero mean Gaussian white noise with
covariance matrix ()Nt.
Table 1. List of states
No. Symbol Error state
1 ,Cx
r
δ
C
X
-axis position error in angle
2 ,Cy
r
δ
C
Y-axis position error in angle
3 ,Cx
v
δ
C
X
-axis velocity error
4 ,Cy
v
δ
C
Y-axis velocity error
5 ,Cx
δ
θ
C
X
-axis tilt error
6 ,Cy
δ
θ
C
Y-axis tilt error
7
γ
sin
δ
α
8
β
cos 1
δ
α
9 C
h
δ
C
Z
-axis altitude error
10 ,Cz
v
δ
C
Z
-axis velocity error
11 ,Cx
b B
X-axis accelerometer bias
12 ,Cy
b B
Y-axis accelerometer bias
13 ,Cz
b B
Z
-axis accelerometer bias
14 ,Cx
d B
X
-axis gyro bias
15 ,Cx
d B
Y -axis gyro bias
16 ,Cx
d B
Z
-axis gyro bias,
17 1,2
,ku
N double differenced ambiguity
By discretizing the state equation (16), we have
(1) ()((),)()
x
kxkfxkktwk
+
=+Δ+ (17)
where ()wk is assumed to be Gaussian white noise with
zero mean and diagonal covariance matrix Q(k), and t
Δ
is a sampling interval of the measurement data.
3.5 Measurement equation
In this paper, the measurements are the double
differenced carrier phase and Doppler data. By ignoring
some errors in the carrier phase data such as the
remaining ionospheric and tropospheric delays and
multipath errors, then the double differenced carrier
phase measurement can be simply modelled by
()( )()
ku g
E
yk hrk
N
λ
ε
=
++ (18)
where
E
r is the position vector in the E frame, the
function h is the nonlinear function that indicates the
distance between satellites and receivers, ku
N is the
ambiguity vector,
λ
is the wave length, and
g
ε
is the
measurement noise.
By linearlizing equation (18) with the first order Taylor
series approximation around the position indicated by
INS, ()
i
E
rk, and applying appropriate transformations of
the coordinate systems, we obtain the measurement
Kubo et al.: Modified Gaussian Sum Filtering Methods for INS/GPS Integration 69
equation of the INS position error in the C frame as
follows.
()()( ())
ˆ()()()
i
E
ku g
C
ykyk hk
r
H
krk k
N
δλ
ε
≡−
=++
(19)
where
ˆ()()() () ()()
EL
L ABC
H
kHkTkT kT kTk≡− (20)
and
()()
(())
() ()
()00 010
00100
cos 001
001
i
EE
E
Ekk
p
p
AB
hr k
Hk rk rr
RhRh
TT
λ
=
⎡⎤
≡,
⎢⎥
⎣⎦
−+
⎡⎤
⎡⎤
⎢⎥
+⎢⎥
⎢⎥
≡,≡
⎢⎥
⎢⎥
⎣⎦
⎢⎥
⎣⎦
where p
R is the earth radius.
The Doppler measurement can be modelled by the
change of the distance between the receiver and the
satellite in the sampling interval tΔ (Misra, 2001). By
using the velocity error vector in the C frame, C
δ
v, and
the appropriate transformations similarly to the above
derivations, the Doppler measurement can be formulated
as
{}
()() ()()
Ep
i
CC d
E
uu u
C
Dk Gkctk
vvv
T
δδ
ε
=−−++1
(21)
where
T
T
TT T
12
T
TT T
()() ()
s
i
E
E
p
pun
uu
uu uE
EEE E
CC C
C
r
Gg
gg grrr
TT T
T
=
⎡⎤
⎢⎥
⎣⎦
⎡⎤
⎡⎤
≡,≡,
⎢⎥
⎣⎦
⎣⎦
and
p
u
r denotes the distance between the receiver u and
the satellite p, p
E
v is the velocity of the satellite p in the E
frame and d
ε
is the measurement noise. Then, we have
the following Doppler measurement equation.
()
1
E
pi
c
E
uuuu
C
C
Ed
uCu
DkDGG
vv
T
v
GTct
δ
ε
δ
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
≡+ −
⎡⎤
=− +
⎣⎦
(22)
Finally, from equations (19) and (22), we have the
measurement equation for the integrated navigation in
general form:
()()() ()zk Hkxkk
ε
=+ (23)
where TTT
() [()()]
u
zkkk
yD
.
4 Nonlinear f iltering
Nonlinear filtering techniques are applied to the
integrated INS/GPS system in order to estimate the state
vector (the errors of INS described above). In this section,
firstly, we briefly review the filter algorithms of the GSF
and the UKF. Then the modified Gaussian Sum filter
(MGSF) algorithm is derived.
4.1 Gaussi an Sum filteri ng
Let k
Z
be the set of the measurement such that
{(1)(1)… ()}
k
Z
zz zk
=
,,, (24)
In the GSF (Alspach, 1972), a posteriori probability
density (() )
k
pxkZ
|
is formed by the convex
combination of the outputs of several Kalman filters
processed in parallel. The a priori density 1
(() )
k
pxkZ
|
is assumed that it is formulated by the sum of several
normal distributions as follows:
11
(() )(1)
(1)( 1)
mj
kj
jj
pxkZk k
NkkPkk
μ
γ
μ
=
|= |−
×|−,|−
(25)
where m is the number of distributions, and
j
γ
is the
weight for the j-th distribution such that
1
(1)1(1)0
mjj
j
kk kk
γγ
=
|
−=, |−≥
And
[
]
N
θ
,
Σ denotes the normal probability density
function with mean
θ
and covariance matrix
Σ
. Then,
by the Bayesian recursion relations, a posteriori density
can be formulated by
1
(())() ()()
mj
jj
kj
pxkZkkNkk Pkk
μ
γμ
=
⎡⎤
|
=| |,|
⎣⎦
(26)
where
()
() (1)
() ()()(1)
jj
j
j
kk kk
Kkzk Hkkk
μ
μμ
μ
|= |−
+
|
()( 1)()()( 1)
jj jj
Pkk PkkKkHkPkk
μμ μμ
|
=|−−|−
T
1
T
()(1) ()
() (1)()
()
jj
j
Kk PkkHk
HkPkkRk
Hk
μμ
μ
=|−
|− +
and the weight ()
jkk
γ
|
is given by
1
(1)()
() {(1) ()}
jj
jmll
l
kk k
kk kk k
γβ
γγβ
=
|−
|= |−
(27)
70 Journal of Global Positioning Systems
where
T
(1)(11)
()( 1)
(1)()()(1)
()(1) ()()
jj
jj
j
j
j
jj
kkkk
kNkk P
kkzk Hkkk
PHkPkkHk Rk
νν
νν μ
γγ
βν
μ
ν
⎡⎤
⎢⎥
⎣⎦
|− =−|−
=|−,
|
−≡ −
|
≡|−+
Therefore, we have the filtered estimator
1
ˆ() ()()
mj
j
j
x
kkkk kk
γμ
=
|
=
||
(28)
The a priori density (( 1))
k
pxkZ+
|
can be rewritten with
the same algorithm as the EKF as follows.
1
(( 1))(1)
(1) (1)
mj
kj
jj
pxkZk k
NkkPkk
μ
γ
μ
=
+| =+|
⎡⎤
×+
|
,+
|
⎣⎦
(29)
where
(1)(( ))
jj
kkfkk
μμ
+| =| (30)
T
(1)()( ) ()()
jjj
j
PkkFkPkkk Qk
F
μμ
+| =|+ (31)
()
()
() j
j
kk
fx
Fk xx
μ
=|
⎡⎤
=⎢⎥
⎣⎦ (32)
(1)( )
jj
kk kk
γγ
+
|
=
|
4.2 Unscented Kalman filter
In the UKF, the predict mean ˆ(1|)
x
kk+ and covariance
(1|)Pk k+ are calculated from a set of samples which is
called the sigma points. This method is called the
unscented transformation (Julier, 2000). Under the
assumption that the system noise is independent and
additive, the predict mean and covariance are computed
as following steps.
Step1: choose the sigma points ()
jkk
χ
| which is
associated with the n-dimensional state vector ()
x
k as
follows.
0
0ˆ
()()
ˆ
()()()()
ˆ
()() ()()
1(12…)
2( )
j
j
jn
j
jjn
kk xkkWn
kk xkknPkk
kk xkknPkk
WWj n
n
κ
χκ
κ
χ
κ
χ
κ
⎛⎞
⎜⎟
⎜⎟
⎝⎠
⎛⎞
⎜⎟
⎜⎟
+⎝⎠
+
|= |,=
+
|= |++|
|= |−+|
==,=,,,
+
Step2: compute a set of transformed samples through the
process model equation (17),
(1)(())
jj
kkfkkk
χχ
+
|= |,
Step3: compute the predicting mean and covariance as
follows
2
0
ˆ(1) (1)
n
jj
j
x
kk Wkk
χ
=
+
|= +|
2T
0
(1) ()
n
jjj
j
Pk kWQk
χχ
=
+| =+

ˆ
where(1 )(1)
jj
kkxkk
χχ
+|−+|
j
W is the weight of the j-th point and
κ
is a scaling
parameter. (()( ))
j
nPkk
κ
+|
is the j-th column of the
matrix square root of ()()nPkk
κ
+
|. Then, once the
observation (1)zk
+
is obtained,ˆ(1)
x
kk+| and (1)Pkk
+
|
are updated to ˆ(1 1)xk k
+
|+ and (1 1)Pk k+| + as
follows.
(1) ()(1)
jj
kkHkkk
Z
χ
+
|= +|
2
0
ˆ(1) (1)
n
j
j
j
zk kWk k
Z
=
+
|= +|
(33)
2T
0
(1) ()
n
jjj
j
Pk kWRk
ZZ
νν
=
+| =+
 (34)
2T
0
(1) n
xj
j
j
j
Pk kW
Z
ν
χ
=
+| =
(35)
ˆ
where(1 )(1 )
j
jkkzkk
Z
Z
+| −+|
1
(1) (1)(1)
x
K
kPkkPkk
ννν
+
=+
|
+
|
(36)
ˆˆ ˆ
(11) (1)()(()(1))
x
k kxk kKkzkzk k
+
|+ =+|+−+| (37)
T
(11) (1)()(1)()Pk kPk kKkPkkKk
νν
+| +=+|−+| (38)
Since the measurement equation (23) is linear in this
navigation problem, above equations (33)-(35) can be
simply expressed by
ˆ
ˆ(1)()(1)zk kHkxk k
+
|= +| (39)
T
(1)()(1)()Pk kHkPk kHkR
νν
+
|=+| + (40)
T
(1) (1)()
x
PkkPkkHk
ν
+| =+ (41)
4.3 Modified Gaussian Sum filter
In the Gaussian Sum filtering algorithm, we can see from
equations (30)-(32) that the linearization technique is
Kubo et al.: Modified Gaussian Sum Filtering Methods for INS/GPS Integration 71
employed similarly to the extended Kalman filter. In this
paper, we propose the modified Gaussian Sum filter by
applying the unscented transformation algorithm to the
time updating algorithm of the GSF, equations (30)-(32).
Step1: similarly to the step 1 of the UKF, for j-
th (12…)jN=,, , density in GSF, choose the sigma points
and weights as follows.
0
0()()
jj
kkkkWn
κ
χμ κ
|= |,=
+
()() ()()
jj j
ll
kkkkn Pkk
μ
κ
χμ
⎛⎞
⎜⎟
⎜⎟
⎝⎠
|=|++|
()() ()()
jj j
ln l
kkkkn Pkk
μ
κ
χμ
⎛⎞
⎜⎟
⎜⎟
+⎝⎠
|=|−+|
1(12…)
2( )
lln
WWl n
n
κ
+
==,=,,,
+
Step2: compute a set of transformed samples through the
process model equation (17),
(1)(( ))
jj
ll
kkfkkk
χχ
+|=|,
Step3: compute the j-th predicting mean and covariance
as follows.
2
0
(1) (1)
n
jj
ll
l
kk Wkk
μχ
=
+| =+|
(42)
2T
0
(1)() ()
njj
jlll
l
Pk kWQk
μ
χχ
=
+| =+
(43)
where(1 )(1 )
jj j
ll
kk kk
χχ μ
≡+|−+|
In the MGSF, the original time updating algorithm of
equations (30) and (31) are substituted by (42) and (43)
respectively.
5 Experimental results
The experiments of the INS/GPS In-Motion Alignment
and navigation algorithms described above were carried
out by using simulated INS and GPS data. In the
experiments, we assume the vehicle runs at a speed of
around 15 [km/h] for about 10 minutes. The speed at the
start point was 0 [km/h], and the initial azimuth angle was
60 [deg]. The test trajectory in the local level horizontal
plane is shown in Fig. 3. The data were obtained by
utilizing the Matlab6.5 and INS Toolbox1.0 (GPSoft
LLC.) at 50 [Hz] rate for IMU and at 1 [Hz] rate for GPS.
Four types of filters, i.e. EKF, GSF, UKF and MGSF are
used in the experiments and compared. The nonlinearity
of the INS usually occurs when there exist large attitude
errors. So in the experiments, the initial state estimates
are set to have large azimuth error. And we assume that
there exist no errors in the other initial estimates.
Therefore, in the EKF and UKF, the initial estimate
ˆ(0 1)
x
|
is set to 0, and (01)P|− and Q are configured
from the nominal equipment specifications in Table 2. In
this case, the states related to the azimuth error, i.e. 7th
and 8th components of the state vector have 60 [deg]
initial estimation error respectively.
Fig. 3. Test trajectory
Table 2. Sensor error sp e c if i ca t i on
In the GSF and MGSF, three normal distributions are
utilized, i.e. 3m
=
, and (01)123
j
Pj
μ
|−,=,, are set to the
same value of the EKF and UKF, i.e.
(0 1)(0 1)
j
PP
μ
|
−=|−. The initial estimates (01)
j
μ
|
−,
123j
=
,, are also set to 0 except for the 7th and 8th
components of the state vector (see Table ),
β
and
γ
,
that represent the azimuth error. They are assumed to
have the initial azimuth error estimates such that
60 060
δ
α
=
−,,+ [deg].
The processing r esults are shown in Fig. (4)-(7). Figs. (4 )
and (5) show the results of the positioning and
comparison of the positioning errors. Table 3 also shows
RMS (Root Mean Square) values of the position errors.
From Fig. (5) we can see that the MGSF shows faster
convergence than the others, and the GSF and MGSF
show better performances than EKF and UKF. Therefore,
the GSF and MGSF can work well when there exist large
Accelerometer Specification
Bias 80
[
μ
G](1 )
σ
Scale factor 150 [ppm](1 )
σ
Random error 0.0003 [m/s]2
Gyroscope Specification
Bias 20
[deg/h](1 )
σ
Scale factor 500 [ppm](1 )
σ
Random error 0.06 [deg/ h]
72 Journal of Global Positioning Systems
azimuth error because they can treat large azimuth error
by assuming multiple initial error distributions. From
Table 3, we can also see that the MGSF achieves the best
performance in this simulation.
Fig. 4. Positioning results
Fig. 5. Positioning errors
Table 3. Root mean square of position errors
North error [m] East error [m]
EKF 0.08 0.61
UKF 0.09 0.10
GSF 0.02 0.13
MGSF 0.02 0.09
Table 4. Root mean square of velocity errors
North error [m/s] East error [m/s]
EKF 0.36 0.25
UKF 0.34 0.19
GSF 0.36 0.22
MGSF 0.34 0.19
Fig. 6 and Table 4 show the velocity errors and their
RMS values respectively. From these figures and tables,
the all filters show almost same performance, whereas th e
UKF and MGSF show slightly better performance than
the EKF and GSF.
Finally, Fig. 7 shows the results of the azimuth errors.
From Fig. 7, we can see that all filters show almost same
results after 200 [sec], but the MGSF shows faster
convergence in its transien t response from 0 to 200 [sec].
Therefore, from these results of the simulation, we can
consider that the UKF and MGSF can achieve better
performance than the EKF and GSF, and the MGSF can
work well when there exist a large initial azimuth error.
Fig. 6. Velocity erro rs
Fig. 7. Azimuth errors
6 Conclusion s
In this paper, the modified Gaussian Sum filtering
algorithm was derived by applying the unscented
transformation to the Gaussian Sum filter, and it was
applied to the GPS/INS integrated system. The algorithm
was tested and compared with the EKF, GSF and UKF by
using simulated data. From the experimental results, it
was found that the derived MGSF show the quick
transient response for azimuth error estimation. Therefore
the MGSF has an ability to improve the navigation
performance, when there are large initial azimuth errors.
Kubo et al.: Modified Gaussian Sum Filtering Methods for INS/GPS Integration 73
References
Alspach D. L. and Sorenson H. W. (1972) Nonlinear Bayesian
Estimation Using Gaussian Sum Approximations, IEEE
Trans. on Automatic Control, Vol. AC-17, No. 4, pp. 439-
448, 1972.
An D. and Liccardo D. (2005) A UKF Based GPS/DR
Positioning System for General Aviation, Proc. of the
Institute of Navigation, ION GNSS 2005, pp. 990-998,
Long Beach, CA, 2005.
Doucet A., Godsill J. and Andrieu C. (2000) On sequential
Monte Carlo sampling methods for Bayesian filtering,
Statistics and Computing, Vol. 3, pp. 197-208, 2000.
Fujioka S., Tanikawara M., Nishiyama M., Kubo Y. and
Sugimoto S. (2005) Comparison of Nonlinear Filtering
Methods for INS/GPS In-Motion Alignment, Proc. of the
Institute of Navigation, ION GNSS 2005, pp. 467-477,
Long Beach, CA, 2005.
Gelb A. (1974) Applied Optimal Estimation, MIT Press,
Massachusetts, 1974.
Grewal M. S., Weill L. R. and Andrews A. P. (2001) Global
Positioning Systems, Inertial Navigation, and
Integration, John Wiley & Sons, New York, 2001.
Julier S., Uhlmann J. and Durrant-Whyte H. F. (2000) A New
Method for the Nonlinear Transformation of Means and
Covariances in Filters and Estimators, IEEE Trans. On
Automatic Control, Vol. 45, No. 3, pp. 477-482, 2000.
Kitagawa G. (1996) Monte Carlo Filter and Smoother for
Non-Gaussian Nonlinear State Space Models, Jounal of
Computational and Graphical Statistics, Vol. 5, pp. 1-25,
1996.
Maybeck P. S. (1979) Stochastic Models, Estimation and
Control (Mathematics in Science and Engineering),
Academic Press, New York, 1979.
Misra P. and Enge P. (2001) Global Positioning System --
Signals, Measurements, and Performance, Ganga-Jamuna
Press, Massachusetts, 2001.
Nishiyama M., Fujioka S., Kubo Y., Sato T. and Sugimoto S.
(2006) Performance Studies of Nonlinear Filtering
Methods in INS/GPS In-Motion Alignment, Proc. of the
Institute of Navigation, ION GNSS 2006, pp. 2733-2742,
Fort Worth, TX, 200 6.
Rogers R. M. (2001) Large Azimuth INS Error Models for In-
Motion Alignment Land-Vehicle Positioning, Proc. of the
Institute of Navigation National Technical Meeting 2001,
pp. 1104-1114, Long Beach, CA, 2001.
Rogers R. M. (2003) Applied Mathematics in Integrated
Navigation Systems, 2nd edition, AIAA, Virginia, 2003.
Shin E.-H. and El-Sheimy N. (2007) Unscented Kalman Filter
and Attitude Errors of Low-Cost Inertial Navigation
Systems, Navigation: Journal of The Institute of
Navigation, Vol. 54, No. 1, pp. 1-9, Spring, 2007.
Siouris G. M. (1993) Aerospace Avionics Systems A Modern
Synthesis, Academic Press, San Diego, 1993.
Sunahara Y. (1970) An Approximate Method of State
Estimation for Nonlinear Dynamical Systems, Journal of
Basic Engineering, Trans. ASME, Series D, Vol. 92, No. 2,
pp. 385-393, 1970.
Tanikawara M., Asaoka N., Oiwa M., Kubo Y. and Sugimoto S.
(2004) Real-Time Nonlinear Filtering Methods for
INS/DGPS In-Motion Alignment, Proc. of the Institute of
Navigation ION GNSS 2004, pp. 1104-1114, Long Beach,
CA, 2004.
Yi Y. and Grejner-Brzezinka D. (2005) Nonlinear Bayesian
Filter: Alternative To The Extended Kalman Filter In
The GPS/INS Fusion Systems, Proc. of the Institute of
Navigation, ION GNSS 2005, pp. 1392-1400, Long Beach,
CA, 2005.