# Real-Time On-Orbit Estimation Method for Microthruster Thrust Based on High-Precision Orbit Determination

Authors: Qinglin Yang; Weijing Zhou; Hao Chang
Journal: International Journal of Aerospace Engineering
 (2021)

---

## 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.


---

## Body

## 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 isx=rvthrustT. Therefore, the dynamic model of the satellite in the process of thrust is
(1)ẋ=ṙv̇thrust·=vacentral+aconservative+anonconservative+athrust−1τ⋅thrust+w,where the reference coordinate system is J2000 geocentric inertial coordinate system, r and v are the position and velocity vectors of the satellite in the inertial coordinate system, thrust is the microthrust generated by the microthruster, acentral is the gravitational acceleration of the central celestial body, aconservative and anonconservative are the perturbation accelerations of conservative force and nonconservative force, respectively, τ is the time parameter of the first-order Markov process, m is the mass of the satellite, and w is the Gaussian noise.Compared with the traditional orbit determination problem, the state variable of this paper isx=rvthrustT. 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].
(2)thrust·t=−1τ⋅thrustt+2σ2τ,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
(3)acentral=−μr3r,where μ is the gravitational constant of the earth, and μ=3.986004415×1014m3/s2. Here, r 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
(4)aconservative=anonspherical+atide+asun&moon+arelativistic,whereanonspherical is the nonspherical perturbation, using the EGM2008 earth gravity field model, and this paper considers 70×70 order; atide is the tidal perturbation, and this paper considers the solid tide perturbation; asun&moon is the gravitational perturbation of the third body, mainly the gravitational perturbation of the sun and the moon; and arelativistic 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]
(5)V=μr∑n=2∞C¯n,0RErnP¯nsinφ+∑n=2∞∑m=1nRErnP¯nmsinφC¯nmcosmλ+S¯nmsinmλ,where φ and λ are the geocentric longitude and latitude of the satellite; C¯n,0, C¯nm, and S¯nm are the coefficients of zonal harmonics and tesseral harmonics which can be obtained from EGM08; and P¯nsinφ and P¯nmsinφ are Legendre functions. In this paper, we consider the earth nonspherical perturbations of order 70×70, that is, the case of n=70. 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
(6)anonspherical=C1C0∂V∂r∂V∂φ∂V∂λ,where C0 is the coordinate transformation matrix from the spherical coordinate system to the rectangular coordinate system [13] and C1 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 n≥3, only the additional potential function caused by the solid tide when n=2 is considered. Therefore, under the action of the sun and the moon, the additional position of the solid tide is ΔC,ΔS, so the potential function caused by the additional potential of the earth tide is ΔV.The perturbation acceleration caused by the combination of the earth nonspherical perturbation and the earth tide perturbation is(7)anonspherical&tide=C1C0∂V+ΔV∂r∂V+ΔV∂φ∂V+ΔV∂λ.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 μsun and μmoon are the gravitational constant of the sun and the gravitational constant of the moon, respectively, rs and rM are the position vectors of the sun and the moon from the center of the earth in the J2000 inertial coordinate system.
(8)asun&moon=−μsun⋅r−rsr−rs3+rsrs3+−μmoon⋅r−rMr−rM3+rMrM3.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]
(9)arelativisitic=3μr3⋅vc2⋅r,where c is the speed of light and c=3×108m/s.


#### 2.1.2. Perturbation Acceleration of Nonconservative Forceanonconservative

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).
(10)anonconservative=asolarpressure+adrag.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]
(11)asolarpressure=−F⋅ρSR⋅AU2⋅CR⋅Sm⋅r−rsr−rs3,where F is the exposure factor; ρSR is the solar pressure constant, and ρSR=4.5605×10‐6N/m2; AU is the distance from the earth to the sun (one astronomical unit), and AU=1.4959787×1011m; CR is the light pressure coefficient, which is a dimensionless parameter and is related to the material of the satellite surface, and CR=1.5; and S/m is the surface quality ratio of the satellite exposed to sunlight. For the sun exposure factor F, the orbit of the satellite is low earth orbit. Considering the cylindrical shadow, the relationship between F and the position of the satellite is shown in Figure 1. The pseudocode for calculating the exposure factor F is shown in Algorithm 1.Algorithm 1: Calculate the exposure factor F.
(1) Input: r, rs Re.(2) D⟵r, rs(3) mo_r⟵r.(4) mo_ rs⟵ rs(5) if D>0(6)  thenF⟵1;(7) elseif D<=0(8)  then H2⟵mo_r2-(D/mo_rs)2;(9)  if H2>=Re2(10)    thenF⟵1;(11)  elseif H2<Re2(12)    thenF⟵1;(13)  end(14) endFigure 1
Relationship between sun exposure factorF and satellite position.For the perturbation acceleration of atmospheric drag, there is [15]
(12)adrag=−12⋅ρ⋅Sm⋅CD⋅ve⋅vewhere ρ is the atmospheric resistance density, and this article uses the NRLMISISE 00 model; S/m is the windward area of the satellite; CD is the atmospheric drag coefficient, which is related to the satellite body shape, and CD=2.2; and ve is the speed of the satellite relative to the atmosphere, and ve is the modulus of ve.(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
(13)athrust=thrustm.


#### 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(14)z=Hx+v,where z is the measured value with noise, which is the position and speed of the satellite in operation,H⋅ is the function relation between the measurement value and state variable, and v 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
(15)F=12mμa3ΔaΔt,where F is the thrust, a is the semimajor axis of the orbit, and Δa is the change of the semimajor axis of the orbit after the thrust, and Δt 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:(16)ẋ=ṙv̇thrust·=vacentral+aconservative+anonconservative+athrust−1τ⋅thrust+w=fx+w,z=Hx+v.It can be seen that the model is continuous and strongly nonlinear, but the filtering algorithm is the estimation fromk to k+1, 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 x=x1,thrustT, wherex1=r,vT. After using Runge-Kutta discretization, there is
(17)k1=fx1,k2=fx1+12Δt⋅k1+12Δt,k3=fx1+12Δt⋅k2+12Δt,k4=fx1+Δt⋅k3+Δt,x1,k+1=x1,k+16Δt⋅k1+2k2+2k3+k4.After the discretization of the state differential equation of the microthruster, there is(18)thrustk+1=e−Δtτ⋅thrustk+ωk,where Δt is the sampling time and ωk is discrete driving white noise. Therefore, the discretized state-space model is
(19)xk+1=Fxk+wk,zk=Hxk+vk,where
(20)xk=x1,kthrustk.


### 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 κ=0 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 isx∧0+, the initial covariance is P0+, the process noise covariance matrix is Qk, and the measurement noise covariance matrix is Rk. It is worth noting that
(21)x=rvthrust=x1thrust,therefore, there is
(22)x0+=r0+v0+thrust0+=x1,0+thrust0+.


#### 2.2.2. Time Updating

The posterior estimated value of the state variable at timek is
(23)x^k+=r^k+v^k+thrustk+∧=x^1,k+thrustk+∧,and the posterior estimated value of the covariance matrix is Pk+.
(1)
Calculating the cubature points(24)x^ki=x^k++Pk+e,−ei,i=1,2,⋯,2n,where
(25)x^ki=x^1,kithrustki∧.n is the dimension of the state variable and i is the number of cubature points, which is twice the dimension of state variable n. In the formula, ei represents the i-th column of
(26)e,−ei=10⋮001⋮0⋯00⋮1−10⋮00−1⋮0⋯00⋮−1n×2ni=1,2,⋯,2n,where e is the n-dimensional unit vector.(2)
Nonlinear propagation of cubature points(27)X1,k+1i=fx^1,ki,Thrustk+1i=e−Δt/τ⋅thrustki,∧So, there is(28)Xk+1i=X1,k+1iThrustk+1i.(3)
Calculating the predicted state(29)x^k+1−=12n∑i=12nXk+1i.(4)
Calculating the predicted error covariance(30)Pk+1−=12n∑i=12nXk+1i−x^k+1−Xk+1i−x∧k+1−T+Qk.


#### 2.2.3. Measurement Updating

(1)
Calculating the cubature points(31)x^k+1i=x^k+1−+Pk+1−e,−ei,i=1,2,⋯,2n.(2)
Calculating the measurement of the propagation cubature points(32)Zk+1i=Hx^k+1i.(3)
Calculating the predicted measurement(33)z^k+1=12n∑i=12nZk+1i.(4)
Calculating the measurement error covariance matrix(34)Pz=12n∑i=12nZk+1i−z^k+1Zk+1i−z∧k+1T+Rk+1.(5)
Calculating the cross-covariance matrix(35)Pxz=12n∑i=12nXk+1i−x^k+1−Zk+1i−z∧k+1T.(6)
Calculating the Kalman gain(36)Kk+1=Pxz,k+1Pz,k+1−1.(7)
Calculating the updated state(37)x^k+1+=x^k+1−+Kk+1zk+1−z^k+1.(8)
Calculating the posterior error covariance matrix(38)Pk+1+=Pk+1−−Kk+1Pz,k+1Kk+1T.The whole process of the cubature Kalman filter is shown in the pseudocode of Algorithm2. 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.Algorithm 2:Pseudocode of the cubature Kalman filter process.
(1) Computeei⟵Ιni(2) Input:x^0+, P0+, zk+1(3)fork=1toz_lengthdo(4)Sk+⟵chol Pk+;(5)fori=1to 2ndo(6) Calculatex^ki using x^k+, Sk+, ei, Pk+;(7)Xk+11⟵fx^ki;(8)end(9)x^ki⟵1/2n∑i=12nXk+1i;(10)Pk+1−⟵1/2n∑i=12nXk+1i−x^k+1−Xk+1i−x∧k+1−T+Qk;(11)Sk+1+⟵cholPk+1+;(12)fori=1to 2ndo(13) Calculatex^k+1i using x^k+1−, Sk+1+, ei, Pk+1+;(14)Zk+1i⟵Hx^k+1i;(15)end(16)z^k+1⟵1/2n∑i=12nZk+1i;(17)Pz⟵1/2n∑i=12nZk+1i−z^k+1Zk+1i−z∧k+1T+Rk+1;(18)Pxz⟵1/2n∑i=12nXk+1i−x^k+1−Zk+1i−z∧k+1T;(19)Kk+1⟵Pxz,k+1Pz,k+1−1;(20)x^k+1+⟵x^k+1−+Kk+1zk+1−z^k+1;(21)Pk+1+⟵Pk+1−−Kk+1Pz,k+1Kk+1T;(22)end(23) Output:x^k+1+, Pk+1+Figure 2
Flowchart of the whole filtering process.


## 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 isx=rvthrustT. Therefore, the dynamic model of the satellite in the process of thrust is
(1)ẋ=ṙv̇thrust·=vacentral+aconservative+anonconservative+athrust−1τ⋅thrust+w,where the reference coordinate system is J2000 geocentric inertial coordinate system, r and v are the position and velocity vectors of the satellite in the inertial coordinate system, thrust is the microthrust generated by the microthruster, acentral is the gravitational acceleration of the central celestial body, aconservative and anonconservative are the perturbation accelerations of conservative force and nonconservative force, respectively, τ is the time parameter of the first-order Markov process, m is the mass of the satellite, and w is the Gaussian noise.Compared with the traditional orbit determination problem, the state variable of this paper isx=rvthrustT. 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].
(2)thrust·t=−1τ⋅thrustt+2σ2τ,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
(3)acentral=−μr3r,where μ is the gravitational constant of the earth, and μ=3.986004415×1014m3/s2. Here, r 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
(4)aconservative=anonspherical+atide+asun&moon+arelativistic,whereanonspherical is the nonspherical perturbation, using the EGM2008 earth gravity field model, and this paper considers 70×70 order; atide is the tidal perturbation, and this paper considers the solid tide perturbation; asun&moon is the gravitational perturbation of the third body, mainly the gravitational perturbation of the sun and the moon; and arelativistic 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]
(5)V=μr∑n=2∞C¯n,0RErnP¯nsinφ+∑n=2∞∑m=1nRErnP¯nmsinφC¯nmcosmλ+S¯nmsinmλ,where φ and λ are the geocentric longitude and latitude of the satellite; C¯n,0, C¯nm, and S¯nm are the coefficients of zonal harmonics and tesseral harmonics which can be obtained from EGM08; and P¯nsinφ and P¯nmsinφ are Legendre functions. In this paper, we consider the earth nonspherical perturbations of order 70×70, that is, the case of n=70. 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
(6)anonspherical=C1C0∂V∂r∂V∂φ∂V∂λ,where C0 is the coordinate transformation matrix from the spherical coordinate system to the rectangular coordinate system [13] and C1 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 n≥3, only the additional potential function caused by the solid tide when n=2 is considered. Therefore, under the action of the sun and the moon, the additional position of the solid tide is ΔC,ΔS, so the potential function caused by the additional potential of the earth tide is ΔV.The perturbation acceleration caused by the combination of the earth nonspherical perturbation and the earth tide perturbation is(7)anonspherical&tide=C1C0∂V+ΔV∂r∂V+ΔV∂φ∂V+ΔV∂λ.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 μsun and μmoon are the gravitational constant of the sun and the gravitational constant of the moon, respectively, rs and rM are the position vectors of the sun and the moon from the center of the earth in the J2000 inertial coordinate system.
(8)asun&moon=−μsun⋅r−rsr−rs3+rsrs3+−μmoon⋅r−rMr−rM3+rMrM3.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]
(9)arelativisitic=3μr3⋅vc2⋅r,where c is the speed of light and c=3×108m/s.


### 2.1.2. Perturbation Acceleration of Nonconservative Forceanonconservative

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).
(10)anonconservative=asolarpressure+adrag.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]
(11)asolarpressure=−F⋅ρSR⋅AU2⋅CR⋅Sm⋅r−rsr−rs3,where F is the exposure factor; ρSR is the solar pressure constant, and ρSR=4.5605×10‐6N/m2; AU is the distance from the earth to the sun (one astronomical unit), and AU=1.4959787×1011m; CR is the light pressure coefficient, which is a dimensionless parameter and is related to the material of the satellite surface, and CR=1.5; and S/m is the surface quality ratio of the satellite exposed to sunlight. For the sun exposure factor F, the orbit of the satellite is low earth orbit. Considering the cylindrical shadow, the relationship between F and the position of the satellite is shown in Figure 1. The pseudocode for calculating the exposure factor F is shown in Algorithm 1.Algorithm 1: Calculate the exposure factor F.
(1) Input: r, rs Re.(2) D⟵r, rs(3) mo_r⟵r.(4) mo_ rs⟵ rs(5) if D>0(6)  thenF⟵1;(7) elseif D<=0(8)  then H2⟵mo_r2-(D/mo_rs)2;(9)  if H2>=Re2(10)    thenF⟵1;(11)  elseif H2<Re2(12)    thenF⟵1;(13)  end(14) endFigure 1
Relationship between sun exposure factorF and satellite position.For the perturbation acceleration of atmospheric drag, there is [15]
(12)adrag=−12⋅ρ⋅Sm⋅CD⋅ve⋅vewhere ρ is the atmospheric resistance density, and this article uses the NRLMISISE 00 model; S/m is the windward area of the satellite; CD is the atmospheric drag coefficient, which is related to the satellite body shape, and CD=2.2; and ve is the speed of the satellite relative to the atmosphere, and ve is the modulus of ve.(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
(13)athrust=thrustm.


### 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(14)z=Hx+v,where z is the measured value with noise, which is the position and speed of the satellite in operation,H⋅ is the function relation between the measurement value and state variable, and v 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
(15)F=12mμa3ΔaΔt,where F is the thrust, a is the semimajor axis of the orbit, and Δa is the change of the semimajor axis of the orbit after the thrust, and Δt 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:(16)ẋ=ṙv̇thrust·=vacentral+aconservative+anonconservative+athrust−1τ⋅thrust+w=fx+w,z=Hx+v.It can be seen that the model is continuous and strongly nonlinear, but the filtering algorithm is the estimation fromk to k+1, 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 x=x1,thrustT, wherex1=r,vT. After using Runge-Kutta discretization, there is
(17)k1=fx1,k2=fx1+12Δt⋅k1+12Δt,k3=fx1+12Δt⋅k2+12Δt,k4=fx1+Δt⋅k3+Δt,x1,k+1=x1,k+16Δt⋅k1+2k2+2k3+k4.After the discretization of the state differential equation of the microthruster, there is(18)thrustk+1=e−Δtτ⋅thrustk+ωk,where Δt is the sampling time and ωk is discrete driving white noise. Therefore, the discretized state-space model is
(19)xk+1=Fxk+wk,zk=Hxk+vk,where
(20)xk=x1,kthrustk.


## 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 isx=rvthrustT. Therefore, the dynamic model of the satellite in the process of thrust is
(1)ẋ=ṙv̇thrust·=vacentral+aconservative+anonconservative+athrust−1τ⋅thrust+w,where the reference coordinate system is J2000 geocentric inertial coordinate system, r and v are the position and velocity vectors of the satellite in the inertial coordinate system, thrust is the microthrust generated by the microthruster, acentral is the gravitational acceleration of the central celestial body, aconservative and anonconservative are the perturbation accelerations of conservative force and nonconservative force, respectively, τ is the time parameter of the first-order Markov process, m is the mass of the satellite, and w is the Gaussian noise.Compared with the traditional orbit determination problem, the state variable of this paper isx=rvthrustT. 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].
(2)thrust·t=−1τ⋅thrustt+2σ2τ,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
(3)acentral=−μr3r,where μ is the gravitational constant of the earth, and μ=3.986004415×1014m3/s2. Here, r 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
(4)aconservative=anonspherical+atide+asun&moon+arelativistic,whereanonspherical is the nonspherical perturbation, using the EGM2008 earth gravity field model, and this paper considers 70×70 order; atide is the tidal perturbation, and this paper considers the solid tide perturbation; asun&moon is the gravitational perturbation of the third body, mainly the gravitational perturbation of the sun and the moon; and arelativistic 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]
(5)V=μr∑n=2∞C¯n,0RErnP¯nsinφ+∑n=2∞∑m=1nRErnP¯nmsinφC¯nmcosmλ+S¯nmsinmλ,where φ and λ are the geocentric longitude and latitude of the satellite; C¯n,0, C¯nm, and S¯nm are the coefficients of zonal harmonics and tesseral harmonics which can be obtained from EGM08; and P¯nsinφ and P¯nmsinφ are Legendre functions. In this paper, we consider the earth nonspherical perturbations of order 70×70, that is, the case of n=70. 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
(6)anonspherical=C1C0∂V∂r∂V∂φ∂V∂λ,where C0 is the coordinate transformation matrix from the spherical coordinate system to the rectangular coordinate system [13] and C1 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 n≥3, only the additional potential function caused by the solid tide when n=2 is considered. Therefore, under the action of the sun and the moon, the additional position of the solid tide is ΔC,ΔS, so the potential function caused by the additional potential of the earth tide is ΔV.The perturbation acceleration caused by the combination of the earth nonspherical perturbation and the earth tide perturbation is(7)anonspherical&tide=C1C0∂V+ΔV∂r∂V+ΔV∂φ∂V+ΔV∂λ.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 μsun and μmoon are the gravitational constant of the sun and the gravitational constant of the moon, respectively, rs and rM are the position vectors of the sun and the moon from the center of the earth in the J2000 inertial coordinate system.
(8)asun&moon=−μsun⋅r−rsr−rs3+rsrs3+−μmoon⋅r−rMr−rM3+rMrM3.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]
(9)arelativisitic=3μr3⋅vc2⋅r,where c is the speed of light and c=3×108m/s.


## 2.1.2. Perturbation Acceleration of Nonconservative Forceanonconservative

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).
(10)anonconservative=asolarpressure+adrag.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]
(11)asolarpressure=−F⋅ρSR⋅AU2⋅CR⋅Sm⋅r−rsr−rs3,where F is the exposure factor; ρSR is the solar pressure constant, and ρSR=4.5605×10‐6N/m2; AU is the distance from the earth to the sun (one astronomical unit), and AU=1.4959787×1011m; CR is the light pressure coefficient, which is a dimensionless parameter and is related to the material of the satellite surface, and CR=1.5; and S/m is the surface quality ratio of the satellite exposed to sunlight. For the sun exposure factor F, the orbit of the satellite is low earth orbit. Considering the cylindrical shadow, the relationship between F and the position of the satellite is shown in Figure 1. The pseudocode for calculating the exposure factor F is shown in Algorithm 1.Algorithm 1: Calculate the exposure factor F.
(1) Input: r, rs Re.(2) D⟵r, rs(3) mo_r⟵r.(4) mo_ rs⟵ rs(5) if D>0(6)  thenF⟵1;(7) elseif D<=0(8)  then H2⟵mo_r2-(D/mo_rs)2;(9)  if H2>=Re2(10)    thenF⟵1;(11)  elseif H2<Re2(12)    thenF⟵1;(13)  end(14) endFigure 1
Relationship between sun exposure factorF and satellite position.For the perturbation acceleration of atmospheric drag, there is [15]
(12)adrag=−12⋅ρ⋅Sm⋅CD⋅ve⋅vewhere ρ is the atmospheric resistance density, and this article uses the NRLMISISE 00 model; S/m is the windward area of the satellite; CD is the atmospheric drag coefficient, which is related to the satellite body shape, and CD=2.2; and ve is the speed of the satellite relative to the atmosphere, and ve is the modulus of ve.(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
(13)athrust=thrustm.


## 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(14)z=Hx+v,where z is the measured value with noise, which is the position and speed of the satellite in operation,H⋅ is the function relation between the measurement value and state variable, and v 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
(15)F=12mμa3ΔaΔt,where F is the thrust, a is the semimajor axis of the orbit, and Δa is the change of the semimajor axis of the orbit after the thrust, and Δt 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:(16)ẋ=ṙv̇thrust·=vacentral+aconservative+anonconservative+athrust−1τ⋅thrust+w=fx+w,z=Hx+v.It can be seen that the model is continuous and strongly nonlinear, but the filtering algorithm is the estimation fromk to k+1, 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 x=x1,thrustT, wherex1=r,vT. After using Runge-Kutta discretization, there is
(17)k1=fx1,k2=fx1+12Δt⋅k1+12Δt,k3=fx1+12Δt⋅k2+12Δt,k4=fx1+Δt⋅k3+Δt,x1,k+1=x1,k+16Δt⋅k1+2k2+2k3+k4.After the discretization of the state differential equation of the microthruster, there is(18)thrustk+1=e−Δtτ⋅thrustk+ωk,where Δt is the sampling time and ωk is discrete driving white noise. Therefore, the discretized state-space model is
(19)xk+1=Fxk+wk,zk=Hxk+vk,where
(20)xk=x1,kthrustk.


## 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 κ=0 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 isx∧0+, the initial covariance is P0+, the process noise covariance matrix is Qk, and the measurement noise covariance matrix is Rk. It is worth noting that
(21)x=rvthrust=x1thrust,therefore, there is
(22)x0+=r0+v0+thrust0+=x1,0+thrust0+.


### 2.2.2. Time Updating

The posterior estimated value of the state variable at timek is
(23)x^k+=r^k+v^k+thrustk+∧=x^1,k+thrustk+∧,and the posterior estimated value of the covariance matrix is Pk+.
(1)
Calculating the cubature points(24)x^ki=x^k++Pk+e,−ei,i=1,2,⋯,2n,where
(25)x^ki=x^1,kithrustki∧.n is the dimension of the state variable and i is the number of cubature points, which is twice the dimension of state variable n. In the formula, ei represents the i-th column of
(26)e,−ei=10⋮001⋮0⋯00⋮1−10⋮00−1⋮0⋯00⋮−1n×2ni=1,2,⋯,2n,where e is the n-dimensional unit vector.(2)
Nonlinear propagation of cubature points(27)X1,k+1i=fx^1,ki,Thrustk+1i=e−Δt/τ⋅thrustki,∧So, there is(28)Xk+1i=X1,k+1iThrustk+1i.(3)
Calculating the predicted state(29)x^k+1−=12n∑i=12nXk+1i.(4)
Calculating the predicted error covariance(30)Pk+1−=12n∑i=12nXk+1i−x^k+1−Xk+1i−x∧k+1−T+Qk.


### 2.2.3. Measurement Updating

(1)
Calculating the cubature points(31)x^k+1i=x^k+1−+Pk+1−e,−ei,i=1,2,⋯,2n.(2)
Calculating the measurement of the propagation cubature points(32)Zk+1i=Hx^k+1i.(3)
Calculating the predicted measurement(33)z^k+1=12n∑i=12nZk+1i.(4)
Calculating the measurement error covariance matrix(34)Pz=12n∑i=12nZk+1i−z^k+1Zk+1i−z∧k+1T+Rk+1.(5)
Calculating the cross-covariance matrix(35)Pxz=12n∑i=12nXk+1i−x^k+1−Zk+1i−z∧k+1T.(6)
Calculating the Kalman gain(36)Kk+1=Pxz,k+1Pz,k+1−1.(7)
Calculating the updated state(37)x^k+1+=x^k+1−+Kk+1zk+1−z^k+1.(8)
Calculating the posterior error covariance matrix(38)Pk+1+=Pk+1−−Kk+1Pz,k+1Kk+1T.The whole process of the cubature Kalman filter is shown in the pseudocode of Algorithm2. 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.Algorithm 2:Pseudocode of the cubature Kalman filter process.
(1) Computeei⟵Ιni(2) Input:x^0+, P0+, zk+1(3)fork=1toz_lengthdo(4)Sk+⟵chol Pk+;(5)fori=1to 2ndo(6) Calculatex^ki using x^k+, Sk+, ei, Pk+;(7)Xk+11⟵fx^ki;(8)end(9)x^ki⟵1/2n∑i=12nXk+1i;(10)Pk+1−⟵1/2n∑i=12nXk+1i−x^k+1−Xk+1i−x∧k+1−T+Qk;(11)Sk+1+⟵cholPk+1+;(12)fori=1to 2ndo(13) Calculatex^k+1i using x^k+1−, Sk+1+, ei, Pk+1+;(14)Zk+1i⟵Hx^k+1i;(15)end(16)z^k+1⟵1/2n∑i=12nZk+1i;(17)Pz⟵1/2n∑i=12nZk+1i−z^k+1Zk+1i−z∧k+1T+Rk+1;(18)Pxz⟵1/2n∑i=12nXk+1i−x^k+1−Zk+1i−z∧k+1T;(19)Kk+1⟵Pxz,k+1Pz,k+1−1;(20)x^k+1+⟵x^k+1−+Kk+1zk+1−z^k+1;(21)Pk+1+⟵Pk+1−−Kk+1Pz,k+1Kk+1T;(22)end(23) Output:x^k+1+, Pk+1+Figure 2
Flowchart of the whole filtering process.


## 2.2.1. Initialization

The initial state isx∧0+, the initial covariance is P0+, the process noise covariance matrix is Qk, and the measurement noise covariance matrix is Rk. It is worth noting that
(21)x=rvthrust=x1thrust,therefore, there is
(22)x0+=r0+v0+thrust0+=x1,0+thrust0+.


## 2.2.2. Time Updating

The posterior estimated value of the state variable at timek is
(23)x^k+=r^k+v^k+thrustk+∧=x^1,k+thrustk+∧,and the posterior estimated value of the covariance matrix is Pk+.
(1)
Calculating the cubature points(24)x^ki=x^k++Pk+e,−ei,i=1,2,⋯,2n,where
(25)x^ki=x^1,kithrustki∧.n is the dimension of the state variable and i is the number of cubature points, which is twice the dimension of state variable n. In the formula, ei represents the i-th column of
(26)e,−ei=10⋮001⋮0⋯00⋮1−10⋮00−1⋮0⋯00⋮−1n×2ni=1,2,⋯,2n,where e is the n-dimensional unit vector.(2)
Nonlinear propagation of cubature points(27)X1,k+1i=fx^1,ki,Thrustk+1i=e−Δt/τ⋅thrustki,∧So, there is(28)Xk+1i=X1,k+1iThrustk+1i.(3)
Calculating the predicted state(29)x^k+1−=12n∑i=12nXk+1i.(4)
Calculating the predicted error covariance(30)Pk+1−=12n∑i=12nXk+1i−x^k+1−Xk+1i−x∧k+1−T+Qk.


## 2.2.3. Measurement Updating

(1)
Calculating the cubature points(31)x^k+1i=x^k+1−+Pk+1−e,−ei,i=1,2,⋯,2n.(2)
Calculating the measurement of the propagation cubature points(32)Zk+1i=Hx^k+1i.(3)
Calculating the predicted measurement(33)z^k+1=12n∑i=12nZk+1i.(4)
Calculating the measurement error covariance matrix(34)Pz=12n∑i=12nZk+1i−z^k+1Zk+1i−z∧k+1T+Rk+1.(5)
Calculating the cross-covariance matrix(35)Pxz=12n∑i=12nXk+1i−x^k+1−Zk+1i−z∧k+1T.(6)
Calculating the Kalman gain(36)Kk+1=Pxz,k+1Pz,k+1−1.(7)
Calculating the updated state(37)x^k+1+=x^k+1−+Kk+1zk+1−z^k+1.(8)
Calculating the posterior error covariance matrix(38)Pk+1+=Pk+1−−Kk+1Pz,k+1Kk+1T.The whole process of the cubature Kalman filter is shown in the pseudocode of Algorithm2. 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.Algorithm 2:Pseudocode of the cubature Kalman filter process.
(1) Computeei⟵Ιni(2) Input:x^0+, P0+, zk+1(3)fork=1toz_lengthdo(4)Sk+⟵chol Pk+;(5)fori=1to 2ndo(6) Calculatex^ki using x^k+, Sk+, ei, Pk+;(7)Xk+11⟵fx^ki;(8)end(9)x^ki⟵1/2n∑i=12nXk+1i;(10)Pk+1−⟵1/2n∑i=12nXk+1i−x^k+1−Xk+1i−x∧k+1−T+Qk;(11)Sk+1+⟵cholPk+1+;(12)fori=1to 2ndo(13) Calculatex^k+1i using x^k+1−, Sk+1+, ei, Pk+1+;(14)Zk+1i⟵Hx^k+1i;(15)end(16)z^k+1⟵1/2n∑i=12nZk+1i;(17)Pz⟵1/2n∑i=12nZk+1i−z^k+1Zk+1i−z∧k+1T+Rk+1;(18)Pxz⟵1/2n∑i=12nXk+1i−x^k+1−Zk+1i−z∧k+1T;(19)Kk+1⟵Pxz,k+1Pz,k+1−1;(20)x^k+1+⟵x^k+1−+Kk+1zk+1−z^k+1;(21)Pk+1+⟵Pk+1−−Kk+1Pz,k+1Kk+1T;(22)end(23) Output:x^k+1+, Pk+1+Figure 2
Flowchart of the whole filtering process.


## 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 Figure3.Figure 3
The specific simulation process.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.Table 1
Simulation parameter setting under 500μN thrust.
ParameterNumerical valueInitial stateInitial state variablex^0+x^0+=6878140,0,−1.433,0,6789.53,3686.41,0,0.000439,0.000239TInitial covariance matrixP0+P0+=0.0520000000000.0520000000000.0520000000000.000520000000000.000520000000000.0005200000000010−1600000000010−1600000000010−16CKFFirst-order Markov process parameterττ=τ1,τ2,τ3T=6.32,6.32,6.32TFirst-order Markov process parameterσσ=σ1,σ2,σ3T=5e−7,5e−7,5e−7TProcess noise covariance matrixQQ=0.0120000000000.0120000000000.0120000000000.00120000000000.00120000000000.00120000000002σ12/τ10000000002σ22/τ20000000002σ32/τ3Measurement noise covariance matrixRR=0.0520000000000.0520000000000.0520000000000.0120000000000.0120000000000.012000000000000000000000000000000The simulation step1 sThe simulation time1800 sFor the initial covariance matrixP0+, it represents the degree of thrust to the initial state variable estimate x∧0+. The closer P0+ 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 Q, it represents the accuracy of model establishment. The smaller Q 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 R, 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.Figure 4
Position of the satellite under the action of 500μN thrust.
(a)(b)(c)(d)Figure 5
Velocity of the satellite under the action of 500μN thrust.
(a)(b)(c)(d)Figure 6
Three-dimensional position of the satellite orbit and partial enlarged view.
(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 Figure7. The specific RMSE expression is shown in equation (39), which, respectively, represents the root mean square error of position and velocity at time k, where N is the run times of Monte Carlo experiments and n is the n-th Monte Carlo experiment. x,y,z and Vx,Vy,Vz are the positions and velocities of satellites obtained by simulating GNSS. x^,y^,z^ and V^x,V^y,V^z are the estimates of the position and velocity of the satellite obtained by the algorithm in this paper.
(39)RMSEposition,k=1N∑n=1Nxn−x^n2+yn−y^n2+zn−z^n2,RMSEvelocity,k=1N∑n=1NVx,n−V^x,n2+Vy,n−V^y,n2+Vz,n−V^z,n2.Figure 7
RMSE of satellite position and velocity under the action of 500μN thrust.
(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 Figure8, in which the average RMSE of the x-axis is 15.92 cm, the average RMSE of y-axis is 15.93 cm, and the average RMSE of z-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.Figure 8
RMSE of a three-axis satellite under the action of 500μN thrust.According to Figure9, 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%.Figure 9
Real-time on-orbit estimation results and errors of 500μN thrust.
(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.Figure 10
Three-dimensional position and partial enlargement of the satellite orbit under different thrusts.
(a)(b)Figure 11
Real-time on-orbit estimates and errors for different thrusts.
(a)(b)Table 2
Three kinds of thrust real-time on-orbit estimation accuracy and required stabilization time.
Thrust of microthruster (N)Real-time on-orbit thrust estimation accuracyThe stability time required for thrust estimation (s)100×10−63.98%250500×10−61.79%605×10−31.43%12


## 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.



---

*Source: 7733495-2021-07-01.xml*