# An Efficient Nonlinear Filter for Spacecraft Attitude Estimation

Authors: Bing Liu; Zhen Chen; Xiangdong Liu; Fan Yang
Journal: International Journal of Aerospace Engineering
 (2014)

---

## Abstract

Increasing the computational efficiency of attitude estimation is a critical problem related to modern spacecraft, especially for those with limited computing resources. In this paper, a computationally efficient nonlinear attitude estimation strategy based on the vector observations is proposed. The Rodrigues parameter is chosen as the local error attitude parameter, to maintain the normalization constraint for the quaternion in the global estimator. The proposed attitude estimator is performed in four stages. First, the local attitude estimation error system is described by a polytopic linear model. Then the local error attitude estimator is designed with constant coefficients based on the robustH
2 filtering algorithm. Subsequently, the attitude predictions and the local error attitude estimations are calculated by a gyro based model and the local error attitude estimator. Finally, the attitude estimations are updated by the predicted attitude with the local error attitude estimations. Since the local error attitude estimator is with constant coefficients, it does not need to calculate the matrix inversion for the filter gain matrix or update the Jacobian matrixes online to obtain the local error attitude estimations. As a result, the computational complexity of the proposed attitude estimator reduces significantly. Simulation results demonstrate the efficiency of the proposed attitude estimation strategy.


---

## Body

## 1. Introduction

Attitude determination is a very important part for a spacecraft to achieve its designed mission. There are various methods for the spacecraft attitude determination. They can be divided into two classes: deterministic method and optimal estimation method [1]. The deterministic method, the TRIAD algorithm, uses a minimal set of date and then solves three possibly nonlinear equations to obtain the attitude [1]. It is simple and elegant; however, it is suboptimal and of limited use because it makes use of only two unit-vector measurements and ignores one piece of information from one of the unit vectors [2]. The optimal estimation method is based on the solutions to the Wahba’s problems, which obtains the optimal attitude estimation by minimizing an appropriate loss function.There are many nonlinear estimation algorithms for the spacecraft attitude estimation since it is essentially a nonlinear problem. The most widely used algorithm for real time attitude estimation is the EKF. The EKF is recursive and easy to implement, but the accuracy can be surprisingly bad in the cases that the dynamic and the measurement models have highly nonlinearities or the system is with large process noise [3, 4]. The poor performance has driven several nonlinear filters for attitude estimation, among which the sigma point filters have attracted much attention, such as the Unscented Kalman Filter (UKF) [5], the Cubature Kalman Filter (CKF) [6], the Gauss-Hermite Quadrature Filter (GHQF) [7], and the Particle Filter (PF) [8]. They deal with the nonlinear functions directly by choosing some points to approximate the probability density function of the nonlinear functions according to certain rules. It is generally believed that the sigma point filters are more accurate than the EKF; nonetheless, the computational cost of the sigma point filter seems high for engineering implementation [9]. Even the implementation of the EKF is also computationally complex [10], because the Jacobian matrixes are required to update online which can be a very cumbersome and error-prone process, and it needs to calculate the matrix inversion for the gain matrix, resulting in heavy burden of the onboard computer especially for the systems with high dimension.Several alternatives on the computational cost of the filter have been developed for attitude estimation. Wei used the optimal-REQUEST to estimate the attitude and the UKF to estimate the gyro drifts [11]. The computational cost of the attitude estimator was reduced by setting the state dimension to three rather than six. Fan and Kiani improved the real-time performance of the attitude estimator by minimizing the number of the required sigma points [9, 12]. Tang et al., Miao et al., and Choukroun et al. presented a reduced quaternion measurement model without losing information to reduce the computational complexity of the attitude estimator [6, 13, 14]. The improved methods mentioned above have reduced the computational cost of the attitude estimator efficiently; nevertheless, the improved sigma point methods remain to have high computational cost for spacecraft which is small and with limited computational abilities. The calculation of the matrix inversion remains in the attitude estimator on the bias of reducing the dimension of the measurement model.In this paper, an efficient nonlinear filtering method is developed for spacecraft attitude estimation. By introducing the polytopic linear differential inclusion (PLDI) theory given by Boyd et al. [15], the local attitude estimation error system is represented by an uncertain polytopic linear model. This leads to the local error attitude estimator designed with constant filter coefficients, without calculating matrix inversion or updating the Jacobian matrixes online. Thus the computational cost is sharply reduced.The rest of the paper is organized as follows: Section2 briefly introduces the attitude kinematics and the sensor models. Section 3 presents the implementation of the efficient nonlinear attitude estimator in detail. Section 4 demonstrates the performance of the attitude estimator and compares the results of this method with the Multiplicative EKF (MEKF) and other filters. Section 5 gives the conclusion remarks.


## 2. Attitude Kinematics and Sensor Models

In this section, the attitude kinematics and the sensor models are briefly introduced.The spacecraft attitude can be described by various parameters, such as the direction cosine matrix, the principal axis and angle, the Euler angels, quaternion, Rodrigues parameters (RPs), and the modified Rodrigues parameters (MRPs). The most widely used attitude parameter is the quaternion because of its nonsingular character for any arbitrary rotation angle and its bilinear kinematic equation.The quaternion is defined as(1)
q
=
[
q
→
q
4
]
=
[
e
→
sin
(
θ
2
)
cos
⁡
(
θ
2
)
]
,
where e
→ is principal rotation axis and θ is the corresponding rotation angle.The attitude kinematic equation in the quaternion form is given by(2)
q
˙
=
1
2
Ω
(
ω
)
q
=
1
2
Ξ
(
q
)
ω
,
where
(3)
Ξ
(
q
)
=
[
q
4
I
3
×
3
+
[
q
→
×
]
-
q
→
T
]
,
Ω
(
ω
)
=
[
-
[
ω
×
]
ω
-
ω
T
0
]
,
and ω is the angular velocity of the spacecraft.The gyro is commonly used to measure the angular velocity of the spacecraft; a general gyro model is given by(4)
u
=
ω
+
b
+
n
g
,
where n
g is the measurement noise and b is the drift rate bias driven by a white noise
(5)
b
˙
=
n
b
,
where n
b is assumed to be white noise.In most practical applications, a typical attitude estimation system for the spacecraft comprises several gyros and vector sensors, such as the sun sensor, star sensor, and magnetometer. Therefore, the vector observation is chosen as the attitude measurement for the most general case. The measurement model for a signal vector observation is described as(6)
r
b
j
=
R
(
q
)
r
i
j
+
v
j
,
where R
(
q
) is the attitude matrix, r
b
j is the observation vector, and r
i
j is the known reference vector. v
j is the measurement noise assumed to be white noise, and the superscript j denotes the index of the observation vector.


## 3. The Computationally Efficient Attitude Estimator

The spacecraft attitude state is given by the attitude quaternion and the gyro drift rate bias(7)
x
=
[
q
T
b
T
]
T
.Then the dynamic model for the attitude estimation system is(8)
x
˙
=
[
q
˙
b
˙
]
=
[
1
2
Ω
(
u
-
b
)
q
0
3
×
1
]
+
[
-
1
2
Ξ
(
q
)
n
g
n
b
]
.Denote the error attitude state as(9)
δ
q
=
q
⊗
q
^
-
1
,
δ
b
=
b
-
b
^
,
where q
^ and b
^ are the predicted attitude quaternion and gyro drift rate bias and q
^
-
1 is the inversion of q
^. The dynamic model for the attitude estimation error system is [16]
(10)
δ
q
˙
=
[
-
[
ω
^
×
]
δ
q
→
0
]
+
1
2
Ω
(
δ
ω
)
δ
q
,
δ
b
˙
=
n
b
,
where ω
^ is the predicted angular velocity; the angular velocity error is expressed as follows:
(11)
δ
ω
=
ω
-
ω
^
=
u
-
b
-
n
g
-
(
u
-
b
^
)
=
-
δ
b
-
n
g
.Substitute the above equation into (10); the dynamic model for the attitude estimation error system can be rewritten as
(12)
δ
q
˙
=
[
-
[
ω
^
×
]
δ
q
→
0
]
-
1
2
[
δ
q
4
δ
b
+
[
δ
q
→
×
]
δ
b
δ
q
→
T
δ
b
]
-
1
2
Ξ
(
δ
q
)
n
g
,
δ
b
˙
=
n
b
.The quaternion must obey a normalization constraint. In order to ensure that the quaternion maintains the normalization constraint, the most common attitude quaternion estimation method uses an unconstrained three-component vector to represent the local attitude error. In this paper, the RPs are chosen as the local attitude error parameter. The error RPs is defined in the terms of the error quaternion by(13)
δ
g
=
δ
q
→
δ
q
4
,and the inverse transformation is given by(14)
δ
q
→
=
δ
g
1
+
δ
g
T
δ
g
,
δ
q
4
=
1
1
+
δ
g
T
δ
g
.If the attitude error is quite small or tends to be zero, (13) can be approximated as
(15)
δ
g
=
δ
q
→
.Therefore, the dynamic model (12) can be easily approximated to the first order form by using the error RPs
(16)
δ
x
˙
=
[
δ
g
˙
δ
b
˙
]
=
F
[
δ
g
δ
b
]
+
n
w
,
where
(17)
F
=
-
[
[
ω
^
×
]
1
2
I
3
0
3
0
3
]
,
n
w
=
[
-
1
2
n
g
n
b
]
.The corresponding measurement equation of the attitude estimation error system is(18)
δ
r
b
=
r
b
-
r
^
b
=
H
δ
x
+
v
,
where
(19)
r
b
=
[
r
b
1
⋮
r
b
j
⋮
r
b
n
]
,


r
^
b
=
[
r
^
b
1
⋮
r
^
b
j
⋮
r
^
b
n
]
,


v
=
[
v
1
⋮
v
j
⋮
v
n
]
,


H
=
[
h
1
⋮
h
j
⋮
h
n
]
,
where h
j
=
[
[
2
(
R
(
q
^
)
r
i
j
)
×
]
0
3
]
,
j
=
1
,
…
,
n, and n is the total number of observation vectors.


### 3.1. The Multiplicative EKF for Attitude Estimation

The most well-known nonlinear filter for spacecraft attitude estimation is the EKF. There are several different implementations of the attitude EKF, depending on both the attitude parameter used in the state vector and the form in which observations are input [3]. The best known and most widely used attitude EKF is MEKF. The attitude MEKF is derived from the following equations.Attitude predictions [5]:
(20)
q
^
k
/
k
-
1
=
Ω
-
k
/
k
-
1
(
ω
^
k
-
1
)
q
^
k
-
1
,
b
^
k
/
k
-
1
=
b
^
k
-
1
.Predicted measurement error:(21)
Δ
r
b
k
=
r
b
k
-
r
^
b
k
,
r
^
b
k
j
=
R
(
q
^
k
/
k
-
1
)
r
i
j
,
r
^
b
k
=
[
(
r
^
b
k
1
)
T
⋯
(
r
^
b
k
j
)
T
⋯
(
r
^
b
k
n
)
T
]
T
.Covariance matrix for the attitude predicted errors:(22)
P
k
/
k
-
1
=
Φ
k
/
k
-
1
P
k
-
1
Φ
k
/
k
-
1
T
+
Q
w
k
.Gain matrix:(23)
K
k
=
P
k
/
k
-
1
H
k
T
[
H
k
P
k
/
k
-
1
H
k
T
+
Q
v
k
]
-
1
.Covariance matrix for the attitude estimation errors:(24)
P
k
=
(
I
6
-
K
k
H
k
)
P
k
/
k
-
1
(
I
6
-
K
k
H
k
)
T
+
K
k
Q
v
k
K
k
T
.Local error attitude estimations:(25)
Δ
x
^
k
=
[
Δ
g
^
k
Δ
b
^
k
]
=
K
k
Δ
r
b
k
.Error quaternion estimation:(26)
Δ
q
^
k
=
1
(
1
+
Δ
g
^
k
T
Δ
g
^
k
)
[
Δ
g
^
k
1
]
.Attitude estimations:(27)
q
^
k
=
Δ
q
^
k
⊗
q
^
k
/
k
-
1
,
b
^
k
=
b
^
k
/
k
-
1
+
Δ
b
^
k
,
where
(28)
Ω
-
k
/
k
-
1
=
[
cos
⁡
(
ψ
k
-
1
)
I
3
-
[
ϕ
k
-
1
×
]
ϕ
k
-
1
-
ϕ
k
-
1
cos
⁡
(
ψ
k
-
1
)
]
,
Φ
k
/
k
-
1
=
I
6
-
[
[
ω
^
k
-
1
×
]
1
2
I
3
0
3
0
3
]
·
Δ
t
,
ϕ
k
-
1
=
sin
(
ψ
k
-
1
)
ω
^
k
-
1
∥
ω
^
k
-
1
∥
,
ψ
k
-
1
=
0.5
∥
ω
^
k
-
1
∥
Δ
t
,
where Q
w
k and Q
v
k are the covariance matrixes for the process noise and the vector observation noise, respectively, and Δ
t is the discretization step size.


### 3.2. The Computationally Efficient Attitude Estimator

The evaluation for the gain matrix in attitude MEKF requires to calculate the inverse of a3
n
×
3
n matrix, resulting in heavy computational burden when n is large, especially for small spacecraft with limited computational source. A computationally efficient attitude MEKF based on the reduced vector observation model (RMEKF) is developed to solve this problem [13, 17].According to Cayley transformation, the attitude matrixR can be mapped to a minimum-element attitude parameterization, expressed by the skew symmetric Rodrigues matrix G [18]:
(29)
R
=
(
I
3
-
G
)
(
I
3
+
G
)
-
1
=
(
I
3
+
G
)
-
1
(
I
3
-
G
)
,
and the inverse transformation is expressed as follows:
(30)
G
=
(
I
3
-
R
)
(
I
3
+
R
)
-
1
=
(
I
3
+
R
)
-
1
(
I
3
-
R
)
,
where G is the skew symmetric matrix generated from the RPs g: G
=
[
g
×
].Equation (6) can be rewritten as
(31)
r
b
j
=
R
(
q
)
r
i
j
+
v
j
=
δ
R
R
(
q
^
)
r
i
j
+
v
j
=
δ
R
r
^
b
j
+
v
j
.Substitute the second term of (29) into the above equation; one can obtain the following equation:
(32)
r
b
j
-
r
^
b
j
=
[
(
r
b
j
+
r
^
b
j
)
×
]
δ
g
+
v
1
j
,
where v
1
j
=
(
I
3
+
[
δ
g
×
]
)
v
j. Then, the observation model can be expressed as
(33)
δ
r
b
=
H
1
δ
g
+
v
1
,
where
(34)
H
1
=
[
h
-
1
⋮
h
-
j
⋮
h
-
n
]
,
v
1
=
[
v
-
1
⋮
v
-
j
⋮
v
-
n
]
,
where h
-
j
=
[
(
r
b
j
+
r
^
b
j
)
×
]
,
j
=
1
,
…
,
n.In order to reduce the computational cost of attitude MEKF, the dimension of the observation model equation (33) can be reduced to 3 by multiplying both sides of the equation by H
1
T [13, 17]. The weighted factor is firstly designed on the bias of the information conservation principle, to ensure the information for each vector observation without losing after the dimension of the observation model reduced to 3. The weighted factor is designed as
(35)
w
j
=
Q
v
t
(
Q
v
j
)
-
1
,
where Q
v
t
=
[
∑
j
=
1
n
‍
(
Q
v
j
)
-
1
]
-
1. Then, (33) can be rewritten in the following form:
(36)
δ
r
b
′
=
H
2
δ
g
+
v
2
,
where
(37)
δ
r
′
b
=
[
w
1
δ
r
b
1
⋮
w
j
δ
r
b
j
⋮
w
n
δ
r
b
n
]
,
H
2
=
[
w
1
h
-
1
⋮
w
j
h
-
j
⋮
w
n
h
-
n
]
,
v
2
=
[
w
1
v
1
1
⋮
w
j
v
1
j
⋮
w
n
v
1
n
]
.Multiply both sides of the above equation byH
2
T, one can get
(38)
δ
Z
=
H
-
δ
x
+
v
-
,
where
(39)
δ
Z
=
H
2
T
δ
r
b
′
,
H
-
=
[
H
-
1
0
3
]
,
H
-
1
=
H
2
T
H
2
,
v
-
=
H
2
T
v
2
.The local attitude estimation error system is composed of (16) and (38). It is obvious that the dimension of the vector observation model described by the above equation is 3. It only requires evaluating a 3
×
3 matrix inversion for the gain matrix K
k in the MEKF, rather than the 3
n
×
3
n matrix inversion. As a result, the computational burden is reduced.According to the MEKF, it is easy to get the following implementation of the RMEKF [13, 17]. The equations for attitude predictions, error attitude estimations, and estimations are the same in the RMEKF and MEKF. For the above reason, only the equations for local error attitude estimations in the RMEKF are shown here.Local error attitude estimations:(40)
Δ
x
^
k
=
K
k
Δ
Z
k
.Gain matrix:(41)
K
k
=
[
p
k
/
k
-
1
11
p
k
/
k
-
1
12
]
(
H
-
1
k
p
k
/
k
-
1
11
+
Q
v
t
)
-
1
.Covariance matrix(42)
P
k
=
(
I
6
-
K
k
H
-
k
)
P
k
/
k
-
1
(
I
6
-
K
k
H
-
k
)
T
+
K
k
H
-
1
k
Q
v
t
K
k
T
,
where p
k
/
k
-
1
11 and p
k
/
k
-
1
12 are the 3
×
3 submatrices of P
k
/
k
-
1, namely,
(43)
P
k
/
k
-
1
=
[
p
k
/
k
-
1
11
p
k
/
k
-
1
12
p
k
/
k
-
1
12
p
k
/
k
-
1
22
]
,
and the matrices Δ
Z
k and H
-
1
k are expressed as follows:
(44)
Δ
Z
k
=
∑
j
=
1
n
‍
T
k
j
×
S
k
j
,
H
-
1
k
=
t
r
(
∑
j
=
1
n
‍
S
k
j
(
S
k
j
)
T
)
I
3
-
∑
j
=
1
n
‍
S
k
j
(
S
k
j
)
T
,
S
k
j
=
w
j
(
r
b
k
j
+
r
^
b
k
j
)
,
T
k
j
=
w
j
(
r
b
k
j
-
r
^
b
k
j
)
.


### 3.3. The Improved Computationally Efficient Attitude Estimator

A 3-dimension vector observation model is used in attitude RMEKF instead of the original3
n-dimension model. Only the calculation of a 3
×
3 matrix inversion is required for the gain matrix, resulting in much less computational cost than that of attitude MEKF. However, the Jacobian matrixes need to update online in attitude RMEKF too, and the process of the matrix inversion remains in the calculation for the gain matrix. In order to further reduce the computational burden of attitude EKF, a new attitude estimation strategy is developed in this section.The local attitude error estimation system composed of (16) and (18) is rewritten as follows:
(45)
δ
x
˙
=
f
1
(
p
)
δ
x
+
B
n
n
,
δ
r
b
=
f
2
(
p
)
δ
x
+
D
n
n
,
where
(46)
f
1
(
p
)
=
-
[
[
ω
^
×
]
1
2
I
3
0
3
0
3
]
,
n
n
=
[
n
w
v
]
,
f
2
(
p
)
=
H
,
p
=
[
ω
^
T
q
^
T
]
T
,
B
=
[
I
6
0
6
×
3
n
]
,
D
=
[
0
3
n
×
6
I
3
n
]
.It is assumed that the process noisen
w and the measurement noise v are uncorrelated white noise. The parameter p is bounded by the 6-dimension space and the value set of p belongs to a compact set, since the angular velocity varies in a finite interval in most practical applications and each element of attitude quaternion takes value on the interval [
0,1
]. That is to say, p
∈
Ω
⊂
ℝ
6 and Ω is a compact set.Denote(47)
F
(
p
)
=
[
f
1
(
p
)
B
f
2
(
p
)
D
]
.Then, the above equation can be approximated as a convex combination of thel constant linear system matrixes F
i
,
i
=
1
,
…
,
l [19], namely,
(48)
F
(
p
)
≈
[
A
B
C
D
]
=
∑
i
=
1
l
‍
λ
i
(
p
)
[
A
i
B
C
i
D
]
,
where λ
i
(
p
) is the bias of the convex combination, which satisfies ∑
i
=
1
l
‍
λ
i
(
p
)
=
1, 0
≤
λ
i
(
p
)
≤
1. Thus, the local attitude estimation error system equation (45) can be described with the following form:
(49)
δ
x
˙
=
A
δ
x
+
B
n
n
,
δ
r
b
=
C
δ
x
+
D
n
n
,
where the system matrixes A, B, C, and D are denoted as (48).The equivalent discrete-time form of (49) can be approximated as
(50)
Δ
x
k
=
A
-
(
t
k
,
t
k
-
1
)
Δ
x
k
-
1
+
B
n
-
n
,
k
,
Δ
r
b
k
=
C
Δ
x
k
+
D
n
-
n
,
k
,
where n
-
n
,
k is the equivalent discrete-time noise
(51)
n
-
n
,
k
=
[
n
-
w
,
k
v
k
]
,
n
-
w
,
k
=
∫
t
k
-
1
t
k
‍
A
-
(
t
k
,
τ
)
n
w
(
τ
)
d
τ
.The state transition matrix can be approximated by(52)
A
-
(
t
k
,
t
k
-
1
)
=
I
6
+
A
·
Δ
t
=
I
6
+
(
∑
i
=
1
l
‍
λ
i
(
p
)
A
i
)
·
Δ
t
.Sinceλ
i
(
p
) satisfies ∑
i
=
1
l
‍
λ
i
(
p
)
=
1, 0
≤
λ
i
(
p
)
≤
1, the above equation can be easily reexpressed as
(53)
A
-
(
t
k
,
t
k
-
1
)
=
∑
i
=
1
l
‍
λ
i
(
p
)
(
I
6
+
A
i
·
Δ
t
)
.For convenience of notation, let(54)
A
-
=
A
-
(
t
k
,
t
k
-
1
)
,
A
-
i
=
I
6
+
A
i
;then (53) can be rewritten as
(55)
A
-
=
∑
i
=
1
l
‍
λ
i
(
p
)
A
-
i
.It is now a straightforward matter to show that the discrete-time form of the local attitude error estimation system can be described by an uncertain discrete linear polytopic model in the form of (50). With this model, the local error attitude estimation problem is converted to a robust linear one, that is to find a stable local error attitude estimator in the form
(56)
Δ
x
^
k
=
C
F
Δ
x
^
m
k
+
D
F
Δ
r
b
k
,
Δ
x
^
m
k
+
1
=
A
F
Δ
x
^
m
k
+
B
F
Δ
r
b
k
,such that the local attitude estimation error variance,E
[
(
Δ
x
k
-
Δ
x
^
k
)
T
(
Δ
x
k
-
Δ
x
^
k
)
], is minimized, where (
A
F
,
B
F
,
C
F
,
D
F
) are constant matrixes to be determined. According to the robust H
2 filtering algorithm given by Duan et al. [19], the matrixes (
A
F
,
B
F
,
C
F
,
D
F
) can be obtained by solving an optimization given in the following Lemma.Lemma 1 (see [19]).
Consider the system (50); a filter of the form (56) that achieves a suboptimal guaranteed filtering error covariance bound can be derived from the following optimization:(57)
min
⁡
G
11
,
G
21
,
G
2
,
F
11
,
F
21
,
S
A
,
S
B
,
S
C
,
S
D
,
P
11
i
,
P
12
i
,
P
22
i
t
r
a
c
e
(
P
x
^
)
s
.
t
.
[
G
11
+
G
11
T
-
P
11
i
G
2
+
G
21
T
-
P
12
i
ψ
1
i
*
G
2
+
G
2
T
-
P
22
i
ψ
2
i
*
*
ψ
3
i
*
*
*
*
*
*
S
A
-
F
21
T
G
11
B
+
S
B
D
S
A
-
α
2
G
2
T
G
21
B
+
S
B
D
ψ
4
i
-
F
11
B
-
α
1
S
B
D
P
22
i
-
α
2
S
A
-
α
2
S
A
T
-
F
21
B
-
α
2
S
B
D
*
I
3
n
+
6
]
>
0
,
[
P
x
^
I
6
-
S
D
C
i
-
S
C
-
S
D
D
*
P
11
i
P
12
i
0
6
×
(
3
n
+
6
)
*
*
P
22
i
0
6
×
(
3
n
+
6
)
*
*
*
I
3
n
+
6
]
>
0
,
i
=
1,2
,
…
,
l
,where α
1 and α
2 are fixed parameters, and
(58)
ψ
1
i
=
G
11
A
-
i
+
S
B
C
i
-
F
11
T
,
ψ
2
i
=
G
21
A
-
i
+
S
B
C
i
-
α
1
G
2
T
,
ψ
3
i
=
P
11
i
-
F
11
A
-
i
-
α
1
S
B
C
i
-
A
-
i
T
F
11
T
-
α
1
C
i
T
S
B
T
,
ψ
4
i
=
P
12
i
-
α
1
S
A
-
A
-
i
T
F
21
T
-
α
2
C
i
T
S
B
T
.
The suboptimal filter is given by(59)
A
F
=
G
2
-
1
S
A
,
B
F
=
G
2
-
1
S
B
,
C
F
=
S
C
,
D
F
=
S
D
.As a consequence, the improved computationally efficient attitude estimator can be presented as shown in Table1.Table 1
The improved computationally efficient attitude estimator.
Initialization
Step  1. Determine the vertexes of the polytopic linear model described as (50).
Step  2. Search for optimal solutions (
G
2
,
S
A
,
S
B
,
S
C
,
S
D
) of LMIs (57); then calculate the constant filter coefficients (
A
F
,
B
F
,
C
F
,
D
F
) for the local attitude estimator by (59).
Estimator (one cycle)
Givenq
^
k
-
1
,
b
^
k
-
1
,
ω
^
k
-
1 and measurements u
k
,
r
b
k, one has the following.
Step  1. Compute the attitude predictions q
^
k
/
k
-
1 and b
^
k
/
k
-
1 by (20).
Step  2. Compute measurement prediction errors Δ
r
b
k by (21).
Step  3. Compute the local error attitude estimation Δ
x
^
k by (50).
Step  4. Update the attitude estimations q
^
k and b
^
k by (26) and (27).
Step  5. Update the angular velocity estimation: ω
^
k
=
u
k
-
b
^
k.The filter coefficients(
A
F
,
B
F
,
C
F
,
D
F
) for the local error attitude estimator can be solved before the recursive process of the attitude estimation. Since they are constant, the matrix inversion is not required for the gain matrix and the updating of the Jacobian matrixes online neither. As a result, the computational complexity of the proposed efficient attitude estimator (PEF) reduces extensively. As can be seen in the above attitude estimation strategy, the key to implement the PEF is to determine the vertexes of the polytopic linear model (49). There are two kinds of methods: parameter boundary method and TP model transformation method [20]. The polytopic linear model for an affine parameter system obtained by the former is quite accurate, but, for other systems, the polytopic linear model obtained by the latter is with less conservativeness. The TP model transformation algorithm is described in Algorithm 1.Algorithm 1:TP model transformation algorithm [20].
Step  1. Determine the bounded domain Ω where the parameter p varying on;
Step  2. Set the transformation space Ω into N uniformly distributed grid: Ω
=
[
a
1
,
b
1
]
×
⋯
×
[
a
N
,
b
N
];
Step  3. Sample the functions F
(
p
) over the hyper rectangular grid and store the sample
matrixes in the tensorF
-;
Step  4. Execute the higher order singular value decomposition (HOSVD) on tensor F
- and
extract the minimal basis, the result of this step isF
-
≈
V
⊗
m
U
m;
Step  5. Normalize the basis matrix U
m, one can obtain F
(
p
)
≈
V
-
⊗
m
U
-
m;
Step  6. Extract the vertexes of the polytopic linear model (49) from V
-.


## 3.1. The Multiplicative EKF for Attitude Estimation

The most well-known nonlinear filter for spacecraft attitude estimation is the EKF. There are several different implementations of the attitude EKF, depending on both the attitude parameter used in the state vector and the form in which observations are input [3]. The best known and most widely used attitude EKF is MEKF. The attitude MEKF is derived from the following equations.Attitude predictions [5]:
(20)
q
^
k
/
k
-
1
=
Ω
-
k
/
k
-
1
(
ω
^
k
-
1
)
q
^
k
-
1
,
b
^
k
/
k
-
1
=
b
^
k
-
1
.Predicted measurement error:(21)
Δ
r
b
k
=
r
b
k
-
r
^
b
k
,
r
^
b
k
j
=
R
(
q
^
k
/
k
-
1
)
r
i
j
,
r
^
b
k
=
[
(
r
^
b
k
1
)
T
⋯
(
r
^
b
k
j
)
T
⋯
(
r
^
b
k
n
)
T
]
T
.Covariance matrix for the attitude predicted errors:(22)
P
k
/
k
-
1
=
Φ
k
/
k
-
1
P
k
-
1
Φ
k
/
k
-
1
T
+
Q
w
k
.Gain matrix:(23)
K
k
=
P
k
/
k
-
1
H
k
T
[
H
k
P
k
/
k
-
1
H
k
T
+
Q
v
k
]
-
1
.Covariance matrix for the attitude estimation errors:(24)
P
k
=
(
I
6
-
K
k
H
k
)
P
k
/
k
-
1
(
I
6
-
K
k
H
k
)
T
+
K
k
Q
v
k
K
k
T
.Local error attitude estimations:(25)
Δ
x
^
k
=
[
Δ
g
^
k
Δ
b
^
k
]
=
K
k
Δ
r
b
k
.Error quaternion estimation:(26)
Δ
q
^
k
=
1
(
1
+
Δ
g
^
k
T
Δ
g
^
k
)
[
Δ
g
^
k
1
]
.Attitude estimations:(27)
q
^
k
=
Δ
q
^
k
⊗
q
^
k
/
k
-
1
,
b
^
k
=
b
^
k
/
k
-
1
+
Δ
b
^
k
,
where
(28)
Ω
-
k
/
k
-
1
=
[
cos
⁡
(
ψ
k
-
1
)
I
3
-
[
ϕ
k
-
1
×
]
ϕ
k
-
1
-
ϕ
k
-
1
cos
⁡
(
ψ
k
-
1
)
]
,
Φ
k
/
k
-
1
=
I
6
-
[
[
ω
^
k
-
1
×
]
1
2
I
3
0
3
0
3
]
·
Δ
t
,
ϕ
k
-
1
=
sin
(
ψ
k
-
1
)
ω
^
k
-
1
∥
ω
^
k
-
1
∥
,
ψ
k
-
1
=
0.5
∥
ω
^
k
-
1
∥
Δ
t
,
where Q
w
k and Q
v
k are the covariance matrixes for the process noise and the vector observation noise, respectively, and Δ
t is the discretization step size.


## 3.2. The Computationally Efficient Attitude Estimator

The evaluation for the gain matrix in attitude MEKF requires to calculate the inverse of a3
n
×
3
n matrix, resulting in heavy computational burden when n is large, especially for small spacecraft with limited computational source. A computationally efficient attitude MEKF based on the reduced vector observation model (RMEKF) is developed to solve this problem [13, 17].According to Cayley transformation, the attitude matrixR can be mapped to a minimum-element attitude parameterization, expressed by the skew symmetric Rodrigues matrix G [18]:
(29)
R
=
(
I
3
-
G
)
(
I
3
+
G
)
-
1
=
(
I
3
+
G
)
-
1
(
I
3
-
G
)
,
and the inverse transformation is expressed as follows:
(30)
G
=
(
I
3
-
R
)
(
I
3
+
R
)
-
1
=
(
I
3
+
R
)
-
1
(
I
3
-
R
)
,
where G is the skew symmetric matrix generated from the RPs g: G
=
[
g
×
].Equation (6) can be rewritten as
(31)
r
b
j
=
R
(
q
)
r
i
j
+
v
j
=
δ
R
R
(
q
^
)
r
i
j
+
v
j
=
δ
R
r
^
b
j
+
v
j
.Substitute the second term of (29) into the above equation; one can obtain the following equation:
(32)
r
b
j
-
r
^
b
j
=
[
(
r
b
j
+
r
^
b
j
)
×
]
δ
g
+
v
1
j
,
where v
1
j
=
(
I
3
+
[
δ
g
×
]
)
v
j. Then, the observation model can be expressed as
(33)
δ
r
b
=
H
1
δ
g
+
v
1
,
where
(34)
H
1
=
[
h
-
1
⋮
h
-
j
⋮
h
-
n
]
,
v
1
=
[
v
-
1
⋮
v
-
j
⋮
v
-
n
]
,
where h
-
j
=
[
(
r
b
j
+
r
^
b
j
)
×
]
,
j
=
1
,
…
,
n.In order to reduce the computational cost of attitude MEKF, the dimension of the observation model equation (33) can be reduced to 3 by multiplying both sides of the equation by H
1
T [13, 17]. The weighted factor is firstly designed on the bias of the information conservation principle, to ensure the information for each vector observation without losing after the dimension of the observation model reduced to 3. The weighted factor is designed as
(35)
w
j
=
Q
v
t
(
Q
v
j
)
-
1
,
where Q
v
t
=
[
∑
j
=
1
n
‍
(
Q
v
j
)
-
1
]
-
1. Then, (33) can be rewritten in the following form:
(36)
δ
r
b
′
=
H
2
δ
g
+
v
2
,
where
(37)
δ
r
′
b
=
[
w
1
δ
r
b
1
⋮
w
j
δ
r
b
j
⋮
w
n
δ
r
b
n
]
,
H
2
=
[
w
1
h
-
1
⋮
w
j
h
-
j
⋮
w
n
h
-
n
]
,
v
2
=
[
w
1
v
1
1
⋮
w
j
v
1
j
⋮
w
n
v
1
n
]
.Multiply both sides of the above equation byH
2
T, one can get
(38)
δ
Z
=
H
-
δ
x
+
v
-
,
where
(39)
δ
Z
=
H
2
T
δ
r
b
′
,
H
-
=
[
H
-
1
0
3
]
,
H
-
1
=
H
2
T
H
2
,
v
-
=
H
2
T
v
2
.The local attitude estimation error system is composed of (16) and (38). It is obvious that the dimension of the vector observation model described by the above equation is 3. It only requires evaluating a 3
×
3 matrix inversion for the gain matrix K
k in the MEKF, rather than the 3
n
×
3
n matrix inversion. As a result, the computational burden is reduced.According to the MEKF, it is easy to get the following implementation of the RMEKF [13, 17]. The equations for attitude predictions, error attitude estimations, and estimations are the same in the RMEKF and MEKF. For the above reason, only the equations for local error attitude estimations in the RMEKF are shown here.Local error attitude estimations:(40)
Δ
x
^
k
=
K
k
Δ
Z
k
.Gain matrix:(41)
K
k
=
[
p
k
/
k
-
1
11
p
k
/
k
-
1
12
]
(
H
-
1
k
p
k
/
k
-
1
11
+
Q
v
t
)
-
1
.Covariance matrix(42)
P
k
=
(
I
6
-
K
k
H
-
k
)
P
k
/
k
-
1
(
I
6
-
K
k
H
-
k
)
T
+
K
k
H
-
1
k
Q
v
t
K
k
T
,
where p
k
/
k
-
1
11 and p
k
/
k
-
1
12 are the 3
×
3 submatrices of P
k
/
k
-
1, namely,
(43)
P
k
/
k
-
1
=
[
p
k
/
k
-
1
11
p
k
/
k
-
1
12
p
k
/
k
-
1
12
p
k
/
k
-
1
22
]
,
and the matrices Δ
Z
k and H
-
1
k are expressed as follows:
(44)
Δ
Z
k
=
∑
j
=
1
n
‍
T
k
j
×
S
k
j
,
H
-
1
k
=
t
r
(
∑
j
=
1
n
‍
S
k
j
(
S
k
j
)
T
)
I
3
-
∑
j
=
1
n
‍
S
k
j
(
S
k
j
)
T
,
S
k
j
=
w
j
(
r
b
k
j
+
r
^
b
k
j
)
,
T
k
j
=
w
j
(
r
b
k
j
-
r
^
b
k
j
)
.


## 3.3. The Improved Computationally Efficient Attitude Estimator

A 3-dimension vector observation model is used in attitude RMEKF instead of the original3
n-dimension model. Only the calculation of a 3
×
3 matrix inversion is required for the gain matrix, resulting in much less computational cost than that of attitude MEKF. However, the Jacobian matrixes need to update online in attitude RMEKF too, and the process of the matrix inversion remains in the calculation for the gain matrix. In order to further reduce the computational burden of attitude EKF, a new attitude estimation strategy is developed in this section.The local attitude error estimation system composed of (16) and (18) is rewritten as follows:
(45)
δ
x
˙
=
f
1
(
p
)
δ
x
+
B
n
n
,
δ
r
b
=
f
2
(
p
)
δ
x
+
D
n
n
,
where
(46)
f
1
(
p
)
=
-
[
[
ω
^
×
]
1
2
I
3
0
3
0
3
]
,
n
n
=
[
n
w
v
]
,
f
2
(
p
)
=
H
,
p
=
[
ω
^
T
q
^
T
]
T
,
B
=
[
I
6
0
6
×
3
n
]
,
D
=
[
0
3
n
×
6
I
3
n
]
.It is assumed that the process noisen
w and the measurement noise v are uncorrelated white noise. The parameter p is bounded by the 6-dimension space and the value set of p belongs to a compact set, since the angular velocity varies in a finite interval in most practical applications and each element of attitude quaternion takes value on the interval [
0,1
]. That is to say, p
∈
Ω
⊂
ℝ
6 and Ω is a compact set.Denote(47)
F
(
p
)
=
[
f
1
(
p
)
B
f
2
(
p
)
D
]
.Then, the above equation can be approximated as a convex combination of thel constant linear system matrixes F
i
,
i
=
1
,
…
,
l [19], namely,
(48)
F
(
p
)
≈
[
A
B
C
D
]
=
∑
i
=
1
l
‍
λ
i
(
p
)
[
A
i
B
C
i
D
]
,
where λ
i
(
p
) is the bias of the convex combination, which satisfies ∑
i
=
1
l
‍
λ
i
(
p
)
=
1, 0
≤
λ
i
(
p
)
≤
1. Thus, the local attitude estimation error system equation (45) can be described with the following form:
(49)
δ
x
˙
=
A
δ
x
+
B
n
n
,
δ
r
b
=
C
δ
x
+
D
n
n
,
where the system matrixes A, B, C, and D are denoted as (48).The equivalent discrete-time form of (49) can be approximated as
(50)
Δ
x
k
=
A
-
(
t
k
,
t
k
-
1
)
Δ
x
k
-
1
+
B
n
-
n
,
k
,
Δ
r
b
k
=
C
Δ
x
k
+
D
n
-
n
,
k
,
where n
-
n
,
k is the equivalent discrete-time noise
(51)
n
-
n
,
k
=
[
n
-
w
,
k
v
k
]
,
n
-
w
,
k
=
∫
t
k
-
1
t
k
‍
A
-
(
t
k
,
τ
)
n
w
(
τ
)
d
τ
.The state transition matrix can be approximated by(52)
A
-
(
t
k
,
t
k
-
1
)
=
I
6
+
A
·
Δ
t
=
I
6
+
(
∑
i
=
1
l
‍
λ
i
(
p
)
A
i
)
·
Δ
t
.Sinceλ
i
(
p
) satisfies ∑
i
=
1
l
‍
λ
i
(
p
)
=
1, 0
≤
λ
i
(
p
)
≤
1, the above equation can be easily reexpressed as
(53)
A
-
(
t
k
,
t
k
-
1
)
=
∑
i
=
1
l
‍
λ
i
(
p
)
(
I
6
+
A
i
·
Δ
t
)
.For convenience of notation, let(54)
A
-
=
A
-
(
t
k
,
t
k
-
1
)
,
A
-
i
=
I
6
+
A
i
;then (53) can be rewritten as
(55)
A
-
=
∑
i
=
1
l
‍
λ
i
(
p
)
A
-
i
.It is now a straightforward matter to show that the discrete-time form of the local attitude error estimation system can be described by an uncertain discrete linear polytopic model in the form of (50). With this model, the local error attitude estimation problem is converted to a robust linear one, that is to find a stable local error attitude estimator in the form
(56)
Δ
x
^
k
=
C
F
Δ
x
^
m
k
+
D
F
Δ
r
b
k
,
Δ
x
^
m
k
+
1
=
A
F
Δ
x
^
m
k
+
B
F
Δ
r
b
k
,such that the local attitude estimation error variance,E
[
(
Δ
x
k
-
Δ
x
^
k
)
T
(
Δ
x
k
-
Δ
x
^
k
)
], is minimized, where (
A
F
,
B
F
,
C
F
,
D
F
) are constant matrixes to be determined. According to the robust H
2 filtering algorithm given by Duan et al. [19], the matrixes (
A
F
,
B
F
,
C
F
,
D
F
) can be obtained by solving an optimization given in the following Lemma.Lemma 1 (see [19]).
Consider the system (50); a filter of the form (56) that achieves a suboptimal guaranteed filtering error covariance bound can be derived from the following optimization:(57)
min
⁡
G
11
,
G
21
,
G
2
,
F
11
,
F
21
,
S
A
,
S
B
,
S
C
,
S
D
,
P
11
i
,
P
12
i
,
P
22
i
t
r
a
c
e
(
P
x
^
)
s
.
t
.
[
G
11
+
G
11
T
-
P
11
i
G
2
+
G
21
T
-
P
12
i
ψ
1
i
*
G
2
+
G
2
T
-
P
22
i
ψ
2
i
*
*
ψ
3
i
*
*
*
*
*
*
S
A
-
F
21
T
G
11
B
+
S
B
D
S
A
-
α
2
G
2
T
G
21
B
+
S
B
D
ψ
4
i
-
F
11
B
-
α
1
S
B
D
P
22
i
-
α
2
S
A
-
α
2
S
A
T
-
F
21
B
-
α
2
S
B
D
*
I
3
n
+
6
]
>
0
,
[
P
x
^
I
6
-
S
D
C
i
-
S
C
-
S
D
D
*
P
11
i
P
12
i
0
6
×
(
3
n
+
6
)
*
*
P
22
i
0
6
×
(
3
n
+
6
)
*
*
*
I
3
n
+
6
]
>
0
,
i
=
1,2
,
…
,
l
,where α
1 and α
2 are fixed parameters, and
(58)
ψ
1
i
=
G
11
A
-
i
+
S
B
C
i
-
F
11
T
,
ψ
2
i
=
G
21
A
-
i
+
S
B
C
i
-
α
1
G
2
T
,
ψ
3
i
=
P
11
i
-
F
11
A
-
i
-
α
1
S
B
C
i
-
A
-
i
T
F
11
T
-
α
1
C
i
T
S
B
T
,
ψ
4
i
=
P
12
i
-
α
1
S
A
-
A
-
i
T
F
21
T
-
α
2
C
i
T
S
B
T
.
The suboptimal filter is given by(59)
A
F
=
G
2
-
1
S
A
,
B
F
=
G
2
-
1
S
B
,
C
F
=
S
C
,
D
F
=
S
D
.As a consequence, the improved computationally efficient attitude estimator can be presented as shown in Table1.Table 1
The improved computationally efficient attitude estimator.
Initialization
Step  1. Determine the vertexes of the polytopic linear model described as (50).
Step  2. Search for optimal solutions (
G
2
,
S
A
,
S
B
,
S
C
,
S
D
) of LMIs (57); then calculate the constant filter coefficients (
A
F
,
B
F
,
C
F
,
D
F
) for the local attitude estimator by (59).
Estimator (one cycle)
Givenq
^
k
-
1
,
b
^
k
-
1
,
ω
^
k
-
1 and measurements u
k
,
r
b
k, one has the following.
Step  1. Compute the attitude predictions q
^
k
/
k
-
1 and b
^
k
/
k
-
1 by (20).
Step  2. Compute measurement prediction errors Δ
r
b
k by (21).
Step  3. Compute the local error attitude estimation Δ
x
^
k by (50).
Step  4. Update the attitude estimations q
^
k and b
^
k by (26) and (27).
Step  5. Update the angular velocity estimation: ω
^
k
=
u
k
-
b
^
k.The filter coefficients(
A
F
,
B
F
,
C
F
,
D
F
) for the local error attitude estimator can be solved before the recursive process of the attitude estimation. Since they are constant, the matrix inversion is not required for the gain matrix and the updating of the Jacobian matrixes online neither. As a result, the computational complexity of the proposed efficient attitude estimator (PEF) reduces extensively. As can be seen in the above attitude estimation strategy, the key to implement the PEF is to determine the vertexes of the polytopic linear model (49). There are two kinds of methods: parameter boundary method and TP model transformation method [20]. The polytopic linear model for an affine parameter system obtained by the former is quite accurate, but, for other systems, the polytopic linear model obtained by the latter is with less conservativeness. The TP model transformation algorithm is described in Algorithm 1.Algorithm 1:TP model transformation algorithm [20].
Step  1. Determine the bounded domain Ω where the parameter p varying on;
Step  2. Set the transformation space Ω into N uniformly distributed grid: Ω
=
[
a
1
,
b
1
]
×
⋯
×
[
a
N
,
b
N
];
Step  3. Sample the functions F
(
p
) over the hyper rectangular grid and store the sample
matrixes in the tensorF
-;
Step  4. Execute the higher order singular value decomposition (HOSVD) on tensor F
- and
extract the minimal basis, the result of this step isF
-
≈
V
⊗
m
U
m;
Step  5. Normalize the basis matrix U
m, one can obtain F
(
p
)
≈
V
-
⊗
m
U
-
m;
Step  6. Extract the vertexes of the polytopic linear model (49) from V
-.


## 4. Simulation Results

In this section, the comparisons of the computational cost and the accuracy between the PEF and other attitude estimators are given.


### 4.1. Computational Complexity

The computer complexity of the attitudes PEF, RMEKF, and MEKF is shown in Table2. Only the computational cost of the local error attitude estimation is given since the equations of the attitude predictions and updatings are the same in the three attitude filters.Table 2
The statistics of the computational cost for three filters.
Operation
PEF
RMEKF
MEKF
Addition
39
n
+
60
16
n
+
726
45
n
3
+
148.5
n
2
+
380.5
n
+
750
Multiplicationa
36
n
+
72
24
n
+
922
45
n
3
+
193.5
n
2
+
449.5
n
+
864
aThe matrix inversions in MEKF and RMEKF are calculated by QR decomposition method.According to the statistics, it is obvious to find that the computational cost of the MEKF is much more than that of the other two filters. The additional cost of the PEF is less than that of the RMEKF ifn
<
28, while the multiplicational cost of the RMEKF is much more than that of the PEF if n
≤
70. It happens that in a realistic situation the observation vector number is less than 28; therefore, one can get the conclusion that the PEF has the least computational cost in the three attitude filters and the MEKF has the most.Denote the computational cost error between the PEF and other attitude filters as(60)
Δ
C
=
C
F
-
C
PEF
,
where C
PEF and C
F are the computational cost of the multiplication or addition for the PEF and other filters (MEKF or RMEKF), respectively.The computational cost error for the MEKF and PEF is shown in Figure1. The increasing rates of the multiplicational cost error and the additional cost error are nearly the same for the two filters. They both increase rapidly as the observation vector number increases, indicating that the computational complexity of the PEF is much less than the MEKF, especially when the observation vector number is large.Figure 1
Error of computational cost.The computational cost errors for the RMEKF and PEF are shown in Figure2. The advantage of less computation complexity is not evident as the observation vector number increases in the PEF, while the decrescent rate of the multiplicational cost error for the two filters reduces more slowly than that of the additional cost error. Since the multiplication is much more complex than the addition, it can be concluded that the computational burden of the PEF reduces much more than the other two filters for the application of spacecraft attitude estimation in the practical engineering.Figure 2
Error of computational cost.


### 4.2. Accuracy Comparison

The initial angular velocity vector and the 3-1-2 Euler angles are given by[
0.05
-
0.1
0.05
]
T deg/s and [
45
60
32
]
T deg, respectively. The initial gyro drift rate is 5 deg/h and the standard covariance of the driven noise is 0.003 deg/h. The standard covariances of the gyro measurement noise and the vector observation noise are 0.5 deg/h and 5 arcsec, respectively. The initial covariance of attitude estimation error is 1
0
-
8
I
6 and the discretization step size is 200 ms. The initial attitude estimations are set as their true values, avoiding the poor performance or divergence of the three estimators caused by lacking a good a priori estimate of attitude.The diagrammatic representation of the attitude estimation system (attitude PEF or the other two estimators) designed in the simulation is shown in Figure3; the attitude estimations are given by the attitude estimator (attitude PEF or the other attitude estimators) based on the gyro and the vector observations. The observation vectors can be given by unit sun, star, and Earth’s magnetic field vectors. Two unit star vectors are chosen in the simulation, because several star vector observations can be obtained from the star tracker at a time. The Monte Carlo simulation is computed over an ensemble of 100 independent runs. In order to make the simulation results more intuitive, the attitude estimation errors are shown by the error attitude 3-1-2 Euler angles instead of quaternion, because the quaternion does not have intuitive physical meanings.Figure 3
The diagrammatic representation of the attitude estimation system.The average estimation errors of the attitude and the gyro drift rate bias with their respective3
σ bounds in the PEF are shown in Figure 4. As can be seen in the figure, the estimation errors of the attitude and gyro drift rate bias in the PEF all converge to within their respective 3
σ, indicating that the PEF performs well for the attitude estimation.Estimation errors with3
σ bounds in the attitude PEF.
(a)
Estimation errors of attitude with3
σ bounds
(b)
Estimation errors of angular velocities with3
σ boundsThe attitude principal rotation angle is used to scale the attitude estimation error, which is expressed in the form of quaternion as follows:(61)
Δ
θ
=
2
arccos
(
Δ
q
4
)
.The average initial attitude estimation errors of the three filters are shown in Figure5. The maximal estimation error of the attitude angle is the smallest in the PEF during the transient process of the filters. The initial attitude estimation errors in the RMEKF and MEKF are nearly identical.Figure 5
Initial attitude angle estimation errors.The steady attitude estimation errors of attitudes PEF and MEKF are given in Figure6. The steady attitude estimation error in attitude RMEKF is not shown here, because it is nearly the same as that in attitude MEKF. From the figure, one can get that the steady attitude error angle in attitude PEF varies around 4.13 arcsec, while it varies around 1.41 arcsec in attitude MEKF. The magnitude order of the accuracy achieved by attitude PEF is 1 arcsec and it is identical to that of attitude MEKF.Figure 6
Steady attitude angle estimation errors.The maximum values of the absolute steady estimation errors for the triaxial angular velocities and gyro drift rate bias are shown in Table3. The magnitude order of the steady estimation errors for the triaxial angular velocities and gyro drift rate bias are 1
0
-
4 deg/s and 1
0
-
6 deg/s in attitude PEF, respectively. They are nearly identical to those of attitude MEKF. Therefore, one can conclude that the proposed attitude PEF performs well for the spacecraft attitude estimation. Since its computational cost is much lower than that of the MEKF and RMEKF, it should be a good choice for the spacecraft attitude estimation application with limited computing resources and low accuracy demand.Table 3
The maximum absolute steady estimation errors.
Δ
ω  (
1
0
-
4
deg/s
)
Δ
b  (
1
0
-
6
deg/s
)
Δ
ω
x
Δ
ω
y
Δ
ω
z
Δ
b
x
Δ
b
y
Δ
b
z
PEF
1.904
1.971
1.789
3.826
4.458
2.858
MEKF
1.897
1.963
1.852
2.320
2.789
2.179


## 4.1. Computational Complexity

The computer complexity of the attitudes PEF, RMEKF, and MEKF is shown in Table2. Only the computational cost of the local error attitude estimation is given since the equations of the attitude predictions and updatings are the same in the three attitude filters.Table 2
The statistics of the computational cost for three filters.
Operation
PEF
RMEKF
MEKF
Addition
39
n
+
60
16
n
+
726
45
n
3
+
148.5
n
2
+
380.5
n
+
750
Multiplicationa
36
n
+
72
24
n
+
922
45
n
3
+
193.5
n
2
+
449.5
n
+
864
aThe matrix inversions in MEKF and RMEKF are calculated by QR decomposition method.According to the statistics, it is obvious to find that the computational cost of the MEKF is much more than that of the other two filters. The additional cost of the PEF is less than that of the RMEKF ifn
<
28, while the multiplicational cost of the RMEKF is much more than that of the PEF if n
≤
70. It happens that in a realistic situation the observation vector number is less than 28; therefore, one can get the conclusion that the PEF has the least computational cost in the three attitude filters and the MEKF has the most.Denote the computational cost error between the PEF and other attitude filters as(60)
Δ
C
=
C
F
-
C
PEF
,
where C
PEF and C
F are the computational cost of the multiplication or addition for the PEF and other filters (MEKF or RMEKF), respectively.The computational cost error for the MEKF and PEF is shown in Figure1. The increasing rates of the multiplicational cost error and the additional cost error are nearly the same for the two filters. They both increase rapidly as the observation vector number increases, indicating that the computational complexity of the PEF is much less than the MEKF, especially when the observation vector number is large.Figure 1
Error of computational cost.The computational cost errors for the RMEKF and PEF are shown in Figure2. The advantage of less computation complexity is not evident as the observation vector number increases in the PEF, while the decrescent rate of the multiplicational cost error for the two filters reduces more slowly than that of the additional cost error. Since the multiplication is much more complex than the addition, it can be concluded that the computational burden of the PEF reduces much more than the other two filters for the application of spacecraft attitude estimation in the practical engineering.Figure 2
Error of computational cost.


## 4.2. Accuracy Comparison

The initial angular velocity vector and the 3-1-2 Euler angles are given by[
0.05
-
0.1
0.05
]
T deg/s and [
45
60
32
]
T deg, respectively. The initial gyro drift rate is 5 deg/h and the standard covariance of the driven noise is 0.003 deg/h. The standard covariances of the gyro measurement noise and the vector observation noise are 0.5 deg/h and 5 arcsec, respectively. The initial covariance of attitude estimation error is 1
0
-
8
I
6 and the discretization step size is 200 ms. The initial attitude estimations are set as their true values, avoiding the poor performance or divergence of the three estimators caused by lacking a good a priori estimate of attitude.The diagrammatic representation of the attitude estimation system (attitude PEF or the other two estimators) designed in the simulation is shown in Figure3; the attitude estimations are given by the attitude estimator (attitude PEF or the other attitude estimators) based on the gyro and the vector observations. The observation vectors can be given by unit sun, star, and Earth’s magnetic field vectors. Two unit star vectors are chosen in the simulation, because several star vector observations can be obtained from the star tracker at a time. The Monte Carlo simulation is computed over an ensemble of 100 independent runs. In order to make the simulation results more intuitive, the attitude estimation errors are shown by the error attitude 3-1-2 Euler angles instead of quaternion, because the quaternion does not have intuitive physical meanings.Figure 3
The diagrammatic representation of the attitude estimation system.The average estimation errors of the attitude and the gyro drift rate bias with their respective3
σ bounds in the PEF are shown in Figure 4. As can be seen in the figure, the estimation errors of the attitude and gyro drift rate bias in the PEF all converge to within their respective 3
σ, indicating that the PEF performs well for the attitude estimation.Estimation errors with3
σ bounds in the attitude PEF.
(a)
Estimation errors of attitude with3
σ bounds
(b)
Estimation errors of angular velocities with3
σ boundsThe attitude principal rotation angle is used to scale the attitude estimation error, which is expressed in the form of quaternion as follows:(61)
Δ
θ
=
2
arccos
(
Δ
q
4
)
.The average initial attitude estimation errors of the three filters are shown in Figure5. The maximal estimation error of the attitude angle is the smallest in the PEF during the transient process of the filters. The initial attitude estimation errors in the RMEKF and MEKF are nearly identical.Figure 5
Initial attitude angle estimation errors.The steady attitude estimation errors of attitudes PEF and MEKF are given in Figure6. The steady attitude estimation error in attitude RMEKF is not shown here, because it is nearly the same as that in attitude MEKF. From the figure, one can get that the steady attitude error angle in attitude PEF varies around 4.13 arcsec, while it varies around 1.41 arcsec in attitude MEKF. The magnitude order of the accuracy achieved by attitude PEF is 1 arcsec and it is identical to that of attitude MEKF.Figure 6
Steady attitude angle estimation errors.The maximum values of the absolute steady estimation errors for the triaxial angular velocities and gyro drift rate bias are shown in Table3. The magnitude order of the steady estimation errors for the triaxial angular velocities and gyro drift rate bias are 1
0
-
4 deg/s and 1
0
-
6 deg/s in attitude PEF, respectively. They are nearly identical to those of attitude MEKF. Therefore, one can conclude that the proposed attitude PEF performs well for the spacecraft attitude estimation. Since its computational cost is much lower than that of the MEKF and RMEKF, it should be a good choice for the spacecraft attitude estimation application with limited computing resources and low accuracy demand.Table 3
The maximum absolute steady estimation errors.
Δ
ω  (
1
0
-
4
deg/s
)
Δ
b  (
1
0
-
6
deg/s
)
Δ
ω
x
Δ
ω
y
Δ
ω
z
Δ
b
x
Δ
b
y
Δ
b
z
PEF
1.904
1.971
1.789
3.826
4.458
2.858
MEKF
1.897
1.963
1.852
2.320
2.789
2.179


## 5. Conclusion

This paper develops a computationally efficient attitude estimation method for the spacecraft based on the gyro measurements and the vector observations. The global attitude is described by the quaternion, while the local attitude error is represented by the Rodrigues parameter to ensure that the quaternion satisfies the normalization constraint. Rather than reducing the computational complexity of attitude estimator by reducing the number of the sigma points or the measurement model dimension, the local error attitude estimator is designed with constant coefficients in the proposed attitude estimator. It does not need to calculate the matrix inversion for gain matrix or update the Jacobian matrixes online during the recursive process of the attitude estimator. As a result, the computational cost reduces extensively and the convergence speed for the attitude is much faster than other attitude filters. Simulation results show that the proposed attitude estimator improves the computational efficiency, though the accuracy for the attitude is a little larger than that of MEKF. It should be preferred over the MEKF in the practical spacecraft attitude estimation application with limited computing resources and low accuracy demand.



---

*Source: 540235-2014-03-12.xml*