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

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.


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][2][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 sat-ellite 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 cm 3 and 7.8 W, respectively, while the size and maximum power supply of our 2U CubeSat are 20 × 10 × 10 cm 3 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 Section 2. 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.

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 where q 0 is the scalar part; q v is the vector part, q v = q 1 q 2 q 3 ½ Τ . The kinematic differential equation represented by quaternion is described as [14] where ω is the angular rate vector; ΩðωÞ is represented as [14] Ω ω ð Þ = 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 where N and K are the angle random walk and rate random walk, respectively; T s is the sampling period; η is the independent Gaussian white noise, η~Nð0, 1Þ. W b is the bias noise, denoted as In general, two STR are mounted at a tilted angle equal to 90°f or 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 sunsynchronous 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.
Therefore, we consider mounting the two STR on C2 at a tilted angle equal to 90°for protecting them from stray light disturbance, as Figure 2 shows.
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:  Figure 3 illustrates the angles between sun vector and orbit plane titled beta angle within one year from the above date.
As Figure 3 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   International Journal of Aerospace Engineering simulation duration, we can obtain α1 and α2 within the above period, as Figure 4 shows.
As Figure 3 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.

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 whereq 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 whereq si and q si are the measured and true attitude quaternion of the geocentric inertial frame relative to the STR body frame, respectively; q bs is the attitude quaternion of the STR body frame relative to the CubeSat body frame, denoted as 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 where Δr and Δp are the measurement errors in the directions of rolling and pointing, respectively. The measurement residual is denoted as whereq bi is the estimated attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame; q oi 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 where δq bi v , Δq bo v , and Δq bi v are the vector part of δq bi , Δq bo , and Δq bi , respectively.

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,  International Journal of Aerospace Engineering Ref. [23] proposes an effective nonlinear weighted average method to calculate the average of two quaternions as follows: where ; q 1 and q 2 are the two quaternions probably from two STR, Particle Filter (PF), multiple-model adaptive estimation, and other occasions [23]; w 1 and w 2 are the weights of q 1 and q 2 , 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. Figure 5 shows the structure diagram of 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 _ andP are the prediction and estimation error covariance matrix, respectively.

Centralized
Layout II. Figure 6 shows the structure diagram of the centralized layout II. q bi 2 is the quaternion from STR 2;q bi w is the nonlinear weighted average ofq bi 1 andq bi 2 . Figure 7 shows the structure diagram of the centralized layout III.

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. Figure 8 shows the structure diagram of 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 andP 2 are the prediction estimation error covariance matrices of respectively;P w is the linear weighted average ofP 1 andP 2 .

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 where Δt is the filter period; b ω bo,k-1 is the estimated angular rate of the CubeSat body frame relative to the orbit frame, denoted as where ω oi,k-1 is the angular rate of the orbit frame relative to the geocentric inertial frame, denoted as where ω o,k-1 is the mean orbit angular rate. TðqÞ is the direction cosine matrix represented by the quaternion, denoted as

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 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; Q k−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 where F k−1 is the Jacobian matrix with respect to state error, denoted as where b ω bi,k-1 is the estimated angular rate of the CubeSat body frame relative to the geocentric inertial frame, b ω bi,k-1 =ω bi,k−1 −b k−1 ; ½ω × is denoted as where G k−1 is the Jacobian matrix with respect to process noise, denoted as Q k−1 is denoted as  Figure 8: The decentralized layout. 6 International Journal of Aerospace Engineering (2) For the centralized layout III, Φ k/k−1 is denoted as where F k−1 is denoted as : ð26Þ where G k−1 is denoted as : ð28Þ Q k−1 is denoted as Step (2). Kalman gain matrix where H k is the state error measurement matrix;R k is the measurement noise variance matrix.
(1) For the centralized layout I, where S 1 is the variance of Δq bi 1 ; Δq bi 1 is the converted measurement noise of STR 1.
(2) For the centralized layout II, where S w 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, where S 2 is the variance of Δq bi 2 ; Δq bi 2 is the converted measurement noise of STR 2.
(4) For the decentralized layout, Step (3). Estimation of state error where δZ k is the integrated measurement residual.
(1) For the centralized layout I, where δq bi 1 v,k is the measurement residual of STR 1.
(2) For the centralized layout II, where δq bi w v,k is the measurement residual of STR 1 combined with STR 2.
(3) For the centralized layout III, where δq bi 2 v,k is the measurement residual of STR 2.
(1) For the centralized layout I and centralized layout II, (2) For the centralized layout III, whereP a,k andP b,k are the estimation error covariance matrices of Δq bo a v,k andΔq bo b v,k , respectively, inP k .
(3) For the decentralized layout,P

Correction of State
(1) For the centralized layout I and centralized layout II, (2) For the centralized layout III,  International Journal of Aerospace Engineering (3) For the decentralized layout, Δq bo ω v,k = ΔX w,k 4 : 6 ð Þ, ð45bÞ Δb w,k = ΔX w,k 1 : 3 ð Þ, ð45dÞ    4.1.6. True State Update. q bo,k+1 used for comparison with estimated attitude quaternion is updated by where ω bo,k is the true angular rate of the CubeSat body frame relative to the orbit frame, denoted as where ω bi,k is denoted as where I is the moment of inertia of CubeSat; its value is ð49Þ h r,k-1 is the angular momentum of reaction wheels; T c,k-1 is the control torque derived from zero-momentum attitude control [28]; T m,k-1 , T a,k-1 , and T g,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   Figure 11: Estimation error of Euler angles (initial phase).

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. Figure 9 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 10 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 11 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. Table 5 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.

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. Figure 12 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 13 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. Table 8 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.

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.

Data Availability
All data included in this study are available upon request by contacting the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.