CubeSat attitude determination with decomposed Kalman filters

CubeSats are the cost-effective entry to space research and applications. As mission requirements increase to carry out more complex tasks, the constraints on the satellite challenge how attitude estimation and control systems are designed. Limited energy, sensors, and computational capacity require compromises. In this paper, we propose a Kalman filter architecture to reduce the computational cost of attitude estimation, leveraging the conditional independence structure of its physical model. Our method decomposes attitude dynamics and kinematics, leading to a linear attitude quaternion and a nonlinear angular velocity filter. As accommodating all vector measurements would require a nonlinear filter, we propose the virtual sensor paradigm that transforms the nonlinear observation model into a linear one, without relying on approximations. Our numerical experiments showcase superior error dynamics and robustness to epistemic uncertainty compared to a nonlinear quaternionic filter, and we also investigate performance against star tracker measurement frequency and sensitivity to the angle between Sun and Earth magnetic field measurements.

This work focuses on orientation estimation with Kalman Filters (KFs), motivated by aerospace applications [2][3][4][5][6][7] and exploits the physical model of a small satellite to improve sensor fusion. We will address the need for estimators with stationary error dynamics, which is essential to prepare CubeSat s for high-accuracy Earth observation missions.
We contribute to the development of efficient orientation estimation algorithms. Our goal is to reduce computational complexity while maintaining attitude stability. First, we use the conditional independence of the physical model to decompose attitude kinematics and dynamics estimation; thus, decreasing computational cost. Second, we introduce the virtual sensor paradigm to transform the nonlinear observation model of the attitude filter into a linear one. This way, a linear attitude filter can be used, further reducing the computational cost. We compare the proposed algorithm to a quaternionic filter [8] and show its robustness and superior error dynamics.

Related works
Wahba's problem [9] describes the general attitude estimation objective as where n stands for the number of observations, w i express measurement confidence, b i denote vectors in the measurement (body), whereas r i those in the reference frame. A r,b maps the attitude from the reference to the body frame. Here the focus is on quaternion methods due to their computational efficiency, compare to, e.g., Euler angles [10].
The eigenproblem formulation of (1) offers diverse solutions by trading off efficiency and robustness [9,11]. Among which, the Estimator of the Optimal Quaternion (ESOQ) [12] and Second Estimator of the Optimal Quaternion (ESOQ2) [13] are the most efficient [14,15]. Moreover, ESOQ2 can be solved in closed-form for two observations [13], which we exploit in Sect. 3.2.
Whereas the Marginalized Square Root Quadrature Kalman Filter (MSQKF) [21] uses a linear KF to marginalize out linear states; nonlinear measurements make this filter inapplicable for orientation estimation (cf. Sect. 3.1). As our goal is to reduce computational cost, we do not consider using KFs for Riemannian manifolds [8], although they can accommodate quaternions.
Since attitude kinematics is time-variant (cf. Sect. 3.1), independence does not hold-as required for the Generalized Compressed Kalman Filter (GCKF) [22], so we rely on conditional independence (cf. Sect. 3.2). Our contributions are: 1. We utilize conditional independences in the spacecraft's kinematic and dynamic models to decompose the estimation problem into two filters; 2. We transform the nonlinear observation model of the attitude filter into a linear one, reducing the computational cost of the algorithm; and 3. We provide empirical evidence based on numerical simulations to show the superiority of our method compared to the Quaternionic Additive Square Root Unscented Kalman Filter (QuAdSRUKF) [8], especially in terms of error dynamics and robustness to epistemic uncertainty.

Proposal
This section describes the physical spacecraft model, then introduces and analyzes a decomposed estimation method, exploiting the equations of angular motion. To accommodate nonlinear vector measurements in a linear KF, we propose the virtual sensor paradigm and reason about its reduced computational cost.

Equations of angular motion
The equations of angular motion describe attitude and angular rate evolution over time. We parametrize attitude with a quaternion due to its advantages (no singularity, computational efficiency) [10]. A quaternion q is an element of the quaternion group ℍ and is composed of a tuple q = (s, w) = cos 2 , t sin 2 , where s denotes the scalar, w the vector part [10]. An equivalent representation uses the half angle ∕2 and a rotation axis t ∈ ℝ 3 . In accordance with the literature [10], we assume unit quaternions.
The four coordinates of unit quaternions are linearly dependent (they form a Riemannian manifold) [8]. Since KFs are designed for Euclidean systems, their additive updates will produce off-manifold samples for quaternions, since Riemannian manifolds are not closed w.r.t. addition [8]. Interestingly, H.M.T. Menegaz has shown [8] that simply normalizing q after filter updates generally suffices in practice, which we exploit in our filter implementation to reduce computational cost.
Rotating a vector r by the quaternion q is nonlinear in q, expressed as q⊗r⊗q, where q is the conjugate quaternion. We propose the virtual sensor paradigm to include vector measurements in the linear KF without using this nonlinear operation (Sect. 3.2).
The merit of our proposal is that it utilizes the linearity of the attitude kinematics [10]: where is the angular velocity, [ ×] the matrix of the cross product.
However, the dependence on requires modeling the attitude dynamics as well, which is done with the Newton-Euler equation [10] where is the resultant torque affecting the body; J the inertia matrix, while L stands for angular momentum (e.g., if reaction wheels are present). In practice, both equations are discretized with time index t.

Decomposed estimation
To prove that the discrete joint probability density estimation problem can be decomposed, the Markov property should hold [23], i.e., (for given states X and observations Y) In this setting, X = {q, } and Y = {q,̃}, where the accent • denotes measurements. For q and , it is not sufficient to condition only on q t or t , as indirectly relates both (3). Moreover, there is also direct dependence between t and q t .
Nonetheless, additionally conditioning on q t and on t solves the problem. This gives fulfilling the Markov property for states. For observations t and q t are sufficient, too.
Exploiting the Markov property, we can factorize the joint state distribution over time: Thus, we can estimate q and separately if we use the most recent estimates of both quantities.

The virtual sensor
The factorized distribution (4) can be estimated more efficiently than the joint distribution [24]. This is independent of whether a linear or nonlinear filter is used-i.e., deploying a linear filter could further reduce computational costs.
To use a linear filter, the measurements need to be linear in the filter states , q. However, this does not hold for all sensors on CubeSat s [2,4]: gyroscope ( ) and Star Tracker (ST) (q) measurements are linear, but magnetometer, and Sun sensor measurements are nonlinear in the states q, (measuring the Earth's magnetic field B and the Sun vector s, respectively). Namely, B and s need to be compared to a reference vector (as in (1)) to extract attitude information. This requires that the observed vectors are rotated into the reference frame, which is nonlinear in q.
The virtual sensor avoids this nonlinear operation via an intermediate step that extracts the attitude information from the vector observations B, s by solving Wahba's problem (1). We use the computationally efficient ESOQ2 [13] algorithm, which has a closed-form solution for two vector observations (i.e., B and s). However, a different algorithm could have been chosen based on the requirements of the estimation problem. We group the measurements as: where the accent • denotes the measured values, B the magnetic field vector, s the Sun vector, ESOQ2 the virtual sensor, and q ESOQ2 = ESOQ2 B ,s . The subscripts ST and ESOQ2 denote the (virtual) sensor producing the quaternions, whereas the KF and UKF subscripts indicate which filter uses the particular measurements. The above formulation enables us to decompose the estimation into a nonlinear (for ) and a linear (for q) filter. The architecture is shown in Fig. 1, where the roundheaded arrow denotes the time-variant parameters of the matrices of the KF (2). Input arrows (w.r.t. the boxes) denote measurements, while output arrows the predictions. It describes one filter step visually, namely: 1. Collect measurements ̃,q,s,B from the physical sensors. 2. Estimate ̂ from ̃ with the UKF; 3. Convert s,B to q ESOQ2 with the virtual sensor (e.g., ESOQ2); and 4. Estimate q from q, q ESOQ2 ,̃.
The proposed decomposition is independent of the application; only the measurements and the virtual sensor needs to be adapted. ESOQ2 requires a reference vector database ( r i in (1)), which might be a bottleneck for other domains. On the other hand, the virtual sensor can influence robustness as the price for computational efficiency [14,15]-we investigate the latter in Sect. 4 w.r.t. the angle between B and s.

Computational complexity
Computational complexity savings are achieved by: (1) decomposing the estimation of the (N + L)-dimensional state vector (N = dim , L = dim q) into two filters; and (2) using the virtual sensor to make the information content (i.e., an attitude quaternion) of the observations B, s linear in q, which admits a linear KF for estimating the attitude. Using a nonlinear UKF for jointly estimating and q is of O [N + L] 3 , since the UKF has a cubic cost in the state dimensionality [25]. By factorizing the joint state distribution (cf. (4)), separate nonlinear filters can be used with a total cost of O N 3 + L 3 . Since the dynamics of is nonlinear (3), it needs to be estimated with an UKF. However, since the gyroscopes measure directly, this UKF will have a linear observation model.
On the other hand, the estimation of q can be made more efficient with the virtual sensor: since q has linear dynamics (2), if the corresponding observation model can be made linear, the attitude estimate can rely on a linear KF. This is achieved with the virtual sensor, which converts s and B into a quaternion-which is linear in the state q, so a linear KF with cost O(L) can be used [25].
Together, the cost is O N 3 + L , since the virtual sensor's cost does not depend on the dimensionality of the filter's state, adding O(1) [13]. The main advantage of the proposed method is that it relies on conditional independences to decompose the estimation problem into two filters; thus, it uses no approximations: it does not linearize the equations of angular motion as in the EKF [18], but solves Wahba's problem with the virtual sensor to make the observation model linear.

Numerical experiments
This section describes the satellite model, estimation algorithms, control loop, environmental models, and the initial parameters considered for the numerical simulations, and discusses the results.

Experimental setup
Our numerical experiments model a 3U Earth observation CubeSat on a polar, Sun-Synchronous Orbit (SSO) with 98.2 • inclination and 10 AM LTDN and we choose the initial conditions and all parameters based on a real CubeSat mission [26]. Our simulator (written in Python) constitutes the decomposed filter, the control loop, and the environmental models. The Earth's magnetic field is described with the IGRF [27] model (we use the implementation in [28]), the gravitational field by the WGS84 [29] model, whereas satellite position is calculated with SGP4 [30], using TLE data as implemented in [31].
We initialized the attitude quaternion q to [1, 0, 0, 0] , the angular velocity to the zero vector for both the proposed and reference filters, the attitude target for the control loop with longitude 19.040236 • and latitude 47.497913 • (i.e., the coordinates for Budapest). The initial angular velocity norm for the detumbling experiments was set to 1.745 s −1 , and detumbling was considered successful if ‖ ‖ < 0.021 s −1 -these parameters reflect a real CubeSat mission [26]. We use the satellite inertia matrix for a 3U CubeSat [26] (the satellite mass was not required for the simulations) The control system consists of two control laws. Detumbling control is used to to dissipate the rotational energy of the satellite after being ejected from the rocket and determines the control torque with the standard law [10] where B 0 is normalized to unit length and k det = 8 ⋅ 10 −3 is a scalar coefficient-the analysis of [32] provides a formula for choosing this parameter. After the angular velocity of the CubeSat is reduced, we use the control law of [33] for tracking the target attitude. The control torque h is based on the error quaternion q err =q ref ⊗q ref ( • denotes the conjugate quaternion) as: where w err denotes the vector part of q err , k h,q is a scalar, while K h, a diagonal matrix. We empirically selected K h, = 8.5 ⋅ 10 −2 × I 3 and k h,q = 3 ⋅ 10 −3 and use a zero target angular velocity ref .

Results
We compared the attitude determination error of our decomposed filter to the quaternionic QuAdSRUKF of [8], and investigated the robustness of our proposal (Table 1). We assessed the effect of varying the ST sampling frequency (as reduced incoming energy might force the human operators who control the satellite from the Earth to reduce ST usage), epistemic uncertainty of the dynamic model via the inertia matrix (which can change due to the vibration during launch), and its sensitivity to the angle between B and s. We parametrize the attitude with a quaternion q, but we report attitude-related quantities in terms of Euler angles for simplicity.
First, we compare the attitude error in Fig. 2, converted to Euler angles , , (the discontinuity in the plot is only an artifact of converting the quaternion into Euler angles, it is not present in the quaternion components-both algortihms had the same initial conditions). The proposed method has a time-invariant error dynamics, while the quaternionic filter's error is time-variant (denoted by the ref subscript in the figures)-a crucial advantage for the proposed method, since high-resolution Earth observation missions require attitude stability. The error vector's norm's variance is slightly higher for the proposed method, but as the component-wise plot shows, it has a stationary error. We hypothesize that the reason for a slightly larger variance is coming from the ESOQ2 algorithm, the sensitivity of which is analyzed below.
To model epistemic uncertainty we perturbed the satellite's inertia matrix: while ensuring positive definiteness, at the start of each simulation, we multiplied all diagonal elements with a random number uniformly drawn from [0.95;1.05] , for off-diagonal elements, we used samples of 10 −6 magnitude. The Monte Carlo simulations (30 experiments) with perturbed inertia matrix show the superior performance of the proposed method. Both the error norm and the standard deviation were lower than the quaternionic filter's values (Fig. 3). Figure 4 shows the distribution of the attitude error difference between the proposed method and the quaternionic QuAdSRUKF [8] for different ST sampling times, simulating the realistic scenario when the CubeSat does not have enough incoming energy to use the ST frequently. When the ST sampling time is increased over 500 ms, the reference Table 1 Filter parameters S is the square root of the state, Q the process noise, R the observation noise covariance. For the of UKF's -points, w 0 is the center weight ( ∀i ≠ 0, w i = (1 − w 0 )∕2N , where N is the number of states), and a scaling factor for the covariance, for details, cf. [8] 1 This filter uses P = S 2 in the implementation.  Comparing the norm of the attitude error (converted to Euler angles) when the inertia matrix is uncertain. The proposed method is more robust than the QuAdSRUKF of [8] ( ref subscript), as the mean error is slightly less, but variance is significantly smaller for the former (a 3 interval is calculated over 30 runs) Fig. 4 The difference between the attitude determination error of the QuAdSRUKF of [8] and that of the proposed method. Negative values indicate that the proposed method is more sensitive to less frequent ST measurements. The orange lines stand for the median, the boxes spread from the lower to the upper quartile of the data, and the whiskers show the data range quaternionic filter [8] had a lower attitude error, until that point the proposed filter was superior. This shows a limitation of the proposed method compared to the QuAdSRUKF. Nonetheless, even for using the ST only every 125 s , the median error is below 0.005 rad , which can still be considered as an acceptable error, e.g., compared to [34][35][36]. An (almost) parallel s and B can decrease the accuracy of ESOQ2. Indeed, calculating the norm and standard deviation of the error quaternion (‖Δq ESOQ2 ‖ , i.e., the error between the true attitude quaternion and the output of ESOQ2) yields 0.0187 ± 0.0146 when the angle between B and s is greater than 10 • , whereas 0.138 ± 0.148 for smaller angles. The blue line in Fig. 5 shows the norm of the error quaternion. The periodic spikes correspond to scenarios when s and B become (approximately) parallel. To confirm this, we fitted a nonlinear curve onto the spikes (shown in orange) that depends on the inner product of these two vector observations-i.e., Fig. 5 shows that we can capture the extreme variability (the spikes) in (‖Δq ESOQ2 ‖ by a nonlinear function of the inner product of s and B. The proposed formula is: where, c = 40 and k = 90 are empirically determined constants to ensure a low attitude error and f s,B ∈ [0;c] . The quantity f s,B is used as a penalty affecting the diagonal elements of the observation covariance matrix R of the linear KF for the output of the virtual sensor ( q ESOQ2 ), i.e., R = diag 1 ⋅ 10 −3 ×I 4 , f s,B + 1 ⋅ 10 −1 ×I 4 . We note that this strategy is specific for ESOQ2. We experimented with where s,B = �⟨s;B⟩� ‖s‖‖B‖ multiplying s,B with (1 − ‖q ESOQ2 ‖), since our empirical observations suggested that smaller ‖q ESOQ2 ‖ increases the error norm, but that did not change the value of f s,B substantially.
We also report the attitude control error, stability, and the required time to detumble the spacecraft (i.e., when ‖ ‖ < 0.021 s −1 ) with the proposed filter in Table 2. Our numerical evaluation suggests competitive performance compared to attitude determination systems in current Cube-Sat s [34][35][36]. However, since our results are simulated, they should be validated in a physical model.

Conclusion
In this work, we studied attitude estimation for small spacecrafts. We analyzed the governing equations of angular motions and utilized conditional independences via the Markov property to decompose the attitude determination problem into two filters: one for estimating angular velocity, and one for estimating the attitude quaternion. Furthermore, we proposed the concept of a virtual sensor, with the help of which we transformed the nonlinear observation model of the attitude filter into a linear one. These two steps were essential for reducing computational cost: (1) the decomposition of the filter enabled the estimation of a factorized probability density; and (2) the virtual sensor made it possible to use a linear filter for attitude estimation.
Our analysis shows reduced computational cost, whereas our numerical experiments demonstrate a more robust filter w.r.t. epistemic uncertainty and better error dynamics than the nonlinear quaternionic QuAdSRUKF [8]. A possible bottleneck of our method is when ST measurements are less frequently available. Given these considerations, we believe that exploring the proposed methods in real-world CubeSat missions could potentially improve attitude determination in small spacecraft. Whereas the idea of the virtual sensor might be leveraged in other applications of attitude estimation.