1. INTRODUCTION
Attitude determination of spacecraft in low earth
orbits using magnetic sensors is attractive as a backup
for attitude determination system using sun sensors,
earth sensors or star sensors not only because magnetic
sensors are inexpensive, small and lightweight but also
they are robust and complete in the sense that sun
sensors and star sensors are vulnerable to occultation
and sole earth sensors cannot provide three-axis
information[1]. Furthermore, magnetic sensors are
applicable to wide-range of body rate. They can provide
rate information as a backup for gyros. For example,
coarse rate information can be derived when the gyro
output becomes saturated[2].
However, magnetic sensors measure magnetic field in
the spacecraft coordinates and donft provide three-axis
information with a single measurement. Therefore, to
determine a three-axis attitude, information from other
attitude sensors such as sun sensors or another
measurement of the magnetic field is required (Fig, 1).
Fig.1 Conventional attitude sensor configuration
(Magnetometer and sun sensor)
In the former case, requiring other attitude sensors is
unfavorable for a backup system. In the latter case, it is
necessary to wait for a while so that the orbital motion
of the spacecraft brings a well different position in the
earth magnetic field. Waiting for a while presumes a
small attitude change during the two measurements.
Then it is not applicable to rotating spacecraft.
We propose to combine magnetic sensors with rate
sensors to determine a three-axis attitude of spacecraft
that may change its attitude quickly, e.g. in case of a
rotating spacecraft, without using other attitude sensors.
This paper is organized as follows. In section 2, an
attitude determination algorithm using only magnetic
sensors is introduced presuming that the attitude change
during the measurement period is very small. In section
3, the algorithm is extended to combine magnetic
sensors with rate sensors, where both the three-axis
attitude and rate sensor biases are estimated. Section 4
gives numerical examples. Section 5 concludes the
paper.
2. STATIONARY ATTITUDE DETERMINATION
BY MAGNETIC SENSORS
We discuss an attitude determination algorithm using
only magnetic sensors presuming that the attitude
change during the measurement period is very small.
2.1 Measurement model
Magnetic sensors measure the local geomagnetic field.
After calibration of their scale factors and biases, their
outputs have a simple relationship with the magnetic
field vector given by a geomagnetic field model such as
Earth
Sun
Magnetic field
the International Geomagnetic Field Model (IGRF) as
B I
BI U r =r (1)
where BI U is a coordinate transformation matrix
from the spacecraft body coordinates to the inertial
coordinates, rB is a measured magnetic field vector
expressed in the body coordinates, and rI is a
reference magnetic field vector computed from the
geomagnetic field model and expressed in the inertial
coordinates[3-4].
One measurement gives two degree-of-freedom
(d.o.f.) information. To determine a three-axis attitude,
it is necessary to take measurements at two points that
are enough separated from each other (Fig. 2).
Fig.2 Measurements at two points
(Stationary case)
To average out the measurement errors and model
errors contained in Eq. (1), we take N measurements.
Then the attitude estimation problem is to find a
coordinates transformation matrix BI U that satisfies
the following relationship. The subscript i stands for i-th
measurement.
B I ( 1, , )
B I i i i N U r=r =?? (2)
2.2 Attitude determination algorithm
Suppose that two measurements by magnetic sensors
provide a magnetic unit vector in the body coordinates
at each orbit position and that earth magnetic field
model provides a magnetic unit vector in the earth fixed
coordinates (which can be easily converted into the
inertial coordinates) at each orbit position. Then a
simple method called algebraic method is known that
constitutes an orthogonal frame from two unit vectors in
each coordinates and computes a coordinates
transformation matrix (equivalently a three-axis
attitude) by combining the two orthogonal frames.
First, three axes of the body coordinates are
constructed from the measured vectors as
1 2
1 2
1 2
1 2
B B
BM
B B
B B
BM
B B
B B B
M M M
= +
+
= ~
~
= ~
x r r
r r
y r r
r r
z x y
(3)
Similarly, three axes of the inertial coordinates are
constructed from the reference vectors as
1 2
1 2
1 2
1 2
I I
IM
I I
I I
IM
I I
I I I
M M M
= +
+
= ~
~
= ~
x r r
r r
y r r
r r
z x y
(4)
Then, the coordinates transformation matrix BI U
is determined as
? , , , , I I I B B B T
BI M M M M M M U =??x y z?????x y z?? (5)
where symbol ?? indicates an estimated value, and the
notation ( )T indicates the transpose of the vector or
the matrix in question.
Because the measurement errors of the magnetic
sensors are comparatively large, just two measurements
may not achieve sufficient accuracy in attitude
determination. Therefore, we take N-times
measurements and solve for a coordinates
transformation matrix that match N vectors from
sensors to N vectors from the model. We define a cost
function J and find BI U that minimizes J.
2
1
min ( )
B I
N
I B
B I i B I i
i
J
=
= ? U
U r U r (6)
We start with the three-axis attitude estimated by the
algebraic method; ? BI U and update the attitude
estimates based on the differences between the
measured vectors and vectors that are transformed near
to the measured ones through the initial estimated
coordinates transformation matrix.
2
1
( ) ?
N
T I B
B I i i
i
J
=
U=U r?Ur (7)
? (new) ? B I B I U U U (8)
Noting that the attitude estimate update U is
small, we apply small angle approximations to U .
[ ] 3 U?I ? Á~ (9)
where 3 I is a 3x3 identity matrix and the notation
[Á~] indicates the 3x3 cross-product-equivalent
matrix associated with the 3-dimensional vector , i.e.
[ ]
3 2
3 1
2 1
0
0
0
? ? ?
~ =?? ? ??
??? ??
(10)
Substituting Eq. (9) into Eq. (7), we obtain
[ ] 2
1
N
B B
i i
i
J
=
?r +Á~r
(11)
where a residual vector B
i r is defined as
B ?T I B
i BIi i r =U r ?r (12)
From Eq. (11), we derive linearized equations in terms
of attitude estimate update.
2
1
N
B
i i
i
J
=
?r ?B (13)
where matrix i B is defined as
B
i i B=??r ~?? (14)
The global minimum of the cost function J can be
computed by solving for the first-order necessary
condition for an optimum:
J =
0
(15)
Then we obtain the attitude update ? as
1
1 1
?
N N
T T B
i i i i
i i
?
= =
? ?
=? ?
? ?
B B B r (16)
We can apply an iteration to improve estimation
accuracy.
2.3 Evaluation of estimation accuracy
The measured magnetic field vector contains a
measurement error. This error is modeled as a small
angle rotation.
[ ] B
i i i w = ˁ~r (17)
The small angle rotation i is supposed to be a
Gaussian, zero-mean random vector.
[ ] 2
3 , T
i i j ij E E =0 ?? ??= I (18)
where E[?] represents an operation of taking an
expectation, is the standard deviation of the small
angle rotations, and ij is Kronecker's delta.
From Eqs. (16) and (17), estimation errors are
computed as
[ ]
1
1 1
?
N N
T T B
i i i i i
i i
?
= =
? ?
= ? =? ? ~
? ?
?? B B B r (19)
From Eqs. (18) and (19), statistical evaluation of
attitude estimation errors is obtained as
E[??] = 0 ,
1
2
1
N
T T
i i
i
E
?
=
?? ?? = ??? ???
???? B B (20)
Noting that the following identity holds for an arbitrary
3x1 unit vector x ,
[ ][ ] 3
x~ x~T=I ?xxT (21)
we rewrite Eq. (20) into
1
2
3
1
T N BBT
i i
i
E
?
=
?? ??= ??? ? ???
???? I r r (22)
Recalling that the measured magnetic vector rB does
not change much in a short period, we can deduce two
conclusions from Eq. (22).
1) The attitude estimation error perpendicular to rB is
approximately evaluated as / N
2) The attitude estimation error parallel to rB depends
on the amount of changes in rB and can be fairly
large.
3. ROTATING ATTITUDE
DETERMINATION BY MAGNETIC
SENSORS AND RATE SENSORS
We extend the attitude determination algorithm
developed in section 2 to combine magnetic sensors
with rate sensors in order to determine a three-axis
attitude of spacecraft that may change its attitude
quickly. The basic idea is to apply body rate estimates
to compensate for the rotational so that the problem of
attitude determination is reduced similar to a stationary
one (Fig. 3).
Fig. 3 Measurement at two points
(Rotating case)
Rotational motion
is compensated by
applying coordinates
transformation.
3.1 Measurement model
The rate sensors measure the body rates of the
spacecraft. Their outputs contain a bias, which can be
regarded as constant during a short time period.
mes
i =i+b (23)
where mes ,
i i and b is measured body rates, their
true values and a measurement bias vector, respectively,
and are expressed in the body coordinates.
In combination with the measured body rates mes
i ,
the relationships between the measured magnetic field
vector and the reference vector become
(i) (1) (i) (1)
B I I
i IBt i BtBt IBt i r =U r=U U r (24)
where matrix
B(t1)B(ti) U is a coordinates
transformation matrix from the initial body coordinates
to the body coordinates at i-th sampling time.
(1) ( ) 1
2 1
( )
( ) ( )
i
mes b
B t B t i
mes b mes b
t
t t
? = ? ~
~ ? ~ ?
U U ??
U U
(25)
When a rate sensor is allocated to each axis, the
unknown variables are the three-axis attitude and three
biases, which means 6 unknowns in total. Therefore, it
is necessary to take measurements at three points that
are enough separated from each other. To average out
the measurement errors and model errors contained in
Eqs. (23) and (24), we take N measurements. Then the
attitude estimation problem is to find an initial
coordinates transformation matrix B(t1)I U and a bias
vector b that satisfy Eqs. (23) and (24) for
i=1??N.
Then we can compute a coordinates transformation
matrix at i-th sampling time as
IB(ti) B(t1)B(ti) IB(t1) U =U U (26)
3.2 Attitude determination algorithm
Similar to subsection 2.2, we define a cost function J
and find B(t1)I U and b that minimizes J.
1
(1)
, ( )
2
( )
1
min ( , )
B t I b
i
B t I b
N
I B
i Bt Ii
i
J
=
= ?
U
U
r U r
(27)
We first assume b is zero and compensate the
rotational motion of the spacecraft by applying the
following coordinates transformation.
2 (1) (2) 2
B TB
B t B t r U r (29)
With this rotation-compensated vector, we compute the
initial estimate of (1)
? B t I U using Eqs. (3), (4) and (5).
Then we update the attitude estimates based on the
differences between the measured vectors and vectors
that are transformed near to the measured ones through
the initial estimated coordinates transformation matrix
combined with measured body rates subtracted by the
estimated bias vector.
(1) (1)
?
B t I B t I U =U U (30)
mes
i i b b = ?? (31)
[ ] 3 U?I ? Á~ (32)
Similar to section 2, we derive linearized equations in
terms of attitude estimate update and a bias vector.
(1)
2
, 1
min ( , )
B t I b
N
B
U b i i i i b i
J t
=
=r ?B?C (33)
where a residual vector B
i r and matrices i B and
i C are defined as
(1) (1) ( )
?
i
B T I T B
i BtI i BtBt i r =U r ?U r (34)
( ) (1) ( i)
T B
i Bt B t i
B=??U r ~?? (35)
( )
( )
( )
1
1 2 2
1 1 1
( ) ( )
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
i
i
i i i
T B
i Bt B t i
T TB
B t B t B t B t i
T TB
B t B t? B t? B t i
?? ~??
+ ?? ~??
+ + ?? ~??
C U r
U U r
U U r
??
??
(36)
The global minimum of the cost function J can be
computed by solving for the first-order necessary
condition for an optimum:
,
b
J J
= =
0 0
(37)
Then we obtain the attitude update ? and a bias update
b as
1
0 1 0
1 2 1
?
b
? ? ? ? ? ? ?
? ?= ? ? ? ?
? ? ? ? ? ?
C C d
C C d
(38)
0 0
1 1
1 0
1 1
2
1
,
,
N N
T T B
i i i i
i i
N N
T T B
i i i i
i i
N
T
i i
i
= =
= =
=
= =
= =
=
C BB d B r
C CB d C r
C CC
(39)
We can apply an iteration to improve estimation
accuracy.
3.3 Evaluation of estimation accuracy
Similarly to section 2.3, statistical evaluation of
attitude and bias estimation errors is obtained as
?
b
E
? ? ? ?
? ?= ? ?
? ? ? ?
0
0
(40)
1
2 0 1
1 2
T
b b
E
?? ?? ? ? ? ??
?? ?? ? ?= ? ?
??? ?? ? ?? ? ?
C C
C C
?? ??
?? ??
(41)
4. NUMERICAL EXAMPLES
4.1 Simulation conditions
To evaluate the proposed algorithm, a slowly rotating
spacecraft is chosen as an numerical example.
The simulation condition is summarized in Table 1.
Initial attitude (quaternion) and rate sensor bias are
randomly chosen.
Table 1 Simulation condition.
Orbit height 1000 km
Eccentricity 0
Inclination 30 degrees
Geomagnetic field model Dipole model
Attitude estimation error 0.05 rad-sigma per axis
Initial attitude -0.247,-0.952,0.072, 0.164
Initial body rate 2, 2, 2 deg/s
Rate sensor bias 0.1, 0.1, 0.1 deg/s
Total measurement period 600 seconds
Sampling period 2 seconds
4.2 Simulation results
Figures 4 and 5 are the results of parametric study
and they show attitude and rate sensor bias estimation
error with respect to magnetic sensor output error for a
spacecraft whose initial body rates are two degrees/s per
each axis.
In the figures, the estimation errors are divided into
three components, i.e. parallel to the magnetic field
vector ( r?? ), parallel to the time derivative of r?? ( r?? ?? ),
and parallel to the cross product of r and r??
((r~r)?? ?? ).
From the figures, it is shown that
1) Both attitude and body rate estimation errors are in
proportion to magnetic sensor output error,
2) Attitude estimation error parallel to the magnetic
field vector is much larger than attitude estimation
errors perpendicular to the magnetic field vector
3) Rate bias estimation errors are on the same order for
each axis.
The third finding suggests specific application of
attitude determination using magnetic sensors, i.e.
estimating body rate information from magnetic
sensors.
5. CONCLUSIONS
We proposed a new attitude determination algorithm
for a spacecraft in low earth orbits using magnetic
sensors and rate sensors, which is applicable to a
rotating spacecraft and does not require additional
attitude sensors. First, an attitude determination
algorithm using only magnetic sensors was introduced
presuming that the attitude change is very small. Then
the algorithm was extended to combine magnetic
sensors with rate sensors, where both the three-axis
attitude and rate sensor biases are estimated. The
estimation error was statistically analyzed. Numerical
examples showed validity of the attitude determination
algorithm and its accuracy evaluation.
REFERENCES
[1] J. Merayo, et al, gThe Spinning Astrid-2 Satellite
Used for Modeling the Earthfs Main Magnetic
Fieldh, IEEE Trans. Geoscience and Remote
Sensing, Vol. 40, No. 4, pp. 898-909, 2002.
[2] M. Challa, et. Al, gMagnetometer-Only Attitude
and Rate Estimates during the Earth Radiation
Budget Satellite 1987 Control Anomalyh,
AIAA-97-3615, 1997.
[3] J. R. Wertz (ed.), Spacecraft Attitude
Determination and Control, pp. 249-254, D. Reidel
Pub. Co., Boston, 1978.
[4] M. Challa and R. Harman, gA New Magnetometer
Calibration Algorithm and Applicationsh, Proc.
AIAA Guidance, Navigation, and Control
Conference, Boston, MA, 1998, pp. 697-707.
r??
r?? ??
(r~r)?? ??
Each sample
}1 sigma line
Magnetic sensor output error [deg-]
Attitude estimation error [deg]
Fig. 4 Attitude estimation error (rotating case)
(A dot represents an each sample. Short lines show +/- 1 sigma values derived from analytical study.
Symbol r?? shows the component parallel to the magnetic field vector, r?? ?? shows the component parallel to the
time derivative of r?? , and (r~r)?? ?? shows the component parallel to the cross product of r and r?? .)
Each sample
r }1 sigma line ??
r?? ??
(r~r)?? ??
Rate sensor bias estimation error [deg/s]
Magnetic sensor output error [deg-]
Fig. 5 Rate sensor bias estimation error (rotating case)
(A dot represents an each sample. Short lines show +/- 1 sigma values derived from analytical study.
Symbol r?? shows the component parallel to the magnetic field vector, r?? ?? shows the component parallel to the time
derivative of r?? , and (r~r)?? ?? shows the component parallel to the cross product of r and r?? .)