Attitude Determination System for a CubeSat Experiencing Eclipse

1. Introduction

Cube satellites (CubeSats) are being developed because they are capable of carrying out missions such as earth observation, astronomical physics etc., in place of large satellites [1,2]. Additionally, CubeSats which are predominant in the low earth orbit (LEO) have been considered to be key enablers in global coverage and connectivity [3]. Furthermore, they use commercial off the shelf (COTS) components, which reduces development and production costs. A CubeSat is a standard cubic satellite with a side of 10 cm and a mass of approximately 1 kg. Large CubeSats can be built using several cubic units, for instance, 2, 3 and 6 unit CubeSat. CubeSats have increased access to space by providing low launch cost since they are a carried as secondary payloads in the launch vehicle. Although they are relatively cheap to produce and launch, they have limitations in size, weight, and power (SWaP) and this limits the size and performance of the subsystems [4].

The attitude determination and control system (ADCS) is one of the satellite subsystems, and its key function is centered around steering the spacecraft to its intended orientation. The ADCS can further be divided into two subsystems; attitude determination and attitude control subsystems. The former uses sensors to compute the rotation of a rigid body about its center of mass (attitude) while the latter uses actuators to reorient the rigid body to the desired rotation. The desired orientation could be nadir pointing to allow ground stations to telecommand, or sun pointing to maximize power generation. It is evident in Figure 1 that the attitude determination system (ADS) serves as a reference for control.

The ADS sensors are broadly categorized into reference and inertial sensors. The reference sensors provide directional measurements of the satellite regarding another celestial object in space, for instance, horizon sensors, sun sensors etc. [4]. Inertial sensors are used to determine the rate of attitude change of an object about an inertial frame of reference [5]. These sensors typically include gyroscopes, which are used to measure angular velocity. Since the ADS uses COTS components, the sensors are usually categorized by high noise levels, which affect the systems’ accuracy. The accuracy of CubeSat attitude estimation can also be influenced by the accumulated attitude error over time, which is caused by the drifts of the inertial sensors used.

The attitude determination or estimation algorithms are used to process sensor measurements to compute the attitude. The computed attitude can take several parameterizations, for instance, it may be described using Euler angles, quaternions, axis angle or directional cosine matrix (DCM) [6,7,8]. However, minimal parameterization such as Euler angles, Rodriguez parameters and modified Rodriguez parameters (MRPs)are frequently avoided in filter designs for the global attitude because they can fall victim of singularities, but they are commonly used to define local error attitudes [9]. The most popular attitude parameterization technique is the quaternion since they encounter no singularity problem.

The attitude determination algorithms or single frame methods use measurement vectors provided by reference sensors to calculate an object’s orientation instantaneously. These deterministic algorithms include triaxial attitude determination (TRIAD), Davenport’s q-method, quaternion estimator (QUEST) and single value decomposition (SVD) [10,11]. These deterministic methods have been well-defined in [12], therefore we will not go into detail here. Attitude estimation algorithms make use of vector measurements, inertial measurements and previous attitude information to provide the current optimal attitude estimate. The estimation process can be divided into filtering of sensor measurements and computation of attitude from sensor observations [13]. As stated in [13,14], the filtering process can be accomplished through various methods, such as particle filtering and derivatives of the Kalman filter.

It is common to couple a magnetometer and a sun sensor to measure attitude in small satellites such as CubeSat because each sensor cannot provide 3-axis attitude knowledge independently [15,16,17,18,19]. When one vector is available, the satellite is free to rotate about that vector, thus, at least two vectors are used [15]. When the satellite is in orbit, it passes through regions where the sunlight is blocked by the earth, referred to as eclipse phase. For low earth orbit satellites, they spend about 35 minutes of the orbital time in this phase [20]. During this region, only one vector is available for computing the satellite’s attitude, resulting in high attitude error. The high attitude error affects attitude control accuracy, and this may result in loss of communication during overpasses. It is possible to use additional auxiliary sensors during this phase, but due to the CubeSat’s limited mass and volume, only minimal components are used.

In an attempt to make the attitude estimation process accurate during both sun and eclipse phase, several methods have been employed. Authors in [21] propose the use of magnetic field derivatives as the second vector in the QUEST aided multiplicative EKF architecture to increase the accuracy of attitude estimation during the eclipse phase. Using only magnetometer measurements, the proposed method can estimate rough attitude and angular rate even during eclipse phases; therefore, it can also be used in case of gyroscope failure. In [22], it was found that the SVD-aided EKF outperformed the traditional EKF during the eclipse phase. This improved performance was due to the capability of the SVD-aided EKF to adapt the measurement error covariance matrix ($\mathbf{R}$) when the ADS lost sensor measurements. During the eclipse phase, when the satellite is in shadow and the reference sensors may not be providing reliable measurements, the SVD-aided EKF utilized an adaptation rule to update the $\mathbf{R}$ matrix. This adaptation allowed the filter to effectively account for the loss of sensor measurements and adjust the covariance matrix accordingly, leading to improved attitude estimation accuracy. This adaptive behavior of the SVD-aided EKF during the eclipse phase was found to be a significant advantage over the traditional EKF, which did not incorporate such adaptation. As a result, the SVD-aided EKF demonstrated superior performance in terms of attitude estimation accuracy during this challenging phase of the satellite’s orbit, as shown in the research findings by [22]. In contrast, the research revealed that the SVD-aided EKF may exhibit poor performance over a longer period of time (e.g., over 1000 seconds) compared to the traditional EKF, particularly during extended eclipse phases. For instance, in the case of low Earth orbit (LEO) CubeSats, the eclipse phase typically lasts for approximately 35 minutes (2100 seconds), as shown in [20]. As a result, the SVD-aided EKF may not be well-suited for attitude estimation during such prolonged eclipse phases. Another paper [23] proposes a prediction approach for estimating satellite attitude during the eclipse phase, and an adaptive scheme that utilizes the SVD-aided EKF to determine the sun phase. However, the study identified that the accuracy of the prediction algorithm was compromised due to the prolonged duration of the eclipse period.

The SVD-aided EKF was developed as an attempt to enhance the performance of the traditional EKF by incorporating measurement noise covariance matrix adaptation. However, it has been observed that the SVD-aided EKF may face limitations in prolonged eclipse phases. In response, further research has been conducted to develop improved approaches that can make the EKF more adaptable to anomalies such as measurement faults during extended eclipse periods. Researchers in [24], the authors present a technique for computing satellite attitude estimation using an EKF with singular value decomposition (SVD) assistance. The algorithm also incorporates simultaneous adjustments of the process and measurement covariance utilizing data obtained from magnetometers and sun sensors. When the process noise increases or the spacecraft is in eclipse, SVD-aided adaptive EKF performs better than SVD aided EKF. Authors in [25], propose a hierarchical and efficient framework for satellite attitude determination that aims to compensate for observation errors in raw attitude data. It includes a simplified adaptive Kalman filtering module, a neural network-based system error compensation module, and a weighted attitude smoothing module. The performance of the proposed framework, trained with full matrix elements, was found to be comparable to the optimal accuracy when compared with conventional algorithms. Another study in [26] introduces a novel algorithm called adaptive iterated extended Kalman filter (AIEKF) for relative position and attitude estimation, considering model uncertainty in a nonlinear stochastic discrete-time system with unknown disturbance. The AIEKF algorithm employs Gauss-Newton iterative optimization steps for maximum a posteriori (MAP) estimation, and incorporates a switch-mode combination technique to achieve adaptability. The mean-square estimation error (MSE) of the state estimate is derived, and it is proven that the AIEKF outperforms traditional extended Kalman filter (EKF) or iterated extended Kalman filter (IEKF) in terms of MSE.

The structure of the paper is as follows: The attitude estimation and static algorithms are presented in Section 2 in detail. Section 3 provides an introduction to the mathematical models utilized for the presented CubeSat. The proposed filter adaptation method is presented in Section 4. In Section 5, the simulation results of the adaptive EKF algorithm for a hypothetical CubeSat are presented. In Section 6, a brief conclusion is given.

2. Attitude Estimation Algorithms

In this section, we look at the singular value decomposition, extended Kalman filter and their hybrid as they have been proved to be robust algorithms.

2.1. The SVD Algorithm

The SVD algorithm is a deterministic technique for computing the attitude matrix that minimizes Wahba’s problem [11,27]. When more than two sets of measurements are available, Wahba defined a problem aiming to minimize the loss between reference and measured unit vectors. The loss function is

$$J\left({\mathbf{A}}_{bi}\right)=\frac{1}{2}\sum _{k=1}^N{w}_k{|{\mathbf{b}}_k-{\mathbf{A}}_{bi}{\mathbf{r}}_k|}^2,$$

(1)

where $w_k$ is the non-negative weight, ${\mathbf{r}}_k$ is a set of reference vectors in the reference frame and ${\mathbf{b}}_k$ is a set of measured vectors in the satellite body frame. An alternative way to simplify the loss function is to use a matrix $\mathbf{B}$ given as,

$$\mathbf{B}=\sum _{k=1}^N{w}_k{\mathbf{b}}_k{\mathbf{r}}_k^{\mathsf{T}}={\mathbf{USV}}^T,$$

(2)

where U and V are orthogonal matrices and S is a diagonal matrix with w entries. The optimal attitude matrix is calculated using $\mathbf{U}$ and $\mathbf{V}$ matrices, the matrix has a maximized trace and $det\left(\mathbf{A}\right)=1$;

$${\mathbf{A}}_{opt}=\mathbf{U}\mathrm{diag}\left[\begin{array}{ccc}1& 1& det\left(\mathbf{U}\right)det\left(\mathbf{V}\right)\end{array}\right]{\mathbf{V}}^{\mathsf{T}}.$$

(3)

To find the rotation angles of the satellite, a transformation matrix should be found from the (3) first with the determinant of one. To examine the rotation angle error, the error covariance matrix is defined as

$${\mathbf{P}}_{svd}=\mathbf{U}\mathrm{diag}\left[{(s_2+s_3)}^{-}1{(s_3+s_1)}^{-}1{(s_1+s_2)}^{-}1\right],$$

(4)

where $s_1=S_{11}$, $s_2=S_{22}$ and $s_3=det\mathbf{U}det\mathbf{V}{S}_{33}$ When the satellite has only two sensors, that is, the sun and magnetometer sensor, the SVD-method will fail when the satellite is in eclipse because it requires at least two vectors to work properly. This draw back of algebraic methods led to the use of extended Kalman filters for attitude estimation.

2.2. The Extended Kalman Filter

Various estimation techniques have been employed for determining the orientation of satellites. Since Extended-Kalman filters (EKF) are appropriate for nonlinear state equations, they can be used as alternative attitude estimators [28]. Nonlinear discrete time systems can be defined as follows:

$${\mathbf{x}}_k=f({\mathbf{x}}_{k-1},w_k,k),$$

(5)

$$z_k=h({\mathbf{x}}_{k+1},k)$$

(6)

where state vector $\mathbf{x}$ is involved in the system and measurement differential functions denoted by f and h, respectively, y is the expected measurement, w and v is the process and measurement noise with known covariance, respectively. Optimal linear state estimation requires linearizing a nonlinear model about an operating point and applying it to the linearized system. The first order extended Kalman filter uses Taylor series expansion to linearize the system about ${\mathbf{x}}_{k-1}$, other terms are truncated after the first order terms [29]. The EKF consists of two steps: prediction and measurement update, that is,

Prediction step

$${\mathbf{x}}_{k+1}={\mathbf{F}}_k{\mathbf{x}}_{k-1},$$

(7)

$${\mathbf{P}}_{k+1}={\mathbf{F}}_k{\mathbf{P}}_{k-1}{\mathbf{F}}_k^{\mathsf{T}}+{\mathbf{Q}}_k.$$

(8)

The prediction step of the recursive extended Kalman filter computes the predicted state vector and predicted error covariance matrix, denoted as ${\mathbf{x}}_{k+1}$ and ${\mathbf{P}}_{k+1}$, respectively. $\mathbf{F}$ is the state transition matrix, which is obtained by computing the Jacobian matrix of the dynamic system nonlinear function f. The state transition can be defined using (15) since it which shows the change of attitude information over time.

Measurement Update

In the measurement update step of the recursive extended Kalman filter, the estimated state vector ${\widehat{\mathbf{x}}}_k$ and error covariance matrix ${\mathbf{P}}_k$ are updated using the Kalman gain $\mathbf{K}k$ and the Jacobian matrix $\mathbf{H}$ of the differential function $h({\mathbf{x}}_{k+1},k)$ in Equation (6), evaluated at ${\mathbf{x}}_{k+1}$. The mathematical formulation of the Kalman gain $\mathbf{K}k$ and ${\widehat{\mathbf{x}}}_k$ are given by

$${\mathbf{K}}_k=\frac{{\mathbf{P}}_{k+1}{\mathbf{H}}_k^{\mathsf{T}}}{{\mathbf{H}}_k{\mathbf{P}}_{k+1}{\mathbf{H}}_k^{\mathsf{T}}+{\mathbf{R}}_k},$$

(9)

$${\widehat{\mathbf{x}}}_k={\mathbf{x}}_{k+1}+{\mathbf{K}}_k({\mathbf{z}}_k-z_k)$$

(10)

$${\mathbf{P}}_k={\mathbf{P}}_{k+1}-{\mathbf{K}}_k{\mathbf{H}}_k{\mathbf{P}}_{k+1}.$$

(11)

The differential function h can be deduced from the relationship between reference vectors $\mathbf{r}$ and the satellite measured vectors $\mathbf{b}$ given by

$$\mathbf{b}=\mathbf{A}\mathbf{r},$$

(12)

where $\mathbf{A}$ is given by (20). The measured vector $\mathbf{b}$ is the expected measurement given a set of references vector; therefore, it is similar to $z_k$ from (6). From here, the Jacobian matrix is computed based on the resulting vector from (12) to get $\mathbf{H}$.

Algorithm 1 provides a summary of the EKF algorithm, in which the variables N and i are utilized to indicate the time window and iteration step, respectively. The inputs of the filter are the initial state ${\mathbf{x}}_{k-1}$, state covariance matrix ${\mathbf{P}}_{k-1}$, sensor measurements ${\mathbf{z}}_k$, measurement, and process covariance matrices ${\mathbf{R}}_k$ and ${\mathbf{Q}}_k$, respectively. All the other inputs except the sensor measurements are initialized based on assumptions, especially in the absence of covariance data. In the Kalman filter, the state equation consists of kinematics and the accuracy of attitude estimation is determined by how realistic the state model is. Due to the numerous possible disturbances in space, establishing a good state model is difficult. Satellites equipped with 3-axis rate gyros are capable of detecting rotational motion without explicitly modeling disturbances. The rate gyros, however, have biases due to temperature and other factors. As a result, it is generally advisable to estimate the gyro drift rate biases with the attitude parameter.    

Algorithm 1: The extended Kalman filter

Data: inputs: ${\mathbf{x}}_{k-1},{\mathbf{P}}_{k-1},{\mathbf{R}}_k,{\mathbf{Q}}_k,{\mathbf{z}}_k$ Result: Optimal state vector ${\mathbf{x}}_k$ Prediction Step: 1 Compute the Jacobian matrix of $f({\mathbf{x}}_{k-1},k)$ to find ${\mathbf{F}}_k^i$ 2 ${\mathbf{x}}_{k+1}={\mathbf{F}}_k^i{\mathbf{x}}_{k-1}$ 3 ${\mathbf{P}}_{k+1}={\mathbf{F}}_k^i{\mathbf{P}}_{k-1}{\left({\mathbf{F}}_k^i\right)}^{\mathsf{T}}+{\mathbf{Q}}_k$ Measurement update 4 Compute the Jacobian matrix of $h({\mathbf{x}}_{k+1},k)$ to find ${\mathbf{H}}_k^i$ 5 $z_k^i={\mathbf{H}}_k^i{\mathbf{x}}_{k+1}$ 6 Update the Kalman gain ${\mathbf{K}}_k^i$ using (9) 7 Update the predicted state ${\mathbf{x}}_{k+1}$ to find the optimal state estimate ${\widehat{\mathbf{x}}}_k^i$ ,using (10) 8 ${\mathbf{P}}_k^i={\mathbf{P}}_{k+1}-{\mathbf{K}}_k^i{\mathbf{H}}_k^i{\mathbf{P}}_{k+1}$

2.3. SVD aided EKF algorithm

The EKF generally performs better than the algebraic methods such as the SVD, however it performs linearization at every stage as shown in Algorithm 1. This computational cost is eliminated by a hybrid algorithm that makes use of both the single frame method and an estimation algorithm such as the SVD aided EKF[24]. The Extended Kalman Filter (EKF) employs the Singular Value Decomposition (SVD) method to obtain a rough estimate of the attitude, such as Euler angles, for state correction. Moreover, the estimation covariance obtained from the SVD is used as the measurement noise covariance matrix in the filter, which inherently makes it adaptive to noise increments in measurements. The SVD acts as a preliminary filter, while the measurement model is linear and the system equations remain nonlinear. In the SVD-aided EKF, the state prediction and measurement model are expressed in discrete time as follows,

$${\mathbf{x}}_{k+1}=f[{\mathbf{x}}_k,k]+{\mathbf{w}}_k$$

(13)

$${\mathbf{z}}_k=\mathbf{G}{\mathbf{x}}_{k+1}+{\mathbf{v}}_k,$$

(14)

where $\mathbf{G}$ is the linear measurement matrix. The algorithm is summarized in Algorithm 1.

Unlike Algorithm 1, Algorithm 2 makes use of the environmental models’ data because the SVD method uses a set of measured vectors and corresponding vectors in the ECI to compute the rough attitude. The best estimate is computed using the attitude information computed using the SVD method and gyroscope measurements instead of unprocessed vector measurements.

Algorithm 2: The SVD aided EKF for attitude estimation

Data: inputs: ${\mathbf{x}}_{k-1},{\mathbf{P}}_{k-1},{\mathbf{Q}}_k,{\mathbf{B}}_b,{\mathbf{B}}_i,{\mathbf{S}}_b,{\mathbf{S}}_i,\omega ,N,t$ Result: Optimal state vector ${\mathbf{x}}_k$ Prediction Step: 1 Similarly as in Algorithm 1. Measurement update 2 Compute $\mathbf{B}$ using ${\mathbf{B}}_b,{\mathbf{B}}_o,{\mathbf{S}}_b,{\mathbf{S}}_o$ in (2). 3 Compute the optimal attitude matrix ${\mathbf{A}}_{opt}$ using (3). 4 Compute $\mathbf{P}$ from (Section 2.1) and use it as ${\mathbf{R}}_k$. 5 Convert ${\mathbf{A}}_{opt}$ to Euler angles using $\varphi =\arctan 2(A_{32},A_{33})$, $\theta =-arcsin{A}_{31}$ and $\psi =arctan2(A_{21},A_{11})$ 6 Define the measurement vector ${\mathbf{z}}_k$ using Euler angles and the angular velocity components. 7 $z_k^i={\mathbf{G}}_k^i{\mathbf{x}}_{k+1}$ 8 Update the Kalman gain ${\mathbf{K}}_k^i$ using (9). 9 Update the predicted state ${\mathbf{x}}_{k+1}$ to find the optimal state estimate ${\widehat{\mathbf{x}}}_k^i$ ,using (10). 10 ${\mathbf{P}}_k^i={\mathbf{P}}_{k+1}-{\mathbf{K}}_k^i{\mathbf{H}}_k^i{\mathbf{P}}_{k+1}$

3. System Model

3.1. System Kinematics

The change in attitude of a satellite is modeled using the kinematics equation, which describes rotational motion despite the cause. The kinematic equation of a satellite using quaternions attitude representation is given by

$$\dot{\mathbf{q}}=\left[\begin{array}{c}\dot{q_0}\\ \dot{q_1}\\ \dot{q_2}\\ \dot{q_3}\end{array}\right]=\frac{1}{2}\left[\begin{array}{ccc}-q_1& -q_2& -q_3\\ q_0& -q_3& q_2\\ q_3& q_0& -q_1\\ -q_2& q_1& q_0\end{array}\right]\left[\begin{array}{c}p\\ q\\ r\end{array}\right]$$

(15)

where $p,q$ and r are the components of the angular velocity $\omega$ provided by the gyroscope [30] and $q_0$ is the scalar component of the quaternion. The kinematic equation is useful when defining the state transition matrix in the extended Kalman filter. In this research, quaternions are employed for attitude computation due to their mathematical advantages, such as singularity avoidance and efficient numerical representation. However, for graphical presentation and visualization, Euler angles are used as they are more commonly used in aerospace engineering and provide a more intuitive representation of attitude. This choice allows for a clear and easy-to-understand presentation of the results in a familiar format. The computed quaternion is converted to Euler angles using;

$$\begin{array}{cc}& \varphi =arctan2(2(q_0{q}_1+q_2{q}_3),1-2(q_1^2+q_2^2)),\\ & \theta =arcsin2\left(2(q_0{q}_2-q_3{q}_1)\right),\\ & \psi =arctan2(2(q_0{q}_3+q_1{q}_2),1-2(q_2^2+q_3^2)).\end{array}$$

(16)

3.2. Environment Models

3.2.1. International Geomagnetic Reference Field (IGRF)

The IGRF denotes the main geomagnetic field generated by internal sources, primarily inside the core of the earth [31]. The absence of electric currents on the surface allows for the computation of the Earth’s geomagnetic field as the negative gradient of a scalar potential, such that the magnetic field can be expressed as $\mathbf{B}=-\delta V$. The potential function is represented by a finite series expansion in terms of Gauss coefficients, ${\mathbf{g}}_n^m$ and ${\mathbf{h}}_n^m$ [32]:

$$\begin{array}{cc}\hfill \mathbf{B}(r,\beta ,\alpha )=V(r,\beta ,\alpha ,t)=& a\sum _{n=1}^N\sum _{m=1}^n{\left(\frac{a}{r}\right)}^{n+1}\hfill \\ \hfill & [{\mathbf{g}}_n^m\left(t\right)cosm\alpha +{\mathbf{h}}_n^m\left(t\right)sinm\alpha ]P_n^m(cos\beta ).\hfill \end{array}$$

(17)

Here, $r,\beta ,\alpha ,t$ is the radial distance from the Earth’s center, geocentric co-latitude, longitude and time respectively. The variable a is the Earth’s radius. More information of the spherical harmonic coefficients can be found in [31]. The latest version in the series, IGRF-13, extends up to thirteen harmonic degrees.

3.2.2. PSA Sun Position Algorithm

The paper utilizes the PSA algorithm to calculate the solar vector in the inertial frame due to its small computational footprint and high accuracy. It has been established that when applied over a period ranging from 2020 to 2050, the maximum error in angular deviation from the actual solar vector is only 35 arc-sec, with the algorithm maintaining its computational structure and simplicity [33]. The PSA algorithm’s inputs consist of the location’s latitude and longitude, along with the time specified in universal time (UT1), encompassing year, month, day, hours, minutes, and seconds. The sun vector is given by

$$\mathbf{S}=\left[\begin{array}{c}sin{\vartheta }_{z}sin\gamma \\ sin{\vartheta }_{z}cos\gamma \\ cos{\vartheta }_z\end{array}\right]$$

(18)

where $\vartheta$ and $\gamma$ is the zenith and azimuth angle. More information on how these angles are obtained can be found in [33].

3.3. Measurement Model

In this model, we assume that only zero-mean Gaussian white noise can corrupt the sensor measurements. The measurement model for the vector observations can be given as

$${\mathbf{b}}_i=\mathbf{A}{\mathbf{r}}_i+{\mu }_i,$$

(19)

where $\mathbf{b}$ is the sun and magnetic field vector observations in the body frame, $\mathbf{r}$ can either be $\mathbf{B}$ or $\mathbf{S}$ modelled in the reference frame, $\mu$ is the Gaussian white noise and $\mathbf{A}$ is given by (20).

$$\mathbf{A}\left(\mathbf{q}\right)=\left[\begin{array}{ccc}{q}_0^2+q_3^2-q_1^2-q_2^2& 2(q_0{q}_1+q_2{q}_3)& 2(q_0{q}_2-q_1{q}_3)\\ 2(q_0{q}_1-q_2{q}_3)& q_1^2+q_3^2-q_0^2-q_2^2& 2(q_1{q}_2+q_0{q}_3)\\ 2(q_0{q}_2+q_1{q}_3)& 2(q_1{q}_2-q_0{q}_3)& q_0^2-q_1^2+q_2^2+q_3^2\end{array}\right]$$

(20)

For gyroscopes, the inherent sensor bias introduces errors throughout the attitude determination process; therefore, these errors are also included in the state vector estimate. The gyroscope mathematical model presented in this paper is

$${\widehat{\omega }}_b={\omega }_b+\mathbf{b}+\mathbf{v},$$

(21)

where gyroscope noise and bias of a satellite are represented by $\mathbf{b}$ and $\mathbf{v}$, respectively, while the true angular velocity of the satellite relative to the inertial frame is denoted by ${\omega }_b$.

3.4. Coordinate Systems

Several reference frames are used to describe the satellite’s attitude. This paper employs three reference frames, each of which is defined by its fixed element or direction, as well as by the location of its center or origin, that is, satellite body coordinate frame, orbit referenced coordinate frame and Earth Centered Inertial (ECI) coordinate frame. These coordinate frames have been well explained in [14].

4. Adaptive EKF

During extended eclipse phases, it was noted that the traditional Extended Kalman Filter (EKF) exhibited superior performance compared to other methods like the Singular Value Decomposition (SVD)-aided EKF. However, if the accuracy of the process and measurement noise covariance matrices is compromised, the EKF may produce substantial estimation errors, despite outperforming other algorithms in prolonged eclipses. As a result, there is a great need for developing an adaptive EKF that is compatible with inaccurate covariance matrices. Firstly, we discuss how each covariance matrix affects the performance of the extended Kalman filter (EKF). The process noise covariance matrix ${\mathbf{Q}}_k$ can cause inaccurate prediction of the error covariance matrix ${\mathbf{P}}_{k+1}$ in (8) which directly affects the computation of the Kalman gain $\mathbf{K}$ as shown in (9). Finally, the optimal estimate ${\mathbf{x}}_k$ and the estimated error covariance matrix are derived from the incorrect ${\mathbf{P}}_{k+1}$ and $\mathbf{K}$ as seen in (10) and (11), respectively. Additionally, it is clear that an imprecise measurement noise covariance matrix ${\mathbf{R}}_k$ can lead to an imprecise Kalman gain, resulting in inadequate state estimation. It is apparent that $\mathbf{P}$ and $\mathbf{R}$ directly affects the state estimate as compared to $\mathbf{Q}$. In this paper, a $\mathbf{P}$ and $\mathbf{R}$ adaptive filter is proposed for attitude estimation. The filter performance will be compared to the SVD aided EKF in the sun and eclipse phases, which performs better than the conventional EKF in attitude estimation.

4.1. Expectation Maximization

There is usually inadequate information for training Bayesian networks; for example, it is difficult to accurately determine noise covariance matrices in a Kalman filter [34]. Such latent parameters can explicitly be estimated using the expectation maximization (EM) algorithm, which is an iterative algorithm used in unsupervised learning. The EM algorithm is an approach to find the maximum likelihood (ML) estimates of the latent variables in a statistical method. It is used to find the maximum likelihood estimate of the parameters, given the observed data. The algorithm consists of two steps: the Expectation step, where the expectation of the log-likelihood function is calculated, given the current estimates of the parameters, and the Maximization step (M-step), where the parameters are updated to maximize the expected log-likelihood [35]. The algorithm iterates between these two steps until convergence, at which point the parameters are considered to be the maximum likelihood estimates. The goal of this paper is to use the joint log-likelihood function to estimate inaccurate noise covariance matrices ${\mathbf{P}}_k$ and ${\mathbf{R}}_k$ based on the available measurements ${\mathbf{z}}_{1:k}$ (i.e., measurements taken from time 1 to k) and the unknown state vector ${\mathbf{x}}_k$, that is,

$$L_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_k)=log{p}_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_k),$$

(22)

where ${\sigma }_k=[{\mathbf{P}}_k,{\mathbf{R}}_k]$. The joint log-likelihood function contains the natural logarithm function log and the probability density function p, which depends on the parameter $\sigma$. To obtain an approximate maximum likelihood (ML) solution for the parameter ${\sigma }_k$, we employ the expectation-maximization (EM) algorithm, which maximizes the joint log-likelihood function, such that

$$\begin{array}{cc}\hfill {\widehat{\sigma }}_k& =argmax{L}_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_k)\hfill \\ & \triangleq argmaxlog{p}_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_k),\hfill \end{array}$$

(23)

where ${\widehat{\sigma }}_k$ is the ML estimate of ${\sigma }_k$. It is impossible to solve the complete data log-likelihood function $L_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_k)$ due to the unavailable ${\mathbf{x}}_k$. The EM algorithm addresses the problem by approximating the joint likelihood function in (23) as its minimum variance estimate $\mathbf{Q}({\sigma }_k,{\sigma }_k^i)$, such that [36,37]

$$\begin{array}{cc}\hfill \mathbf{Q}({\sigma }_k,{\sigma }_k^i)\triangleq & {\mathbf{E}}_{{\sigma }_{k^i}}\left[L_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_k)\right|{\mathbf{z}}_k]\hfill \\ \hfill =& \int L_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_k)p_{{\sigma }_k}\left({\mathbf{x}}_k\right|{\mathbf{z}}_k)d{\mathbf{x}}_k,\hfill \end{array}$$

(24)

where ${\sigma }_k^i$ is an approximation of ${\widehat{\sigma }}_k$ at the $i^{th}$ step.

The EM algorithm takes advantage of this property to generate a series of values ${\sigma }_k$, where $k=1,2,\cdots$, with the aim of successively improving the accuracy of the maximum likelihood (ML) estimate. The algorithm is summarized in Algorithm 3. During the E-step, the complete data log-likelihood’s expectation is evaluated, which relies on the current estimates ${\sigma }_k^i$ and measurements ${\mathbf{z}}_{1:k}$. Following this, the M-step maximizes the computed $\mathbf{Q}({\sigma }_k,{\sigma }_k^i)$ utilizing an arithmetic technique.

Algorithm 3: The Expectation Maximization Algorithm

1. Initialization: ${\sigma }_k$ 2. E-step: Compute $\mathbf{Q}({\sigma }_k,{\sigma }_k^i)$ 3. Maximization: Calculate, ${\sigma }_{k+1}=argmax\mathbf{Q}({\sigma }_k,{\sigma }_k^i)$ 4. If convergence has not been reached, proceed to update the variable k by incrementing it by 1, and return to step 2.

4.1.1. E Step

The EM algorithm approximates the joint likelihood function as it minimum variance estimate as shown in (24). To determine the minimum variance estimate, we compute the joint log likelihood $logp{\sigma }_k({\mathbf{x}}_k,{\mathbf{z}}_{1:k})$ and the second posterior PDF $p_{{\sigma }_k}\left({\mathbf{x}}_k\right|{\mathbf{z}}_k)$. Firstly, the joint log likelihood in (22) is written as the product of the conditional likelihood, such that,

$$\begin{array}{cc}\hfill p_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_{1:k})& =p_{{\sigma }_k}\left(z_k\right|{\mathbf{x}}_k,{\mathbf{z}}_{1:k})p_{{\sigma }_k}\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})p\left({\mathbf{z}}_{1:k}\right)\hfill \\ & =p_{{\sigma }_k}\left(z_k\right|{\mathbf{x}}_k)p_{{\sigma }_k}\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})p\left({\mathbf{z}}_{1:k}\right)\hfill \end{array}$$

(25)

where $p_{{\sigma }_k}\left(z_k\right|{\mathbf{x}}_k,{\mathbf{z}}_{1:k})=p_{{\sigma }_k}\left(z_k\right|{\mathbf{x}}_k)$ since expected measurement ${\mathbf{z}}_k$ is only dependent on state ${\mathbf{x}}_k$ as shown in (6). It is worth noting that the probability distribution $p\left({\mathbf{z}}_{1:k}\right)$ of the sensor measurements is typically not dependent on the current state error and measurement noise covariance matrix ${\sigma }_k$. This is because the measurements are generally assumed to be independent of both the state error covariance matrix and the measurement noise covariance matrix, which together determine the joint distribution of the state estimates and the measurements [37]. In the EKF model, the predicted PDF $p_{{\sigma }_k}\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})$ is approximated as Gaussian, that is,

$$p_{{\sigma }_k}\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})=N({\mathbf{x}}_k;{\mathbf{x}}_{k+1},{\mathbf{P}}_{k-1}).$$

(26)

The symbol $N(.;\mu ,\sum )$ represents the Gaussian probability density function (PDF) with mean vector $\mu$ and covariance matrix ∑, and ${\mathbf{x}}_{k+1}$ is obtained using Equation (7). At this stage ${\mathbf{P}}_{k-1}$ is considered inaccurate because of ${\mathbf{Q}}_k$. The likelihood PDF of the measurement model in (8) is given by

$$p_{{\sigma }_k}\left(z_k\right|{\mathbf{x}}_k)=N(z_k;h({\mathbf{x}}_{k+1},k),{\mathbf{R}}_k).$$

(27)

The joint log-likelihood function can be computed using (25)-(27). The logarithm of the normal distribution $N(z_k;h({\mathbf{x}}_{k+1},k),{\mathbf{R}}_k)$ is:

$$logN(z_k;h({\mathbf{x}}_{k+1},k),{\mathbf{R}}_k)=-0.5log\left|{\mathbf{R}}_k\right|-0.5{[z_k-h({\mathbf{x}}_{k+1},k)]}^{\mathsf{T}}{\mathbf{R}}^{-1}[z_k-h({\mathbf{x}}_{k+1},k)].$$

(28)

The logarithm of the normal distribution $N({\mathbf{x}}_k;{\mathbf{x}}_{k+1},{\mathbf{P}}_{k-1})$ is:

$$logN({\mathbf{x}}_k;{\mathbf{x}}_{k+1},{\mathbf{P}}_{k-1})=-0.5log\left|{\mathbf{P}}_{k+1}\right|-0.5{({\mathbf{x}}_k-{\mathbf{x}}_{k+1})}^{\mathsf{T}}{\mathbf{P}}_{k+1}^{-1}({\mathbf{x}}_k-{\mathbf{x}}_{k+1}).$$

(29)

Putting everything together:

$$\begin{array}{cc}\hfill log{p}_{{\sigma }_k}({\mathbf{x}}_k,{\mathbf{z}}_{1:k})=& -0.5|{\mathbf{R}}_k|-0.5{[z_k-h({\mathbf{x}}_{k+1},k)]}^{\mathsf{T}}{\mathbf{R}}^{-1}[z_k-h({\mathbf{x}}_{k+1},k)]\hfill \\ & -0.5|{\mathbf{P}}_{k+1}|-0.5{({\mathbf{x}}_k-{\mathbf{x}}_{k+1})}^{\mathsf{T}}{\mathbf{P}}_{k+1}^{-1}({\mathbf{x}}_k-{\mathbf{x}}_{k+1})+d_{{\sigma }_k},\hfill \end{array}$$

(30)

where |.| is the determinant operation of a matrix. After determining the joint likelihood in (30), the posterior PDF is computed. At the $i+1^{th}$ step, the state vector ${\mathbf{x}}_k^i$ estimate and corresponding state estimation error covariance matrix ${\mathbf{P}}_k^i$ have been calculated. From here, the nonlinear measurement is function is linearized using the intermediate state estimate ${\mathbf{x}}_k^i$, that is,

$$h({\mathbf{x}}_k,k)=h({\widehat{\mathbf{x}}}_k^i,k)+{\mathbf{H}}_k^i({\mathbf{x}}_{k+1}-{\widehat{\mathbf{x}}}_k^i),$$

(31)

where ${\mathbf{H}}_k^i$ is the Jacobian matrix of the measurement function at ${\mathbf{x}}_k^i$. The second posterior PDF $p_{{\sigma }_k^i}\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})$ is approximated as Gaussian by using measurement update of the EKF, such that,

$$p_{{\sigma }_k^i}\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})=N({\mathbf{x}}_k;{\mathbf{x}}_k^{i+1},{\mathbf{P}}_k^{i+1})$$

(32)

where ${\mathbf{P}}_k^{i+1}$ and ${\mathbf{x}}_k^{i+1}$ are obtained using ${\mathbf{P}}_k^i$ and ${\mathbf{x}}_k^i$, respectively given as,

$${\mathbf{K}}_k^{i+1}={\mathbf{P}}_{k+1}^i{\left({\mathbf{H}}_k^i\right)}^{\mathsf{T}}{[{\mathbf{H}}_k^i{\mathbf{P}}_{k+1}^i{\left({\mathbf{H}}_k^i\right)}^{\mathsf{T}}+R_k^i]}^{-1}$$

(33)

$${\mathbf{x}}_k^{i+1}={\mathbf{x}}_{k+1}+{\mathbf{K}}_k^{i+1}[z_k-(h({\widehat{\mathbf{x}}}_k^i,k)+{\mathbf{H}}_k^i({\mathbf{x}}_{k+1}-{\widehat{\mathbf{x}}}_k^i))]$$

(34)

$${\mathbf{P}}_k^{i+1}={\mathbf{P}}_{k+1}^i-{\mathbf{K}}_k^{i+1}{\mathbf{H}}_k^i{\mathbf{P}}_{k+1}^i$$

(35)

To obtain the minimum variance estimate, Equations (30) and (32) are substituted into Equation (24). The resulting expression is

$$\begin{array}{cc}\hfill & \mathbf{Q}({\sigma }_k,{\sigma }_k^i)=-0.5\left|{\mathbf{R}}_k\right|-0.5\mathtt{tr}\left({\mathbf{A}}_k{\mathbf{R}}_k^{-1}\right)\hfill \\ & -0.5|{\mathbf{P}}_{k+1}|-0.5\mathtt{tr}\left({\mathbf{B}}_k{\mathbf{P}}_{k+1}^{-1}\right)+d_{{\sigma }_k}\hfill \end{array}$$

(36)

where $\mathtt{tr}$ is the trace of a matrix. The matrices ${\mathbf{A}}_k$ and ${\mathbf{B}}_k$ are given by

$${\mathbf{A}}_k=\int [z_k-h\left({\mathbf{x}}_{k+1}\right)]{[z_k-h\left({\mathbf{x}}_{k+1}\right)]}^{\mathsf{T}}N({\mathbf{x}}_k;{\mathbf{x}}_k^{i+1},{\mathbf{P}}_k^{i+1})d{\mathbf{x}}_k$$

(37)

$${\mathbf{B}}_k=\int ({\mathbf{x}}_k-{\mathbf{x}}_{k+1}){({\mathbf{x}}_k-{\mathbf{x}}_{k+1})}^{\mathsf{T}}N({\mathbf{x}}_k;{\mathbf{x}}_k^{i+1},{\mathbf{P}}_k^{i+1})d{\mathbf{x}}_k.$$

(38)

${\mathbf{A}}_k$ represents the squared Mahalanobis distance between the prediction and the observation, and ${\mathbf{B}}_k$ can be seen as an intermediate step in the calculation of the expected value. The measurement is further linearized at ${\mathbf{x}}_k^{i+1}$ such that,

$$h({\mathbf{x}}_k,k)=h({\widehat{\mathbf{x}}}_k^{i+1},k)+{\mathbf{H}}_k^{i+1}({\mathbf{x}}_{k+1}-{\widehat{\mathbf{x}}}_k^{i+1}),$$

(39)

Equation (39) is substituted in (37) thus ${\mathbf{A}}_k$ is computed as

$${\mathbf{A}}_k=[{\mathbf{z}}_k-h({\widehat{\mathbf{x}}}_k^{i+1},k)][{\mathbf{z}}_k-h({\widehat{\mathbf{x}}}_k^{i+1},k)){]}^{\mathsf{T}}+{\mathbf{H}}_k^{i+1}{\mathbf{P}}_k^{i+1}{\left({\mathbf{H}}_k^{i+1}\right)}^{\mathsf{T}}.$$

(40)

From there ${\mathbf{B}}_k$ can be computed as

$${\mathbf{B}}_k={\mathbf{P}}_k^{i+1}+({\widehat{\mathbf{x}}}_k^{i+1}-{\mathbf{x}}_{k+1}){({\widehat{\mathbf{x}}}_k^{i+1}-{\mathbf{x}}_{k+1})}^{\mathsf{T}}$$

(41)

4.1.2. M Step

The M-step entails maximizing $\mathbf{Q}({\sigma }_k,{\sigma }_k^i)$ regarding ${\sigma }_k$, such that

$${\sigma }_k^{i+1}\approx \mathtt{argmax}\mathbf{Q}({\sigma }_k,{\sigma }_k^i).$$

(42)

At maximum point, ${\sigma }_k^{i+1}$ satisfies

$$\frac{\delta \mathbf{Q}({\sigma }_k,{\sigma }_k^i)}{\delta {\sigma }_k}{|}_{{\sigma }_k={\sigma }_k^{i+1}}=0.$$

(43)

Since ${\sigma }_k=\left[{\mathbf{P}}_{k-1}\right],{\mathbf{R}}_k]$, (43) becomes a partial derivative problem, that is, $\frac{\delta \mathbf{Q}({\sigma }_k,{\sigma }_k^i)}{\delta {\mathbf{P}}_{k-1}}$ and $\frac{\delta \mathbf{Q}({\sigma }_k,{\sigma }_k^i)}{\delta {\mathbf{R}}_k}$. The partial derivative can be calculated by exploiting (36) as

$$\frac{\delta \mathbf{Q}({\sigma }_k,{\sigma }_k^i)}{\delta {\mathbf{P}}_{k-1}}=-0.5{\mathbf{P}}_{k-1}^{-1}+0.5{\mathbf{P}}_{k-1}^{-1}{\mathbf{B}}_k{\mathbf{P}}_{k-1}^{-1}$$

(44)

$$\frac{\delta \mathbf{Q}({\sigma }_k,{\sigma }_k^i)}{\delta {\mathbf{R}}_k}=-0.5{\mathbf{R}}_k{-}^{-1}+0.5{\mathbf{R}}_k^{-1}{\mathbf{B}}_k{\mathbf{R}}_k^{-1}$$

(45)

where ${\mathbf{A}}_k$ and ${\mathbf{B}}_k$ are given by (37) and (38), respectively. Substituting (44) and (45) in (43) yields

$$-0.5{\left({\mathbf{P}}_{k-1}^{i+1}\right)}^{-1}+{\left({\mathbf{P}}_{k-1}^{i+1}\right)}^{-1}{\mathbf{B}}_k{\left({\mathbf{P}}_{k-1}^{i+1}\right)}^{-1}=0$$

(46)

$$-0.5{\left({\mathbf{R}}_k^{i+1}\right)}^{-1}+{\left({\mathbf{R}}_k^{i+1}\right)}^{-1}{\mathbf{A}}_k{\left({\mathbf{R}}_k^{i+1}\right)}^{-1}=0,$$

(47)

therefore, the error and measurement noise covariance matrices are defined as ${\mathbf{P}}_{k-1}^{i+1}={\mathbf{B}}_k$ and ${\mathbf{R}}_k^{i+1}={\mathbf{A}}_k$, respectively.

4.2. $\mathbf{P}$ and $\mathbf{R}$ adaptive EKF

The proposed adaptive Extended Kalman Filter (EKF) algorithm involves two primary steps: the time update and the iterated measurement update. In the time update, the predicted state vector ${\mathbf{x}}_{k+1}$ and the nominal predicted state error covariance matrix ${\mathbf{P}}_{k+1}$ are calculated using Equations (7) and (8), respectively. The algorithm requires several inputs, including the initial state vector ${\mathbf{x}}_{k-1}$, state error covariance matrix ${\mathbf{P}}_{k-1}$, measurement covariance matrix ${\mathbf{R}}_k$, process noise covariance matrix ${\mathbf{Q}}_k$, and the number of iterated measurements N. The prediction step only computes the state vector and error covariance matrix, considering the system model kinematics and the previous state error covariance matrix. The nominal predicted error covariance matrix ${\mathbf{P}}_{k+1}$ serves as an appropriate initial value for ${\mathbf{P}}_{k+1}^0$, as it includes the information on the state transition matrix $\mathbf{F}$, process noise covariance matrix $\mathbf{Q}$, and the previous estimation error covariance matrix ${\mathbf{P}}_{k-1}$. The for loop estimates the optimal state vector ${\mathbf{x}}_k$, predicted error covariance matrix ${\mathbf{P}}_k$, and measurement noise covariance matrix ${\mathbf{R}}_k$ in an adaptive manner.

Algorithm 4: $\mathbf{P}$ and $\mathbf{R}$ adaptive EKF for attitude estimation

Inputs: ${\mathbf{x}}_{k-1},{\mathbf{P}}_{k-1},{\tilde{\mathbf{R}}}_k,{\mathbf{Q}}_k,{\mathbf{B}}_b,{\mathbf{S}}_b,\omega ,N,t$ Prediction Step: 1. ${\mathbf{x}}_{k+1}={\mathbf{F}}_k{\mathbf{x}}_{k-1}+\mathbf{B}\omega$ 2. ${\mathbf{P}}_{k+1}={\mathbf{F}}_k{\mathbf{P}}_{k-1}{\mathbf{F}}_k^{\mathsf{T}}+{\mathbf{Q}}_k$ Measurement update 3. Initialize at $i=0,{\mathbf{P}}_{k+1}^{\left(0\right)}={\mathbf{P}}_{k+1},$ ${\mathbf{R}}_k^0={\tilde{\mathbf{R}}}_k,{\mathbf{x}}_{k+1}^0={\mathbf{x}}_{k+1}$ for$i=0:N-1$ 4. Compute the Jacobian matrix of $h({\mathbf{x}}_{k+1},k)$ to find ${\mathbf{H}}_k^i$ 5. ${\mathbf{K}}_k^i=\frac{{\mathbf{P}}_{k+1}^i{\mathbf{H}}_k^{iT}}{{\mathbf{H}}_k^i{\mathbf{P}}_{k+1}^i{\left({\mathbf{H}}_k^i\right)}^{\mathsf{T}}+{\mathbf{R}}_k^i}$ 6. $\left[\begin{array}{cc}{\mathbf{B}}_b& {\mathbf{S}}_b\end{array}\right]$ 7. ${\mathbf{x}}_k^i={\mathbf{x}}_{k+1}+{\mathbf{K}}_k^i({\mathbf{z}}_k-h\left(\mathbf{x}\right)-{\mathbf{H}}_k^i{\mathbf{x}}_{k+1})$ 8. ${\mathbf{P}}_k^i={\mathbf{P}}_{k+1}-{\mathbf{K}}_k^i{\mathbf{H}}_k^i{\mathbf{P}}_{k+1}$ 9. Compute Jacobian matrix ${\mathbf{H}}_k^{i=1}$ with updated states estimates ${\mathbf{x}}_k^i$ 10. ${\mathbf{A}}_k=[{\mathbf{z}}_k-h\left(\mathbf{x}\right)]{[{\mathbf{z}}_k-{\mathbf{H}}_k^i{\mathbf{x}}_{k+1}]}^{\mathsf{T}}+{\mathbf{H}}_k^{i+1}{\mathbf{P}}_k^i{\left({\mathbf{H}}_k^{i+1}\right)}^{\mathsf{T}}$ 11. ${\mathbf{B}}_k={\mathbf{P}}_k^i+({\mathbf{x}}_k^i-{\mathbf{x}}_{k+1}){({\mathbf{x}}_k^i-{\mathbf{x}}_{k+1})}^{\mathsf{T}}$ 12. ${\mathbf{P}}_k^{i+1}={\mathbf{B}}_k,{\mathbf{R}}_k^{i+1}={\mathbf{A}}_k$ end for 13. ${\widehat{\mathbf{P}}}_k={\mathbf{P}}_k^{i+1},{\widehat{\mathbf{R}}}_k={\mathbf{R}}_k^{i+1}$ Outputs:${\mathbf{x}}_k,{\widehat{\mathbf{P}}}_k,{\widehat{\mathbf{R}}}_k$

5. Results and Analysis

In the simulations, a CubeSat in a low Earth orbit is considered, with a principal moment of inertia of $\mathbf{J}=\mathtt{diag}\left[\begin{array}{ccc}0.0071& 0.0035& 0.0035\end{array}\right]$kg.m${}^{\mathbf{2}}$. The spacecraft is assumed to be tumbling, and its initial state is arbitrarily chosen as ${\mathbf{x}}_0={[0.03\mathrm{rad}-0.001\mathrm{rad}0.002\mathrm{rad}0.08\mathrm{rad}/\mathrm{s}-0.02\mathrm{rad}/\mathrm{s}0.03\mathrm{rad}/\mathrm{s}]}^{\mathsf{T}}$. The orbit parameters, satellite kinematics, sensor and environmental models are used to validate the proposed algorithm for attitude estimation. In this study, sun sensor, magnetometer, and gyroscope are used as attitude sensors. The local magnetic vector is computed using the international geomagnetic reference field model, while a sun position algorithm is used to extract the local reference sun vector. The orbit used in this paper is circular with an inclination of 0${}^{\circ }$ at an altitude of 600 km. The simulation epoch is 2021.01.01, UTC 00:00 with a time-step of 1 second.

The first part of this section is devoted to analyzing the performance of SVD-EKF and proposed AEKF during sun phase on the orbit. The second part analyses the performance of the two algorithms during the eclipse phase. The simultaneous $\mathbf{P}$ and $\mathbf{R}$ adaptive extended Kalman filter is referred to as AEKF in this study. The evaluation of the estimation accuracy for the measurement noise covariance matrix and predicted error covariance matrix was carried out using the SRNFN and ASRNFN metrics, which were chosen as the error measures. These metrics provide information about the magnitude of the errors in the estimated matrices by measuring the difference between the estimated and reference values, and are defined as

$$\mathtt{SRNFN}={\left(\frac{1}{n^2}{{\mathbf{P}}_{k+1}-{\mathbf{P}}_{r,k+1}}^2\right)}^{\frac{1}{4}},$$

(48)

$$\mathtt{ASRNFN}=\frac{1}{T}\sum _{k=1}^T{\left(\frac{1}{n^2}{{\mathbf{P}}_{k+1}-{\mathbf{P}}_{r,k+1}}^2\right)}^{\frac{1}{4}}.$$

(49)

The predicted error covariance at the time k is represented by ${\mathbf{P}}_{k+1}$, while ${\mathbf{P}}_{r,k+1}$ refers to the reference value of the predicted error covariance matrix at that same time. Similarly, the formulas used for calculating SRNFN and ASRNFN for the predicted error covariance matrix are also applied to the measurement noise covariance matrix. The SRNFN is a metric that measures the similarity between two matrices. The SRNFN is widely used in the field of computer vision and image processing as a measure of the similarity between two images [38]. The Averaged Square Root of Normalized Frobenius Norm (ASRNFN) is an extension of the SRNFN, which is used to measure the similarity between multiple matrices. The ASRNFN is calculated as the average of the SRNFN between each pair of matrices [39].

5.1. Sun Phase

Firstly, we compare the performance of the SVD-EKF and the proposed AEKF during sun phases. During the sun phase or normal operation mode, the filter depends on sun sensor and magnetometer measurements for measurement update. The proposed algorithm is affected by the number of iterated measurement N, the more you iterate the, the more accurate the state estimates, however this means increased computational time. For these simulations $N=5$. Table 1 shows the root-mean-square error of the attitude angles for both filters considered in this study. It is evident that both filters perform reasonably well under the sun mode, as all sensor measurements are available.

Table 2 shows the single step run time for the two algorithms, it shows that the proposed algorithm has more computational time as compared to the SVD-EKF. Even though, the computational load is more by a small margin, the accuracy of the proposed algorithm is significant. To demonstrate the efficacy of the proposed algorithm, we calculated the SRNFN and ASRNFN values of the predicted error and measurement noise covariance matrices, which are useful in analyzing the convergence of these matrices. Analysis of the results presented in Figure 2 and Table 3 indicate that the proposed algorithm outperforms the SVD-EKF algorithm, as evidenced by its lower SRNFN and ASRNFN values.

5.2. eclipse phase

According to a study in [20], satellites in low earth orbit are in the eclipse phase for approximately 28 minutes. During this phase, the attitude determination system depends solely on magnetometer measurements. Table 4 displays the root-mean-square error of the SVD-EKF method and AEKF. While the proposed algorithm’s error has increased slightly, it still performs significantly better than the SVD-aided EKF. Table 4 confirms that the gyroscope bias has increased during the eclipse phase, as compared to the values in Table 1. This increase in bias is likely due to the temperature changes that occur during the eclipse phase. Although modern gyroscopes are designed to minimize temperature-related errors and often include temperature compensation features, the effects of temperature changes cannot be eliminated. Therefore, it is crucial to consider the potential impact of temperature changes on gyroscope bias when designing and implementing attitude determination systems for satellites.

The inferior estimation accuracy of the SVD-EKF algorithm during the eclipse phase is likely due to the sensitivity of the SVD algorithm, which requires at least two vector measurements to compute attitude. Despite relying solely on magnetometer measurements, the proposed algorithm performs reasonably well. Furthermore, Figure 3 and Table 5 show that the proposed algorithm has a smaller SRNFN and ASRNFN than the SVD -EKF even in eclipse phase. The noise during the eclipse phase affects the measurements used to estimate the measurement noise and state covariance matrices; therefore, the resulting matrices are less accurate or less precise.

6. Conclusions

This paper presents the design and numerical analysis of an error and measurement noise covariance matrices adaptive extended Kalman filter algorithm based on expectation maximization method. A novel attitude estimation algorithm was created utilizing an adaptive filter approach. The performance of this newly developed algorithm, as well as a state-of-the-art existing algorithm, was evaluated through application to CubeSat attitude estimation in both the eclipse and sun phases of the orbit. During both the phases, the proposed adaptive EKF outperforms the existing algorithm, being the SVD-EKF. The SVD algorithm requires a minimum of two vector measurements and associated reference models to function optimally. However, during the eclipse phase, only magnetometer measurements are available, which causes the performance of the SVD-aided algorithm to deteriorate. Based on the simulation findings outlined in the paper, the proposed algorithm demonstrated significantly improved accuracy in attitude estimation compared to SVD aided EKF algorithms. However, it was also observed that the proposed algorithm has marginally higher computational complexity.

To ensure the accuracy and reliability of the proposed algorithm in various conditions and applications, it is necessary to conduct rigorous validation and verification through both simulations and real-world tests. This will help to confirm the effectiveness of the algorithm and identify any potential limitations or areas for improvement. Furthermore, the proposed algorithm can be optimized to improve its performance and reduce computational complexity. This could involve refining the algorithm’s underlying mathematical models and algorithms, as well as optimizing its implementation and hardware requirements. By doing so, the algorithm could potentially be made more efficient and practical for use in a wider range of systems and applications.