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?? .)