Three-axis high-accuracy spacecraft attitude estimation via sequential extended Kalman filtering of single-axis magnetometer measurements

Magnetometer is a highly advantageous sensor for determining a spacecraft’s attitude. This article provides a solution to the problem of spacecraft attitude estimation using magnetometer measurements only. To ensure full observability of spacecraft attitude states, it is necessary to use at least two types of sensors. Consequently, utilizing a single sensor, such as the magnetometer, poses a significant challenge for any attitude estimation algorithm, including the extended Kalman filter (EKF). Moreover, implementing the EKF algorithm, or any other attitude estimation algorithm, is computationally intensive. To address these issues, an algorithm has been developed that estimates spacecraft attitude angles and attitude rates using a sequential extended Kalman filter (SEKF). This algorithm offers numerous benefits over those found in the literature such as high accuracy, low computational resource requirements, the ability to converge even with large initial attitude and angular velocity estimation errors, and the ability to function even if two of the three measurement channels of the magnetometer are not functioning. With these benefits, the developed SEKF algorithm is capable of operating in all spacecraft operational modes, delivering accurate performance and computation time. In spite of measurements with large noise values, the high accuracy achieved by the SEKF algorithm enables the magnetometer to serve as the sole source of attitude information, even if one or two magnetometer measurement channels are not functioning.


Introduction
The estimation process of spacecraft attitude relies on the measurement of specific physical quantities within the spacecraft body axes, such as the Earth's magnetic field. The most commonly utilized sensors for this purpose are magnetometers, with almost all spacecraft attitude and orbit control systems (AOCS) integrating a three-axis magnetometer as the primary sensor. Since magnetometers have no moving parts, they typically have a long lifetime and higher reliability compared to other sensors. Additionally, magnetometers are relatively inexpensive, which further contributes to their widespread application in AOCS systems. These advantages are significant reasons of why magnetometers are preferred and commonly used as sensors in AOCS systems.
Magnetometers can be used in conjunction with other spacecraft attitude sensors, as in the case of single frame attitude determination methods [1]. However, such an approach introduces some difficulties. For instance, during the eclipse period at each revolution of the spacecraft around the Earth, sun sensor measurements are unavailable. Gyros are also unable to run for long periods due to the accumulation of errors. Star sensors, which typically provide high accuracy attitude estimates, consume a significant amount of electrical power. Since the available power is usually limited, star sensors may not be the ideal solution for prolonged operational time periods. In such scenarios, the only sensor that may be available is the magnetometer. These elongated operational time periods often occur during the spacecraft's standby mode, which is the phase when the spacecraft spends most of its lifetime. When the spacecraft is in standby operation mode, it must maintain its attitude angles at a coarse level to ensure that it can switch to the high accuracy operational mode at any time. Magnetometer is an ideal sensor for such objective. As such, in addition to the previously mentioned benefits of magnetometers, we can also identify the following advantages: 1. Long operational periods 2. Non-intermittent measurements 3. Low power requirements for operation.
With all of these benefits in mind, a considerable amount of research papers and theses focused on the topic of spacecraft attitude estimation using magnetometer measurements only. Unfortunately, the proposed solutions suffer from several challenges. This is due to several reasons. The first reason is that the spacecraft attitude determination problem becomes fully observable only if the spacecraft is equipped with two or more types of sensors. Thus, utilizing single type of sensors, such as a magnetometer, is considered a challenge because the system loses its full observability conditions. The system observability problem becomes more severe as one of the three magnetometer channels becomes malfunctioning. In addition, observability reaches its worst condition if two magnetometer channels become malfunctioning. The second reason for the encountered difficulties is due to the selection process of system representing states. Most of the attitude representing states found in the literature suffer from singularity at certain spacecraft attitude angles. This inhibits the general usage of many developed algorithms in any arbitrary spacecraft attitude angle and may result in a total algorithm blowup. The third reason is due to the assumption of any small angle. This assumption limits greatly the capability of any developed algorithm to only small angles. Thus, the estimator could diverge if the initial attitude estimation error becomes large. The fourth reason is that the computational load of any developed algorithm should be kept minimum as much as possible in order to be feasible for real-time application.
As an example of the aforementioned difficulties, in Ref. [2], a three-axis magnetometer (TAM) is employed to estimate both spacecraft orbit and attitude. However, the developed algorithms were only capable of handling initial attitude estimation errors of up to 16°. In Ref. [3], spacecraft attitude was estimated solely based on magnetometer measurements. Although the developed algorithms were able to address initial attitude estimation errors of up to 45°, they had limitations in this regard. In Ref. [4], a number of spacecraft attitude estimation models were described. In Ref. [5], there were certain attitude determination methods that are not commonly employed due to their high computational demands. Such demands can often restrict or even preclude the use of these algorithms in real-time conditions aboard the spacecraft. It is worth noting that most commercially operational spacecraft still use processors that are 15-20 years old [6]. This is because the procedure of qualifying a processor for use in space is lengthy and requires multiple missions to be flown with the processor onboard. Algorithms such as those described in Ref. [7] may not be easily implemented on a 20-year-old processor. Similarly, the standard Kalman filter (KF) or extended Kalman filter (EKF) may not be the optimal choice, since they involve matrix inversion operations that are computationally expensive for the spacecraft's onboard computer (OBC). In Ref. [8], various methods for estimating the attitude of gyro-less spacecraft using sun sensors and magnetometers were reviewed. In Ref. [9], the authors presented algorithms for spacecraft attitude estimation based on magnetometer and sun sensor measurements. However, the maximum initial estimation error that these algorithms could handle is only 60°. In Ref. [10], the authors utilized TAM to estimate spacecraft attitude, but the filter had a steady state error of approximately 14°in the yaw axis only, which is relatively high. In Ref. [11], the authors employed a magnetometer for single-axis attitude estimation, making it impossible to obtain three-axis attitude estimates. In Ref. [12], an algorithm to provide linear estimates of spacecraft attitude was developed based on magnetometer-only with a steady-state error of 3°, which is considered a high value. In Ref. [13], attitude estimates of 0.3°were achieved using magnetometer measurements via the unscented Kalman filter (UKF). However, it should be noted that the computational load of UKF is much higher than that of EKF due to the need to propagate a number of sigma points. The UKF algorithm developed in Ref. [22], suffered also from singularity problems due to selection of spacecraft attitude representing states.
The primary aim and contribution of this research article is to propose a spacecraft attitude estimation algorithm utilizing only magnetometer measurements. The proposed algorithm could deal with all of the aforementioned challenges, existing in the above literature survey, simultaneously, and effectively. The algorithm introduced in this research utilizes sequential extended Kalman filter (SEKF), which is designed to tackle the matrix inversion problem associated with standard KF or EKF. By doing so, the algorithm enhances computational speed and performance in real-time applications. Moreover, the proposed algorithm is compared against a verified EKF algorithm that was implemented to estimate the attitude of the former Egyptian Satellite EGYPTSAT-1, ensuring the algorithm's effectiveness.
The organization of the present research article is as follows: in Sect. 2, nonlinear spacecraft motion models are introduced. In Sect. 3, the necessary estimation algorithms for EKF and SEKF are developed. Simulation results and discussions are presented in Sect. 4. Finally, in Sect. 5, conclusions are summarized, and future work is proposed.

Nonlinear spacecraft motion models
To proceed further, some definitions of the coordinate systems utilized need to be given. The Earth Centered Inertial (ECI) coordinate system has its X-axis pointing to the Vernal Equinox, Z-axis pointing to the direction of Earth's rotation, and the Y -axis completing the right hand rule. There is also a reference coordinate system. The X-axis of the reference coordinate system points in the direction of spacecraft velocity around the Earth, the Z-axis points to Nadir, and the Y -axis completes the right hand rule. The spacecraft body coordinate system (BCS) coincides with the reference coordinate system when all of the spacecraft attitude angles are nullified.
Cowell's formulation finds utility in explicating the model of a spacecraft's translational motion [21]. The formulation can be represented as follows [4]: In ECI coordinate system, the spacecraft position vector is denoted by R I , and the double dot over it represents the second time derivative of it. This vector specifies the position of the spacecraft with respect to the center of the Earth. The gravitational constant of the Earth is represented by μ E , and a a is the perturbing acceleration arising from the non-spherical shape of the Earth (see [15] for its calculation details).
The relationship that describes the motion of a spacecraft's attitude can be explained by the Kinematic equations given as [16,23] where the dot over the symbol represents the first time derivative and the transformation of a vector from the ECI to BCS is carried out using a quaternion vector, denoted by , which comprises of a real part, q 4 , and imaginary parts, q 1 ,q 2 , and q 3 . The relation which defines the matrix, denoted as , can be expressed as follows where ω = ω x ω y ω z T is the vector of spacecraft inertial angular velocities. The spacecraft rotational motion dynamics is represented by the relation [16] where T g , T a , T r , T s , are torques due to gravity gradient, aerodynamics, residual magnetic, solar pressure respectively, and, H ω is the angular momentum vector of any moving parts, w k is a zero-mean white Gaussian noise associated with the process noise. Details of computing disturbance torques could be found in Ref. [15].
[β×] is the cross product equivalent matrix for the vector β = β x β y β z T given as and the spacecraft inertia matrix, J , is expressed as Choosing the state vector that represents the attitude motion of a spacecraft is extremely important, and can be viewed as a critical factor for any attitude estimation algorithm. Therefore, selecting the appropriate spacecraft attitude motion state vector as follows is of utmost significance. Thus, the state vector is selected as The state vector chosen in Eq. (7) possesses the advantageous quality of being capable of representing all attitude maneuvers without encountering any singularity, thereby preventing the entire algorithm from collapsing due to singularities. Furthermore, it facilitates the representation of maneuvers that involve significant angles. Additionally, the selection of the orbital state vector, as represented by Eq. (1), which consists of the inertial position and inertial velocity, prevents any singularities from arising in the simulation for certain types of orbits.

Spacecraft attitude estimation algorithms
The primary responsibility of any algorithm for estimating spacecraft attitude is to estimate the spacecraft's attitude as accurately as possible, using a group of measurements that may contain errors. Algorithms which rely on singular value decomposition (SVD) are not desirable because of their high computational cost [5], especially considering the limited computational resources available onboard a spacecraft while it is in orbit. Therefore, algorithms that are based on extended Kalman filter (EKF) have been widely used due to their versatility in various applications. By leveraging these EKF-based algorithms, an efficient real-time algorithm can be developed using the following steps:

Extended Kalman filter
For nonlinear systems which include nonlinearities in either of the plant or the measurement processes, the extended Kalman filter (EKF) is utilized [2]. The state projection ahead of time is accomplished through the relation where f is a nonlinear function describing system behavior.
In the case at hand, f , is computed by concatenating Eqs. (2), and (4). Thus, A priori estimation error covariance is computed by the relation and the measurement update stage is given by whereX (+) k , P (+) k , and K k represent a posteriori state estimate, estimation error covariance, and EKF gain matrix, respectively. Due to its numerical stability, the Joseph form of the estimation error covariance matrix is used in Eq. (13). The state transition and the measurement matrices could be evaluated from the relations with T is defined as the sampling period, v k is a zeromean white Gaussian noise associated with magnetometer measurements. In addition, the measurements in the current article are provided by a magnetometer. Thus, the measurement vector, h k (X k ),is given by where B B , is the Earth's magnetic field vector measured in BCS, B I , is Earth's magnetic field vector expressed in the ECI, and D I →B is the transformation matrix from ECI to BCS. The transformation matrix, D I →B , could be computed from [16] Taking into consideration that the matrix representing the measurement noise covariance in its discrete form, R k , is related to its continuous form, R(t), through the relation [18] and similarly, the matrix representing the process noise covariance in its discrete form, Q k , is related to its continuous form, Q(t), through the relation [19] where η is an intermediate variable representing time?

Sequential extended Kalman filter
In certain cases, the conventional extended Kalman filter (EKF) described in Sect. 3.1 is referred to as the "batch EKF" because it incorporates the measurements at each time step as a batch. However, embedded systems may encounter issues with the batch Kalman filter due to the requirement of inverting a matrix with dimensions equal to the dimension of the measurement noise covariance matrix. To address this challenge, a solution is proposed in the form of the SEKF as outlined in Ref. [17]. With SEKF, the measurement vector is comprised of individual measurement elements, which are introduced one at a time to mitigate the aforementioned problem. Thus, where Z i, k : is the measurement element number i (i = 1, 2, . . . , M), and M is the measurement vector length. H i, k is the ith row of the measurement matrix, H k . Filter representing equations are now given as follows: The above formulation of SKF assumes that the measurement matrix, R k , is diagonal.

Observability
Floquet theory can be utilized to assure observability of the discrete filter system through computing the relation [3], All of the eigenvalues of the matrix, , should have a magnitude less than unity in the filter steady-state region in order for the estimator to be stable. The smallness of the maximum eigenvalue reflects the rate of convergence.

Simulation results and discussions
The developed SEKF in the current study has been compared to a previously verified EKF algorithm developed for the Egyptian spacecraft EGYPTSAT-1, which was launched in April of 2007. Figure 1 shows EGYPTSAT-1 spacecraft and its layout which was owned by the National Authority for Remote Sensing and Space Sciences. EGYPTSAT-1 was the first Egyptian Remote sensing low Earth orbit spacecraft.
Parameters of the EGYPTSAT-1 spacecraft are drawn from the design simulations during different spacecraft design phases, and were given in Ref. [14]. In addition, a verified EKF algorithm was developed and tested in Refs. [4,14]. The spacecraft had a semi-major axis of 7,039,200 m, eccentricity of 0, orbit inclination angle of 98.085°, argument of perigee angle of 69°, true anomaly angle of 0°, and right ascension of ascending node of 337.5°. The epoch time for the spacecraft was 17 April 2007 at 00 h:00 m:00 s. The actual initial attitude angles of the spacecraft were − 165°, 170°, and 85°for the yaw, roll, and pitch angles, respectively. Additionally, the initial angular rates for the attitude of the spacecraft were 0.7°per second, 0.8°/s, and − 0.2°/s for the yaw, roll, and pitch motions, respectively. Both EKF and SEKF algorithms commence with an initial condition of zero. This means that no initial information regarding the spacecraft's attitude is known, which is typically the case when the spacecraft is in detumbling mode immediately after launch. The spacecraft's inertial characteristics are defined by six values: J x , J y , J z , J xy , J xz , and J yz . The numerical values for each of these characteristics are as follows: J x equals 11.2 kg m 2 , J y equals 11.4 kg m 2 , J z equals 9.2 kg m 2 , J xy equals 0.02 kg m 2 , J xz equals negative 0.08 kg m 2 , and J yz equals 0.2 kg m 2 . The time step for the algorithm is 4 s. The magnetometer error standard deviation is 200 nT for each axis of the three axes of the magnetometer and this is considered a very high value. Additionally, the maximum permitted angular error is about ± 7°for low accuracy operational modes and ± 0.5°for high accuracy operational modes.

Case 1 (measurements of three-axis magnetometer)
In the first scenario, we are examining measurements from a three-axis magnetometer. The time history of attitude estimation errors can be seen in Fig. 2  which reaches 170°and 85°that are typically encountered during the detumbling mode. As a result, the algorithms that have been developed are versatile and can be implemented during any operational mode of a spacecraft, irrespective of whether it is a high accuracy or low accuracy mode.
Calculations are performed on a personal computer equipped with a Core(TM) i7 CPU 2.6 GHz and 8 GB RAM, using the MATLAB R2022b environment. Table 1 demonstrates the execution time and accuracy of the SEKF algorithm that has been developed. We can observe from Table 1 that both EKF and SEKF have a remarkably low error standard deviation in comparison to the maximum permitted angular error, even though the magnetometer measurement error is considerably larger than that used in other research articles (which typically utilized a magnetometer standard deviation of 50 nT or a maximum of 100 nT). Based on the achieved accuracy levels, the magnetometer can serve as the sole source of required attitude information during all spacecraft operational modes, without the need for any other sensor. Table 1 also shows that the EKF and the developed SEKF have comparable levels of accuracy. The average execution time depicted in Table 1 has nearly the same values, but the maximum execution time is extremely different. The maximum execution time is considered of utmost importance for real-time application on-board spacecraft. This is due to the fact that at each calculation step, the processor puts a hard time limit to execute any function so as to prevent processor hanging up. Thus, if this limit is exceeded, the processor terminates the execution and may shutdown or restart all of the spacecraft subsystems. This could represent a serious problem for the spacecraft if it has not enough electric power to restart, and may lead to whole mission failure. The maximum execution time of the EKF is approximately twice that of the SEKF. Thus, the proposed SEKF has a significant improvement.

Fig. 4 SEKF maximum eigenvalue time history
A study published in Ref. [20] detailed the development of a filter called the αβ filter and its comparison to the wellestablished EKF. The performance of these two filters was tested on a personal computer (PC), where the αβ filter demonstrated a faster execution time of 0.39 s, while the EKF's execution time was measured at 0.74 s. By comparing the current measured value of 0.0537 s to the values reported in the previous study [20], it is possible to notice the improved performance of the developed algorithm in the current research article, considering the variations in hardware and software configurations. Furthermore, it is worth noting that the current research article does not make any small angle approximations, which enabled the SEKF to effectively handle large initial attitude estimation errors and converge towards accurate results. Figure 4 displays the time history of the maximum eigenvalue magnitude of the matrix . As observed in the figure, the estimator achieves its steady state mode after a single orbital period. Subsequently, the maximum eigenvalue magnitude of the matrix, , decreases below unity and continues to decrease with the passage of time during SEKF operation. This trend signifies a superior convergence rate, and assures system observability despite using a single sensor for spacecraft attitude estimation. Based on the obtained results in Figs. 2, 3 and 4, and Table 1 The developed SEKF has the following contribution compared to the aforementioned literature: 1. Convergence despite large initial attitude estimation error. 2. Lower levels of attitude estimation errors (which reaches a standard deviation of 0.14°) compared to the maximum permissible error which is predefined at (± 0.5°). 3. Lower execution time compared to EKF, with a maximum execution time of nearly half that of EKF.

Case 2 (measurements of a two-axis magnetometer)
For Case 2, we use the same parameters for the spacecraft as in Case 1. We assume that the z-channel of the magnetometer is malfunctioning, which means that we have only measurements available for the x and y components of the Earth's magnetic field in BCS. Figure 5 displays the time history of the attitude estimation error in the case of a malfunctioning z-channel. As depicted in Fig. 5, the estimation error of the EKF, and SEKF are superimposed. This indicates identical accuracy performance for both filters. Moreover, Both EKF and SEKF reach steady state performance after 0.6 orbits, which is a bit longer compared to Case 1. This is because we have only two measurements available at a time instead of three. As a result, both EKF and SEKF require more time and more measurements to converge. Table 2 presents the primary performance indicators for both EKF and SEKF in the scenario where the magnetometer's z-channel malfunctions. In case 2, the steady state errors of the estimators are slightly greater than those in case 1. This is attributed to the fact that there are only two measurements available at any given time, as opposed to three. Furthermore, despite the SEKF having lower computational requirements, the average execution time of the SEKF is slightly less than that of the EKF. However, the maximum execution time of the EKF is 10 times higher than that of SEKF. This indicates superior performance for the SEKF. Consequently, in Case 2 another contribution is clear in addition to the contributions obtained in Case 1. The contribution is the ability of the developed SEKF algorithm to provide the desired accuracy levels even with a malfunctioning measurement channel.

Case 3 (measurements of a single-axis magnetometer)
For Case 3, we used the same parameters for the spacecraft as in Case 1. However, in this scenario, we assumed that the y and z-channels of the magnetometer are malfunctioning, which means we could only obtain measurements of the xcomponent of the Earth's magnetic field in BCS. Figure 6 depicts the time history of attitude estimation error when the y and z-channels are malfunctioning. As depicted in Fig. 6, the estimation error of the EKF, and SEKF are superimposed. This indicates again an identical accuracy performance for both filters. Both EKF and SEKF achieve steady state performance after three orbits, which is longer than in the case of the first and second scenarios. This can be attributed to the fact that only one measurement is available at a time, as opposed to three or two measurements. As a result, the estimators require more time and additional measurements from the functioning channel to converge. Table 3 presents the key performance indices for both EKF and SEKF in the case of malfunctioning y and z-channels of the magnetometer. The steady state errors of the estimators in the third scenario are higher than those in the first and second scenarios, which is again due to the fact that only one measurement is available at a time instead of three. It is worth noting that the steady state error falls well below the predefined range of ± 7°in the low accuracy operational modes. Thus, the developed algorithm can operate efficiently even when two of the three magnetometer channels malfunction. Additionally, the average execution run time of the SEKF is identical to that of the EKF because, in the third scenario, the measurement is a scalar quantity. As a result, both algorithms require the same computational resources. Thus, the added contribution in Case 3 is the ability of the developed SEKF algorithm to converge and provide the required standby mode accuracy levels despite having only one working measurement channel.

Conclusion and future work
In the current article, an algorithm that estimated spacecraft attitude angles and angular velocities was developed using only a three-axis magnetometer. The developed algorithm used SEKF formulation and was compared to a verified EKF. The SEKF algorithm was just as accurate as EKF but had a faster execution time. It provided attitude accuracy of at least 0.14°in each axis, meaning it could be used without any other sensors. No small angle approximation was utilized, and the quaternion vector was used to represent spacecraft attitude, which avoided singularity problems at large spacecraft attitude angles. The SEKF algorithm could converge even with initial attitude estimation errors of up to 170°. SEKF performance was assessed based on Floquet theory, and it showed good observability and convergence rate. The developed algorithm could still provide high accuracy attitude estimates even if the z-channel of the magnetometer malfunctioned. In case of malfunctioning y and z measurement channels, it could still provide suitable attitude estimates for low accuracy spacecraft operational modes. A future work is planned to apply the developed algorithms with in hardware in the loop simulator (HIL). Afterwards, applying the developed algorithms on a real spacecraft mission as reserve algorithms, and then as a main standalone algorithm.