# High-Accuracy and Low-Cost Attitude Measurement Unit of the CubeSat

Authors: Haining Ma; Zhengliang Lu; Xiang Zhang; Wenhe Liao; Klaus Briess
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

This paper proposes high-accuracy and reliable attitude measurement methods exclusive for CubeSat with restrictions of low cost, limited space, and low power consumption. The attitude measurement unit is equipped with Commercial Off-The-Shelf (COTS) components including Micro-Electro-Mechanical System (MEMS) gyro and two simultaneously operating star trackers (STR) to enhance the measurement accuracy. The Multiplicative Extended Kalman Filter (MEKF) is used to estimate the attitude of CubeSat, and four kinds of attitude estimation layouts are put forward according to the idea of weighted average of two quaternions from two STR and different architectures of information fusion. Using the proposed methods, the attitude measurement unit can continuously provide accurate and reliable attitude knowledge for attitude control unit when the CubeSat is running in orbit. Numerical simulation is performed to verify the effectiveness of the proposed methods, and it offers a reference for CubeSat developers from the perspective of engineering application.


---

## Body

## 1. Introduction

CubeSat has become an interesting innovation in the space frontier due to advantages such as low cost, good flexibility, short development period, and high functional density [1–3]. The CubeSats used for imaging and videoing are mainstream among those that have been launched [4]. High-accuracy measurement and control are key to achieving the above tasks.Gyro and STR are widely used in satellite with a high demand of attitude control accuracy. The small satellite Flying Laptop adopts one STR and four fiber-optic gyros (FOG) to accomplish inertial attitude measurement when the satellite is working in earth observation mode [5]; a novel STR named NSTT with high-accuracy attitude measurement capacity (<4″) and high acquisition rate (>3°/s) is utilized in Super Low Altitude Test Satellite (SLATS) as a technology demonstration payload [6]; the nanosatellite TechnoSat uses one FOG to integrate angular rate to generate attitude knowledge when the satellite is working in attitude maneuver mode [7]; there are two STR vertical to each other in the small satellite Bispectral InfraRed Optical System (BIROS), one of them works as cold backup [8]. The above FOG and STR are mainly applied in nanosatellites or larger satellite platforms, although they have high enough measurement accuracy; however, they are not suitable for CubaSat platforms because of large size and high power consumption. For example, the size and power consumption of the STR used in Flying Laptop are 10×10×10 cm3 and 7.8 W, respectively, while the size and maximum power supply of our 2U CubeSat are 20×10×10 cm3 and 4.6 W [9, 10].In recent years, many scholars have carried out studies on the attitude measurement of satellites with gyro and STR integrated. Ref. [11] proposes Moving Horizon Filter (MHF) to estimate the attitude and gyro calibration parameters by jointly using one STR and one gyro and draws a conclusion that the proposed method results in an increased accuracy on nonlinear systems with respect to Extended Kalman Filter (EKF); in Ref. [12], attitude fusion data obtained by Complementary Filter (CF) according to different noise frequency characteristics of one gyro and one STR is introduced as the observed values of the Unscented Kalman Filter (UKF) system to the process of measurement update, finally improving the accuracy of attitude estimation compared with only UKF used; in Ref. [13], Cubature Kalman Filter (CKF) is utilized to estimate satellite attitude aiming at the satellite star sensor/gyro attitude measurement system; the results indicate that CKF is better than EKF and UKF in accuracy and stability when dealing with strongly nonlinear attitude kinematics and dynamics equations; Ref. [14] proposes Robust Adaptive Unscented Kalman Filter (RAUKF) to achieve precise spacecraft attitude estimation in the presence of measurement faults by using one STR and one gyro and proves that the proposed algorithm outperforms the standard EKF and UKF under measurement faults through error covariance adaptation. The above literature mainly focuses on the research of attitude measurement algorithms, although the results show that all the proposed algorithms perform better than EKF in theory; however, they cannot perform as well as we expect in practical application. Firstly, the kinematics equation of CubeSat is not so strongly nonlinear, so filters applied to strongly nonlinear systems such as UKF, CKF, and their variants cannot provide more accurate attitude knowledge than EKF actually in this case. Secondly, UKF, CKF, and their variants involve complicated mathematical calculation such as Cholesky and QR decomposition [15]; although we can invoke relevant functions in MATLAB directly just for numerical simulation, we have to rewrite the complex codes in other programming languages suitable for running in an on-board computer of CubeSat; moreover, a large amount of loop iteration and complicated function calculation in UKF, CKF, and their variants lead to more execution time than EKF; for example, the execution time of UKF is almost twice as much as EKF under the same simulation condition according to our comparison, which greatly increases computation burden of an on-board computer of CubeSat with limited capacity.Physically, the rotation angle about the line-of-sight axis of STR is less accurate than that about the other two axes; hence, the total attitude measurement accuracy is influenced when only one STR is used [16]. Therefore, some scholars have carried out research on attitude measurement by using multiple STR. Singular Value Decomposition (SVD) and Quaternion Estimator (QUEST) algorithms are proposed in Ref. [17] to estimate the attitude of satellite based on stellar vectors from multiple (usually more than two) simultaneously operating STR; Ref. [18] reviews a kind of gyroless attitude measurement method which utilizes the TRIAD algorithm to estimate the attitude angles of satellite based on stellar vectors from two simultaneously operating STR then calculates the angular rate of a satellite by a first-order difference; Federated Extended Kalman Filter (FEKF) and Federated Unscented Kalman Filter (FUKF) are, respectively, proposed in Ref. [19] and Ref. [20]; the former contains three independent local filters using measurements from one gyro and three simultaneously operating STR and one master filter while the latter contains two independent local filters using measurements from one gyro and two simultaneously operating STR and one master filter, and their respective estimated results are fused in the master filter to finally generate the attitude knowledge of satellite. The above SVD, QUEST, and TRIAD belong to single-frame methods; they do not deal with measurement errors of attitude sensors and use any knowledge about kinematics and dynamics of satellites compared with Kalman filter-based methods such as EKF, UKF, and CKF; hence, the Kalman filter-based methods generally provide more accurate estimates than the single-frame methods [21, 22]. For the large satellite equipped with high-accuracy attitude sensors, these algorithms can provide sufficiently accurate attitude knowledge; however, for the CubeSat equipped with relatively coarse attitude sensors, they are not good choices actually.Inspired by the idea of federated filters proposed in Ref. [19] and Ref. [20] and a weighted average of two quaternions proposed in Ref. [23], we propose a novel attitude measurement method by using one gyro and two simultaneously operating STR in this paper, as part of the development task of the CubeSat NJUST-2. On the hardware level, we adopt COTS components with the properties of small size, low power dissipation, and low cost as attitude sensors, which include one MEMS gyro and two simultaneously operating STR; on the algorithm level, we develop four kinds of attitude estimation layouts based on EKF according to the architecture of information fusion; they are called centralized layout I, centralized layout II, centralized layout III, and decentralized layout.The paper is organized as follows. Attitude kinematics of CubeSat and attitude sensor modeling are described in Section2. Four kinds of attitude estimation layouts and attitude estimation processes based on MEKF are given in Section 3. Numerical simulation, comparison, and analysis are conducted in Section 4. Finally, conclusions are presented in Section 5.


## 2. Attitude Kinematics and Sensor Modeling

### 2.1. Quaternion Kinematics

Quaternion is used to describe the attitude of CubeSat in our study due to its nonsingularity. The quaternion is defined as(1)q=q0qvΤΤ,where q0 is the scalar part; qv is the vector part, qv=q1q2q3Τ.The kinematic differential equation represented by quaternion is described as [14]
(2)q̇=12Ωωq,where ω is the angular rate vector; Ωω is represented as [14]
(3)Ωω=0−ωx−ωy−ωzωx0ωz‐ωyωy−ωz0ωxωzωy−ωx0.


### 2.2. Gyro Modeling

The gyro body axes are aligned with principal axes of inertia of the CubeSat body; it can measure the angular rate of the CubeSat body frame relative to the geocentric inertial frame. Considering the presence of bias, measurement noise, and random drift, the discrete mathematical model of the gyro is built as [24]
(4a)ω~bi=ω¯bi+12b¯k+b¯k‐1+Wω,(4b)b¯k+1=b¯k+Wb,where ω~bi and ω¯bi are the measured and true angular rate of the CubeSat body frame relative to the geocentric inertial frame, respectively; b¯k+1, b¯k, and b¯k‐1 are the true biases at time k+1, k, and k−1, respectively; Wω is the measurement noise, denoted as
(5)Wω=N2Ts+112K2Ts1/2η,where N and K are the angle random walk and rate random walk, respectively; Ts is the sampling period; η is the independent Gaussian white noise, η~N0,1.Wb is the bias noise, denoted as
(6)Wb=KTs1/2η.


### 2.3. STR Modeling

#### 2.3.1. Mounting and Layout

Two pico STR titled PST3 manufactured by TY-Space Technology are used in the attitude measurement unit of NJUST-2; its parameters are listed in Table1.Table 1
PST3 pico STR performance.
Dimension3.3×3.3×3.7 cm3Field of view (FOV)15°×12° (effective)Measurement accuracy<10″ (3σ, pointing)<70″ (3σ, rolling)Stray light suppression angle>35°Power consumption<0.6 WOutput formatQuaternionIn general, two STR are mounted at a tilted angle equal to 90° for the highest measurement accuracy and complete observability [25]. Figure 1 shows the on-orbit flight model of the CubeSat built in Systems Tool Kit (STK). The CubeSat will run in a sun-synchronous orbit with an inclination of 97.5°; A1, B1, B2, and C1 are the primary surfaces illuminated by sunlight; C2 is the secondary surface illuminated by sunlight; A2 is the surface illuminated by reflected light form atmosphere.Figure 1
On-orbit flight model of CubeSat in STK.Therefore, we consider mounting the two STR on C2 at a tilted angle equal to 90° for protecting them from stray light disturbance, as Figure2 shows.Figure 2
Mounting diagram of two STR.To verify the effectiveness of the above mounting method, relevant simulation and analysis are carried out by STK. Since the two STR are not exposed to sunlight during the eclipse area, we only consider the sunlight area. Orbit parameters in our study are from Two-Line Element (TLE) sets of our last CubeSat NJUST-1 whose orbit is similar to NJUST-2, published by North American Aerospace Defense Command (NORAD). The first available TLE is as follows:(1)
1 43156U 18008B 18039.44439993 .00000182 00000-0 15203-4 0 9998(2)
2 43156 97.5404 116.0398 0014251 219.2554 201.5169 15.09764844 3059The initial orbit epoch is 8 Feb 2018 10:39:56.154 Coordinated Universal Time Greenwich (UTCG). Figure3 illustrates the angles between sun vector and orbit plane titled beta angle within one year from the above date.Figure 3
Beta angles within one year.As Figure3 shows, the beta angle is minimal at noon on June 1st. The angle between boresight 1 and sun vector titled α1 and the angle between boresight 2 and sun vector titled α2 are also minimal at that time, which means sunlight has the greatest influence on the two STR. Therefore, we only need to consider the above worst condition. Selecting 1 Jun 2018 12:10:56.997 UTCG~1 Jun 2018 13:11:20.285 UTCG as simulation duration, we can obtain α1 and α2 within the above period, as Figure 4 shows.Figure 4
Angles between STR boresight and sun vector.As Figure3 shows, the minimums of α1 and α2 are both larger than the stray light suppression angle of the STR, which means that the mounting method of the two STR can effectively protect them from stray light disturbance during the CubeSat’s lifetime.


#### 2.3.2. Mathematical Model

The STR can measure attitude quaternion of the geocentric inertial frame relative to the STR body frame. Considering that the STR body axes are not aligned with principal axes of inertia of the CubeSat body and the presence of measurement noise, the discrete mathematical model of STR is built as(7)q~bi=q¯bi⊗Δq~bi,where q~bi and q¯bi are the measured and true attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame, respectively; Δq~bi is the converted measurement noise, denoted as
(8)Δq~bi=q¯bi−1⊗q~bi=q¯si⊗qbs−1⊗q~si⊗qbs=qbs−1⊗q¯si−1⊗q~si⊗qbs=qbs−1⊗Δq~si⊗qbs,where q~si and q¯si are the measured and true attitude quaternion of the geocentric inertial frame relative to the STR body frame, respectively; qbs is the attitude quaternion of the STR body frame relative to the CubeSat body frame, denoted as
(9a)qbs0=cosθbs2cosφbs2cosψbs2−sinθbs2sinφbs2sinψbs2,(9b)qbs1=cosθbs2sinφbs2cosψbs2−sinθbs2cosφbs2sinψbs2,(9c)qbs2=sinθbs2cosφbs2cosψbs2+cosθbs2sinφbs2sinψbs2,(9d)qbs3=sinθbs2sinφbs2cosψbs2+cosθbs2cosφbs2sinψbs2,where ψbs, φbs, and θbs are the Euler angles of the STR body frame relative to the CubeSat body frame.Δq~si is the original measurement noise, under the condition of small angle; it can be denoted as
(10)Δq~si=1−Δr2η+Δp2η+Δp2η2Δr2ηΔp2ηΔp2ηΤ,where Δr and Δp are the measurement errors in the directions of rolling and pointing, respectively.The measurement residual is denoted as(11)δq~bi=q∧bi−1⊗q~bi=q∧bi−1⊗q¯bi⊗Δq~bi=qoi⊗q∧bo−1⊗qoi⊗q¯bo⊗Δq~bi=q∧bo−1⊗qoi−1⊗qoi⊗q¯bo⊗Δq~bi=q∧bo−1⊗q¯bo⊗Δq~bi=Δq^bo⊗Δq~bi,where q^bi is the estimated attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame; qoi is the attitude quaternion of the geocentric inertial frame relative to the orbit frame; q^bo and q¯bo are the estimated and true attitude quaternion of the orbit frame relative to the CubeSat body frame, respectively; Δq^bo is the estimation error.Extracting the vector part ofδq~bi, the linear equation of measurement residual can be denoted as
(12)δq~bi_v=Δq^bo_v+Δq~bi_v,where δq~bi_v, Δq^bo_v, and Δq~bi_v are the vector part of δq~bi, Δq^bo, and Δq~bi, respectively.


## 2.1. Quaternion Kinematics

Quaternion is used to describe the attitude of CubeSat in our study due to its nonsingularity. The quaternion is defined as(1)q=q0qvΤΤ,where q0 is the scalar part; qv is the vector part, qv=q1q2q3Τ.The kinematic differential equation represented by quaternion is described as [14]
(2)q̇=12Ωωq,where ω is the angular rate vector; Ωω is represented as [14]
(3)Ωω=0−ωx−ωy−ωzωx0ωz‐ωyωy−ωz0ωxωzωy−ωx0.


## 2.2. Gyro Modeling

The gyro body axes are aligned with principal axes of inertia of the CubeSat body; it can measure the angular rate of the CubeSat body frame relative to the geocentric inertial frame. Considering the presence of bias, measurement noise, and random drift, the discrete mathematical model of the gyro is built as [24]
(4a)ω~bi=ω¯bi+12b¯k+b¯k‐1+Wω,(4b)b¯k+1=b¯k+Wb,where ω~bi and ω¯bi are the measured and true angular rate of the CubeSat body frame relative to the geocentric inertial frame, respectively; b¯k+1, b¯k, and b¯k‐1 are the true biases at time k+1, k, and k−1, respectively; Wω is the measurement noise, denoted as
(5)Wω=N2Ts+112K2Ts1/2η,where N and K are the angle random walk and rate random walk, respectively; Ts is the sampling period; η is the independent Gaussian white noise, η~N0,1.Wb is the bias noise, denoted as
(6)Wb=KTs1/2η.


## 2.3. STR Modeling

### 2.3.1. Mounting and Layout

Two pico STR titled PST3 manufactured by TY-Space Technology are used in the attitude measurement unit of NJUST-2; its parameters are listed in Table1.Table 1
PST3 pico STR performance.
Dimension3.3×3.3×3.7 cm3Field of view (FOV)15°×12° (effective)Measurement accuracy<10″ (3σ, pointing)<70″ (3σ, rolling)Stray light suppression angle>35°Power consumption<0.6 WOutput formatQuaternionIn general, two STR are mounted at a tilted angle equal to 90° for the highest measurement accuracy and complete observability [25]. Figure 1 shows the on-orbit flight model of the CubeSat built in Systems Tool Kit (STK). The CubeSat will run in a sun-synchronous orbit with an inclination of 97.5°; A1, B1, B2, and C1 are the primary surfaces illuminated by sunlight; C2 is the secondary surface illuminated by sunlight; A2 is the surface illuminated by reflected light form atmosphere.Figure 1
On-orbit flight model of CubeSat in STK.Therefore, we consider mounting the two STR on C2 at a tilted angle equal to 90° for protecting them from stray light disturbance, as Figure2 shows.Figure 2
Mounting diagram of two STR.To verify the effectiveness of the above mounting method, relevant simulation and analysis are carried out by STK. Since the two STR are not exposed to sunlight during the eclipse area, we only consider the sunlight area. Orbit parameters in our study are from Two-Line Element (TLE) sets of our last CubeSat NJUST-1 whose orbit is similar to NJUST-2, published by North American Aerospace Defense Command (NORAD). The first available TLE is as follows:(1)
1 43156U 18008B 18039.44439993 .00000182 00000-0 15203-4 0 9998(2)
2 43156 97.5404 116.0398 0014251 219.2554 201.5169 15.09764844 3059The initial orbit epoch is 8 Feb 2018 10:39:56.154 Coordinated Universal Time Greenwich (UTCG). Figure3 illustrates the angles between sun vector and orbit plane titled beta angle within one year from the above date.Figure 3
Beta angles within one year.As Figure3 shows, the beta angle is minimal at noon on June 1st. The angle between boresight 1 and sun vector titled α1 and the angle between boresight 2 and sun vector titled α2 are also minimal at that time, which means sunlight has the greatest influence on the two STR. Therefore, we only need to consider the above worst condition. Selecting 1 Jun 2018 12:10:56.997 UTCG~1 Jun 2018 13:11:20.285 UTCG as simulation duration, we can obtain α1 and α2 within the above period, as Figure 4 shows.Figure 4
Angles between STR boresight and sun vector.As Figure3 shows, the minimums of α1 and α2 are both larger than the stray light suppression angle of the STR, which means that the mounting method of the two STR can effectively protect them from stray light disturbance during the CubeSat’s lifetime.


### 2.3.2. Mathematical Model

The STR can measure attitude quaternion of the geocentric inertial frame relative to the STR body frame. Considering that the STR body axes are not aligned with principal axes of inertia of the CubeSat body and the presence of measurement noise, the discrete mathematical model of STR is built as(7)q~bi=q¯bi⊗Δq~bi,where q~bi and q¯bi are the measured and true attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame, respectively; Δq~bi is the converted measurement noise, denoted as
(8)Δq~bi=q¯bi−1⊗q~bi=q¯si⊗qbs−1⊗q~si⊗qbs=qbs−1⊗q¯si−1⊗q~si⊗qbs=qbs−1⊗Δq~si⊗qbs,where q~si and q¯si are the measured and true attitude quaternion of the geocentric inertial frame relative to the STR body frame, respectively; qbs is the attitude quaternion of the STR body frame relative to the CubeSat body frame, denoted as
(9a)qbs0=cosθbs2cosφbs2cosψbs2−sinθbs2sinφbs2sinψbs2,(9b)qbs1=cosθbs2sinφbs2cosψbs2−sinθbs2cosφbs2sinψbs2,(9c)qbs2=sinθbs2cosφbs2cosψbs2+cosθbs2sinφbs2sinψbs2,(9d)qbs3=sinθbs2sinφbs2cosψbs2+cosθbs2cosφbs2sinψbs2,where ψbs, φbs, and θbs are the Euler angles of the STR body frame relative to the CubeSat body frame.Δq~si is the original measurement noise, under the condition of small angle; it can be denoted as
(10)Δq~si=1−Δr2η+Δp2η+Δp2η2Δr2ηΔp2ηΔp2ηΤ,where Δr and Δp are the measurement errors in the directions of rolling and pointing, respectively.The measurement residual is denoted as(11)δq~bi=q∧bi−1⊗q~bi=q∧bi−1⊗q¯bi⊗Δq~bi=qoi⊗q∧bo−1⊗qoi⊗q¯bo⊗Δq~bi=q∧bo−1⊗qoi−1⊗qoi⊗q¯bo⊗Δq~bi=q∧bo−1⊗q¯bo⊗Δq~bi=Δq^bo⊗Δq~bi,where q^bi is the estimated attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame; qoi is the attitude quaternion of the geocentric inertial frame relative to the orbit frame; q^bo and q¯bo are the estimated and true attitude quaternion of the orbit frame relative to the CubeSat body frame, respectively; Δq^bo is the estimation error.Extracting the vector part ofδq~bi, the linear equation of measurement residual can be denoted as
(12)δq~bi_v=Δq^bo_v+Δq~bi_v,where δq~bi_v, Δq^bo_v, and Δq~bi_v are the vector part of δq~bi, Δq^bo, and Δq~bi, respectively.


## 2.3.1. Mounting and Layout

Two pico STR titled PST3 manufactured by TY-Space Technology are used in the attitude measurement unit of NJUST-2; its parameters are listed in Table1.Table 1
PST3 pico STR performance.
Dimension3.3×3.3×3.7 cm3Field of view (FOV)15°×12° (effective)Measurement accuracy<10″ (3σ, pointing)<70″ (3σ, rolling)Stray light suppression angle>35°Power consumption<0.6 WOutput formatQuaternionIn general, two STR are mounted at a tilted angle equal to 90° for the highest measurement accuracy and complete observability [25]. Figure 1 shows the on-orbit flight model of the CubeSat built in Systems Tool Kit (STK). The CubeSat will run in a sun-synchronous orbit with an inclination of 97.5°; A1, B1, B2, and C1 are the primary surfaces illuminated by sunlight; C2 is the secondary surface illuminated by sunlight; A2 is the surface illuminated by reflected light form atmosphere.Figure 1
On-orbit flight model of CubeSat in STK.Therefore, we consider mounting the two STR on C2 at a tilted angle equal to 90° for protecting them from stray light disturbance, as Figure2 shows.Figure 2
Mounting diagram of two STR.To verify the effectiveness of the above mounting method, relevant simulation and analysis are carried out by STK. Since the two STR are not exposed to sunlight during the eclipse area, we only consider the sunlight area. Orbit parameters in our study are from Two-Line Element (TLE) sets of our last CubeSat NJUST-1 whose orbit is similar to NJUST-2, published by North American Aerospace Defense Command (NORAD). The first available TLE is as follows:(1)
1 43156U 18008B 18039.44439993 .00000182 00000-0 15203-4 0 9998(2)
2 43156 97.5404 116.0398 0014251 219.2554 201.5169 15.09764844 3059The initial orbit epoch is 8 Feb 2018 10:39:56.154 Coordinated Universal Time Greenwich (UTCG). Figure3 illustrates the angles between sun vector and orbit plane titled beta angle within one year from the above date.Figure 3
Beta angles within one year.As Figure3 shows, the beta angle is minimal at noon on June 1st. The angle between boresight 1 and sun vector titled α1 and the angle between boresight 2 and sun vector titled α2 are also minimal at that time, which means sunlight has the greatest influence on the two STR. Therefore, we only need to consider the above worst condition. Selecting 1 Jun 2018 12:10:56.997 UTCG~1 Jun 2018 13:11:20.285 UTCG as simulation duration, we can obtain α1 and α2 within the above period, as Figure 4 shows.Figure 4
Angles between STR boresight and sun vector.As Figure3 shows, the minimums of α1 and α2 are both larger than the stray light suppression angle of the STR, which means that the mounting method of the two STR can effectively protect them from stray light disturbance during the CubeSat’s lifetime.


## 2.3.2. Mathematical Model

The STR can measure attitude quaternion of the geocentric inertial frame relative to the STR body frame. Considering that the STR body axes are not aligned with principal axes of inertia of the CubeSat body and the presence of measurement noise, the discrete mathematical model of STR is built as(7)q~bi=q¯bi⊗Δq~bi,where q~bi and q¯bi are the measured and true attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame, respectively; Δq~bi is the converted measurement noise, denoted as
(8)Δq~bi=q¯bi−1⊗q~bi=q¯si⊗qbs−1⊗q~si⊗qbs=qbs−1⊗q¯si−1⊗q~si⊗qbs=qbs−1⊗Δq~si⊗qbs,where q~si and q¯si are the measured and true attitude quaternion of the geocentric inertial frame relative to the STR body frame, respectively; qbs is the attitude quaternion of the STR body frame relative to the CubeSat body frame, denoted as
(9a)qbs0=cosθbs2cosφbs2cosψbs2−sinθbs2sinφbs2sinψbs2,(9b)qbs1=cosθbs2sinφbs2cosψbs2−sinθbs2cosφbs2sinψbs2,(9c)qbs2=sinθbs2cosφbs2cosψbs2+cosθbs2sinφbs2sinψbs2,(9d)qbs3=sinθbs2sinφbs2cosψbs2+cosθbs2cosφbs2sinψbs2,where ψbs, φbs, and θbs are the Euler angles of the STR body frame relative to the CubeSat body frame.Δq~si is the original measurement noise, under the condition of small angle; it can be denoted as
(10)Δq~si=1−Δr2η+Δp2η+Δp2η2Δr2ηΔp2ηΔp2ηΤ,where Δr and Δp are the measurement errors in the directions of rolling and pointing, respectively.The measurement residual is denoted as(11)δq~bi=q∧bi−1⊗q~bi=q∧bi−1⊗q¯bi⊗Δq~bi=qoi⊗q∧bo−1⊗qoi⊗q¯bo⊗Δq~bi=q∧bo−1⊗qoi−1⊗qoi⊗q¯bo⊗Δq~bi=q∧bo−1⊗q¯bo⊗Δq~bi=Δq^bo⊗Δq~bi,where q^bi is the estimated attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame; qoi is the attitude quaternion of the geocentric inertial frame relative to the orbit frame; q^bo and q¯bo are the estimated and true attitude quaternion of the orbit frame relative to the CubeSat body frame, respectively; Δq^bo is the estimation error.Extracting the vector part ofδq~bi, the linear equation of measurement residual can be denoted as
(12)δq~bi_v=Δq^bo_v+Δq~bi_v,where δq~bi_v, Δq^bo_v, and Δq~bi_v are the vector part of δq~bi, Δq^bo, and Δq~bi, respectively.


## 3. Attitude Estimation Algorithm

### 3.1. Attitude Estimation Layout

Due to the nonadditive characteristic of the quaternion, the common linear weighted average method is unable to directly calculate the average of multiple quaternions because it has the feature of nonuniqueness and cannot guarantee the normalization constraint of the quaternion. Aiming at the above problem, Ref. [23] proposes an effective nonlinear weighted average method to calculate the average of two quaternions as follows:
(13)qw=w1w1−w2+zzw1+w2+zq1+signq1Τq2w2w2−w1+zzw1+w2+zq2,where z=w1−w22+4w1w2q1Τq22; q1 and q2 are the two quaternions probably from two STR, Particle Filter (PF), multiple-model adaptive estimation, and other occasions [23]; w1 and w2 are the weights of q1 and q2, respectively.In terms of the architecture of information fusion, nonlinear filters are divided into centralized and decentralized types [26]. According to the above ideas, we develop four kinds of attitude estimation layouts based on EKF to estimate the attitude of CubeSat by using one gyro and two STR. They are centralized layout I, centralized layout II, centralized layout III, and decentralized layout, respectively.


#### 3.1.1. Centralized Layout I

Figure5 shows the structure diagram of the centralized layout I.Figure 5
The centralized layout I.q~bi_1 is the quaternion from STR 1; b^, b⌢, and Δb^ are the estimation, prediction, and estimation error of gyro bias, respectively; q⌢bo is the predicted attitude quaternion of the orbit frame relative to the CubeSat body frame; P⌢ and P^ are the prediction and estimation error covariance matrix, respectively.


#### 3.1.2. Centralized Layout II

Figure6 shows the structure diagram of the centralized layout II.Figure 6
The centralized layout II.q~bi_2 is the quaternion from STR 2; q~bi_w is the nonlinear weighted average of q~bi_1 and q~bi_2.


#### 3.1.3. Centralized Layout III

Figure7 shows the structure diagram of the centralized layout III.Figure 7
The centralized layout III.Δq^bo_a_v and Δq^bo_b_v are the vector parts of two estimation error quaternions from MEKF; Δq^bo_w_v is the linear weighted average of Δq^bo_a_v and Δq^bo_b_v.


#### 3.1.4. Decentralized Layout

Figure8 shows the structure diagram of the decentralized layout.Figure 8
The decentralized layout.Δb^1 and Δb^2 are the two bias estimation errors from MEKF 1 and MEKF 2, respectively; Δb^w is the linear weighted average of Δb^1 and Δb^2; Δq^bo_1_v and Δq^bo_2_v are the vector parts of two estimation error quaternions from MEKF 1 and MEKF 2, respectively; P^1 and P^2 are the prediction estimation error covariance matrices of Δb∧1TΔq∧bo_1_vTT and Δb∧2TΔq∧bo_2_vTT, respectively; P^w is the linear weighted average of P^1 andP^2.


### 3.2. Attitude Estimation Process

#### 3.2.1. Prediction of State

According to the attitude kinematics model (2), gyro bias model (4b), and estimation of state at time k−1, the prediction of state at time k can be obtained by
(14a)q⌢bo,k=q^bo,k‐1+12Ωω^bo,k‐1q^bo,k‐1Δt,(14b)b⌢k=b^k‐1,where Δt is the filter period; ω^bo,k‐1 is the estimated angular rate of the CubeSat body frame relative to the orbit frame, denoted as
(15)ω^bo,k‐1=ω~bi,k‐1−b^k−1−Tq^bo,k‐1ωoi,k‐1,where ωoi,k‐1 is the angular rate of the orbit frame relative to the geocentric inertial frame, denoted as
(16)ωoi,k‐1=0−ωo,k‐10T,where ωo,k‐1 is the mean orbit angular rate.Tq is the direction cosine matrix represented by the quaternion, denoted as
(17)Tq=2q02+q12−12q1q2+q0q32q1q3−q0q22q1q2−q0q32q02+q22−12q2q3+q0q12q1q3+q0q22q2q3−q0q12q02+q32−1.


#### 3.2.2. MEKF Equation

Since the estimation of state error is identically equal to zero at the initial time of each iteration [24], MEKF is simplified as follows.Step (1).
Prediction error covariance matrix(18)P⌢k=Φk/k‐1P^k‐1ΦTk/k‐1+Γk/k‐1Qk‐1ΓTk/k‐1,where Φk/k−1 is the state error transition matrix from time k−1 to time k; Γk/k−1 is the process noise input matrix; Qk−1 is the process noise variance matrix.
(1)
For the centralized layout I, centralized layout II, and decentralized layout,Φk/k−1 is denoted as(19)Φk/k‐1=I6×6+Fk−1Δt+Fk−12Δt22!+Fk−13Δt33!,where Fk−1 is the Jacobian matrix with respect to state error, denoted as
(20)Fk−1=03×303×3−12I3×3−ω^bi,k‐1×,where ω^bi,k‐1 is the estimated angular rate of the CubeSat body frame relative to the geocentric inertial frame, ω^bi,k‐1=ω~bi,k−1−b^k−1; ω× is denoted as
(21)ω×=0−ωzωyωz0−ωx−ωyωx0.
Γk/k−1 is denoted as
(22)Γk/k−1=I6×6Δt+Fk−112Δt2+Fk−1216Δt3+Fk−13124Δt4Gk−1,where Gk−1 is the Jacobian matrix with respect to process noise, denoted as
(23)Gk‐1=I3×303×303×3−12I3×3.
Qk−1 is denoted as
(24)Qk−1=K2TsI3×3−14K2Ts2I3×3−14K2Ts2I3×314N2Ts+112K2Ts3I3×3.(2)
For the centralized layout III,Φk/k−1 is denoted as(25)Φk/k‐1=I9×9+Fk−1Δt+Fk−12Δt22!+Fk−13Δt33!,where Fk−1 is denoted as
(26)Fk−1=03×303×303×3−12I3×3−ω^bi,k‐1×03×3−12I3×303×3−ω^bi,k‐1×.
Γk/k−1 is denoted as
(27)Γk/k−1=I9×9Δt+Fk−112Δt2+Fk−1216Δt3+Fk−13124Δt4Gk−1,where Gk−1 is denoted as
(28)Gk‐1=I3×303×303×303×3−12I3×303×303×303×3−12I3×3.
Qk−1 is denoted as
(29)Qk−1=K2TsI3×3−14K2Ts2I3×3−14K2Ts2I3×3−14K2Ts2I3×314N2Ts+112K2Ts3I3×3112K2Ts3I3×3−14K2Ts2I3×3112K2Ts3I3×314N2Ts+112K2Ts3I3×3.Step (2).
Kalman gain matrix(30)Kk=P⌢kHkTHkP⌢kHkT+Rk−1,where Hk is the state error measurement matrix;Rkis the measurement noise variance matrix.
(1)
For the centralized layout I,(31a)Hk=03×3I3×3,(31b)Rk=S1I3×3,where S1 is the variance of Δq~bi_1; Δq~bi_1 is the converted measurement noise of STR 1.
(2)
For the centralized layout II,(32a)Hk=03×3I3×3,(32b)Rk=SwI3×3,where Sw is the variance of Δq~bi_w; Δq~bi_w is the converted measurement noise of STR 1 combined with STR 2.
(3)
For the centralized layout III,(33a)Hk=03×3I3×303×303×303×3I3×3,(33b)Rk=S1I3×303×303×3S2I3×3,where S2 is the variance of Δq~bi_2; Δq~bi_2 is the converted measurement noise of STR 2.
(4)
For the decentralized layout,(34a)H1,k=03×3I3×3,(34b)H2,k=03×3I3×3,(34c)R1,k=S1I3×3,(34d)R2,k=S2I3×3.Step (3).
Estimation of state error(35)ΔX^k=KkδZk,where δZk is the integrated measurement residual.
(1)
For the centralized layout I,(36a)ΔX^k=Δb∧kTΔq∧bo_v,kTT,(36b)δZk=δq~bi_1_v,k,where δq~bi_1_v,k is the measurement residual of STR 1.
(2)
For the centralized layout II,(37a)ΔX^k=Δb∧kTΔq∧bo_v,kTT,(37b)δZk=δq~bi_w_v,k,where δq~bi_w_v,k is the measurement residual of STR 1 combined with STR 2.
(3)
For the centralized layout III,(38a)ΔX^k=Δb∧kTΔq∧bo_a_v,kTΔq∧bo_b_v,kTT,(38b)δZk=δq~bi_1_v,kδq~bi_2_v,kT,where δq~bi_2_v,k is the measurement residual of STR 2.
(4)
For the decentralized layout,(39a)ΔX^1,k=Δb∧1,kTΔq∧bo_1_v,kTT,(39b)ΔX^2,k=Δb∧2,kTΔq∧bo_2_v,kTT,(39c)δZ1,k=δq~bi_1_v,k,(39d)δZ2,k=δq~bi_2_v,k.Step (4).
Estimation error covariance matrix.(1)
For the centralized layout I and centralized layout II,(40)P^k=I6×6−KkHkP⌢kI6×6−KkHkT+KkRkKkT.(2)
For the centralized layout III,(41a)P^k=I9×9−KkHkP⌢kI9×9−KkHkT+KkRkKkT,(41b)P^w,k=P∧a,k−1+P∧b,k−1−1,where P^a,k and P^b,k are the estimation error covariance matrices of Δq^bo_a_v,k andΔq^bo_b_v,k, respectively, in P^k.
(3)
For the decentralized layout,(42a)P^1,k=I6×6−K1,kH1,kP⌢kI6×6−K1,kH1,kT+K1,kR1,kK1,kT,(42b)P^2,k=I6×6−K2,kH2,kP⌢kI6×6−K2,kH2,kT+K2,kR2,kK2,kT,(42c)P^w,k=P∧1,k−1+P∧2,k−1−1.


#### 3.2.3. Correction of State

(1)
For the centralized layout I and centralized layout II,(43a)q^bo,k=q⌢bo,k⊗1−Δq∧bo_v,k2Δq∧bo_v,kTT,(43b)b^k=b⌢k+Δb^k.(2)
For the centralized layout III,(44a)Δq^bo_ω_v,k=P^w,kP∧a,k−1Δq^bo_a_v,k+P∧b,k−1Δq^bo_b_v,k,(44b)q^bo,k=q⌢bo,k⊗1−Δq∧bo_ω_v,k2Δq∧bo_ω_v,kTT,(44c)b^k=b⌢k+Δb^k.(3)
For the decentralized layout,(45a)ΔX^w,k=P^w,kP∧1,k−1ΔX^1,k+P∧2,k−1ΔX^2,k(45b)Δq^bo_ω_v,k=ΔX^w,k4:6,(45c)q^bo,k=q⌢bo,k⊗1−Δq∧bo_ω_v,k2Δq∧bo_ω_v,kTT,(45d)Δb^w,k=ΔX^w,k1:3,(45e)b^k=b⌢k+Δb^w,k.


## 3.1. Attitude Estimation Layout

Due to the nonadditive characteristic of the quaternion, the common linear weighted average method is unable to directly calculate the average of multiple quaternions because it has the feature of nonuniqueness and cannot guarantee the normalization constraint of the quaternion. Aiming at the above problem, Ref. [23] proposes an effective nonlinear weighted average method to calculate the average of two quaternions as follows:
(13)qw=w1w1−w2+zzw1+w2+zq1+signq1Τq2w2w2−w1+zzw1+w2+zq2,where z=w1−w22+4w1w2q1Τq22; q1 and q2 are the two quaternions probably from two STR, Particle Filter (PF), multiple-model adaptive estimation, and other occasions [23]; w1 and w2 are the weights of q1 and q2, respectively.In terms of the architecture of information fusion, nonlinear filters are divided into centralized and decentralized types [26]. According to the above ideas, we develop four kinds of attitude estimation layouts based on EKF to estimate the attitude of CubeSat by using one gyro and two STR. They are centralized layout I, centralized layout II, centralized layout III, and decentralized layout, respectively.


### 3.1.1. Centralized Layout I

Figure5 shows the structure diagram of the centralized layout I.Figure 5
The centralized layout I.q~bi_1 is the quaternion from STR 1; b^, b⌢, and Δb^ are the estimation, prediction, and estimation error of gyro bias, respectively; q⌢bo is the predicted attitude quaternion of the orbit frame relative to the CubeSat body frame; P⌢ and P^ are the prediction and estimation error covariance matrix, respectively.


### 3.1.2. Centralized Layout II

Figure6 shows the structure diagram of the centralized layout II.Figure 6
The centralized layout II.q~bi_2 is the quaternion from STR 2; q~bi_w is the nonlinear weighted average of q~bi_1 and q~bi_2.


### 3.1.3. Centralized Layout III

Figure7 shows the structure diagram of the centralized layout III.Figure 7
The centralized layout III.Δq^bo_a_v and Δq^bo_b_v are the vector parts of two estimation error quaternions from MEKF; Δq^bo_w_v is the linear weighted average of Δq^bo_a_v and Δq^bo_b_v.


### 3.1.4. Decentralized Layout

Figure8 shows the structure diagram of the decentralized layout.Figure 8
The decentralized layout.Δb^1 and Δb^2 are the two bias estimation errors from MEKF 1 and MEKF 2, respectively; Δb^w is the linear weighted average of Δb^1 and Δb^2; Δq^bo_1_v and Δq^bo_2_v are the vector parts of two estimation error quaternions from MEKF 1 and MEKF 2, respectively; P^1 and P^2 are the prediction estimation error covariance matrices of Δb∧1TΔq∧bo_1_vTT and Δb∧2TΔq∧bo_2_vTT, respectively; P^w is the linear weighted average of P^1 andP^2.


## 3.1.1. Centralized Layout I

Figure5 shows the structure diagram of the centralized layout I.Figure 5
The centralized layout I.q~bi_1 is the quaternion from STR 1; b^, b⌢, and Δb^ are the estimation, prediction, and estimation error of gyro bias, respectively; q⌢bo is the predicted attitude quaternion of the orbit frame relative to the CubeSat body frame; P⌢ and P^ are the prediction and estimation error covariance matrix, respectively.


## 3.1.2. Centralized Layout II

Figure6 shows the structure diagram of the centralized layout II.Figure 6
The centralized layout II.q~bi_2 is the quaternion from STR 2; q~bi_w is the nonlinear weighted average of q~bi_1 and q~bi_2.


## 3.1.3. Centralized Layout III

Figure7 shows the structure diagram of the centralized layout III.Figure 7
The centralized layout III.Δq^bo_a_v and Δq^bo_b_v are the vector parts of two estimation error quaternions from MEKF; Δq^bo_w_v is the linear weighted average of Δq^bo_a_v and Δq^bo_b_v.


## 3.1.4. Decentralized Layout

Figure8 shows the structure diagram of the decentralized layout.Figure 8
The decentralized layout.Δb^1 and Δb^2 are the two bias estimation errors from MEKF 1 and MEKF 2, respectively; Δb^w is the linear weighted average of Δb^1 and Δb^2; Δq^bo_1_v and Δq^bo_2_v are the vector parts of two estimation error quaternions from MEKF 1 and MEKF 2, respectively; P^1 and P^2 are the prediction estimation error covariance matrices of Δb∧1TΔq∧bo_1_vTT and Δb∧2TΔq∧bo_2_vTT, respectively; P^w is the linear weighted average of P^1 andP^2.


## 3.2. Attitude Estimation Process

### 3.2.1. Prediction of State

According to the attitude kinematics model (2), gyro bias model (4b), and estimation of state at time k−1, the prediction of state at time k can be obtained by
(14a)q⌢bo,k=q^bo,k‐1+12Ωω^bo,k‐1q^bo,k‐1Δt,(14b)b⌢k=b^k‐1,where Δt is the filter period; ω^bo,k‐1 is the estimated angular rate of the CubeSat body frame relative to the orbit frame, denoted as
(15)ω^bo,k‐1=ω~bi,k‐1−b^k−1−Tq^bo,k‐1ωoi,k‐1,where ωoi,k‐1 is the angular rate of the orbit frame relative to the geocentric inertial frame, denoted as
(16)ωoi,k‐1=0−ωo,k‐10T,where ωo,k‐1 is the mean orbit angular rate.Tq is the direction cosine matrix represented by the quaternion, denoted as
(17)Tq=2q02+q12−12q1q2+q0q32q1q3−q0q22q1q2−q0q32q02+q22−12q2q3+q0q12q1q3+q0q22q2q3−q0q12q02+q32−1.


### 3.2.2. MEKF Equation

Since the estimation of state error is identically equal to zero at the initial time of each iteration [24], MEKF is simplified as follows.Step (1).
Prediction error covariance matrix(18)P⌢k=Φk/k‐1P^k‐1ΦTk/k‐1+Γk/k‐1Qk‐1ΓTk/k‐1,where Φk/k−1 is the state error transition matrix from time k−1 to time k; Γk/k−1 is the process noise input matrix; Qk−1 is the process noise variance matrix.
(1)
For the centralized layout I, centralized layout II, and decentralized layout,Φk/k−1 is denoted as(19)Φk/k‐1=I6×6+Fk−1Δt+Fk−12Δt22!+Fk−13Δt33!,where Fk−1 is the Jacobian matrix with respect to state error, denoted as
(20)Fk−1=03×303×3−12I3×3−ω^bi,k‐1×,where ω^bi,k‐1 is the estimated angular rate of the CubeSat body frame relative to the geocentric inertial frame, ω^bi,k‐1=ω~bi,k−1−b^k−1; ω× is denoted as
(21)ω×=0−ωzωyωz0−ωx−ωyωx0.
Γk/k−1 is denoted as
(22)Γk/k−1=I6×6Δt+Fk−112Δt2+Fk−1216Δt3+Fk−13124Δt4Gk−1,where Gk−1 is the Jacobian matrix with respect to process noise, denoted as
(23)Gk‐1=I3×303×303×3−12I3×3.
Qk−1 is denoted as
(24)Qk−1=K2TsI3×3−14K2Ts2I3×3−14K2Ts2I3×314N2Ts+112K2Ts3I3×3.(2)
For the centralized layout III,Φk/k−1 is denoted as(25)Φk/k‐1=I9×9+Fk−1Δt+Fk−12Δt22!+Fk−13Δt33!,where Fk−1 is denoted as
(26)Fk−1=03×303×303×3−12I3×3−ω^bi,k‐1×03×3−12I3×303×3−ω^bi,k‐1×.
Γk/k−1 is denoted as
(27)Γk/k−1=I9×9Δt+Fk−112Δt2+Fk−1216Δt3+Fk−13124Δt4Gk−1,where Gk−1 is denoted as
(28)Gk‐1=I3×303×303×303×3−12I3×303×303×303×3−12I3×3.
Qk−1 is denoted as
(29)Qk−1=K2TsI3×3−14K2Ts2I3×3−14K2Ts2I3×3−14K2Ts2I3×314N2Ts+112K2Ts3I3×3112K2Ts3I3×3−14K2Ts2I3×3112K2Ts3I3×314N2Ts+112K2Ts3I3×3.Step (2).
Kalman gain matrix(30)Kk=P⌢kHkTHkP⌢kHkT+Rk−1,where Hk is the state error measurement matrix;Rkis the measurement noise variance matrix.
(1)
For the centralized layout I,(31a)Hk=03×3I3×3,(31b)Rk=S1I3×3,where S1 is the variance of Δq~bi_1; Δq~bi_1 is the converted measurement noise of STR 1.
(2)
For the centralized layout II,(32a)Hk=03×3I3×3,(32b)Rk=SwI3×3,where Sw is the variance of Δq~bi_w; Δq~bi_w is the converted measurement noise of STR 1 combined with STR 2.
(3)
For the centralized layout III,(33a)Hk=03×3I3×303×303×303×3I3×3,(33b)Rk=S1I3×303×303×3S2I3×3,where S2 is the variance of Δq~bi_2; Δq~bi_2 is the converted measurement noise of STR 2.
(4)
For the decentralized layout,(34a)H1,k=03×3I3×3,(34b)H2,k=03×3I3×3,(34c)R1,k=S1I3×3,(34d)R2,k=S2I3×3.Step (3).
Estimation of state error(35)ΔX^k=KkδZk,where δZk is the integrated measurement residual.
(1)
For the centralized layout I,(36a)ΔX^k=Δb∧kTΔq∧bo_v,kTT,(36b)δZk=δq~bi_1_v,k,where δq~bi_1_v,k is the measurement residual of STR 1.
(2)
For the centralized layout II,(37a)ΔX^k=Δb∧kTΔq∧bo_v,kTT,(37b)δZk=δq~bi_w_v,k,where δq~bi_w_v,k is the measurement residual of STR 1 combined with STR 2.
(3)
For the centralized layout III,(38a)ΔX^k=Δb∧kTΔq∧bo_a_v,kTΔq∧bo_b_v,kTT,(38b)δZk=δq~bi_1_v,kδq~bi_2_v,kT,where δq~bi_2_v,k is the measurement residual of STR 2.
(4)
For the decentralized layout,(39a)ΔX^1,k=Δb∧1,kTΔq∧bo_1_v,kTT,(39b)ΔX^2,k=Δb∧2,kTΔq∧bo_2_v,kTT,(39c)δZ1,k=δq~bi_1_v,k,(39d)δZ2,k=δq~bi_2_v,k.Step (4).
Estimation error covariance matrix.(1)
For the centralized layout I and centralized layout II,(40)P^k=I6×6−KkHkP⌢kI6×6−KkHkT+KkRkKkT.(2)
For the centralized layout III,(41a)P^k=I9×9−KkHkP⌢kI9×9−KkHkT+KkRkKkT,(41b)P^w,k=P∧a,k−1+P∧b,k−1−1,where P^a,k and P^b,k are the estimation error covariance matrices of Δq^bo_a_v,k andΔq^bo_b_v,k, respectively, in P^k.
(3)
For the decentralized layout,(42a)P^1,k=I6×6−K1,kH1,kP⌢kI6×6−K1,kH1,kT+K1,kR1,kK1,kT,(42b)P^2,k=I6×6−K2,kH2,kP⌢kI6×6−K2,kH2,kT+K2,kR2,kK2,kT,(42c)P^w,k=P∧1,k−1+P∧2,k−1−1.


### 3.2.3. Correction of State

(1)
For the centralized layout I and centralized layout II,(43a)q^bo,k=q⌢bo,k⊗1−Δq∧bo_v,k2Δq∧bo_v,kTT,(43b)b^k=b⌢k+Δb^k.(2)
For the centralized layout III,(44a)Δq^bo_ω_v,k=P^w,kP∧a,k−1Δq^bo_a_v,k+P∧b,k−1Δq^bo_b_v,k,(44b)q^bo,k=q⌢bo,k⊗1−Δq∧bo_ω_v,k2Δq∧bo_ω_v,kTT,(44c)b^k=b⌢k+Δb^k.(3)
For the decentralized layout,(45a)ΔX^w,k=P^w,kP∧1,k−1ΔX^1,k+P∧2,k−1ΔX^2,k(45b)Δq^bo_ω_v,k=ΔX^w,k4:6,(45c)q^bo,k=q⌢bo,k⊗1−Δq∧bo_ω_v,k2Δq∧bo_ω_v,kTT,(45d)Δb^w,k=ΔX^w,k1:3,(45e)b^k=b⌢k+Δb^w,k.


## 3.2.1. Prediction of State

According to the attitude kinematics model (2), gyro bias model (4b), and estimation of state at time k−1, the prediction of state at time k can be obtained by
(14a)q⌢bo,k=q^bo,k‐1+12Ωω^bo,k‐1q^bo,k‐1Δt,(14b)b⌢k=b^k‐1,where Δt is the filter period; ω^bo,k‐1 is the estimated angular rate of the CubeSat body frame relative to the orbit frame, denoted as
(15)ω^bo,k‐1=ω~bi,k‐1−b^k−1−Tq^bo,k‐1ωoi,k‐1,where ωoi,k‐1 is the angular rate of the orbit frame relative to the geocentric inertial frame, denoted as
(16)ωoi,k‐1=0−ωo,k‐10T,where ωo,k‐1 is the mean orbit angular rate.Tq is the direction cosine matrix represented by the quaternion, denoted as
(17)Tq=2q02+q12−12q1q2+q0q32q1q3−q0q22q1q2−q0q32q02+q22−12q2q3+q0q12q1q3+q0q22q2q3−q0q12q02+q32−1.


## 3.2.2. MEKF Equation

Since the estimation of state error is identically equal to zero at the initial time of each iteration [24], MEKF is simplified as follows.Step (1).
Prediction error covariance matrix(18)P⌢k=Φk/k‐1P^k‐1ΦTk/k‐1+Γk/k‐1Qk‐1ΓTk/k‐1,where Φk/k−1 is the state error transition matrix from time k−1 to time k; Γk/k−1 is the process noise input matrix; Qk−1 is the process noise variance matrix.
(1)
For the centralized layout I, centralized layout II, and decentralized layout,Φk/k−1 is denoted as(19)Φk/k‐1=I6×6+Fk−1Δt+Fk−12Δt22!+Fk−13Δt33!,where Fk−1 is the Jacobian matrix with respect to state error, denoted as
(20)Fk−1=03×303×3−12I3×3−ω^bi,k‐1×,where ω^bi,k‐1 is the estimated angular rate of the CubeSat body frame relative to the geocentric inertial frame, ω^bi,k‐1=ω~bi,k−1−b^k−1; ω× is denoted as
(21)ω×=0−ωzωyωz0−ωx−ωyωx0.
Γk/k−1 is denoted as
(22)Γk/k−1=I6×6Δt+Fk−112Δt2+Fk−1216Δt3+Fk−13124Δt4Gk−1,where Gk−1 is the Jacobian matrix with respect to process noise, denoted as
(23)Gk‐1=I3×303×303×3−12I3×3.
Qk−1 is denoted as
(24)Qk−1=K2TsI3×3−14K2Ts2I3×3−14K2Ts2I3×314N2Ts+112K2Ts3I3×3.(2)
For the centralized layout III,Φk/k−1 is denoted as(25)Φk/k‐1=I9×9+Fk−1Δt+Fk−12Δt22!+Fk−13Δt33!,where Fk−1 is denoted as
(26)Fk−1=03×303×303×3−12I3×3−ω^bi,k‐1×03×3−12I3×303×3−ω^bi,k‐1×.
Γk/k−1 is denoted as
(27)Γk/k−1=I9×9Δt+Fk−112Δt2+Fk−1216Δt3+Fk−13124Δt4Gk−1,where Gk−1 is denoted as
(28)Gk‐1=I3×303×303×303×3−12I3×303×303×303×3−12I3×3.
Qk−1 is denoted as
(29)Qk−1=K2TsI3×3−14K2Ts2I3×3−14K2Ts2I3×3−14K2Ts2I3×314N2Ts+112K2Ts3I3×3112K2Ts3I3×3−14K2Ts2I3×3112K2Ts3I3×314N2Ts+112K2Ts3I3×3.Step (2).
Kalman gain matrix(30)Kk=P⌢kHkTHkP⌢kHkT+Rk−1,where Hk is the state error measurement matrix;Rkis the measurement noise variance matrix.
(1)
For the centralized layout I,(31a)Hk=03×3I3×3,(31b)Rk=S1I3×3,where S1 is the variance of Δq~bi_1; Δq~bi_1 is the converted measurement noise of STR 1.
(2)
For the centralized layout II,(32a)Hk=03×3I3×3,(32b)Rk=SwI3×3,where Sw is the variance of Δq~bi_w; Δq~bi_w is the converted measurement noise of STR 1 combined with STR 2.
(3)
For the centralized layout III,(33a)Hk=03×3I3×303×303×303×3I3×3,(33b)Rk=S1I3×303×303×3S2I3×3,where S2 is the variance of Δq~bi_2; Δq~bi_2 is the converted measurement noise of STR 2.
(4)
For the decentralized layout,(34a)H1,k=03×3I3×3,(34b)H2,k=03×3I3×3,(34c)R1,k=S1I3×3,(34d)R2,k=S2I3×3.Step (3).
Estimation of state error(35)ΔX^k=KkδZk,where δZk is the integrated measurement residual.
(1)
For the centralized layout I,(36a)ΔX^k=Δb∧kTΔq∧bo_v,kTT,(36b)δZk=δq~bi_1_v,k,where δq~bi_1_v,k is the measurement residual of STR 1.
(2)
For the centralized layout II,(37a)ΔX^k=Δb∧kTΔq∧bo_v,kTT,(37b)δZk=δq~bi_w_v,k,where δq~bi_w_v,k is the measurement residual of STR 1 combined with STR 2.
(3)
For the centralized layout III,(38a)ΔX^k=Δb∧kTΔq∧bo_a_v,kTΔq∧bo_b_v,kTT,(38b)δZk=δq~bi_1_v,kδq~bi_2_v,kT,where δq~bi_2_v,k is the measurement residual of STR 2.
(4)
For the decentralized layout,(39a)ΔX^1,k=Δb∧1,kTΔq∧bo_1_v,kTT,(39b)ΔX^2,k=Δb∧2,kTΔq∧bo_2_v,kTT,(39c)δZ1,k=δq~bi_1_v,k,(39d)δZ2,k=δq~bi_2_v,k.Step (4).
Estimation error covariance matrix.(1)
For the centralized layout I and centralized layout II,(40)P^k=I6×6−KkHkP⌢kI6×6−KkHkT+KkRkKkT.(2)
For the centralized layout III,(41a)P^k=I9×9−KkHkP⌢kI9×9−KkHkT+KkRkKkT,(41b)P^w,k=P∧a,k−1+P∧b,k−1−1,where P^a,k and P^b,k are the estimation error covariance matrices of Δq^bo_a_v,k andΔq^bo_b_v,k, respectively, in P^k.
(3)
For the decentralized layout,(42a)P^1,k=I6×6−K1,kH1,kP⌢kI6×6−K1,kH1,kT+K1,kR1,kK1,kT,(42b)P^2,k=I6×6−K2,kH2,kP⌢kI6×6−K2,kH2,kT+K2,kR2,kK2,kT,(42c)P^w,k=P∧1,k−1+P∧2,k−1−1.


## 3.2.3. Correction of State

(1)
For the centralized layout I and centralized layout II,(43a)q^bo,k=q⌢bo,k⊗1−Δq∧bo_v,k2Δq∧bo_v,kTT,(43b)b^k=b⌢k+Δb^k.(2)
For the centralized layout III,(44a)Δq^bo_ω_v,k=P^w,kP∧a,k−1Δq^bo_a_v,k+P∧b,k−1Δq^bo_b_v,k,(44b)q^bo,k=q⌢bo,k⊗1−Δq∧bo_ω_v,k2Δq∧bo_ω_v,kTT,(44c)b^k=b⌢k+Δb^k.(3)
For the decentralized layout,(45a)ΔX^w,k=P^w,kP∧1,k−1ΔX^1,k+P∧2,k−1ΔX^2,k(45b)Δq^bo_ω_v,k=ΔX^w,k4:6,(45c)q^bo,k=q⌢bo,k⊗1−Δq∧bo_ω_v,k2Δq∧bo_ω_v,kTT,(45d)Δb^w,k=ΔX^w,k1:3,(45e)b^k=b⌢k+Δb^w,k.


## 4. Simulation and Analysis

### 4.1. Simulation Parameters

#### 4.1.1. Navigation Parameters

The TLE for numerical simulation is as follows.(1)
43156U 18008B 18116.45379836 .00001071 00000-0 69869-4 0 9999(2)
43156 97.5309 191.7767 0014159 325.0023 91.3114 15.09925864 14673The position and velocity of CubeSat at any moment are from the Simplified General Perturbations 4 (SGP4) propagator in MATLAB based on the above TLE [27].


#### 4.1.2. Gyro Parameters

Initial true biasb¯0=101010T°/h; initial estimated bias b^0=000T°/h;N=0.48°/h1/2; K=120.34°/h3/2; Ts=1s.


#### 4.1.3. STR Parameters

Δr1=70″, Δp1=10″, ψbs_1=−90°, φbs_1=45°, θbs_1=−90°, w1=0.5, S1=1.44×10‐8;Δr2=70″, Δp2=10″, ψbs_2=90°, φbs_2=45°, θbs_2=90°, w2=0.5, S2=1.44×10‐8;Sw=7.4×10‐9.


#### 4.1.4. Initial Attitude Knowledge

True Euler anglea¯bo,0=5‐55°; estimated Euler angle a^bo,0=000°;q¯bo,0=0.9972−0.04540.04160.0416T; q^bo,0=1000T; true Euler angular rate r¯bo,0=−111°/s; estimated Euler angular rate r^bo,0=000°/s.


#### 4.1.5. Initial Estimation Error Covariance Matrix

For the centralized layout I, centralized layout II, and decentralized layout,P^0=I6×6; for the centralized layout III, P^0=I9×9.


#### 4.1.6. True State Update

q¯bo,k+1 used for comparison with estimated attitude quaternion is updated by
(46)q¯bo,k+1=q¯bo,k+12Ωω¯bo,kq¯bo,kΔt,where ω¯bo,k is the true angular rate of the CubeSat body frame relative to the orbit frame, denoted as
(47)ω¯bo,k=ω¯bi,k−Tq¯bo,kωoi,k,where ω¯bi,k is denoted as
(48)ω¯bi,k=ω¯bi,k‐1+I−1ω¯bi,k‐1×Iω¯bi,k‐1+hr,k‐1+Tc,k‐1+Tm,k‐1+Ta,k‐1+Tg,k‐1,where I is the moment of inertia of CubeSat; its value is
(49)I=0.003735‐0.000001755‐0.00004026‐0.0000017550.012820.00001191‐0.000040260.000011910.01273.hr,k‐1 is the angular momentum of reaction wheels; Tc,k‐1 is the control torque derived from zero-momentum attitude control [28]; Tm,k‐1, Ta,k‐1, and Tg,k‐1 are the disturbance torques derived from residual magnetism, atmospheric drag, and gravity gradient, respectively.b¯k+1 used for comparison with estimated gyro bias is updated by
(50)b¯k+1=b¯k+KTs1/2η.


#### 4.1.7. Simulation Time

Start timet0=26Apr201810:53:28.178UTCG; duration td=5000s; Δt=1s.


#### 4.1.8. Simulation Condition

Intel Core i5-4590 3.3 GHz, dual-core, RAM 8GB; MATLAB R2015a.


### 4.2. Simulation Results 1

The simulation and performance comparison of the above four kinds of attitude estimation layouts are carried out in terms of accuracy, convergence, and execution time.


#### 4.2.1. Estimation Accuracy of State

Figure9 shows the estimation error of Euler angles in a stable phase, and Table 2 lists the RMSE and MAE of the above data. We can see that the RMSE and MAE of centralized layout II are all smallest among the four layouts, which means that the centralized layout II gives the most accurate attitude knowledge with an accuracy of 0.01989° (3σ), 0.02409° (3σ), and 0.02472° (3σ) in the direction of yaw, roll, and pitch, respectively.Figure 9
Estimation error of Euler angles (stable phase).Table 2
RMSE and MAE of estimation error of Euler angles.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutEvaluation indexYawRMSE (°)0.007810.006330.006950.02085MAE (°)0.006240.005060.005530.01677RollRMSE (°)0.010240.008030.008670.02266MAE (°)0.008160.006410.006890.01809PitchRMSE (°)0.010340.008240.008520.02120MAE (°)0.008270.006570.006800.01684Figure10 shows the estimation error of gyro biases in a stable phase, and Table 3 lists the RMSE and MAE of the above data. We can see that the RMSE and MAE of centralized layout II are all smallest among the four layouts, which means that the centralized layout II gives the most accurate bias knowledge with an accuracy of 0.00711°/s (3σ), 0.00705°/s (3σ), and 0.00693°/s (3σ) in the direction of X, Y, and Z, respectively.Figure 10
Estimation error of gyro biases (stable phase).Table 3
RMSE and MAE of estimation error of gyro biases.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutEvaluation indexXRMSE (°/s)0.002510.002370.002400.00321MAE (°/s)0.002030.001890.001910.00249YRMSE (°/s)0.002470.002350.002420.00283MAE (°/s)0.001970.001900.001930.00226ZRMSE (°/s)0.002410.002310.002370.00284MAE (°/s)0.001900.001860.001870.00234


#### 4.2.2. Convergence Speed

Figure11 shows the estimation error of Euler angles in the initial phase, and Table 4 lists the convergence time. We can see that the convergence time of four kinds of attitude estimation layouts is almost the same and the convergence speed is fast enough.Figure 11
Estimation error of Euler angles (initial phase).Table 4
Convergence time.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutConvergence time (s)1111


#### 4.2.3. Execution Time

Table5 lists the execution time of 5000 iterations under the same simulation condition. We can see that the execution time of centralized layout III is the longest among the four layouts.Table 5
Execution time.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutExecution time (s)22.6623.1129.3823.91


### 4.3. Simulation Results 2

Furthermore, we also conduct the simulation under the same conditions for attitude measurements based on FEKF by using three STR presented in Ref. [19] and FUKF by using two STR presented in Ref. [20], respectively; their results are compared with the decentralized layout.


#### 4.3.1. Estimation Accuracy of State

Figure12 shows the estimation error of Euler angles in a stable phase, and Table 6 lists the RMSE and MAE of the above data. We can see that the attitude estimation accuracy of the decentralized layout is a little lower than FUKF with two STR but higher than FEKF with three STR.Figure 12
Estimation error of Euler angles (stable phase).Table 6
RMSE and MAE of estimation error of Euler angles.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STREvaluation indexYawRMSE (°)0.020850.019650.03805MAE (°)0.016770.015960.02970RollRMSE (°)0.022660.021820.04446MAE (°)0.018090.016920.03585PitchRMSE (°)0.021200.020230.04373MAE (°)0.016840.016440.03505Figure13 shows the estimation error of gyro biases in a stable phase, and Table 7 lists the RMSE and MAE of the above data. We can see that the bias estimation accuracy of the decentralized layout is a little lower than FUKF with two STR but higher than FEKF with three STR.Figure 13
Estimation error of gyro biases (stable phase).Table 7
RMSE and MAE of estimation error of gyro biases.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STREvaluation indexXRMSE (°/s)0.003210.002950.00389MAE (°/s)0.002490.002170.00315YRMSE (°/s)0.002830.002770.00396MAE (°/s)0.002260.002160.00323ZRMSE (°/s)0.002840.002840.00329MAE (°/s)0.002340.002320.00259


#### 4.3.2. Execution Time

Table8 lists the execution time of 5000 iterations under the same simulation condition. We can see that the execution time of FUKF with two STR is the longest among the three layouts, and it is more than twice that of the decentralized layout.Table 8
Execution time.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STRExecution time (s)23.9150.9930.28


## 4.1. Simulation Parameters

### 4.1.1. Navigation Parameters

The TLE for numerical simulation is as follows.(1)
43156U 18008B 18116.45379836 .00001071 00000-0 69869-4 0 9999(2)
43156 97.5309 191.7767 0014159 325.0023 91.3114 15.09925864 14673The position and velocity of CubeSat at any moment are from the Simplified General Perturbations 4 (SGP4) propagator in MATLAB based on the above TLE [27].


### 4.1.2. Gyro Parameters

Initial true biasb¯0=101010T°/h; initial estimated bias b^0=000T°/h;N=0.48°/h1/2; K=120.34°/h3/2; Ts=1s.


### 4.1.3. STR Parameters

Δr1=70″, Δp1=10″, ψbs_1=−90°, φbs_1=45°, θbs_1=−90°, w1=0.5, S1=1.44×10‐8;Δr2=70″, Δp2=10″, ψbs_2=90°, φbs_2=45°, θbs_2=90°, w2=0.5, S2=1.44×10‐8;Sw=7.4×10‐9.


### 4.1.4. Initial Attitude Knowledge

True Euler anglea¯bo,0=5‐55°; estimated Euler angle a^bo,0=000°;q¯bo,0=0.9972−0.04540.04160.0416T; q^bo,0=1000T; true Euler angular rate r¯bo,0=−111°/s; estimated Euler angular rate r^bo,0=000°/s.


### 4.1.5. Initial Estimation Error Covariance Matrix

For the centralized layout I, centralized layout II, and decentralized layout,P^0=I6×6; for the centralized layout III, P^0=I9×9.


### 4.1.6. True State Update

q¯bo,k+1 used for comparison with estimated attitude quaternion is updated by
(46)q¯bo,k+1=q¯bo,k+12Ωω¯bo,kq¯bo,kΔt,where ω¯bo,k is the true angular rate of the CubeSat body frame relative to the orbit frame, denoted as
(47)ω¯bo,k=ω¯bi,k−Tq¯bo,kωoi,k,where ω¯bi,k is denoted as
(48)ω¯bi,k=ω¯bi,k‐1+I−1ω¯bi,k‐1×Iω¯bi,k‐1+hr,k‐1+Tc,k‐1+Tm,k‐1+Ta,k‐1+Tg,k‐1,where I is the moment of inertia of CubeSat; its value is
(49)I=0.003735‐0.000001755‐0.00004026‐0.0000017550.012820.00001191‐0.000040260.000011910.01273.hr,k‐1 is the angular momentum of reaction wheels; Tc,k‐1 is the control torque derived from zero-momentum attitude control [28]; Tm,k‐1, Ta,k‐1, and Tg,k‐1 are the disturbance torques derived from residual magnetism, atmospheric drag, and gravity gradient, respectively.b¯k+1 used for comparison with estimated gyro bias is updated by
(50)b¯k+1=b¯k+KTs1/2η.


### 4.1.7. Simulation Time

Start timet0=26Apr201810:53:28.178UTCG; duration td=5000s; Δt=1s.


### 4.1.8. Simulation Condition

Intel Core i5-4590 3.3 GHz, dual-core, RAM 8GB; MATLAB R2015a.


## 4.1.1. Navigation Parameters

The TLE for numerical simulation is as follows.(1)
43156U 18008B 18116.45379836 .00001071 00000-0 69869-4 0 9999(2)
43156 97.5309 191.7767 0014159 325.0023 91.3114 15.09925864 14673The position and velocity of CubeSat at any moment are from the Simplified General Perturbations 4 (SGP4) propagator in MATLAB based on the above TLE [27].


## 4.1.2. Gyro Parameters

Initial true biasb¯0=101010T°/h; initial estimated bias b^0=000T°/h;N=0.48°/h1/2; K=120.34°/h3/2; Ts=1s.


## 4.1.3. STR Parameters

Δr1=70″, Δp1=10″, ψbs_1=−90°, φbs_1=45°, θbs_1=−90°, w1=0.5, S1=1.44×10‐8;Δr2=70″, Δp2=10″, ψbs_2=90°, φbs_2=45°, θbs_2=90°, w2=0.5, S2=1.44×10‐8;Sw=7.4×10‐9.


## 4.1.4. Initial Attitude Knowledge

True Euler anglea¯bo,0=5‐55°; estimated Euler angle a^bo,0=000°;q¯bo,0=0.9972−0.04540.04160.0416T; q^bo,0=1000T; true Euler angular rate r¯bo,0=−111°/s; estimated Euler angular rate r^bo,0=000°/s.


## 4.1.5. Initial Estimation Error Covariance Matrix

For the centralized layout I, centralized layout II, and decentralized layout,P^0=I6×6; for the centralized layout III, P^0=I9×9.


## 4.1.6. True State Update

q¯bo,k+1 used for comparison with estimated attitude quaternion is updated by
(46)q¯bo,k+1=q¯bo,k+12Ωω¯bo,kq¯bo,kΔt,where ω¯bo,k is the true angular rate of the CubeSat body frame relative to the orbit frame, denoted as
(47)ω¯bo,k=ω¯bi,k−Tq¯bo,kωoi,k,where ω¯bi,k is denoted as
(48)ω¯bi,k=ω¯bi,k‐1+I−1ω¯bi,k‐1×Iω¯bi,k‐1+hr,k‐1+Tc,k‐1+Tm,k‐1+Ta,k‐1+Tg,k‐1,where I is the moment of inertia of CubeSat; its value is
(49)I=0.003735‐0.000001755‐0.00004026‐0.0000017550.012820.00001191‐0.000040260.000011910.01273.hr,k‐1 is the angular momentum of reaction wheels; Tc,k‐1 is the control torque derived from zero-momentum attitude control [28]; Tm,k‐1, Ta,k‐1, and Tg,k‐1 are the disturbance torques derived from residual magnetism, atmospheric drag, and gravity gradient, respectively.b¯k+1 used for comparison with estimated gyro bias is updated by
(50)b¯k+1=b¯k+KTs1/2η.


## 4.1.7. Simulation Time

Start timet0=26Apr201810:53:28.178UTCG; duration td=5000s; Δt=1s.


## 4.1.8. Simulation Condition

Intel Core i5-4590 3.3 GHz, dual-core, RAM 8GB; MATLAB R2015a.


## 4.2. Simulation Results 1

The simulation and performance comparison of the above four kinds of attitude estimation layouts are carried out in terms of accuracy, convergence, and execution time.


### 4.2.1. Estimation Accuracy of State

Figure9 shows the estimation error of Euler angles in a stable phase, and Table 2 lists the RMSE and MAE of the above data. We can see that the RMSE and MAE of centralized layout II are all smallest among the four layouts, which means that the centralized layout II gives the most accurate attitude knowledge with an accuracy of 0.01989° (3σ), 0.02409° (3σ), and 0.02472° (3σ) in the direction of yaw, roll, and pitch, respectively.Figure 9
Estimation error of Euler angles (stable phase).Table 2
RMSE and MAE of estimation error of Euler angles.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutEvaluation indexYawRMSE (°)0.007810.006330.006950.02085MAE (°)0.006240.005060.005530.01677RollRMSE (°)0.010240.008030.008670.02266MAE (°)0.008160.006410.006890.01809PitchRMSE (°)0.010340.008240.008520.02120MAE (°)0.008270.006570.006800.01684Figure10 shows the estimation error of gyro biases in a stable phase, and Table 3 lists the RMSE and MAE of the above data. We can see that the RMSE and MAE of centralized layout II are all smallest among the four layouts, which means that the centralized layout II gives the most accurate bias knowledge with an accuracy of 0.00711°/s (3σ), 0.00705°/s (3σ), and 0.00693°/s (3σ) in the direction of X, Y, and Z, respectively.Figure 10
Estimation error of gyro biases (stable phase).Table 3
RMSE and MAE of estimation error of gyro biases.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutEvaluation indexXRMSE (°/s)0.002510.002370.002400.00321MAE (°/s)0.002030.001890.001910.00249YRMSE (°/s)0.002470.002350.002420.00283MAE (°/s)0.001970.001900.001930.00226ZRMSE (°/s)0.002410.002310.002370.00284MAE (°/s)0.001900.001860.001870.00234


### 4.2.2. Convergence Speed

Figure11 shows the estimation error of Euler angles in the initial phase, and Table 4 lists the convergence time. We can see that the convergence time of four kinds of attitude estimation layouts is almost the same and the convergence speed is fast enough.Figure 11
Estimation error of Euler angles (initial phase).Table 4
Convergence time.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutConvergence time (s)1111


### 4.2.3. Execution Time

Table5 lists the execution time of 5000 iterations under the same simulation condition. We can see that the execution time of centralized layout III is the longest among the four layouts.Table 5
Execution time.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutExecution time (s)22.6623.1129.3823.91


## 4.2.1. Estimation Accuracy of State

Figure9 shows the estimation error of Euler angles in a stable phase, and Table 2 lists the RMSE and MAE of the above data. We can see that the RMSE and MAE of centralized layout II are all smallest among the four layouts, which means that the centralized layout II gives the most accurate attitude knowledge with an accuracy of 0.01989° (3σ), 0.02409° (3σ), and 0.02472° (3σ) in the direction of yaw, roll, and pitch, respectively.Figure 9
Estimation error of Euler angles (stable phase).Table 2
RMSE and MAE of estimation error of Euler angles.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutEvaluation indexYawRMSE (°)0.007810.006330.006950.02085MAE (°)0.006240.005060.005530.01677RollRMSE (°)0.010240.008030.008670.02266MAE (°)0.008160.006410.006890.01809PitchRMSE (°)0.010340.008240.008520.02120MAE (°)0.008270.006570.006800.01684Figure10 shows the estimation error of gyro biases in a stable phase, and Table 3 lists the RMSE and MAE of the above data. We can see that the RMSE and MAE of centralized layout II are all smallest among the four layouts, which means that the centralized layout II gives the most accurate bias knowledge with an accuracy of 0.00711°/s (3σ), 0.00705°/s (3σ), and 0.00693°/s (3σ) in the direction of X, Y, and Z, respectively.Figure 10
Estimation error of gyro biases (stable phase).Table 3
RMSE and MAE of estimation error of gyro biases.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutEvaluation indexXRMSE (°/s)0.002510.002370.002400.00321MAE (°/s)0.002030.001890.001910.00249YRMSE (°/s)0.002470.002350.002420.00283MAE (°/s)0.001970.001900.001930.00226ZRMSE (°/s)0.002410.002310.002370.00284MAE (°/s)0.001900.001860.001870.00234


## 4.2.2. Convergence Speed

Figure11 shows the estimation error of Euler angles in the initial phase, and Table 4 lists the convergence time. We can see that the convergence time of four kinds of attitude estimation layouts is almost the same and the convergence speed is fast enough.Figure 11
Estimation error of Euler angles (initial phase).Table 4
Convergence time.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutConvergence time (s)1111


## 4.2.3. Execution Time

Table5 lists the execution time of 5000 iterations under the same simulation condition. We can see that the execution time of centralized layout III is the longest among the four layouts.Table 5
Execution time.
Attitude estimation layoutsCentralized layout ICentralized layout IICentralized layout IIIDecentralized layoutExecution time (s)22.6623.1129.3823.91


## 4.3. Simulation Results 2

Furthermore, we also conduct the simulation under the same conditions for attitude measurements based on FEKF by using three STR presented in Ref. [19] and FUKF by using two STR presented in Ref. [20], respectively; their results are compared with the decentralized layout.


### 4.3.1. Estimation Accuracy of State

Figure12 shows the estimation error of Euler angles in a stable phase, and Table 6 lists the RMSE and MAE of the above data. We can see that the attitude estimation accuracy of the decentralized layout is a little lower than FUKF with two STR but higher than FEKF with three STR.Figure 12
Estimation error of Euler angles (stable phase).Table 6
RMSE and MAE of estimation error of Euler angles.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STREvaluation indexYawRMSE (°)0.020850.019650.03805MAE (°)0.016770.015960.02970RollRMSE (°)0.022660.021820.04446MAE (°)0.018090.016920.03585PitchRMSE (°)0.021200.020230.04373MAE (°)0.016840.016440.03505Figure13 shows the estimation error of gyro biases in a stable phase, and Table 7 lists the RMSE and MAE of the above data. We can see that the bias estimation accuracy of the decentralized layout is a little lower than FUKF with two STR but higher than FEKF with three STR.Figure 13
Estimation error of gyro biases (stable phase).Table 7
RMSE and MAE of estimation error of gyro biases.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STREvaluation indexXRMSE (°/s)0.003210.002950.00389MAE (°/s)0.002490.002170.00315YRMSE (°/s)0.002830.002770.00396MAE (°/s)0.002260.002160.00323ZRMSE (°/s)0.002840.002840.00329MAE (°/s)0.002340.002320.00259


### 4.3.2. Execution Time

Table8 lists the execution time of 5000 iterations under the same simulation condition. We can see that the execution time of FUKF with two STR is the longest among the three layouts, and it is more than twice that of the decentralized layout.Table 8
Execution time.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STRExecution time (s)23.9150.9930.28


## 4.3.1. Estimation Accuracy of State

Figure12 shows the estimation error of Euler angles in a stable phase, and Table 6 lists the RMSE and MAE of the above data. We can see that the attitude estimation accuracy of the decentralized layout is a little lower than FUKF with two STR but higher than FEKF with three STR.Figure 12
Estimation error of Euler angles (stable phase).Table 6
RMSE and MAE of estimation error of Euler angles.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STREvaluation indexYawRMSE (°)0.020850.019650.03805MAE (°)0.016770.015960.02970RollRMSE (°)0.022660.021820.04446MAE (°)0.018090.016920.03585PitchRMSE (°)0.021200.020230.04373MAE (°)0.016840.016440.03505Figure13 shows the estimation error of gyro biases in a stable phase, and Table 7 lists the RMSE and MAE of the above data. We can see that the bias estimation accuracy of the decentralized layout is a little lower than FUKF with two STR but higher than FEKF with three STR.Figure 13
Estimation error of gyro biases (stable phase).Table 7
RMSE and MAE of estimation error of gyro biases.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STREvaluation indexXRMSE (°/s)0.003210.002950.00389MAE (°/s)0.002490.002170.00315YRMSE (°/s)0.002830.002770.00396MAE (°/s)0.002260.002160.00323ZRMSE (°/s)0.002840.002840.00329MAE (°/s)0.002340.002320.00259


## 4.3.2. Execution Time

Table8 lists the execution time of 5000 iterations under the same simulation condition. We can see that the execution time of FUKF with two STR is the longest among the three layouts, and it is more than twice that of the decentralized layout.Table 8
Execution time.
Attitude estimation layoutsDecentralized layoutFUKF with two STRFEKF with three STRExecution time (s)23.9150.9930.28


## 5. Conclusion

This paper proposes four kinds of attitude estimation layouts with one MEMS gyro and two simultaneously operating STR as part of the development task of the CubeSat NJUST-2. Simulation results show that the centralized layout II gives the highest estimation accuracy of Euler angles and gyro biases compared with the other three layouts and has high execution speed and fast convergence at the same time. Moreover, the comparisons between decentralized layout, FUKF with two STR, and FEKF with three STR indicate that the FUKF with two STR performs only a little better than the decentralized layout in terms of estimation accuracy but the execution time of the former is much longer than that of the latter; the attitude estimation layout with three STR cannot provide more accurate attitude and bias knowledge than that with two STR when using the same filter, and utilizing three STR will increase cost and take up more space in practical engineering application. Therefore, the centralized layout II is the best choice for the attitude measurement unit of the CubeSat with requirements of high accuracy, low cost, and limited computing capacity.



---

*Source: 4973970-2020-08-28.xml*