Abstract
In order to enable the micro-nanosatellites equipped with microthrusters to better complete various space applications, it is necessary to estimate the thrust performance of the microthrusters in real-time on orbit. This paper proposes a real-time on-orbit estimation method for microthrust based on high-precision orbit determination. By establishing a high-precision orbit dynamic model, the microthrust generated by a microthruster is modeled as a first-order Markov model, combined with a high-precision GNSS measuring device, and the satellite position is obtained through the cubature Kalman filter algorithm, velocity, and thrust real-time on-orbit estimates. For a thrust of 100 μN, the error accuracy of the on-orbit estimation is 3.98%; for a thrust of 500 μN, the error accuracy is 1.79%; for a thrust of 5 mN, the error accuracy can be reduced to 1.43%; and when the thrust is 500 μN, the accuracy of orbit determination is 16 cm. This method solves the problem that the traditional on-orbit thrust estimation method cannot perform real-time on-orbit estimation of microthrust on the order of hundreds of μN. The real-time on-orbit estimation of microthrust of micro-nanosatellites equipped with microthrusters of the order of hundreds of micronewtons or even several mN to tens of mN has certain reference value.
1. Introduction
At present, with the development of miniaturization technology and highly integrated technology, the development of micro-nanosatellites has become very hot. The development of micro-nanosatellites is inseparable from the micropropulsion system. Through the micropropulsion system, the micro-nanosatellite can complete orbital applications such as orbit maintenance and orbit maneuvering, thereby completing autonomous orbit control. The micropropulsion system is inseparable from the evaluation of its propulsion performance. Although there is a relatively mature microthrust measurement system on the ground, it is difficult to accurately simulate the ground experimental environment due to the complex space environment [1]. And in order to better accomplish space missions of micro-nanosatellites, in recent years, some scholars have proposed methods for estimating thrust in orbit [2]. The main thrust estimation methods of microthruster are orbit estimation method [3–6] and attitude estimation method [2, 7–11].
However, either the orbit estimation method or the attitude estimation method, most of them are on-orbit estimation of the thrust of the microthruster with the magnitude of tens of mN [6, 7, 9]. Although the attitude estimation method can be used to estimate the thrust of the microthruster in the order of a few μN, the accuracy of the on-orbit estimation is low, and the on-orbit estimation of the thrust of the microthruster in the order of hundreds of μN is even less. There are some applications, that is, when the microthruster is only used for satellite orbit control, the microthruster is not installed eccentrically. Therefore, the microthrust generated by the microthruster passes the center of mass of the satellite. In this case, orbit estimation method should be adopted for the on-orbit estimation of microthrust in this situation. The orbital estimation method is mainly based on the function relationship between the change of the semimajor axis of the orbit and the thrust force to estimate the magnitude of the microthrust [3, 6]. However, this orbital estimation method based on the change of the orbital semimajor axis to determine the thrust is based on the postorbit determination method to determine the change of the orbital semimajor axis before and after the thrust action, so it belongs to the postorbit estimation method [6] and does not have real-time characteristics.
Aimed at the problems existing in the on-orbit estimation of microthruster thrust, a real-time on-orbit estimation method of microthruster thrust based on high-precision orbit determination is proposed in this paper, which solves the problem of real-time on-orbit estimation of microthruster with the magnitude of hundreds of μN. By establishing a high-precision orbit dynamic model considering various perturbation factors, the thrust generated by the microthruster is modeled as a first-order Markov process, and the state equation is obtained. In this paper, the high-precision GNSS measurement equipment is used to obtain the measurement value with noise, and then, the cubature Kalman filter algorithm is used to filter the state-space model to obtain the on-orbit estimation value including satellite position, velocity, and thrust of the microthruster.
2. Real-Time On-Orbit Thrust Estimation Method for a Microthruster
For the real-time on-orbit estimation of the microthruster thrust, this paper considers the case that the thrust of the microthruster passes through the mass center of satellite and only acts on the tangential direction of satellite, low earth circular orbit is considered, and the microthrust thrust is considered to be in the order of hundreds of micronewtons (10-4 N). Through the establishment of a high-precision orbital dynamic model including thruster thrust and various perturbation accelerations, combined with GNSS with high-precision orbit measurement accuracy, the cubature Kalman filter algorithm is used to obtain the real-time estimation of satellite velocity and position and the thrust of the microthruster.
2.1. Real-Time On-Orbit Estimation State-Space Model of Microthrust Based on High-Precision Orbit Determination
Due to the fact that the thrust of the microthruster is small, when the thrust action time is short, the change of orbit height is small. Therefore, it is necessary to establish a high-precision state-space model, including dynamic model and observation model. In this paper, a high-precision orbit dynamic model is used to consider the conservative and nonconservative forces on the LEO (low earth orbit) satellite, and the thrust of the microthruster is modeled into the dynamic model as a state variable to be estimated. Combined with the measurement model with measurement data provided by high-precision measurement device, the thrust in orbit estimation is obtained simultaneously in the process of real-time orbit determination through the filtering process.
2.1.1. Orbital Dynamic Model including Microthruster Thrust
The traditional method of determining the thrust of the thruster based on the relationship between the change of the semimajor axis of the orbit after the thrust is applied, and the thrust is a method of postthrust estimation, which is not real-time. In order to estimate the magnitude of thrust in real time, the thrust generated by the thruster is modeled into the dynamic model as a state variable, so the state variable is . Therefore, the dynamic model of the satellite in the process of thrust is where the reference coordinate system is J2000 geocentric inertial coordinate system, and are the position and velocity vectors of the satellite in the inertial coordinate system, is the microthrust generated by the microthruster, is the gravitational acceleration of the central celestial body, and are the perturbation accelerations of conservative force and nonconservative force, respectively, is the time parameter of the first-order Markov process, is the mass of the satellite, and is the Gaussian noise.
Compared with the traditional orbit determination problem, the state variable of this paper is . The microthrust generated by the microthruster is modeled into the dynamic model as the state variable to be estimated, and the first-order Markov process is used to express the thrust of the microthruster, so the first-order Markov process satisfies the Langevin differential equation [12]. where and are the time parameters and the standard deviation of the first-order Markov process, which determine the degree of change in the thrust estimate on orbit. Therefore, a reasonable selection of these two parameters will help to improve the accuracy of the microthruster thrust on-orbit estimation.
Many typical microthrusters, such as vacuum cathode microarc microthrusters, pulsed plasma microthrusters, MEMS microthrusters, and laser microthrusters, have thrust magnitude of 10-4 N, so their influence on orbit altitude cannot be ignored. Therefore, in the problem of high-precision orbit determination, the microthrust must be modeled in the dynamic model of the system. The first-order Markov stochastic process is used to represent the microthrust, because the first-order Markov process is a process that changes slowly with time, which can well simulate the microthrust characteristics generated by the microthruster on orbit in real time. By adjusting the parameters and , the accuracy of the model can be improved, so as to improve the accuracy of microthrust real-time on-orbit estimation.
(1) Gravitational Acceleration of Central Celestial Body acentral. For the acceleration model, the gravitational acceleration, the conservative force acceleration, the nonconservative force acceleration, and the thrust acceleration are considered. The acceleration caused by the central celestial body is where is the gravitational constant of the earth, and . Here, is the vector, which represents the position vector of the satellite in the J2000 inertial coordinate system.
(2) Perturbation Acceleration of Conservative Force aconservative. In this paper, for the perturbation acceleration of conservative force, we consider the nonspherical perturbation of the earth, the tidal perturbation, the third body gravitational perturbation, and the relativistic effect perturbation, so where is the nonspherical perturbation, using the EGM2008 earth gravity field model, and this paper considers order; is the tidal perturbation, and this paper considers the solid tide perturbation; is the gravitational perturbation of the third body, mainly the gravitational perturbation of the sun and the moon; and is the relativistic gravitational perturbation.
The aspherical perturbation and tidal perturbation of the earth mainly cause the change of the earth’s position function. The potential function caused by the nonspherical perturbation of the earth is [13] where and are the geocentric longitude and latitude of the satellite; , , and are the coefficients of zonal harmonics and tesseral harmonics which can be obtained from EGM08; and and are Legendre functions. In this paper, we consider the earth nonspherical perturbations of order , that is, the case of . For high-precision orbit determination, the accuracy of order 70 has been far satisfied. Therefore, in the J2000 geocentric inertial coordinate system, the perturbation acceleration caused by the nonspherical perturbation of the earth is where is the coordinate transformation matrix from the spherical coordinate system to the rectangular coordinate system [13] and is the coordinate transformation matrix from the earth fixed coordinate system to the J2000 geocentric inertial coordinate system [13].
Tides are mainly caused by the gravity of the moon and the sun on the earth, so the tidal perturbation mainly causes the change of the gravitational potential function of the earth. According to the research in this paper, only the perturbation acceleration caused by the earth tide is considered. Since the magnitude of the gravitational potential function is less than 10-10 when , only the additional potential function caused by the solid tide when is considered. Therefore, under the action of the sun and the moon, the additional position of the solid tide is , so the potential function caused by the additional potential of the earth tide is .
The perturbation acceleration caused by the combination of the earth nonspherical perturbation and the earth tide perturbation is
The gravitational perturbation of the third body mainly considers the gravitational perturbation of the sun and the moon [14], so as shown in equation (8), the gravitational perturbation acceleration of the third body, where and are the gravitational constant of the sun and the gravitational constant of the moon, respectively, and are the position vectors of the sun and the moon from the center of the earth in the J2000 inertial coordinate system.
Traditional celestial body motions are all studied under the framework of Newtonian gravity, but for the high-precision orbit determination of low-orbit satellites, the influence of general relativity must be considered. Therefore, in the framework of Newton’s gravity, the perturbation acceleration is caused by the effect of general relativity, this paper only considers the Schwarzschild term, and there is [13] where is the speed of light and .
2.1.2. Perturbation Acceleration of Nonconservative Force
For the perturbed acceleration caused by nonconservative forces, the solar pressure perturbation and atmospheric drag perturbation are considered in this paper, as shown in equation (12).
The solar pressure perturbation and atmospheric drag perturbation are not only related to the orbit height of the satellite but also related to the mass, shape, and volume of the satellite. For the perturbed acceleration of solar pressure, there is [13] where is the exposure factor; is the solar pressure constant, and ; is the distance from the earth to the sun (one astronomical unit), and ; is the light pressure coefficient, which is a dimensionless parameter and is related to the material of the satellite surface, and ; and is the surface quality ratio of the satellite exposed to sunlight. For the sun exposure factor , the orbit of the satellite is low earth orbit. Considering the cylindrical shadow, the relationship between and the position of the satellite is shown in Figure 1. The pseudocode for calculating the exposure factor is shown in Algorithm 1.
| ||||||||||||||||||||||||||||
For the perturbation acceleration of atmospheric drag, there is [15] where is the atmospheric resistance density, and this article uses the NRLMISISE 00 model; is the windward area of the satellite; is the atmospheric drag coefficient, which is related to the satellite body shape, and ; and is the speed of the satellite relative to the atmosphere, and is the modulus of .
(1) Thrust Acceleration athrust. As for the thrust acceleration generated by the microthruster, when the thrust action time is short, the mass consumption rate per second is considered to be zero; that is, the mass of the satellite remains unchanged, so the thrust acceleration generated by the microthruster is
2.1.3. Measurement Model
The measurement device used in this paper is GNSS (Global Navigation Satellite System), which does not depend on the ground station. The data measured by GNSS is converted into the position and velocity information of the satellite in the J2000 coordinate system as the measurement value with noise. Therefore, the measurement model is where is the measured value with noise, which is the position and speed of the satellite in operation, is the function relation between the measurement value and state variable, and is the measurement noise and is Gaussian white noise.
It is worth noting that due to the fact that the thrust generated by the microthruster is relatively small, when the thrust action time of the microthruster in the order of hundreds of micronewtons is short, the change in the orbital height caused by the thrust is at the level of m or even dm.
According to the relationship between the change of the semimajor axis of the orbit and the thrust after the thrust is applied in the tangential direction of the circular orbit [3, 6], there is where is the thrust, is the semimajor axis of the orbit, and is the change of the semimajor axis of the orbit after the thrust, and is the thrust action time. For a satellite with a mass of 25 kg, a thrust of 100 μN acts on the tangential direction of 1800 s, and the change in orbital height is 13.01 m. Therefore, on the basis of establishing a relatively accurate dynamic model, in order to make the estimation of state variables more accurate, the accuracy of the required orbit measurement device should be as high as possible.
2.1.4. Discretization of the State-Space Model
Combining the first two sections, the state-space model for the low earth orbit microthrust orbit maneuvering process of the micro-nanosatellite studied in this article, there is a state-space model:
It can be seen that the model is continuous and strongly nonlinear, but the filtering algorithm is the estimation from to , so it is necessary to discretize the continuous and strongly nonlinear state-space model.
The fourth-order Runge-Kutta is used to discretize the continuous differential equations of position and velocity in the state equation, while the microthrust modeled by the first-order Markov process is discretized according to the characteristics of the first-order Markov stochastic process [16]. Therefore, the state variables are divided into two parts , where. After using Runge-Kutta discretization, there is
After the discretization of the state differential equation of the microthruster, there is where is the sampling time and is discrete driving white noise. Therefore, the discretized state-space model is where
2.2. Real-Time On-Orbit Estimation of Microthrust Based on the Cubature Kalman Filter
In order to estimate the state variables of the satellite in real time, a filtering algorithm is needed. The two main factors that affect the accuracy of state variable estimation are the accuracy of the state-space model and the accuracy of the filtering algorithm. On the basis of establishing a relatively high-precision model, the selection of a suitable filtering algorithm determines the accuracy of the state variable estimation value.
In this paper, the cubature Kalman filter [17] is used. Compared with the extended Kalman filter, the cubature Kalman filter has higher accuracy. And for the strongly nonlinear state-space model, linearization is not required, but a set of cubature points is used to approximate the state mean and covariance of the nonlinear system with additional Gaussian noise, which is the closest to Bayesian filtering in theory. The approximate algorithm is a powerful tool for solving nonlinear system state estimation. Compared with the unscented Kalman filter, the cubature Kalman filter is in the unscented Kalman filter, but compared with the unscented Kalman filter, the cubature Kalman filter has strict mathematical proof and theoretical support, and the amount of calculation is less than that of the unscented Kalman filter. The state-space model in this paper is strongly nonlinear and has many dimensions of state variables. Therefore, the cubature Kalman filter is selected to filter the state-space model with process noise and measurement noise, so as to obtain the estimated values of position, velocity, and microthrust of the satellite in the low thrust orbit maneuver process.
Like the classical Kalman Filter, the cubature Kalman filter is divided into three processes: initialization, time updating, and measurement updating.
2.2.1. Initialization
The initial state is , the initial covariance is , the process noise covariance matrix is , and the measurement noise covariance matrix is . It is worth noting that therefore, there is
2.2.2. Time Updating
The posterior estimated value of the state variable at time is and the posterior estimated value of the covariance matrix is . (1)Calculating the cubature points where is the dimension of the state variable and is the number of cubature points, which is twice the dimension of state variable . In the formula, represents the -th column of where is the -dimensional unit vector.(2)Nonlinear propagation of cubature points
So, there is (3)Calculating the predicted state (4)Calculating the predicted error covariance
2.2.3. Measurement Updating
(1)Calculating the cubature points (2)Calculating the measurement of the propagation cubature points (3)Calculating the predicted measurement (4)Calculating the measurement error covariance matrix (5)Calculating the cross-covariance matrix (6)Calculating the Kalman gain (7)Calculating the updated state (8)Calculating the posterior error covariance matrix
The whole process of the cubature Kalman filter is shown in the pseudocode of Algorithm 2. On the basis of the continuous strong nonlinear state-space model, the flowchart of the whole filtering process is shown in Figure 2, in which the specific flow of time updating and measurement updating is shown by the pseudocode.
| ||||||||||||||||||||||||||||||||||||||||||||||
3. Simulation
The simulation platform of this study is MATLAB, and the simulation software is developed by the laboratory. The simulation software is used to simulate the satellite position and velocity information measured by GNSS. The satellite orbit dynamic model and the cubature Kalman filter algorithm are established in MATLAB. By filtering the position and velocity information of the satellite with noise, the estimated value of the state variable is estimated; that is, the on-orbit estimated value of the satellite position, velocity, and microthruster thrust is obtained. The specific simulation process is shown in Figure 3.
In this paper, the tangential microthrust applied to a micro-nanosatellite with a mass of 25 kg and an orbit altitude of 500 km is simulated. It is assumed that the initial attitude of the satellite is a three-axis earth orientation mode, so the body coordinate system of the satellite coincides with the orbit coordinate system. When the thrust of the microthruster is 500 μN, the direction is the same as the tangential direction of the orbit and the velocity, and considering that the thrust output of the microthruster is continuous, the thrust action time is the whole simulation process. Specific parameters of the simulation process are shown in Table 1.
For the initial covariance matrix , it represents the degree of thrust to the initial state variable estimate . The closer is to 0, the closer the initial estimate is to the real value. Therefore, in this paper, we choose the initial covariance matrix value which is relatively small. For process noise covariance matrix , it represents the accuracy of model establishment. The smaller is, the higher the thrust to the established model is and the higher the proportion of the predicted state value in the estimated value will be. In this paper, the process noise covariance of position and velocity represents the accuracy of satellite orbit dynamic model establishment. Therefore, by referring to the selected values of covariance in [18, 19] and the relatively accurate orbital dynamic model in this paper, the initial values of process noise covariance are set as 0.012 and 0.0012. For the microthruster thrust which is modeled as a first-order Markov process, the process noise covariance is determined by the properties of the first-order Markov process and is related to the time parameter and process parameter of the first-order Markov process. For the measurement noise covariance matrix , it depends on the accuracy of the measurement device. In this paper, the measurement accuracy of the simulated GNSS is 5 cm and the velocity accuracy is 1 cm/s, so the measurement noise covariance matrix is shown in Table 1.
The simulation results of satellite position and velocity under 500 μN thrust are shown in Figures 4(a) and 5(a). The position and velocity of the satellite are the position and velocity in the J2000 geocentric inertial coordinate system. In order to explore the position and velocity information of the satellite under 500 μN thrust, the position and velocity information of the satellite with and without thrust are compared, as shown in Figures 4 and 5. Figures 4(b)–4(d) are partial enlarged images of the satellite three-axis position. Figures 5(b)–5(d) are partial enlarged images of satellite three-axis velocity. Figure 6 is a three-dimensional simulation of the satellite orbit, and Figure 6(b) is a three-dimensional partial enlarged view of the orbit position. According to the simulation results, it can be found that applying tangential microthrust can improve the orbital altitude and velocity of the satellite.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
The simulation result of the RMSE (root mean square error) of orbit determination based on the microthrust real-time on-orbit estimation method based on high-precision orbit determination is shown in Figure 7. The specific RMSE expression is shown in equation (39), which, respectively, represents the root mean square error of position and velocity at time , where is the run times of Monte Carlo experiments and is the -th Monte Carlo experiment. and are the positions and velocities of satellites obtained by simulating GNSS. and are the estimates of the position and velocity of the satellite obtained by the algorithm in this paper.
(a)
(b)
According to the simulation results, the average RMSE of the position is 27.58 cm and the average RMSE of the velocity is 53.73 cm/s. The RMSE simulation results of satellite three-axis position are shown in Figure 8, in which the average RMSE of the -axis is 15.92 cm, the average RMSE of -axis is 15.93 cm, and the average RMSE of -axis is 15.92 cm. According to the experimental results, it can be found that the orbit determination accuracy of this method can reach decimeter level, which shows that the algorithm has a high accuracy.
According to Figure 9, under the action of 500 μN thrust, the real-time on-orbit estimated value of thrust stabilized after about 60 s. After stabilization, the real-time estimation error of orbital thrust is about 1.7% and the average thrust estimation error is 1.79%.
(a)
(b)
In order to verify the accuracy of the real-time on-orbit estimation method of microthrust, a comparative experiment was carried out on the basis of 500 μN thrust, and two thrusts of 100 μN and 5 mN were added. Among them, the satellite orbit position under different thrusts is shown in Figure 10, and Figure 10(b) is a partial enlarged view of the orbit position. The real-time on-orbit estimation of the three thrusts and the real-time on-orbit estimation error of the thrust are shown in Figure 11. According to the simulation results, it can be found that the smaller the microthrust, the longer it takes for the thrust estimation to stabilize, and the accuracy of the thrust estimation increases with the increase in the thrust. The larger the microthrust, the higher the accuracy of real-time on-orbit thrust estimation. See Table 2 for the specific simulation results.
(a)
(b)
(a)
(b)
4. Conclusion
In this paper, through the establishment of the real-time on-orbit estimation model and simulation process of the micro-nanosatellite equipped with a microthruster, the following conclusions can be drawn: (1)In this paper, a high-precision orbit dynamic model is established, including various perturbed accelerations, and the thrust generated by the microthruster is modeled as a first-order Markov process. The position, velocity, and microthrust of the satellite are estimated in real time by using the cubature Kalman filter combined with the measured data of the simulated high-precision GNSS measurement device. For a microthrust of 500 μN, the accuracy of real-time on-orbit estimation error is 1.79%, and the accuracy of orbit determination can reach 16 cm. Through comparative experiments, the real-time on-orbit estimation error of 100 μN microthrust is 3.98%, and the real-time on-orbit estimation error of 5 mN microthrust can be reduced to 1.43%(2)The method of real-time on-orbit estimation of microthrust based on high-precision orbit determination proposed in this paper can estimate the size of microthrust well and has real-time performance. This method solves the problem of real-time on-orbit estimation of microthrust with the order of hundreds of micronewtons and has a certain reference value for the on-orbit estimation of the thrust of micro-nanosatellites equipped with microthrusters of the order of hundreds of μN to tens of mN, so that micro-nanosatellites can better complete various space application tasks.(3)Since the development of the microthruster is still in the final test stage, this simulation study provides the feasible optimal scheme of on-orbit thrust estimation for the satellite on-orbit flight experiment. Therefore, the selection of the actual parameters after the satellite is launched may be different from the simulation parameters. The research in this paper lays a foundation for the actual parameter selection of the satellite in orbit flight. Subsequent work will compare the actual flight data and simulation data after the satellite was launched to better verify the effectiveness of the method.
Data Availability
The on-orbit experimental data involved in the research process of this paper are obtained by the simulation software independently developed by our laboratory. The corresponding data can be obtained by contacting the first author (yql2314272020@163.com).
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.