Attitude determination for nano-satellites – II. Dead reckoning with a multiplicative extended Kalman filter

This paper is the second part of a series of studies discussing a novel attitude determination method for nano-satellites. Our approach is based on the utilization of thermal imaging sensors to determine the direction of the Sun and the nadir with respect to the satellite with sub-degree accuracy. The proposed method is planned to be applied during the Cubesats Applied for MEasuring and LOcalising Transients (CAMELOT) mission aimed at detecting and localizing gamma-ray bursts with an efficiency and accuracy comparable to large gamma-ray space observatories. In our previous work we determined the spherical projection function of the MLX90640 infrasensors planned to be used for this purpose. We showed that with the known projection function the direction of the Sun can be located with an overall accuracy of ∼40′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 40^\prime$$\end{document}. In this paper we introduce a simulation model aimed at testing the applicability of our attitude determination approach. Its first part simulates the orbit and rotation of a satellite with arbitrary initial conditions while its second part applies our attitude determination algorithm which is based on a multiplicative extended Kalman filter. The simulated satellite is assumed to be equipped with a GPS system, MEMS gyroscopes and the infrasensors. These instruments provide the required data input for the Kalman filter. We demonstrate the applicability of our attitude determination algorithm by simulating the motion of a nano-satellite on Low Earth Orbit. Our results show that the attitude determination may have a 1σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} error of ∼30′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 30'$$\end{document} even with a large gyroscope drift during the orbital periods when the infrasensors provide both the direction of the Sun and the Earth (the nadir). This accuracy is an improvement on the point source detection accuracy of the infrasensors. However, the attitude determination error can get as high as 25∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\circ }$$\end{document} during periods when the Sun is occulted by the Earth. We show that following an occultation period the attitude information is immediately recovered by the Kalman filter once the Sun is observed again.


Introduction
Owing to the enormous funding requirements, satellite missions were only conducted by the economically most powerful countries of the world in the first few decades of the space era. However, as a consequence of the explosive technological development, small satellite missions with substantially lower funding requirement -for instance, nano-satellites including as CubeSats -became a viable alternative for traditional largesize and high-cost satellites, which made space an achievable goal also for countries/ organizations with less financial resources. The last few decades brought along a lot of such missions with more and more scientific aims being targeted by them.
One of the technological difficulties that needs to be handled in connection with small satellite missions is the accurate determination of the satellite's actual orientation, i.e. its attitude. While on large-size satellites this information is usually provided by costly, large-size star trackers, which determine the attitude based on the angular distribution of bright stars in their field of view, these systems do not fit the very restricted size and power budget criteria of nano-satellites. In our recent paper [7] we proposed a new, cost-efficient approach to this problem which is based on the utilization of thermal imaging infrasensors. For this purpose we chose the MLX90640 infrasensor of [10], which is a small-size, low-cost sensor having 32× 24 pixels and a relatively large, 110× 75 degree field of view. This coverage by a single sensor implies that six of these sensors, placed on the six sides of a cube, could cover the full sphere, see Fig. 4 of [7]. This technology might be suitable to even smaller satellites in similar missions like GRBAlpha [12]. As the spherical projection function of MLX90640 infrasensors (to be used for this purpose) is now known with an overall accuracy of ∼ 40 � [7] we now turn to the next step and introduce a simulation model for testing the applicability of our attitude determination approach. As we outline in our recent paper [4], our method is based on a multiplicative extended Kalman filter that uses the information provided by the infrasensors (direction of the Sun and the nadir in the satellite's coordinate frame), the GPS system (the location of the satellite in the Earth centered coordinate frame) and the MEMS gyroscopes (angular velocity of the satellite) carried by the satellite.
The layout of the paper is the following. In Section 2 we describe the simulation of the satellite dynamics and introduce the multiplicative extended Kalman filter method. We demonstrate our results in Section 3, and we summarize our conclusions in Section 4. A detailed description of the equations used for the Kalman filter can be found in Appendix 1.

Simulation model for testing the on-board attitude determination algorithm
The attitude determination algorithm we developed is aimed to run on-board and therefore it needs to be tested for the different situations possible during a space mission. For this purpose we built a simulation model where all parts of the attitude determination process can be tested independently and as a whole as well. The first part of the code simulates the dynamics of the Sun-Earth-satellite system while its second part determines the attitude of the satellite by applying a multiplicative extended Kalman filter to the simulated data provided by the first part of the code. The goal of this process is to see how the recovered attitudes compare to the 'real' ones.

Simulation of the satellite dynamics
This part of the code calculates the position and the attitude of the satellite in the Earth centered J2000 reference system (where the X and Z axes point towards the former positions of the vernal equinox and the Earth's rotation axis in January 1, 2000 at 12:00 TT, respectively) as well as the position of the Sun and the Earth in the satellite's coordinate system (where the origin of the system is fixed to the center of mass of the cubesat and the axes are parallel to its edges, see Fig. 1). In the code time is expressed in Julian dates and the GPS to J2000 coordinate transformations are implemented as well. The code also determines the 'night' part of the orbit, i.e. where the Sun is occulted by the Earth. The orbit of a satellite can be characterized by five orbital elements ( Ω , longitude of ascending node; i, inclination; , argument of periapsis; e, eccentricity; h, altitude of satellite orbit at perigeum; see Fig 1.) which remain constant when assuming a spherical Earth. The actual position of the satellite on this elliptical orbit is given via the true anomaly ( ) which is obtained from the Kepler equation. However, the oblateness of the Earth introduces perturbations, from which the J 2 perturbation has the largest magnitude, and therefore it has to be taken into account while simulating the satellite dynamics. The main effect of the J 2 perturbation is on Ω and . The time derivatives of these orbital elements are the following [8]: where R e is the radius of the Earth, p = a(1 − e 2 ) , n = √ GM e ∕a 3 and J 2 ≃ 1.082629 ⋅ 10 −3 .
Next we discuss the description of the rotation of the satellite. As the variation of the gravitational field is negligible in the range of the satellite's dimensions and we did not apply any kind of attitude control in our model, the net angular momentum transfer due to external torques is negligible during one orbit. We also assumed that the satellite frame coincides with the principal frame, however, in case it does not, an additional constant rotation needs to be taken into account between the two frames. The time evolution of the attitude then can be determined using Euler's rotation equation, which describes the evolution of the angular velocity of a rotating rigid body represented in its principal frame (see e.g. [2]): where is the angular velocity vector of the satellite represented in the principal frame (which coincides with the satellite frame), and I 1 , I 2 and I 3 are the moments of inertia corresponding to the x, y and z axes, respectively. Note that Eq. (3) could imply complex motions, such as tumbling 1 , which we can readily reproduce within our model (see supplementary material).
The attitude can then be calculated by an additional integration over the angular velocities. In this work we use unit quaternions to represent the attitude of the satellite as q s = [nsin( ∕2), cos( ∕2)] T , where n is the axis of the rotation that transforms the J2000 coordinate frame to the satellite frame, is the rotation angle and the superscript T denotes transposition (quantites represented in the satellite frame are denoted by the superscript 's'). The quaternion kinematics is given by the following equation (see e.g., [1,3]): where ( s ) is the matrix representation of the cross product ( s ×):

Attitude determination with Kalman filter
The Kalman filter is an algorithm that provides an efficient way to estimate the state of a dynamic system by a series of measurements with inaccuracies over time [6]. The estimates produced by the algorithm are more accurate than those based on a single measurement alone since the joint probability distribution of the variables is estimated for each discrete time-step of the process. This leads to the minimization of the mean of the squared error of the estimates. The Kalman filter works in a two-step process with a prediction step (time update) and a measurement step. In the prediction step the filter propagates the estimates of state and uncertainties from current to the next time step. In the measurement step the state of the system is measured with some error and the estimate is updated by the weighted average of the estimate and the measurement, where the weights are determined by the respective uncertainties.
In our specific objective the Kalman filter serves to combine the infrasensor data with the angular velocity information provided by the MEMS gyroscopes. Although the 3 axis MEMS gyroscopes yield an accurate attitude information on a short time period, due to the error accumulation effect known as gyroscope drift an absolute attitude information is required as well, which is provided by the thermal imaging infrasensors in our case. In earlier works this kind of absolute attitude information was usually gathered by a 3-axis magnetometer and an optical Sun sensor (see e.g., [1,5,11,14]). The use of infrasensors is more convenient in the sense that as opposed to magnetometers they can be built-in parts of the satellite and do not need an external boom to be mounted on. The infrasensors determine the direction of the Sun and the nadir in the satellite's coordinate frame, which vectors are known in the Earth centered reference frame as well owing to the location information provided by the onboard GPS. The rotation that transforms the reference frame to the satellite frame (which is indeed equivalent to the attitude of the satellite) is unequivocally determined by the pair of these two vectors in the two coordinate frames. Hence we get a prediction of the system's state from the gyroscope which can be corrected by the absolute attitude information provided by the infrasensors.
Since the quaternion kinematics, described by Eq. (4), is nonlinear in the variables s and q s the utilization of an extended Kalman filter is necessary. We use the Multiplicative Extended Kalman Filter (MEKF) method [9], where a multiplicative error state q is introduced, which represents a small rotation from the predicted attitude -which contains measurement errors -to the actual attitude (from now on we omit superscripts, since everything is understood to be represented in the satellite frame, unless noted otherwise) 2 : where the circumflex ' ∧ ' denotes the expected (or predicted) value of a quantity. With this multiplicative approach the conservation of the unit quaternion length is guaranteed and the problem of singular covariance matrices, encountered in the additive approach, is avoided as well.
To describe gyroscope measurements we use the model of [9], where in addition to a zero mean Gaussian error a time dependent bias vector is also introduced, the motion of which is determined by a random walk. Hence, the measured angular velocity m is given by where the Gaussian processes and have covariance matrices 2 3×3 and 2 3×3 , respectively. Therefore the estimated value of the angular momentum is The state space model is then given by the equations: The predicted values of the angular momentum and bias vectors are updated through the Kalman filter using the position of the Sun and the nadir both as seen by the satellite and as calculated in the inertial frame. A vector in the inertial frame is transformed to the satellite frame by It can then be shown that a small multiplicative quaternion error q creates a small r s deviation detemined by the following equation [3]: where q 3 is a three component vector containing the imaginary part of the multiplicative error state q . This determines the so called sensitivity matrix, which is then used by the Kalman filter to calculate the multiplicative error state from the r s quantities.
Since measurements are made at discrete time-steps, to implement these equations one must first discretize the kinematical equations. The Kalman filter then can be applied after each time-step to refine the attitude information predicted by the kinematical equations. The detailed discretized equations can be found in Appendix 1.

Simulation results
The applicability of our attitude determination approach is demonstrated by simulating a satellite on Low Earth Orbit (LEO) described by the following parameters: where L 0 denotes the initial angular momentum of the satellite in the satellite frame. The chosen I 1 , I 2 and I 3 values correspond to a 3U CubeSat with a size of 10x10x30 cm and total mass of 3.3 kg. The initial attitude was selected randomly. Figure 2 shows how the position and the attitude of the satellite changes during 6 hours on such an orbit. The attitude is represented in the form of right ascension ( ), declination ( ) and roll ( ). The conversion rule between these angles and the quaternion representation of the attitude is given by the following formulas: Since MEMS gyroscopes are available with various precision we investigated three different cases for attitude determination characterized by three different values for gyroscope drifts. For the largest error case we used = 4.89 × 10 −3 rad/s 1∕2 and = 3.14 × 10 −4 rad/s 3∕2 as proposed by [1], while for our standard and low-error case we used errors 0.1 and 0.3 times those of the high-error case, respectively. The initial parameters of the MEKF were chosen as follows: q 0 was selected randomly and the standard deviation of the measured vectors had been set to 0.012 rad ( ∼ 40 � , in accordance with our previous result on the pointing accuracy of MLX90640 infrasensors) and had been added to the input vectors. Sensor data were sampled at 1 Hz.
The recovery of the attitude on the orbit shown in Fig. 2 for our standard choice of gyroscope error is presented in Fig. 3. The left panels of Fig. 3 show that the attitude elements are well recovered when the MEKF works with two input vectors ('day'), i.e. when the infrasensors provide both the direction of the Sun and

Fig. 3
The real (orange) and the MEKF recovered (green) attitude for the simulation shown in Fig. 2 for our standard choice of gyroscope error (left), and the error of the attitude determination, i.e. the difference between the real and the recovered attitude elements for the same orbital configuration (right). Shaded areas represent parts of the orbit where the Sun is occulted by the Earth the Earth (the nadir), while the accuracy breaks down significantly when there is only one input vector (only the nadir direction) available for the MEKF ('night'). This behavior is not surprising, since we are lacking the minimum of two linearly independent vectors necessary to gain information about the absolute attitude of the satellite, and since the bias instability of MEMS gyroscopes is relatively high. The right panels of Fig. 3 show that the difference between the real and the recovered attitude elements may reach 25 • during the 'night' phase in our standard case.
Even though the accuracy of recovering the independent attitude parameters breaks down during 'night', the errors are correlated even in this case due to the information gained from observing the horizon. This is shown in the left panel of Fig. 4, where we plot the y and z components of the quaternion error states ( q y and q z ). As the information about the horizon determines the orientation of the satellite with respect to the orbital plane, the error of the quaternion component that describes the rotation within this plane ( q O ) does not increase during the 'night' periods. q O can be produced as a linear combination of q y and q z in our example. The right panel of Fig. 4 shows a short time period around a 'night' to 'day' transition and how the attitude information is immediately recovered once the Sun is visible again.
We also investigated the statistical behavior of the measurement errors. To do so we initialized our simulation with the same parameters except for the direction of angular momentum vector, which we picked randomly. By starting the simulation from several different initial conditions we collected statistical data about the first 'day' and 'night' phases. Figure 5 shows the distribution of the right ascension's measurement error for the 'day' case with different gyroscope precisions (the other two attitude parameters have similar error distributions). The results show that the recovery has a 1 error of ∼ 22 � in our standard case, while this error is ∼ 18 � and ∼ 32 � in the low-and high-error cases, respectively. This is an improvement on the ∼ 40 � error of the MLX infrasensor's point source detection accuracy [7], which shows the power of the MEKF method.
In Figure 6 we show the same errors for different parts of the 'night' phase. We divided the 'night' period to four equal-length segments to investigate the evolution of the errors and to avoid creating statistics from time periods with qualitatively different behavior. We see that the distribution of errors gets smeared as time progresses, and also with larger gyroscope errors. In the last segment of the higherror case the distribution is completely smeared so that not much information is retained about the real attitude. This is in accordance with the results of [1].

Summary
In the present paper we described a simulation model for testing our new concept aimed at determining the attitude of nano-satellites. The attitude was represented by unit quaternions and a MEKF approach was applied to estimate the most probable  6 Probability distributions of the right ascension's measurement error during 'night'. The left, middle and right columns correspond to the cases with low, standard and high gyroscope error, respectively, while the different rows represent different equal-length time segments with the top row being the first and the bottom row the last quarter of the 'night' phase state (attitude) of the system. In our model the prediction step of the Kalman filter utilizes gyroscope measurements while its measurement step is based on infrasensor measurements and GPS location information which provide the direction of the Sun and the nadir in the satellite and in the J2000 reference frames, respectively.
The results of our simulations show that an attitude accuracy of 22 ′ is achievable using combined measurements of the infrasensors and MEMS gyroscopes having a conservative drift. This is an improvement on the accuracy of point source detection with the MLX infrasensors [ ∼ 40 � , see 7]. This accuracy is gradually lost when the Sun is occulted by the Earth whereupon it can reach values of ∼ 15 − 25 • . The attitude information, however, is recovered within a short time once the Sun is observed again.
During the actual mission the satellites will not have infrasensors on all of their six sides and hence not being able to observe the Sun will be more regular. However, these time periods will be relatively short and the gyroscope drift is expected to be manageable during these intervals.
In this work we simulated a satellite on LEO with an inclination of 60 • . This is a reasonable choice for a particle detector experiment like a GRB detector because on this orbit the satellite evades high latitudes with increased noise contamination from the polar regions but is able to cover a large area of the sky. However, on such an orbit the illumination conditions may change substantially on the timescales of a few months due to the motion of the Earth around the Sun, as well as due to the orbital precession caused by J 2 . However, we consider our simulations to represent the average conditions on such an orbit sufficiently well.
The attitude determination method described in this paper is planned to be used in the CAMELOT mission where the attitude data will also serve as additional information for localizing gamma-ray bursts besides triangulation. An in-orbit demonstration of our experiment is planned to be scheduled for the end of 2022.
The covariance matrices are propagated using with being the state transition matrix: where The and matrices determining the process noise matrix are given by

Measurement update
In the measurement update step the MEKF first estimates the quaternion error state using the sensitivity matrix determined by Eq. (14) and then updates the attitude utilizing Eq. (7).
Supposing there are n vectors measured by the satellite, the quaternion error state and the bias vector error can be obtained using the following formula: (20) , where s i denotes a vector measured by the satellite, while ̂ s i = (̂ − ) i i is the predicted value of that vector. k is the Kalman gain defined the usual way: with k being the sensitivity matrix: and k the measurement covariance matrix: The quaternion state, the bias vector and the covariance matrix are then updated by where the quaternion error state is obtained from its imaginary part using the normalization constraint: In our setup the number of measured vectors is n = 2 when the Sun and the horizon is visible at the same time, while it is n = 1 when the Sun is occulted by the Earth.
Development and Innovation Fund together with the European Union. This work has received financial support of the Hungarian National Research, Development and Innovation Office NKFIH Grant K-138962.
Funding Open access funding provided by ELKH Wigner Research Centre for Physics.

Conflicts of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.