- Email: [email protected]
Attitude estimation from magnetometer and earth-albedo-corrected coarse sun sensor measurements
Acta Astronautica 56 (2005) 115 – 126 www.elsevier.com/locate/actaastro
Attitude estimation from magnetometer and earth-albedo-corrected coarse sun sensor measurements Pontus Appel∗ Center of Applied Space Technology and Microgravity, ZARM, University of Bremen, Am Fallturm, 28359 Bremen, Germany Available online 28 October 2004
Abstract For full 3-axes attitude determination the magnetic field vector and the Sun vector can be used. A Coarse Sun Sensor consisting of six solar cells placed on each of the six outer surfaces of the satellite is used for Sun vector determination. This robust and low cost setup is sensitive to surrounding light sources as it sees the whole sky. To compensate for the largest error source, the Earth, an albedo model is developed. The total albedo light vector has contributions from the Earth surface which is illuminated by the Sun and visible from the satellite. Depending on the reflectivity of the Earth surface, the satellite’s position and the Sun’s position the albedo light changes. This cannot be calculated analytically and hence a numerical model is developed. For on-board computer use the Earth albedo model consisting of data tables is transferred into polynomial functions in order to save memory space. For an absolute worst case the attitude determination error can be held below 2◦ . In a nominal case it is better than 1◦ . © 2004 Elsevier Ltd. All rights reserved.
1. Introduction The attitude control system is an important satellite sub-system. On one hand it must achieve a certain accuracy. On the other hand it must be very robust to ensure the mission success. Furthermore, it shall be equipped with low-cost hardware. The solution is the utilization of simple, robust sensors in combination with intelligent algorithms. One very reliable and simple possibility is the combination of a magnetometer and a Coarse Sun Sensor
∗ Tel.: +49 (0) 421 218 3635; fax: +49 (0) 421 218 2526.
E-mail address: [email protected] (P. Appel). 0094-5765/$ - see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.actaastro.2004.09.001
(CSS). This will provide two vectors, and hence a full 3-axis attitude determination. The CSS consists of six solar cells that are placed on the surface on each side of the satellite, see Fig. 1. Using the cosine law, the total light vector, having contributions from both Sun light as well as Earth albedo, is obtained. In order to get the correct vector a measurement model of the albedo must be included. This paper covers the development of an attitude determination system using a CSS and a magnetometer. In the first part an albedo model which is applied for sensor simulation and for attitude determination is developed. In the second part an algorithm for attitude determination using an extended Kalman filter for state estimation is presented.
116
P. Appel / Acta Astronautica 56 (2005) 115 – 126
Fig. 1. Solar cell.
2. Albedo model The CSS used for Sun vector determination consists of six solar cells placed on each of the six outer surfaces of a cube shaped satellite. As the CSS, unlike a fine Sun sensor, sees the whole sky it is much more sensitive to other light sources than the Sun. The largest of these sources for a LEO is the light reflected by the Earth, the Earth albedo. To be able to get a good estimation of the Sun vector the Earth albedo has to be compensated for. This is made by developing an numerical albedo model and including it in the Kalman filter for attitude estimation. In order to calculate the amount of light received by a satellite in orbit a model for the light reflected by Earth is needed. The incoming light is a function of the satellite position, the Sun’s position and the reflectivity of Earth. The Earth is assumed to emit the radiation like a black body and the main part of the light is assumed to be reflected by Earth diffusely. The radiation from a black body is isotropic, that is, it has no preferred direction. For such a body the Lambertian cosine law is applicable. As the incoming sunlight hits the Earth it is diffusely reflected by the atmosphere and Earth surface in all directions. The amount of emitted light is determined by the reflectivity of the element and the angle of the incoming light. The amount of emitted light in the direction of the satellite is decided by the angle
Fig. 2. Description of the path of light geometry.
between the surface element normal and the satellite vector. The voltage delivered by the solar cell receiving the light is also calculated with the cosine law with the satellite vector and solar cell surface normal as components. The path of the light and the vectors are described in Fig. 2. The amount of incoming light to the satellite is equal to the integral of the three vector products over the illuminated and visible area of Earth. The total received energy is ˙ r = Is Acell Q visible and illuminated sphere area
(R e,0 • S 0 ) · (S 0 • n0 ) · (R e,0 • r 0 ) × dA, (1) S2 where Is is the solar intensity at the Earth, Acell the solar cell surface area, the reflectivity of each surface element of Earth, R e,0 the surface element normal, S 0 the vector from the surface element towards the satellite, n0 the solar cell surface normal, and r 0 the Sun vector. The solution of Eq. (1) is needed for every point in the orbit of a satellite. Because of its complexity with the non-trivial limits of integration it cannot be solved analytically. Together with the fact that the reflectivity of Earth is far from being uniform this implies that a numerical solution is the best way to solve the problem.
P. Appel / Acta Astronautica 56 (2005) 115 – 126
117
Fig. 3. Daily reflectivity of the Earth. Fig. 4. Latitude dependency of Earth’s reflectivity.
2.1. Reflectivity model First, a numerical model of the Earth is created by dividing a sphere into a finite number of surface elements. A reflectivity model is developed to describe the reflectivity of every single element. 2.2. Reflectivity data To get a realistic reflectivity distribution reflectivity data from existing measurements is analyzed. The data is taken from the NASA project Total Ozone Mapping Spectrometer [1]. When plotted the data shows Earth’s reflectivity as a function of the latitude and longitude. Fig. 3 shows the reflectivity reading from 1 day. The bright areas represent high reflectivity and the dark low. From the plot the Polar Regions and Greenland (lat. 170◦ , long. 150◦ ) can be recognized as well as several cloud formations. It can be seen that a lot of the reflected light origins from these cloud formations and the Polar Regions. From these plots a pattern could be recognized. It was obvious that the light reflected by Earth had a large dependency on the latitude, but only a small on longitude. Fig. 4 shows an example of the reflectivity of Earth as a function of its latitude. The solid curve is the average reflectivity and the dashed the standard deviation as a function of the latitude of Earth. The reflectivity has a maximum in the Polar Regions where the deviations are smallest. This
is due to the constant ice or snow coverage. For lower latitudes the reflectivity decreases and the deviations increase, respectively. The deviations are due to local differences in reflectivity of the surface but also due to tidal changes of, for example, cloud formations which have a high reflectivity. Knowing this, a statistical model of the reflectivity can be developed. From the data averages and standard deviations are derived for each latitude. The reflectivity model is created by randomizing new reflectivities based on these averages and deviations for each latitude. The data is merged in order to generate a model of the whole surface of the Earth. This reflectivity pattern is applied to the numerical model of Earth. 2.3. Received light As the light reflected by Earth is emitted diffusely the cosine law is applicable to calculate the received light by the satellite in any position above the Earth surface. It will suffice to calculate under which angle each surface element is illuminated by the Sun and from which angle it is seen by the satellite. The vector towards Sun is known and so is also the position of the satellite. For each element a normal is calculated as is also the vector pointing towards the satellite. The vector products between the Sun vector and the surface normal as well as between the surface normal and the satellite vector are calculated, respectively. Both
118
P. Appel / Acta Astronautica 56 (2005) 115 – 126
2.5. Recording data The above procedure is run through several times, each time with a new randomized reflectivity applied to the numerical Earth and for all satellite positions. The mean received albedo intensity of the three components is recorded as well as the standard deviations for every run. This results in three matrices describing the mean of the albedo light vector in the three components and three matrices describing their deviations. 2.6. Analyzing data
Fig. 5. Albedo light vector at satellite position.
products are multiplied with each other and the reflectivity of the element. This is the numerical solution to Eq. (1) and the solution can now be calculated for every position. 2.4. The grid A shell describing a grid of satellite locations at a certain altitude above the Earth surface is created. In every grid point the elements that are visible from the satellite and are illuminated by the Sun are picked out. For the elements surviving the test the emitted albedo light vector is calculated using the cosine law as above. Now the sum of the contributions from all elements is calculated for each satellite point in the grid, see Fig. 5, where Itotal =
n
Ik .
(2)
k=1
This total albedo light vector is split up into three components, north, east and down. This coordinate system is chosen to make the curve fitting to come easier, that is, to avoid sine- and cosine patterns from the transformation to occur in the plots. With knowledge of the attitude of the satellite having its solar cells in this position this would be the numerical representation of Eq. (1). The attitude is not taken into calculation at this moment.
Fig. 5 shows a result from these simulations. It describes the down component of the albedo light vector. The intensity on the z-axis is normalized w.r.t. the solar intensity. The zero longitude plotted is defined as the noon–midnight meridian. The zero latitude is at the north pole; hence the plots show the latitudes beginning from the shadow side of earth, sweeping over the north pole, the south pole and back to the shadow side. These six data sets are only valid for this altitude and for this Sun position. For other setups new data sets have to be developed. For a satellite that stays a long time in orbit and/or has a eccentric orbit this means a lot of data sets and so a lot of memory capacity needed. Because of this the data sets were transferred into functions instead of being saved as look-up tables. Functions are very memory saving as they only need to include the functions themselves and a corresponding parameter vector for the current setup. 2.7. Polynomial derivation From the six three-dimensional plots an attempt to fit polynomials is made. The polynomials are used to describe the albedo intensity and the standard deviation, received by the satellite for an arbitrary position in the orbit. The mean parts of the polynomials are used in the Kalman filter as expected measurements as they deliver the solution to Eq. (1) for all positions. All polynomials together are used to simulate measurements when testing the filter as described later. 2.8. Curve fitting The three-dimensional plots describing the incoming light are fitted with two-dimensional polynomials
P. Appel / Acta Astronautica 56 (2005) 115 – 126
Intensity [percent]
Intensity [percent]
0.2 0.15 0.1 0.05 0 200 150 100
150
50 Longitude [degrees]
0 0 -100 -50
200
250
300
0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 200 150 100
100 50 Latitude [degrees]
50 Longitude [degrees]
Fig. 6. The data representing the down component.
0
-100
-50
0
50
100 150
200 250
300
Latitude [degrees]
Fig. 7. The fitted model representing the down component.
with the structure I = (Af (latitude) + B)(Cg(longitude) + D),
119
3. The sun sensor (3)
where I is the intensity received by the satellite , f (latitude) is a function describing the dependency on the latitude and g(longitude) a function describing the dependency on longitude. A, B, C and D are scalars that are used to fit the plots to the data. The polynomials for the mean values of the vector components are composed of three sine or cosine functions to describe the latitude dependency and one sine function to describe the longitude dependency. The polynomials for the mean values have the form y = (P1 (P5 sin(P6 ) + P7 sin(P8 ) + P9 cos(P10 + P13 )) + P2 )(P3 (P11 sin(P12 )) + P4 ), (4) where P1 –P13 are elements in a parameter vector, the latitude of the satellite and the longitude. P1 –P4 represent the scalars A–D in Eq. (4), the others represent amplitude, phase and phase displacement of the sine functions. The upper row describes the latitude dependency and the lower the longitude dependency. Each vector component has its own parameter vector. The result of the curve fitting can be seen in Figs. 6 and 7. The plot is the fitted model to Fig. 6. The parameter vector creating this plot is saved and later used in the Kalman filter.
The CSS used to measure the Sun vector consists of six solar cells placed on each of the outside surfaces of the satellite. Using the cosine law the amount of incoming light can be calculated for each cell and so the angle towards the Sun. 3.1. The solar cells The sensor head is developed at ZARM [2]. Fig. 1 shows such using two 20 × 20 mm solar cells. 3.2. Cell placement Fig. 8 describes one way to place the cells on a cubeshaped satellite. The cells should be placed where the risk of being shaded by other objects is small. The coordinate system used for determination of the direction of the solar cell surface normal is a body-fixed one. 3.3. Implementation The solar cell delivers a voltage that is proportional to the incoming light. The amplitude of the voltage is decided by the cosine law. With tests the characteristics of the cell were examined. Fig. 9 shows some results from these tests. The figure shows the resulting normalized current through a resistor connected to the cell as a function of the angle of the incoming light. The dark/lower curve represents the test results and the brighter/upper one
120
P. Appel / Acta Astronautica 56 (2005) 115 – 126
Fig. 10. Kalman filter work order.
Fig. 8. Solar cell orientation.
The principle of the filter setup used in this case is shown in Fig. 10, where the lower parts describe the model of the real case with an update part to make corrections when necessary. To work properly the Kalman filter needs information about the system dynamics and its partial derivatives with respect to the state vector [3]. The state vector is composed of the rotation rates and the quaternions describing the spacecraft attitude w.r.t. the inertial frame x(t) = (x , y , z , q1 , q2 , q3 , q4 )T ,
(5)
where is the rotation rate vector and q the quaternion vector. The time derivative of the state vector is also needed: Fig. 9. Solar cell characteristics.
is a cosine curve. The results show that the modelling of the delivered voltage from the cell can be done with a simple cosine function.
4. Kalman filter for attitude estimation
˙ x, ˙ y, ˙ z , q˙1 , q˙2 , q˙3 , q˙4 )T = f (x). x(t) ˙ = (
(6)
The state error covariance propagation part of the filter also needs the partial derivative of the state vector x w.r.t. itself. This is the F-matrix *f (x(t), ˆ t) F (x(t), ˆ t) = (7) *x(t) ˆ x(t)=x(t) ˆ
or An extended Kalman filter is used to estimate the state of the satellite. The filter must have knowledge about the dynamics of the systems. It needs a statistical description of the noises and errors in the system and measurements- and attitude dynamics models. It also uses any available information about the initial condition of the variables.
* ˙ F=
*f *x T
* = *q˙
*
˙ * *q . *q˙
*q
The resulting matrix F is a 7 × 7 element matrix.
(8)
P. Appel / Acta Astronautica 56 (2005) 115 – 126
(9)
in inertial frame, and eb# the vector normal to the solar cell surface. Also
0, f < 0, g{f } = (14) f, f 0.
(10)
Together they deliver the expected measurement vector from the magnetometer and the CSS to the filter
4.1. Attitude dynamics model The attitude dynamics equations of motion for a rigid body are ˙ = I −1 [T dist + T control − × (I )], q˙ = 21 q,
0 −3 = 2 −1
3 0 −1 −2
−2 1 0 −3
1 2 , 3 0
Y = h(x) = [uCSS , umag ]. (11)
where I is the inertia matrix of the spacecraft. The vector T dist includes all external disturbance torques acting on the satellite and the vector T control all control torques. The external torques taken into account are those from the gravity gradient, magnetic and aerodynamic drag and solar pressure. 4.2. Measurement model The most important inputs to the filter are the measurements from the magnetometer and the CSS. These have to be manipulated to be interpretable by the filter. The Kalman filter uses the magnetic field vector measured by the magnetometer and the scaled light intensity on each of the solar cells of the CSS for the attitude determination. These measurements have to be modelled for the Kalman filter to know what to expect. The measurement model for the magnetometer is umag = Abi (q)B i ,
The observation matrix also needed by the Kalman filter for estimation of the expected measurements is the measurement vector y = h(x) partially differentiated w.r.t. the state vector *hk (x(tk )) H k (xˆ k (−)) = . (16) *x T (tk ) x(tk )=xˆ k (−) 4.3. Filter simulation All data needed for the simulations are acquired from the program SATSIM [5]. It is a simulation tool developed by the Center of Applied Space Technology and Microgravity (ZARM) to provide a test bench for the development of small satellites. The environment- and dynamics simulation tool of the SATSIM program simulates the orbital- and attitude dynamics of the satellite. 5. Results 5.1. Scenario definition
uCSS# = g{eb# · Abi (q) · P isun } T
T
(15)
(12)
where Abi is the transformation matrix from inertial to body fixed frame, umag is the magnetometer measurement and B i is the magnetic field vector in inertial frame from the IGRF model [4]. The measurement model for the #:th solar cell of the Sun sensor is
+ g{eb# · Abi (q) · P iEarth }
121
(13)
where # = 1, 2, . . . , 6 and uCSS# is the intensity measurement, P isun the sunlight vector in inertial frame, P iEarth the albedo light vector from polynomials,
The orbit used in the simulations presented in this paper is a polar one with its normal perpendicular to the Sun vector. The Sun is assumed to be located in the equatorial plane. An orbit altitude of 630 km is used. The simulation is started at the equator on the illuminated side of Earth and runs for about three orbits. The satellite is of cubic form with a side length of 1 m and a mass of 300 kg. The moment of inertia matrix is fully occupied because of the displacement of the center of gravity with respect to the geometrical center of the satellite. This displacement is unknown to the Kalman filter. The center of gravity is displaced 0.1 m in the x-direction, 0.05 m in the y-direction and
122
P. Appel / Acta Astronautica 56 (2005) 115 – 126
0.04 m in the z-direction. The moment of inertia matrix is
53 I = 1.5 1.2
1.5 50.75 0.6
1.2 0.6 kg m2 . 50.48
(17)
5.5. Total torque The total torque N Total acting on the satellite is the sum of the torques originating from aerodynamic drag, solar torque, gravity gradient and magnetic torque. N Total = N Aero + N solar + N GG + N B .
5.2. Initial conditions To test the performance of the filter in the initial phase the state vector is assumed not to be known well. The satellite is assumed to have completed the detumbling phase and to have a total rotation rate of 0.01 rad/s in the simulations. This corresponds to about 10 revolutions per orbit. This quite fast rotation has large influence on the attitude determination if the moment of inertia matrix is poorly known. 5.3. Satellite environment To make the simulations more realistic some assumptions on forces and torques acting on the satellite and parameter errors that can occur have been made. In addition two parameter errors have been introduced; one of the moment of inertia matrix elements has been increased by 5%. This is a reasonable value as small satellites with limited budget may have uncertainties in the mass configuration of the spacecraft. Also the difference in propellant mass, if present and not modelled, over time affects the moment of inertia of the satellite. It can furthermore be assumed that the position of the satellite cannot be precisely determined. To simulate this a position error has been introduced. The Kalman filter is given a position that differs up to 1◦ in latitude and longitude from the true value. 5.4. Disturbances The external torques assumed influencing the satellite are those from the gravity gradient, magnetic- and aerodynamic drag and solar pressure [6]. They are all acquired from SATSIM. The only disturbance assumed known to the filter is the gravity gradient. That is, the filter is able to calculate the gravity gradient force and compensate for the errors it could produce.
(18)
It varies between 1×10−5 and 5×10−7 N m depending on the position and attitude of the satellite. 5.6. Filter tuning In the filter the predicted measurements computed from the predicted state are compared with the actual ones. This is done so that the correction of the state vector can be specified in order to minimize the prediction error. There are two major problems when trying to get a good estimate of the attitude; the first are the noises and disturbances in the measurements where the albedo deviation is the most important and in the system dynamics. The second is the lack of observability that occurs when there is a small angle between the two measured vectors. Both problems can be compensated for by changing the amplitude of the measurement covariance matrix. To optimize the Kalman filter performance it must be tuned. This is done by changing the amplitude of the elements in the covariance matrices. The following sections explain where it is necessary and how the changing of the covariance matrix amplitude can enhance the performance of the filter. 5.7. State covariance matrix During the simulations the state covariance matrix is kept constant. That is, the external forces are assumed not changing that much over time. This is of course not true, but the filter is tuned to be able to handle the largest disturbances. The differences in external disturbances can also be compensated for, but that is not done in this phase. In the state covariance matrix the three first elements in the diagonal are representing the rotation rates of the satellite and the four last the quaternions. The rotation rate values are smaller than those of the quaternions; they therefore are represented by smaller values in the matrix.
P. Appel / Acta Astronautica 56 (2005) 115 – 126
5.8. Measurement covariance matrix
123
The components of the measurement covariance matrix used are of different amplitude. The first six elements in the diagonal represent the solar cell measurements and the three last the magnetometer measurements. The noise in the light vector measurements is in order of ten percent. It is composed of the albedo model noise, and sensor noise. The assumed noise in the magnetic field measurements are in the order of one percent and as the covariance matrix is the square of the deviations the elements representing the solar cells are set to a 100 times larger value than those representing the magnetic field vector. This relationship has also been tried and verified with simulations.
the measurement covariance matrix (R-matrix) amplitude is set higher to make the filter rely on the system dynamics and therefore get a steadier and more correct estimation where the noises in the measurements are high. This is only true to a certain point. A greater measurement covariance says that the measurements are not that trustworthy and the filter should rely more on the equations of motions to estimate the state. But when external forces are present the filter will be unable to predict the movements and must use information from the measurements to not drift away. Using the polynomials for the standard deviations of the expected albedo radiation, the R-matrix can be changed continuously to always be adjusted to the measurement noise expected.
5.9. Initial state covariance
5.11. Observability issues
The initial state covariance matrix describes to what extent the initial state can be trusted. A small value indicates that the given parameter is a good approximation and vice versa. The elements are arranged as in the state covariance matrix. The initial state covariance has been set to zero for the rotation rates. This means that the filter is told that the given initial rotation rates are correct although they are not. The initial state covariance of the quaternions is set to 100. This combination makes the filter believe in the rotation rates and therefore it starts correcting the quaternions and makes initially only small attempts correcting the rotation rates. After the right attitude has been found the rotation rates slowly converge also. This can be seen in Fig. 13 which shows results from the same simulation as in Fig. 12. The rotation rates converge after about 100 s and the quaternions after less than 5 s. In Fig. 13 it can also be seen the greater sensitivity to external forces for the rotation rates up to 100 s than after. This is due to the low R-matrix amplitude set during the first 100 s. If the rotation rate covariance is set higher the filter may have problems converging as the filter tries to correct the rotation rates without the reference that a known attitude provides.
The Kalman filter uses the information from the sensors to create two vectors. Using these vectors the filter is able to determine the attitude of the satellite. This is only possible when the two vectors are both present and not pointing in the same direction. The following section explains where it occurs and how, if possible, it can be compensated for.
5.10. Steady state phase During the steady state phase a stable and accurate estimation of the state is wanted. After the initial phase
5.12. Parallel measurement vectors For a high inclined orbit whose normal is close to perpendicular to the Sun vector a problem concerning the observability can occur. This happens four times per orbit as the magnetic field vector is close to parallel to the Sun vector. In the simulations this is visible as larger deviations in the attitude determination than in a nominal case. 5.13. Sensor simulation To test the performance of the filter, several simulations have to be made. Measurements have to be simulated as input to the filter. The magnetic field vector is computed by a function based on the IGRF model. The model is distorted by white noise with a 200 nT amplitude [6] run through a low pass filter to simulate the local and tidal changes. Another white noise with amplitude 2 nT is added to simulate measurement noise. The magnetic field vector delivered to the filter is described in an inertial frame. The received light is
124
P. Appel / Acta Astronautica 56 (2005) 115 – 126
composed of one part coming from the sun, and another part that has been reflected by Earth. The latter part is created using the polynomials for the mean light vector as a base. A random number with a position dependant standard deviation is added. This number is derived from those polynomials of the albedo model describing the uncertainty of the reflected light. The result is put through a low pass filter to simulate the low frequencies at which the albedo measurements can be assumed to change. The light vectors are added and split up into six non-negative scalars representing the solar cells on all sides of the satellite. To these scalars a white, unfiltered noise with amplitude 0.5% of the measured intensity is added before they are given to the filter as measurements. This reflects the electronics noise. 5.14. Filter performance Fig. 11 shows some results from the absolute worst case simulations including the following error sources: • Disturbance torque unknown to Kalman filter with amplitude 10−7 –10−5 Nm. • Measurement noise, about 0.5% of measurement amplitude. • Albedo noise. • Uncertainty in moment of inertia, 5%. • Position error up to 1◦ in latitude and longitude. Picture A describes the intensity of the reflected light, normalized w.r.t. the solar intensity in the NED frame. The curve with the largest amplitude is the down component, the second greatest the north component and the weakest is the east component of the incoming vector. This intensity distribution is natural as the satellite in this case is passing right over the most radiating part of Earth, that is, the part receiving the most sunlight. The B plot shows the modelled standard deviations for the Earth’s albedo vector. They are also expressed in the NED-frame and are normalized w.r.t. the sunlight intensity. It can be seen that the largest deviations are present in the Down component. This is expected as this component is the largest. The deviations are up to 10% of the albedo radiation. The two other components are similar in size and partially also larger
Fig. 11. Albedo model simulation results.
than the mean values they correspond to. This has the consequence that the mean values are to be trusted poorly. The C plot shows the angular deviation between the estimated quaternion vector and the real one over time expressed in degrees. The coordinate system is the body-fixed one. The shaded stripes represent the satellite entering the shadow phase of its orbit. The D plot shows the observability of the filter. This is presented calculating the sine value of the angle between magnetic field vector and total light vector. A high value indicates good observability, that is, the two vectors used are close to being perpendicular to each other. A low value indicates poor observability. This happens when the angle between the vectors decreases. The result of this can be seen in most plots where the
P. Appel / Acta Astronautica 56 (2005) 115 – 126 x y z
50
x10-3
0
30
Rotation rates [radians/sec]
Angular deviation [degrees]
40
1
125
20 10 0 -10 -20 -30 -40
x y z
-1 -2 -3 -4 -5 -6 -7 -8
-50
-9
0
1
2
3
4
5 6 Time[s]
7
8
9
10
-10 0
50
100
150
Time [seconds]
Fig. 12. The Quaternions during the transient phase. Fig. 13. The rotation rates during the transient phase.
deviations are increasing at the same places where the observability shows a low value. From the C plot in every set some general conclusions can be made. The filter performance is best close to the Polar Regions where the deviations are small. At lower latitudes the albedo variance rises and the error rises with it. The peaks in the illuminated phases are a result of poor observability. In the Sun phase they coincide with the lower points in the observability plots. The shaded stripes represent the satellite entering Earth’s shadow. When the satellite enters the Earth shadow the filter can only use the magnetic field vector as measurement and it will therefore be more sensitive to dynamic disturbances. This is visible as larger deviations in this phase. The overall performance of the filter is good. The errors in this worst case scenario are not greater than 2◦ in the illuminated phase, and for a nominal case the errors can be kept below 1◦ . 5.15. Transient phase Figs. 12 and 13 show the transient phase of the same simulation described above. The filter converges quickly. The transient time is shortened down to under 10 s using smaller values for the measurement covariance matrix during this phase. The filter converges and is stable when started in the illuminated phase, and can be kept stable during the shadow phase.
5.16. Requirements To be able to perform correct attitude estimation the on-board computer must have access to the following inputs • satellite position (orbit propagator), • Sun position, • albedo model (polynomial function and parameter vector), • magnetic field model (IGRF), • satellite preferences (moment of inertia matrix, etc.). 6. Summary/Conclusions 6.1. Summary An Extended Kalman filter utilizing a CSS in combination with a model of the Earth’s albedo can drastically improve attitude determination. The absolute worst case accuracy is 2◦ and in the nominal case its better than 1◦ with the same disturbances. The accuracy of the direction determination of the Sun vector is good and the costs much lower than when using a fine Sun sensor. The sensors are also reliable because of their simple construction. As a complement or substitute to the small sensors the ordinary solar cells of a satellite can be used.
126
P. Appel / Acta Astronautica 56 (2005) 115 – 126
A combination of CSS and magnetometer is simple, robust and accurate enough for acquisition and safe modes. It is also sufficient for nominal modes of most small satellite applications. 6.2. Conclusion The tasks presented in this paper can be summarized in the following table: 1. Albedo model for CSS simulation and error correction developed. 2. Polynomial functions for on-board computer use developed. 3. Extended Kalman filter for attitude estimation developed. 4. 2◦ absolute worst case accuracy achieved. 1◦ nominal case accuracy.
References [1] Total ozone mapping spectrometer, http://toms.gsfc.nasa.gov/ [2] ZARM—University of Bremen, http://www.zarm.unibremen.de [3] M.S. Grewal, A.P. Andrews, Kalman Filtering—Theory and Practice, Prentice-Hall, Englewood Cliffs, NJ, 1993. [4] International Geomagnetic Reference Field, http://nssdc.gsfc. nasa.gov/space/model/models/igrf.html [5] S. Theil, M. Wiegand, H.J. Rath, SATSIM—A generic tool for attitude control system development, simulation and testing, 51st International Astronautical Congress, Rio de Janerio, Brazil, 2–6 Oct 2000, IAF-00-U.2.05. [6] M. Wiegand, Autonomous satellite navigation via kalman filtering of magnetometer data, 46th Congress of the International Astronautical Federation, Oslo, Norway, October 2–6, 1995, ST-95-W.1.04.
Attitude estimation from magnetometer and earth-albedo-corrected coarse sun sensor measurements Pontus Appel∗ Center of Applied Space Technology and Microgravity, ZARM, University of Bremen, Am Fallturm, 28359 Bremen, Germany Available online 28 October 2004
Abstract For full 3-axes attitude determination the magnetic field vector and the Sun vector can be used. A Coarse Sun Sensor consisting of six solar cells placed on each of the six outer surfaces of the satellite is used for Sun vector determination. This robust and low cost setup is sensitive to surrounding light sources as it sees the whole sky. To compensate for the largest error source, the Earth, an albedo model is developed. The total albedo light vector has contributions from the Earth surface which is illuminated by the Sun and visible from the satellite. Depending on the reflectivity of the Earth surface, the satellite’s position and the Sun’s position the albedo light changes. This cannot be calculated analytically and hence a numerical model is developed. For on-board computer use the Earth albedo model consisting of data tables is transferred into polynomial functions in order to save memory space. For an absolute worst case the attitude determination error can be held below 2◦ . In a nominal case it is better than 1◦ . © 2004 Elsevier Ltd. All rights reserved.
1. Introduction The attitude control system is an important satellite sub-system. On one hand it must achieve a certain accuracy. On the other hand it must be very robust to ensure the mission success. Furthermore, it shall be equipped with low-cost hardware. The solution is the utilization of simple, robust sensors in combination with intelligent algorithms. One very reliable and simple possibility is the combination of a magnetometer and a Coarse Sun Sensor
∗ Tel.: +49 (0) 421 218 3635; fax: +49 (0) 421 218 2526.
E-mail address: [email protected] (P. Appel). 0094-5765/$ - see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.actaastro.2004.09.001
(CSS). This will provide two vectors, and hence a full 3-axis attitude determination. The CSS consists of six solar cells that are placed on the surface on each side of the satellite, see Fig. 1. Using the cosine law, the total light vector, having contributions from both Sun light as well as Earth albedo, is obtained. In order to get the correct vector a measurement model of the albedo must be included. This paper covers the development of an attitude determination system using a CSS and a magnetometer. In the first part an albedo model which is applied for sensor simulation and for attitude determination is developed. In the second part an algorithm for attitude determination using an extended Kalman filter for state estimation is presented.
116
P. Appel / Acta Astronautica 56 (2005) 115 – 126
Fig. 1. Solar cell.
2. Albedo model The CSS used for Sun vector determination consists of six solar cells placed on each of the six outer surfaces of a cube shaped satellite. As the CSS, unlike a fine Sun sensor, sees the whole sky it is much more sensitive to other light sources than the Sun. The largest of these sources for a LEO is the light reflected by the Earth, the Earth albedo. To be able to get a good estimation of the Sun vector the Earth albedo has to be compensated for. This is made by developing an numerical albedo model and including it in the Kalman filter for attitude estimation. In order to calculate the amount of light received by a satellite in orbit a model for the light reflected by Earth is needed. The incoming light is a function of the satellite position, the Sun’s position and the reflectivity of Earth. The Earth is assumed to emit the radiation like a black body and the main part of the light is assumed to be reflected by Earth diffusely. The radiation from a black body is isotropic, that is, it has no preferred direction. For such a body the Lambertian cosine law is applicable. As the incoming sunlight hits the Earth it is diffusely reflected by the atmosphere and Earth surface in all directions. The amount of emitted light is determined by the reflectivity of the element and the angle of the incoming light. The amount of emitted light in the direction of the satellite is decided by the angle
Fig. 2. Description of the path of light geometry.
between the surface element normal and the satellite vector. The voltage delivered by the solar cell receiving the light is also calculated with the cosine law with the satellite vector and solar cell surface normal as components. The path of the light and the vectors are described in Fig. 2. The amount of incoming light to the satellite is equal to the integral of the three vector products over the illuminated and visible area of Earth. The total received energy is ˙ r = Is Acell Q visible and illuminated sphere area
(R e,0 • S 0 ) · (S 0 • n0 ) · (R e,0 • r 0 ) × dA, (1) S2 where Is is the solar intensity at the Earth, Acell the solar cell surface area, the reflectivity of each surface element of Earth, R e,0 the surface element normal, S 0 the vector from the surface element towards the satellite, n0 the solar cell surface normal, and r 0 the Sun vector. The solution of Eq. (1) is needed for every point in the orbit of a satellite. Because of its complexity with the non-trivial limits of integration it cannot be solved analytically. Together with the fact that the reflectivity of Earth is far from being uniform this implies that a numerical solution is the best way to solve the problem.
P. Appel / Acta Astronautica 56 (2005) 115 – 126
117
Fig. 3. Daily reflectivity of the Earth. Fig. 4. Latitude dependency of Earth’s reflectivity.
2.1. Reflectivity model First, a numerical model of the Earth is created by dividing a sphere into a finite number of surface elements. A reflectivity model is developed to describe the reflectivity of every single element. 2.2. Reflectivity data To get a realistic reflectivity distribution reflectivity data from existing measurements is analyzed. The data is taken from the NASA project Total Ozone Mapping Spectrometer [1]. When plotted the data shows Earth’s reflectivity as a function of the latitude and longitude. Fig. 3 shows the reflectivity reading from 1 day. The bright areas represent high reflectivity and the dark low. From the plot the Polar Regions and Greenland (lat. 170◦ , long. 150◦ ) can be recognized as well as several cloud formations. It can be seen that a lot of the reflected light origins from these cloud formations and the Polar Regions. From these plots a pattern could be recognized. It was obvious that the light reflected by Earth had a large dependency on the latitude, but only a small on longitude. Fig. 4 shows an example of the reflectivity of Earth as a function of its latitude. The solid curve is the average reflectivity and the dashed the standard deviation as a function of the latitude of Earth. The reflectivity has a maximum in the Polar Regions where the deviations are smallest. This
is due to the constant ice or snow coverage. For lower latitudes the reflectivity decreases and the deviations increase, respectively. The deviations are due to local differences in reflectivity of the surface but also due to tidal changes of, for example, cloud formations which have a high reflectivity. Knowing this, a statistical model of the reflectivity can be developed. From the data averages and standard deviations are derived for each latitude. The reflectivity model is created by randomizing new reflectivities based on these averages and deviations for each latitude. The data is merged in order to generate a model of the whole surface of the Earth. This reflectivity pattern is applied to the numerical model of Earth. 2.3. Received light As the light reflected by Earth is emitted diffusely the cosine law is applicable to calculate the received light by the satellite in any position above the Earth surface. It will suffice to calculate under which angle each surface element is illuminated by the Sun and from which angle it is seen by the satellite. The vector towards Sun is known and so is also the position of the satellite. For each element a normal is calculated as is also the vector pointing towards the satellite. The vector products between the Sun vector and the surface normal as well as between the surface normal and the satellite vector are calculated, respectively. Both
118
P. Appel / Acta Astronautica 56 (2005) 115 – 126
2.5. Recording data The above procedure is run through several times, each time with a new randomized reflectivity applied to the numerical Earth and for all satellite positions. The mean received albedo intensity of the three components is recorded as well as the standard deviations for every run. This results in three matrices describing the mean of the albedo light vector in the three components and three matrices describing their deviations. 2.6. Analyzing data
Fig. 5. Albedo light vector at satellite position.
products are multiplied with each other and the reflectivity of the element. This is the numerical solution to Eq. (1) and the solution can now be calculated for every position. 2.4. The grid A shell describing a grid of satellite locations at a certain altitude above the Earth surface is created. In every grid point the elements that are visible from the satellite and are illuminated by the Sun are picked out. For the elements surviving the test the emitted albedo light vector is calculated using the cosine law as above. Now the sum of the contributions from all elements is calculated for each satellite point in the grid, see Fig. 5, where Itotal =
n
Ik .
(2)
k=1
This total albedo light vector is split up into three components, north, east and down. This coordinate system is chosen to make the curve fitting to come easier, that is, to avoid sine- and cosine patterns from the transformation to occur in the plots. With knowledge of the attitude of the satellite having its solar cells in this position this would be the numerical representation of Eq. (1). The attitude is not taken into calculation at this moment.
Fig. 5 shows a result from these simulations. It describes the down component of the albedo light vector. The intensity on the z-axis is normalized w.r.t. the solar intensity. The zero longitude plotted is defined as the noon–midnight meridian. The zero latitude is at the north pole; hence the plots show the latitudes beginning from the shadow side of earth, sweeping over the north pole, the south pole and back to the shadow side. These six data sets are only valid for this altitude and for this Sun position. For other setups new data sets have to be developed. For a satellite that stays a long time in orbit and/or has a eccentric orbit this means a lot of data sets and so a lot of memory capacity needed. Because of this the data sets were transferred into functions instead of being saved as look-up tables. Functions are very memory saving as they only need to include the functions themselves and a corresponding parameter vector for the current setup. 2.7. Polynomial derivation From the six three-dimensional plots an attempt to fit polynomials is made. The polynomials are used to describe the albedo intensity and the standard deviation, received by the satellite for an arbitrary position in the orbit. The mean parts of the polynomials are used in the Kalman filter as expected measurements as they deliver the solution to Eq. (1) for all positions. All polynomials together are used to simulate measurements when testing the filter as described later. 2.8. Curve fitting The three-dimensional plots describing the incoming light are fitted with two-dimensional polynomials
P. Appel / Acta Astronautica 56 (2005) 115 – 126
Intensity [percent]
Intensity [percent]
0.2 0.15 0.1 0.05 0 200 150 100
150
50 Longitude [degrees]
0 0 -100 -50
200
250
300
0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 200 150 100
100 50 Latitude [degrees]
50 Longitude [degrees]
Fig. 6. The data representing the down component.
0
-100
-50
0
50
100 150
200 250
300
Latitude [degrees]
Fig. 7. The fitted model representing the down component.
with the structure I = (Af (latitude) + B)(Cg(longitude) + D),
119
3. The sun sensor (3)
where I is the intensity received by the satellite , f (latitude) is a function describing the dependency on the latitude and g(longitude) a function describing the dependency on longitude. A, B, C and D are scalars that are used to fit the plots to the data. The polynomials for the mean values of the vector components are composed of three sine or cosine functions to describe the latitude dependency and one sine function to describe the longitude dependency. The polynomials for the mean values have the form y = (P1 (P5 sin(P6 ) + P7 sin(P8 ) + P9 cos(P10 + P13 )) + P2 )(P3 (P11 sin(P12 )) + P4 ), (4) where P1 –P13 are elements in a parameter vector, the latitude of the satellite and the longitude. P1 –P4 represent the scalars A–D in Eq. (4), the others represent amplitude, phase and phase displacement of the sine functions. The upper row describes the latitude dependency and the lower the longitude dependency. Each vector component has its own parameter vector. The result of the curve fitting can be seen in Figs. 6 and 7. The plot is the fitted model to Fig. 6. The parameter vector creating this plot is saved and later used in the Kalman filter.
The CSS used to measure the Sun vector consists of six solar cells placed on each of the outside surfaces of the satellite. Using the cosine law the amount of incoming light can be calculated for each cell and so the angle towards the Sun. 3.1. The solar cells The sensor head is developed at ZARM [2]. Fig. 1 shows such using two 20 × 20 mm solar cells. 3.2. Cell placement Fig. 8 describes one way to place the cells on a cubeshaped satellite. The cells should be placed where the risk of being shaded by other objects is small. The coordinate system used for determination of the direction of the solar cell surface normal is a body-fixed one. 3.3. Implementation The solar cell delivers a voltage that is proportional to the incoming light. The amplitude of the voltage is decided by the cosine law. With tests the characteristics of the cell were examined. Fig. 9 shows some results from these tests. The figure shows the resulting normalized current through a resistor connected to the cell as a function of the angle of the incoming light. The dark/lower curve represents the test results and the brighter/upper one
120
P. Appel / Acta Astronautica 56 (2005) 115 – 126
Fig. 10. Kalman filter work order.
Fig. 8. Solar cell orientation.
The principle of the filter setup used in this case is shown in Fig. 10, where the lower parts describe the model of the real case with an update part to make corrections when necessary. To work properly the Kalman filter needs information about the system dynamics and its partial derivatives with respect to the state vector [3]. The state vector is composed of the rotation rates and the quaternions describing the spacecraft attitude w.r.t. the inertial frame x(t) = (x , y , z , q1 , q2 , q3 , q4 )T ,
(5)
where is the rotation rate vector and q the quaternion vector. The time derivative of the state vector is also needed: Fig. 9. Solar cell characteristics.
is a cosine curve. The results show that the modelling of the delivered voltage from the cell can be done with a simple cosine function.
4. Kalman filter for attitude estimation
˙ x, ˙ y, ˙ z , q˙1 , q˙2 , q˙3 , q˙4 )T = f (x). x(t) ˙ = (
(6)
The state error covariance propagation part of the filter also needs the partial derivative of the state vector x w.r.t. itself. This is the F-matrix *f (x(t), ˆ t) F (x(t), ˆ t) = (7) *x(t) ˆ x(t)=x(t) ˆ
or An extended Kalman filter is used to estimate the state of the satellite. The filter must have knowledge about the dynamics of the systems. It needs a statistical description of the noises and errors in the system and measurements- and attitude dynamics models. It also uses any available information about the initial condition of the variables.
* ˙ F=
*f *x T
* = *q˙
*
˙ * *q . *q˙
*q
The resulting matrix F is a 7 × 7 element matrix.
(8)
P. Appel / Acta Astronautica 56 (2005) 115 – 126
(9)
in inertial frame, and eb# the vector normal to the solar cell surface. Also
0, f < 0, g{f } = (14) f, f 0.
(10)
Together they deliver the expected measurement vector from the magnetometer and the CSS to the filter
4.1. Attitude dynamics model The attitude dynamics equations of motion for a rigid body are ˙ = I −1 [T dist + T control − × (I )], q˙ = 21 q,
0 −3 = 2 −1
3 0 −1 −2
−2 1 0 −3
1 2 , 3 0
Y = h(x) = [uCSS , umag ]. (11)
where I is the inertia matrix of the spacecraft. The vector T dist includes all external disturbance torques acting on the satellite and the vector T control all control torques. The external torques taken into account are those from the gravity gradient, magnetic and aerodynamic drag and solar pressure. 4.2. Measurement model The most important inputs to the filter are the measurements from the magnetometer and the CSS. These have to be manipulated to be interpretable by the filter. The Kalman filter uses the magnetic field vector measured by the magnetometer and the scaled light intensity on each of the solar cells of the CSS for the attitude determination. These measurements have to be modelled for the Kalman filter to know what to expect. The measurement model for the magnetometer is umag = Abi (q)B i ,
The observation matrix also needed by the Kalman filter for estimation of the expected measurements is the measurement vector y = h(x) partially differentiated w.r.t. the state vector *hk (x(tk )) H k (xˆ k (−)) = . (16) *x T (tk ) x(tk )=xˆ k (−) 4.3. Filter simulation All data needed for the simulations are acquired from the program SATSIM [5]. It is a simulation tool developed by the Center of Applied Space Technology and Microgravity (ZARM) to provide a test bench for the development of small satellites. The environment- and dynamics simulation tool of the SATSIM program simulates the orbital- and attitude dynamics of the satellite. 5. Results 5.1. Scenario definition
uCSS# = g{eb# · Abi (q) · P isun } T
T
(15)
(12)
where Abi is the transformation matrix from inertial to body fixed frame, umag is the magnetometer measurement and B i is the magnetic field vector in inertial frame from the IGRF model [4]. The measurement model for the #:th solar cell of the Sun sensor is
+ g{eb# · Abi (q) · P iEarth }
121
(13)
where # = 1, 2, . . . , 6 and uCSS# is the intensity measurement, P isun the sunlight vector in inertial frame, P iEarth the albedo light vector from polynomials,
The orbit used in the simulations presented in this paper is a polar one with its normal perpendicular to the Sun vector. The Sun is assumed to be located in the equatorial plane. An orbit altitude of 630 km is used. The simulation is started at the equator on the illuminated side of Earth and runs for about three orbits. The satellite is of cubic form with a side length of 1 m and a mass of 300 kg. The moment of inertia matrix is fully occupied because of the displacement of the center of gravity with respect to the geometrical center of the satellite. This displacement is unknown to the Kalman filter. The center of gravity is displaced 0.1 m in the x-direction, 0.05 m in the y-direction and
122
P. Appel / Acta Astronautica 56 (2005) 115 – 126
0.04 m in the z-direction. The moment of inertia matrix is
53 I = 1.5 1.2
1.5 50.75 0.6
1.2 0.6 kg m2 . 50.48
(17)
5.5. Total torque The total torque N Total acting on the satellite is the sum of the torques originating from aerodynamic drag, solar torque, gravity gradient and magnetic torque. N Total = N Aero + N solar + N GG + N B .
5.2. Initial conditions To test the performance of the filter in the initial phase the state vector is assumed not to be known well. The satellite is assumed to have completed the detumbling phase and to have a total rotation rate of 0.01 rad/s in the simulations. This corresponds to about 10 revolutions per orbit. This quite fast rotation has large influence on the attitude determination if the moment of inertia matrix is poorly known. 5.3. Satellite environment To make the simulations more realistic some assumptions on forces and torques acting on the satellite and parameter errors that can occur have been made. In addition two parameter errors have been introduced; one of the moment of inertia matrix elements has been increased by 5%. This is a reasonable value as small satellites with limited budget may have uncertainties in the mass configuration of the spacecraft. Also the difference in propellant mass, if present and not modelled, over time affects the moment of inertia of the satellite. It can furthermore be assumed that the position of the satellite cannot be precisely determined. To simulate this a position error has been introduced. The Kalman filter is given a position that differs up to 1◦ in latitude and longitude from the true value. 5.4. Disturbances The external torques assumed influencing the satellite are those from the gravity gradient, magnetic- and aerodynamic drag and solar pressure [6]. They are all acquired from SATSIM. The only disturbance assumed known to the filter is the gravity gradient. That is, the filter is able to calculate the gravity gradient force and compensate for the errors it could produce.
(18)
It varies between 1×10−5 and 5×10−7 N m depending on the position and attitude of the satellite. 5.6. Filter tuning In the filter the predicted measurements computed from the predicted state are compared with the actual ones. This is done so that the correction of the state vector can be specified in order to minimize the prediction error. There are two major problems when trying to get a good estimate of the attitude; the first are the noises and disturbances in the measurements where the albedo deviation is the most important and in the system dynamics. The second is the lack of observability that occurs when there is a small angle between the two measured vectors. Both problems can be compensated for by changing the amplitude of the measurement covariance matrix. To optimize the Kalman filter performance it must be tuned. This is done by changing the amplitude of the elements in the covariance matrices. The following sections explain where it is necessary and how the changing of the covariance matrix amplitude can enhance the performance of the filter. 5.7. State covariance matrix During the simulations the state covariance matrix is kept constant. That is, the external forces are assumed not changing that much over time. This is of course not true, but the filter is tuned to be able to handle the largest disturbances. The differences in external disturbances can also be compensated for, but that is not done in this phase. In the state covariance matrix the three first elements in the diagonal are representing the rotation rates of the satellite and the four last the quaternions. The rotation rate values are smaller than those of the quaternions; they therefore are represented by smaller values in the matrix.
P. Appel / Acta Astronautica 56 (2005) 115 – 126
5.8. Measurement covariance matrix
123
The components of the measurement covariance matrix used are of different amplitude. The first six elements in the diagonal represent the solar cell measurements and the three last the magnetometer measurements. The noise in the light vector measurements is in order of ten percent. It is composed of the albedo model noise, and sensor noise. The assumed noise in the magnetic field measurements are in the order of one percent and as the covariance matrix is the square of the deviations the elements representing the solar cells are set to a 100 times larger value than those representing the magnetic field vector. This relationship has also been tried and verified with simulations.
the measurement covariance matrix (R-matrix) amplitude is set higher to make the filter rely on the system dynamics and therefore get a steadier and more correct estimation where the noises in the measurements are high. This is only true to a certain point. A greater measurement covariance says that the measurements are not that trustworthy and the filter should rely more on the equations of motions to estimate the state. But when external forces are present the filter will be unable to predict the movements and must use information from the measurements to not drift away. Using the polynomials for the standard deviations of the expected albedo radiation, the R-matrix can be changed continuously to always be adjusted to the measurement noise expected.
5.9. Initial state covariance
5.11. Observability issues
The initial state covariance matrix describes to what extent the initial state can be trusted. A small value indicates that the given parameter is a good approximation and vice versa. The elements are arranged as in the state covariance matrix. The initial state covariance has been set to zero for the rotation rates. This means that the filter is told that the given initial rotation rates are correct although they are not. The initial state covariance of the quaternions is set to 100. This combination makes the filter believe in the rotation rates and therefore it starts correcting the quaternions and makes initially only small attempts correcting the rotation rates. After the right attitude has been found the rotation rates slowly converge also. This can be seen in Fig. 13 which shows results from the same simulation as in Fig. 12. The rotation rates converge after about 100 s and the quaternions after less than 5 s. In Fig. 13 it can also be seen the greater sensitivity to external forces for the rotation rates up to 100 s than after. This is due to the low R-matrix amplitude set during the first 100 s. If the rotation rate covariance is set higher the filter may have problems converging as the filter tries to correct the rotation rates without the reference that a known attitude provides.
The Kalman filter uses the information from the sensors to create two vectors. Using these vectors the filter is able to determine the attitude of the satellite. This is only possible when the two vectors are both present and not pointing in the same direction. The following section explains where it occurs and how, if possible, it can be compensated for.
5.10. Steady state phase During the steady state phase a stable and accurate estimation of the state is wanted. After the initial phase
5.12. Parallel measurement vectors For a high inclined orbit whose normal is close to perpendicular to the Sun vector a problem concerning the observability can occur. This happens four times per orbit as the magnetic field vector is close to parallel to the Sun vector. In the simulations this is visible as larger deviations in the attitude determination than in a nominal case. 5.13. Sensor simulation To test the performance of the filter, several simulations have to be made. Measurements have to be simulated as input to the filter. The magnetic field vector is computed by a function based on the IGRF model. The model is distorted by white noise with a 200 nT amplitude [6] run through a low pass filter to simulate the local and tidal changes. Another white noise with amplitude 2 nT is added to simulate measurement noise. The magnetic field vector delivered to the filter is described in an inertial frame. The received light is
124
P. Appel / Acta Astronautica 56 (2005) 115 – 126
composed of one part coming from the sun, and another part that has been reflected by Earth. The latter part is created using the polynomials for the mean light vector as a base. A random number with a position dependant standard deviation is added. This number is derived from those polynomials of the albedo model describing the uncertainty of the reflected light. The result is put through a low pass filter to simulate the low frequencies at which the albedo measurements can be assumed to change. The light vectors are added and split up into six non-negative scalars representing the solar cells on all sides of the satellite. To these scalars a white, unfiltered noise with amplitude 0.5% of the measured intensity is added before they are given to the filter as measurements. This reflects the electronics noise. 5.14. Filter performance Fig. 11 shows some results from the absolute worst case simulations including the following error sources: • Disturbance torque unknown to Kalman filter with amplitude 10−7 –10−5 Nm. • Measurement noise, about 0.5% of measurement amplitude. • Albedo noise. • Uncertainty in moment of inertia, 5%. • Position error up to 1◦ in latitude and longitude. Picture A describes the intensity of the reflected light, normalized w.r.t. the solar intensity in the NED frame. The curve with the largest amplitude is the down component, the second greatest the north component and the weakest is the east component of the incoming vector. This intensity distribution is natural as the satellite in this case is passing right over the most radiating part of Earth, that is, the part receiving the most sunlight. The B plot shows the modelled standard deviations for the Earth’s albedo vector. They are also expressed in the NED-frame and are normalized w.r.t. the sunlight intensity. It can be seen that the largest deviations are present in the Down component. This is expected as this component is the largest. The deviations are up to 10% of the albedo radiation. The two other components are similar in size and partially also larger
Fig. 11. Albedo model simulation results.
than the mean values they correspond to. This has the consequence that the mean values are to be trusted poorly. The C plot shows the angular deviation between the estimated quaternion vector and the real one over time expressed in degrees. The coordinate system is the body-fixed one. The shaded stripes represent the satellite entering the shadow phase of its orbit. The D plot shows the observability of the filter. This is presented calculating the sine value of the angle between magnetic field vector and total light vector. A high value indicates good observability, that is, the two vectors used are close to being perpendicular to each other. A low value indicates poor observability. This happens when the angle between the vectors decreases. The result of this can be seen in most plots where the
P. Appel / Acta Astronautica 56 (2005) 115 – 126 x y z
50
x10-3
0
30
Rotation rates [radians/sec]
Angular deviation [degrees]
40
1
125
20 10 0 -10 -20 -30 -40
x y z
-1 -2 -3 -4 -5 -6 -7 -8
-50
-9
0
1
2
3
4
5 6 Time[s]
7
8
9
10
-10 0
50
100
150
Time [seconds]
Fig. 12. The Quaternions during the transient phase. Fig. 13. The rotation rates during the transient phase.
deviations are increasing at the same places where the observability shows a low value. From the C plot in every set some general conclusions can be made. The filter performance is best close to the Polar Regions where the deviations are small. At lower latitudes the albedo variance rises and the error rises with it. The peaks in the illuminated phases are a result of poor observability. In the Sun phase they coincide with the lower points in the observability plots. The shaded stripes represent the satellite entering Earth’s shadow. When the satellite enters the Earth shadow the filter can only use the magnetic field vector as measurement and it will therefore be more sensitive to dynamic disturbances. This is visible as larger deviations in this phase. The overall performance of the filter is good. The errors in this worst case scenario are not greater than 2◦ in the illuminated phase, and for a nominal case the errors can be kept below 1◦ . 5.15. Transient phase Figs. 12 and 13 show the transient phase of the same simulation described above. The filter converges quickly. The transient time is shortened down to under 10 s using smaller values for the measurement covariance matrix during this phase. The filter converges and is stable when started in the illuminated phase, and can be kept stable during the shadow phase.
5.16. Requirements To be able to perform correct attitude estimation the on-board computer must have access to the following inputs • satellite position (orbit propagator), • Sun position, • albedo model (polynomial function and parameter vector), • magnetic field model (IGRF), • satellite preferences (moment of inertia matrix, etc.). 6. Summary/Conclusions 6.1. Summary An Extended Kalman filter utilizing a CSS in combination with a model of the Earth’s albedo can drastically improve attitude determination. The absolute worst case accuracy is 2◦ and in the nominal case its better than 1◦ with the same disturbances. The accuracy of the direction determination of the Sun vector is good and the costs much lower than when using a fine Sun sensor. The sensors are also reliable because of their simple construction. As a complement or substitute to the small sensors the ordinary solar cells of a satellite can be used.
126
P. Appel / Acta Astronautica 56 (2005) 115 – 126
A combination of CSS and magnetometer is simple, robust and accurate enough for acquisition and safe modes. It is also sufficient for nominal modes of most small satellite applications. 6.2. Conclusion The tasks presented in this paper can be summarized in the following table: 1. Albedo model for CSS simulation and error correction developed. 2. Polynomial functions for on-board computer use developed. 3. Extended Kalman filter for attitude estimation developed. 4. 2◦ absolute worst case accuracy achieved. 1◦ nominal case accuracy.
References [1] Total ozone mapping spectrometer, http://toms.gsfc.nasa.gov/ [2] ZARM—University of Bremen, http://www.zarm.unibremen.de [3] M.S. Grewal, A.P. Andrews, Kalman Filtering—Theory and Practice, Prentice-Hall, Englewood Cliffs, NJ, 1993. [4] International Geomagnetic Reference Field, http://nssdc.gsfc. nasa.gov/space/model/models/igrf.html [5] S. Theil, M. Wiegand, H.J. Rath, SATSIM—A generic tool for attitude control system development, simulation and testing, 51st International Astronautical Congress, Rio de Janerio, Brazil, 2–6 Oct 2000, IAF-00-U.2.05. [6] M. Wiegand, Autonomous satellite navigation via kalman filtering of magnetometer data, 46th Congress of the International Astronautical Federation, Oslo, Norway, October 2–6, 1995, ST-95-W.1.04.