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.
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 cm3 and 7.8 W, respectively, while the size and maximum power supply of our 2U CubeSat are 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 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.
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 where is the scalar part; is the vector part, .
The kinematic differential equation represented by quaternion is described as [14] where is the angular rate vector; is represented as [14]
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] where and are the measured and true angular rate of the CubeSat body frame relative to the geocentric inertial frame, respectively; , , and are the true biases at time , , and , respectively; is the measurement noise, denoted as where and are the angle random walk and rate random walk, respectively; is the sampling period; is the independent Gaussian white noise, .
is the bias noise, denoted as
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 Table 1.
In 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.
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: (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 3059
The initial orbit epoch is 8 Feb 2018 10:39:56.154 Coordinated Universal Time Greenwich (UTCG). 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 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.
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 where and are the measured and true attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame, respectively; is the converted measurement noise, denoted as where and are the measured and true attitude quaternion of the geocentric inertial frame relative to the STR body frame, respectively; is the attitude quaternion of the STR body frame relative to the CubeSat body frame, denoted as where , , and are the Euler angles of the STR body frame relative to the CubeSat body frame.
is the original measurement noise, under the condition of small angle; it can be denoted as where and are the measurement errors in the directions of rolling and pointing, respectively.
The measurement residual is denoted as where is the estimated attitude quaternion of the geocentric inertial frame relative to the CubeSat body frame; is the attitude quaternion of the geocentric inertial frame relative to the orbit frame; and are the estimated and true attitude quaternion of the orbit frame relative to the CubeSat body frame, respectively; is the estimation error.
Extracting the vector part of , the linear equation of measurement residual can be denoted as where , , and are the vector part of , , and , 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: where ; and are the two quaternions probably from two STR, Particle Filter (PF), multiple-model adaptive estimation, and other occasions [23]; and are the weights of and , 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.
is the quaternion from STR 1; , , and are the estimation, prediction, and estimation error of gyro bias, respectively; is the predicted attitude quaternion of the orbit frame relative to the CubeSat body frame; and are the prediction and estimation error covariance matrix, respectively.
3.1.2. Centralized Layout II
Figure 6 shows the structure diagram of the centralized layout II.
is the quaternion from STR 2; is the nonlinear weighted average of and .
3.1.3. Centralized Layout III
Figure 7 shows the structure diagram of the centralized layout III.
and are the vector parts of two estimation error quaternions from MEKF; is the linear weighted average of and .
3.1.4. Decentralized Layout
Figure 8 shows the structure diagram of the decentralized layout.
and are the two bias estimation errors from MEKF 1 and MEKF 2, respectively; is the linear weighted average of and ; and are the vector parts of two estimation error quaternions from MEKF 1 and MEKF 2, respectively; and are the prediction estimation error covariance matrices of and , respectively; is the linear weighted average of and.
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 , the prediction of state at time can be obtained by where is the filter period; is the estimated angular rate of the CubeSat body frame relative to the orbit frame, denoted as where is the angular rate of the orbit frame relative to the geocentric inertial frame, denoted as where is the mean orbit angular rate.
is the direction cosine matrix represented by the quaternion, denoted as
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 where is the state error transition matrix from time to time ; is the process noise input matrix; is the process noise variance matrix. (1)For the centralized layout I, centralized layout II, and decentralized layout, is denoted aswhere is the Jacobian matrix with respect to state error, denoted as where is the estimated angular rate of the CubeSat body frame relative to the geocentric inertial frame, ; is denoted as is denoted as where is the Jacobian matrix with respect to process noise, denoted as is denoted as (2)For the centralized layout III, is denoted aswhere is denoted as is denoted as where is denoted as is denoted as
Step (2). Kalman gain matrix where is the state error measurement matrix;is the measurement noise variance matrix. (1)For the centralized layout I,where is the variance of ; is the converted measurement noise of STR 1. (2)For the centralized layout II,where is the variance of ; is the converted measurement noise of STR 1 combined with STR 2. (3)For the centralized layout III,where is the variance of ; is the converted measurement noise of STR 2. (4)For the decentralized layout,
Step (3). Estimation of state error where is the integrated measurement residual. (1)For the centralized layout I,where is the measurement residual of STR 1. (2)For the centralized layout II,where is the measurement residual of STR 1 combined with STR 2. (3)For the centralized layout III,where is the measurement residual of STR 2. (4)For the decentralized layout,
Step (4). Estimation error covariance matrix. (1)For the centralized layout I and centralized layout II,(2)For the centralized layout III,where and are the estimation error covariance matrices of and, respectively, in . (3)For the decentralized layout,
3.2.3. Correction of State
(1)For the centralized layout I and centralized layout II,(2)For the centralized layout III,(3)For the decentralized layout,
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 14673
The 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 bias ; initial estimated bias ;; ; .
4.1.3. STR Parameters
, , , , , , ;, , , , , , ;.
4.1.4. Initial Attitude Knowledge
True Euler angle ; estimated Euler angle ;; ; true Euler angular rate ; estimated Euler angular rate .
4.1.5. Initial Estimation Error Covariance Matrix
For the centralized layout I, centralized layout II, and decentralized layout, ; for the centralized layout III, .
4.1.6. True State Update
used for comparison with estimated attitude quaternion is updated by where is the true angular rate of the CubeSat body frame relative to the orbit frame, denoted as where is denoted as where is the moment of inertia of CubeSat; its value is
is the angular momentum of reaction wheels; is the control torque derived from zero-momentum attitude control [28]; , , and are the disturbance torques derived from residual magnetism, atmospheric drag, and gravity gradient, respectively.
used for comparison with estimated gyro bias is updated by
4.1.7. Simulation Time
Start time ; duration ; .
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
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 , , and , respectively.
4.2.2. Convergence Speed
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.
4.2.3. Execution Time
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.
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
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.
4.3.2. Execution Time
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.
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.
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.
Acknowledgments
This work was supported in part by the National Natural Science Foundation of China (No. 61803204), the Natural Science Foundation of Jiangsu Province (No. BK20180465), and the China Scholarship Council (CSC).