Improvement of Orbit Prediction Algorithm for Spacecraft Through Simplified Precession-Nutation Model Using Cubic Spline Interpolation Method

: For the on-orbit flight missions, the model of orbit prediction is critical for the tasks with high accuracy requirement and limited computing resources of spacecraft. The precession-nutation model, as the main part of extended orbit prediction, affects the efficiency and accuracy of on-board operation. In this paper, the previous research about the conversion between the Geocentric Celestial Reference System and International Terrestrial Reference System is briefly summarized, and a practical concise precession-nutation model is proposed for coordinate transformation computation based on Celestial Intermediate Pole (CIP). The idea that simplifying the CIP-based model with interpolation method is driven by characteristics of precession-nutation parameters changing with time. A cubic spline interpolation algorithm is applied to obtain the required CIP coordinates and Celestial Intermediate Origin locator. The complete precession nutation model containing more than 4000 parameters is simplified to the calculation of a cubic polynomial, which greatly reduces the computational load. In addition, for evaluating the actual performance, an orbit propagator is built with the proposed simplified precession-nutation model. Compared with the orbit prediction results obtained by the truncated series of IAU2000/2006 precession-nutation model, the simplified precession-nutation model with cubic spline interpolation can significantly improve the accuracy of orbit prediction, which implicates great practical application value in further on-orbit missions of spacecraft.


Introduction
The broadly shared interest of concise coordinate transformation between the inertial system and Earth-fixed system is motivated by the need to increase the efficiency of extended orbit prediction, devoted to realize the applications as autonomous constellation maintenance and onboard orbit determination of spacecraft. The balance between computational complexity and accuracy affects the choice of coordinate conversion model. Among all the factors that affect coordinate conversion, precise computation of the bias-precession-nutation parameters is the most demanding on computation time and memory [1]. But the full model is not necessary for every application, such as real-time, on-board orbit determination. In order to simplify the calculation of precession-nutation model, interpolation algorithm constitutes an appealing option with low computing load and high accuracy in this paper. The calculation of the bias-precession-nutation parameters in orbit prediction is regarded as the sampling of the Celestial Intermediate Pole (CIP) coordinates and the Celestial Intermediate Origin (CIO) locator. These variables are calculated at a certain time interval and stored in advance. The values in the required time sequence are obtained by interpolation algorithm. The results of simulation show that the simplified model proposed in this paper is suitable for the extended orbit prediction with high accuracy demand under the condition of limited resources.
Since 1988, several recommended models have been published by the International Earth Rotation Service (IERS) on the transformation between the Geocentric Celestial Reference System (GCRS) and International Terrestrial Reference System (ITRS) [2]. The latest IAU2000A/2006 model provides two implementing methods: traditional equinox-based method and new CIP-based method [3]. There is a tiny difference between the results obtained by these two methods [4]. The CIP-based method is strongly advocated by IAU because it requires only 3 angles instead of 4. And the bias-precession-nutation matrix in CIP-based model is formed directly, thus reducing potential errors [5].
The largest computational burden for CIP-based transformation is the computation of the CIP coordinates and CIP locator. The full accuracy method consists of around 4000 coefficients to compute the three angles in CIP-based matrix [6]. The large number of coefficients and trigonometric evaluations make the parameters evaluation burdensome, both on computation time and memory. Wallace et al. [7] had made remarkable contribution to the simplification of the precession-nutation model. The series are truncated to speed up computation at the expense of accuracy. The research lists 3 typical concise model with 1 mas, 16 mas, and 380 mas accuracy respectively, which is 15, 138, and 890 times quicker than the full model [8]. These concise forms are widely accepted and become the most commonly method.
In Global Navigation Satellite System (GNSS) related applications, sufficient research has been performed to increase the sampling rate of data such as satellite position and clock bias with interpolation algorithm. Wang et al. [9] applied the Lagrange interpolation, Newton interpolation, Hermite interpolation and spline interpolation to the calculation of satellite position respectively and compared the results. Chang et al. [10] used sliding method to improve the accuracy of interpolation. Yuan et al. [11] applied these common methods to extrapolate and predict satellite clock bias, and the research showed that the prediction results obtained by cubic spline interpolation have highest accuracy and smoothness. Given that the orbit prediction is performed with a numerical integration method, the precession-nutation matric is calculated at a set of discrete time, which can be considered as the sampling of the precession and nutation. The simulation shows that interpolation method, which has been widely used in GNSS applications, is applicable to the calculation of the bias-precession-nutation parameters. Coppola et al. [12] proposed a tabulationinterpolation method to reduce the computation load of CIP-based precession-nutation model. Bradley [13] proposed another method based on polynomial fitting, which can greatly reduce the calculation. In this paper, the application of interpolation method in orbit prediction is further analyzed. The proposed method can dramatically decrease the computational complexity caused by precession-nutation model. The simulation shows that the result is extremely close to that obtained by the full model. CMES, 2020, vol.125, no.2 867 The structure of this paper is organized as follows: Section 2 lists the basic mathematical models in orbit propagation and coordinate transformation. In Section 3, the concise precessionnutation model in the orbit prediction algorithm are presented. Afterwards, the simulation is presented in Section 3, analyzing the impact of the selection of the interpolation nodes and comparing the results of interpolated series and truncated series. Conclusions are drawn in Section 4.

Mechanical Models and Numerical Integrator for Orbit Prediction
In the inertial coordinate system, the differential kinematic equation of the spacecraft can be described as: where r is the positon of spacecraft in the inertial coordinate system, a earth is the acceleration caused by the gravitation of the Earth, a others is the acceleration caused by other perturbations. For GNSS satellites, the orbit prediction model in this paper only takes the third-body perturbations (i.e., acceleration caused by the gravitation of the Sun and the Moon) and solar radiation pressure (SRP).
For the gravitational model, in the Earth-center Earth-fixed (ECEF) frame, the geopotential field V at the point (r, θ, λ) can be calculated by spherical harmonics with degree N: where r is the distance from satellite to the centroid of the Earth, θ and λ is the longitude and latitude respectively. According to EGM2008 model, the values of GM and a e are 398,600.4415 km 3 /s 2 and 6,378,136.3 m respectively. C nm and S nm are the normalized geopotential coefficients and P nm are the normalized associated Legendre functions [14]. The EGM2008 model provides the potential coefficients of up to 2159 degree. For GNSS orbits, IERS recommends a 12-order truncating model. Denoting the transformation matrix from the ECEF frame to the inertial frame as T, the gravitational acceleration of the Earth in the inertial frame can be obtained: For spacecraft in LEO, the gravitational perturbation from the Sun and the Moon is relatively small. However, as the orbital altitude increases, the magnitude of the third-body perturbation becomes larger. For navigation satellites, the influence of the third-body gravitational perturbation is obvious [15]. This paper regards the Sun and the Moon as mass points, and the acceleration caused by the third body can be obtained in the inertial coordinate system: The SRP is another practical consideration especially for the satellites with larger area-mass ratio. Considering the SRP is determined by a combination of many factors such as the intensity of the light, the effective area of the satellite, and the ability of the material to absorb and reflect light, it is hard to obtain an accurate physical model. Thus, this paper uses a two-parameter empirical model for navigation satellites [17]: where, λ ∈ [0, 1] is a function describing the positional relationship between the satellite and the Earth's shadow. λ = 0 indicates that the satellite is completely in the Earth's shadowed area, and λ = 1 indicates that the satellite is completely out of the Earth's shadowed area. α 1 and α 2 are the empirical parameters estimated by fitting the calculated orbit data to the precise ephemeris.

Coordinate Transformation and CIP Based Precession-Nutation Model
Since the kinematic equations of spacecraft are based on Newtonian mechanics, the spacecraft dynamic model needs to be built in an inertial coordinate system. One of the most commonly used inertial coordinate systems for orbit prediction is the GCRS prescribed by the International Astronomical Union (IAU) [18]. Nevertheless, the related applications of the satellite navigation and positioning system are established in the ECEF frame, such as the WGS84 coordinate system established for GPS system. The ECEF coordinate systems are based on the ITRS [19]. Therefore, the coordinate transformation between GCRS and ITRS is one of indispensable parts of satellite navigation applications. From another aspect, the required Earth gravity model is established in the ITRS. In each step of the integral calculation, the GCRS position coordinates obtained in the previous step need to be converted into the ITRS to calculate the geopotential [20]. The obtained gravitational vector needs to be converted into the GCRS for the integral calculation [21]. Obviously, the coordinate transformation between GCRS and ITRS largely determines the accuracy and speed of the orbit prediction calculation.
The origins of GCRS and ITRS are both located in the centroid of the Earth. Hence, it only needs three rotations to achieve coordinate transformation, which are accordingly precessionnutation, earth rotation and pole motion. Since the coordinate axes of the GCRS and the ITRS are artificially defined and have no observable physical characteristics, there is no way to directly measure the angles between the two coordinate systems. To this end, two intermediate coordinate systems, the Terrestrial Intermediate Reference System (TIRS) and the Celestial Intermediate Reference System (CIRS), are introduced to realize the transformation scheme from the ITRS to the GCRS as shown in Fig. 1.
where, Q (t) is the rotation matrix generated by the movement of CIP in GCRS (i.e., precession and nutation), R (t) is the rotation matrix due to the rotation of the earth, and W (t) is the rotation matrix due to the movement of CIP in ITRS (i.e., polar motion).
In the CIP based precession-nutation model, the coordinate transformation matrix between CIRS and GCRS is determined by two factors: the Celestial Intermediate Origin (CIO) locator denoted as s, and the coordinates of CIP in GCRS denoted as (X , Y ). The precession-nutation matrix is calculated as follows: where X , Y and s are functions of time defined as follows: Among the three factors (i.e., CIO locator, CIP coordinates, and transformation matrix) having effect on the calculation of the precession-nutation model, the calculation of CIP coordinates has the greatest influence on the accuracy and speed, followed by the calculation of CIO locator. Therefore, the main purpose of simplifying the precession-nutation model is to simplify the calculation of X , Y and s.

Simplification of the Precession-Nutation Model with Cubic Spline Interpolation
As for the precession-nutation, IAU2000/2006 precession-nutation model has good long-term prediction accuracy. The difference between the calculated and the measured value from the International Earth Rotation and Reference Systems Service (IERS) can be guaranteed within 0.2 mas. The error of GNSS orbit caused by this difference is within 2.6 cm [23]. Nevertheless, in order to achieve this level of performance, the model is required to contain more than 6000 amplitude coefficients and more than 1500 angles [24]. In mobile devices or embedded systems with limited computing and storage resources, such as onboard computer on the Cube-Sat, excessive parameters will make a large impact on the computational efficiency. Thus it is necessary to simplify the precession-nutation model according to the accuracy requirements of the specific application.
The simplification of the calculation of CIP coordinates and CIO locator in the IAU2000 precession-nutation model is usually implemented by neglecting the terms with smaller amplitude. Drawing lessons from the ephemeris interpolation method commonly used in GNSS applications, this paper applies interpolation method to simplify the precession-nutation model. Fig. 2 shows the calculated values of CIP coordinates and CIO locator from March 25 to April 24, 2019. It is obvious that the parameters change smoothly, indicating that it is feasible to obtain the parameters by interpolation. The CIP coordinates and the CIO locator is calculated with a certain time interval previously and stored. During the integral process, the pre-stored parameters are interpolated to obtain the values of the CIP coordinates and the CIO locator at every step of integration. There are some common interpolation algorithms including Lagrange interpolation, Hermite interpolation, cubic spline interpolation, etc. Lagrange interpolation constructs a polynomial by selecting appropriate nodes, but there is a Runge phenomenon leading to poor convergence. Hermite interpolation requires not only the same function value, but also the same derivative value. Cubic spline interpolation requires a second-order continuous derivative function, having better convergence and smoothness [25]. According to the characteristics of the curve of X , Y and s, the constructed interpolation function needs excellent smoothness. Thus, the cubic spline interpolation method is chosen to achieve the simplification.
It is supposed that there is an set of discontinuities {x 1 , x 2 , . . . , x n } on domain [a, b] meeting a = x 1 < x 2 < . . . < x n = b, and the corresponding function values are denoted as {y 1 , y 2 , . . . , y n }, the cubic spline interpolation function needs to satisfy the following conditions: S x j = y j , j = 1, 2, 3, . . ., n; the function is a cubic polynomial on domain x j , x j+1 ; the function has a second-order continuous derivative.
According to the definition, S (x) is a linear function. The expression of S x j can be obtained with twice integrals. The integral constant term can be obtained in terms of S x j = y j . Denoting the second-order derivative value at interpolation nodes as S x j = M j , S x j can be expressed as: where h j = x j+1 −x j is the interpolation step. Since the first-order derivative function is continuous, the following equations can be obtained: Increasing boundary conditions that S ( another two equations can be obtained: Above all, Eqs. (10)-(12) constitute a linear equation set with n elements. Since the coefficient matrix of this liner equation set is a three-diagonal matrix, the Tridiagonal Matrix Algorithm (TDMA) can be applied to solve the value of M j . Substitute the obtained value to Eq. (9) and finish the construction of the interpolation function.
The selection of the interpolation nodes determines the accuracy and speed of the interpolation algorithm. More nodes will slow down the calculation, but not necessarily leading to explicit accuracy improvement, even leading to a slight decrease in accuracy. Thus, the appropriate selection of nodes is a pivotal issue to establish an efficient interpolation algorithm. Note that the interpolation function has large error at both ends of the sample area, the interpolated point should locate in the middle of the area. Therefore, the idea of sliding interpolation is applied to make sure that there are the equal number of nodes on both sides of the interpolated point. During the construction of interpolation function, n nodes before and after the calculated point in the timeline are selected (i.e., the interpolation function is constructed by using 2n interpolation nodes).
Moreover, when calculating the points between two nodes, the selected nodes to construct the interpolation function are the same. Thus, as the process shown in Fig. 3, the construction of the interpolation function only needs to be performed once a day, which can greatly reduce the computing load.

The Selection of Interpolation Nodes
In order to determine the optimal number of interpolation nodes, it needs to compare the errors of X , Y and s caused by different number of nodes. X , Y and s of 24:00 every day are calculated previously and stored, and the interpolation functions are built with 4, 6, 8, 12, and 14 nodes respectively. The calculation of full model is based on Standards of Fundamental Astronomy (SOFA) [26]. The results from April 5 to April 15, 2019 are shown in Fig. 4 and Tab. 1.
The best result is attained by the concise model with error closest to 0. The comparison shows that the minimum error is generated by interpolation function built with 6 nodes. If there are more than 6 nodes, the increasing computational complexity cannot bring increasing accuracy. Orbit prediction with the precession-nutation model simplified by constructing interpolation function with 6 nodes will be more cost efficient. Therefore, the pre-storage nodes of 3 days before and after are chosen to establish the simplified model.

Comparison of the Values of CIP Coordinates and CIO Locator Obtained by Interpolation Method and Truncated Series
The extremely concise formulation proposed by Capitaine (i.e., the CPN d in [8]), whose accuracy is acceptable for many simple applications, retains only six parameters. CPN d can meet the requirements in most applications and can attain the most saving. Taking into account a higher accuracy requirement, the interpolating model is compared with the truncated series retaining 42 parameters (i.e., the CPN c in [8]) in this paper.  The values of CIP coordinates and CIO locator from April 5 to April 6, 2019 from are calculated by full accuracy IAU2000 model, CPN c and interpolation respectively, and results are shown as follow.
It is shown in Fig. 5 that compared to truncating simplified model, interpolating simplified model is obviously more consistent with full accuracy model. Interpolation function can fill points between daily tabulated X , Y and s very well. This is due to the fact that the bias-precessionnutation parameter within 6 days is very close to a polynomial function, and even close to a linear function most of the time. Tab. 2 also lists the mean errors caused by two methods. The error generated by the interpolation method is significantly smaller than that generated by the truncation method.   The results in Tab. 3 are calculated with a typical desktop computer. The main computation load of interpolation method is the construction of the interpolation function. Once the interpolation function is constructed, calculations at the remaining points time in the same interval are only cubic polynomials. This can explain that when the step size is larger, the intepolation method takes much more time cost than the truncated series. For the application with a step size less than 10 s, the time cost of truncated series and interpolation method is almost the same. Therefore, in application scenarios it requires small step size, the advantage of interpolation method in computing speed is more prominent.

Comparison of the Values of CIP Coordinates and CIO Locator Obtained by Interpolation Method and Truncated Series
In order to intuitively demonstrate the practical value of interpolation method, the orbit errors caused by different simplified methods are presented to evaluate the simplified model proposed in this paper.
The broadcast ephemeris data of PRN 14 GPS satellite is selected as the initial value of the orbit prediction [27]. It runs on an orbit with altitude of 20,182.642 km, inclination of 55 • , and eccentricity of 0.0095. The initial time of simulation is 03:30 on April 5, 2019. The orbit in the next 24 h is calculated with the orbit prediction algorithm established in this paper. The initial conditions and related parameters in this simulation are listed in Tab. 4.   For the coordinate transformation, the complete IAU2000 precession-nutation model, the interpolating simplified model and the truncating simplified model are applied to the orbit prediction respectively and the orbit prediction data obtained by using the full accuracy model is regarded as the true value. In the interpolating simplified model, the sampling interval of the previously stored nodes is 1 day. Recalling the 1-min step size in integral, it needs to calculate the values of CIP coordinates and CIO locator at every 1 min with cubic spline interpolation method. The relative errors of the interpolating simplified model are compared with that of the truncating simplified model. The normalized computing time of propagator with the full model, interpolation method and truncated series is about 17.3, 1.25 and 1 respectively. Thanks to the way it is used in the procedure, the interpolation method shows similar performance to CPN c in resources saving. It can be seen from Fig. 7 that the accuracy of the simplified model using interpolation is significantly better than that using truncated series. Tab. 5 lists the mean value and standard deviation of orbit error. The difference of relative errors caused by the two methods is up to 2 orders of magnitude. Compared with the CPN c , the interpolation model can greatly improve the accuracy of the orbit prediction. The error caused by interpolation method is small enough to ignore, which shows that the interpolation algorithm is appropriate to save the computing cost especially when a higher accuracy is required.

Conclusion
In order to meet the requirement of low computational load deriving from the computing and storage capacity of the mobile devices and embedded systems, an innovative method to simplify the IAU2000 CIP based precession-nutation model is proposed, and an orbit prediction propagator for GNSS satellites is established. The characteristics of the variation of the CIP coordinates and CIO locator shows that the sliding cubic spline interpolation with 6 nodes is most appropriate for simplifying the related parameters. Compared with the used truncated series, the accuracy of orbit prediction with interpolation method is greatly improved. The study demonstrates that the cubic spline interpolation can dramatically decrease the computing complexity and satisfy the accuracy demand of the coordinate transformation in modern navigation satellite application. The proposed method not only can be implemented for the orbit prediction of GNSS satellites, and also suitable for LEO satellites, serving orbit tasks such as autonomous safety monitoring during the reconfiguration of satellite formation.
Funding Statement: The authors would like to express gratitude for supporting funding from the Natural Science Foundation of China (No. 51905272).