Re ﬁ nement of GRACE-FO Kinematic Orbit by Combining Least Squares Method and Orbital Dynamic Model

,


Introduction
The Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) mission, successfully launched on May 22, 2018, is a joint NASA-GFZ project to continue the objectives of the original GRACE (2002-2017) mission and provide continuity for the GRACE dataset. GRACE-FO satellites are equipped with several instruments, including geodetic quality GPS receivers, accelerometers, star trackers, K-band ranging system, laser-ranging interferometer (LRI), and laser retroreflectors [1,2]. GRACE-FO satellites were injected into the 500 AE 10 km altitude, where orbit eccentricities are less than 0.005 and orbit inclinations are about 89° [3]. During the 5-year period of GRACE-FO, the satellites have performed global measurements of monthly gravity and surface mass changes, for which orbit accuracy is crucial in research. GRACE-FO orbit accuracy is mainly determined by the quality of on-board GPS data, so it is necessary to study the methods of precise orbit determination (POD) [4,5].
In the process of solving the POD of low-Earth orbit (LEO) satellites, researchers have made contributions in certain aspects [6][7][8][9][10][11][12][13][14][15][16][17][18]. Jäggi et al. [7] used the least squares (LS) method to estimate GPS-based GRACE orbits, processing GPS data in the form of satellite-ground double-differenced (DD) carrier phase. Satellite laser range (SLR) residuals indicate that the accuracy of GRACE orbits is about 2-2.5 cm. Roh and Hwang [15] research the POD of KOMPSAT-5 and GRACE by analyzing the characteristics of GPS observation. They verify that the GRACE orbital position accuracy achieves 2-3 cm. Similarly, Kang et al. [14] used LS to estimate GRACE-FO orbits, with SLR residuals showing orbit accuracy of about 2-5 cm. Guo et al. [16,17] using Bernese 5.2 to analysis the high-order ionospheric delay correction for precise reduced-dynamic determination of LEO satellite orbits. The research results indicate that the influence of third-order ionospheric delay correction on LEO satellite orbit determination is less than 0.5 mm, so it can be ignored. However, the influence of second-order ionospheric delay correction can reach the millimeter level, which means it has an important impact on LEO POD. GPS DD data can be used for GRACE-FO POD, but the accuracy is affected by ground stations and the data processing is complicated. Therefore, in this contribution, the Kalman filter is used to estimate GPSbased GRACE-FO kinematic orbits, in which the GPS data comprise undifferenced (UD) code and phase measurements from on-board GRACE-FO GPS receivers. The kinematic orbits are affected by GPS receivers, which means that missing GPS data or lower GPS data quality will decrease orbital accuracy. To improve kinematic orbit accuracy, we propose a dynamical method. First, we consider that most kinematic orbits are stable and have high accuracy (Section 3.1). Second, we use the LS method to estimate the initial orbital state and orbital dynamical parameters. Some orbit coordinates with outliers are deleted during the data processing. Final, the initial state and dynamical parameters are used to generate the dynamical orbits with high accuracy.
The paper is arranged as follows: the related model, theory, and processing strategies are described in Section 2, case studies are described in Section 3, and conclusions are provided in Section 4.

Model, Theory, and Strategy
We used the Kalman filter to process GRACE-FO on-board GPS data and estimate the orbit. To improve GRACE-FO orbit accuracy, a dynamical method is proposed. This section is arranged as follows: Kalman filter theory is described in Section 2.1, the kinematic orbit estimation strategy is described in Section 2.2, and a dynamical method for GRACE-FO is proposed in Section 2.3.

Kalman Filtering Theory
Applied to Orbit Determination of LEO Satellites. LEO on-board GPS data include codes (P1 and P2) and carrier phase (L1 and L2) observations. To eliminate the ionosphere delay, ionosphere-free (IF) combinations (code P s r;IF and carrier phase L s r;IF ) are adopted. The unknown parameters of LEO satellites include orbit coordinates, receiver clock errors, and carrier phase ambiguities. Assuming the on-board receiver observes n GPS satellites, the unknown parameters are expressed as follows: where x; y; z are LEO satellite coordinates, cdt is receiver clock error, and λ IF N s 1 r;IF ⋯ λ IF N s n r;IF are ambiguity parameters. The Kalman filter formula is expressed as follows: where D 0 is the prior variance matrix of X, B is a normal matrix, R is a covariance matrix of measurements, J is the filter-enhanced giant matrix, and X and X 0 are estimated and a priori parameters. l is the observed value minus the computed value (O-C), D is a covariance matrix of the estimated parameters, and B is a linearized matrix, which is expressed as follows: ; ð3Þ where x s ; y s ; z s are GPS satellite coordinates. l is expressed as follows: l s L;r and l s P;r are O-C of carrier phase and code (unit: m). The initial D 0 is expressed as follows: where D x ð Þ; D y ð Þ; D z ð Þ; D λ IF N s 1 r;IF À Á and D λ IF N s n r;IF À Á are initialized to 60 2 (unit: m 2 ). R is expressed as follows: ; ð7Þ 2 Journal of Sensors where R s L and R s P are covariance matrices of carrier phase and code (unit: m 2 ), el s is the elevation angle of GPS satellites, a L ¼ 0:003 and b L ¼ 0:003 are the error factors of phase observation (unit: m), and a P ¼ 0:3 and b P ¼ 0:3 are the error factors of code. In this work, pseudorange single-point positioning (SPP) is used to generate a priori satellite coordinates and clock error in X 0 , in which prior ambiguity parameters adopt floating-point values that were estimated in a previous epoch.

GRACE-FO Kinematic Orbit Determination Method
and Data Processing Strategy. LEO satellite kinematic POD includes the following steps: use pseudorange SPP to initialize satellite coordinates; detect and mark cycle slips, build normal and covariance matrices and O-C vectors, and initialize ambiguity parameters.
(1) Use SPP to initialize satellite coordinates and receiver clock error. (2) Detect and mark the cycle slip.
Carrier-phase data with cycle slips will decrease orbit accuracy. Hence, it is necessary to detect and mark the cycle slips [19][20][21][22][23][24]. We used the TurboEdit method to detect the cycle slips, which includes Melbourne-Wübbena (MW) and geometry free (GF) combinations. MW and GF are expressed as follows: The threshold of cycle slip detection between epochs (i and i + 1) is set as: (3) Initialize parameters and covariance matrices.
Parameters include satellite coordinates, receiver clock error, and ambiguities. SPP is used to initialize the satellite coordinates and receiver clock error. If there are no cycle slips, ambiguity parameters estimated in the previous epoch will be directly adopted. If there are cycle slips, ambiguity parameters are initialized as follows: (4) Use Formulas (3)-(5) to build a normal matrix and O-C vector. (5) Use Formula (2) to estimate parameters and build a covariance matrix of parameters. (6) Use the O-C vector to eliminate measurements with outliers, in which the threshold of outliers is set as: where v P is code O-C and v L is carrier phase O-C.
(7) Repeat steps 4-6 until v P and v L conform to Formula (13). Figure 1 shows how the Kalman filter is used to estimate LEO satellite kinematic orbits. SPP is used to initialize LEO satellite coordinates and receiver clock errors, so initial covariances for these coordinates and clock error are set to 60 2 per epoch.

Dynamic Smooth Orbit
Method. The speed of LEO satellites is about 6-7 km/s, and the time for an on-board GPS receiver to observe one GPS satellite is short, usually 10-30 min. Thus, initialization of carrier phase ambiguity is frequent, which leads to worse accuracy in some orbital arcs. To solve this problem, we propose a new method. First, most kinematic orbits are stable and have high accuracy. We use the LS method to process the kinematic orbits and estimate the initial orbital state and dynamical parameters of the LEO satellite [25]. Second, kinematic orbit coordinates containing outliers are eliminated when processing the data. Finally, initial orbital parameters and dynamical parameters are used to determine high-accuracy dynamical orbits.
Assume that X LEO ¼ṙ r c Estimating the dynamical orbit is a process of improving kinematic orbits and eliminating outliers in the overall orbit of the LEO satellite. The dynamical orbital observation equation V X ¼ ΦX LEO;0 − l X in matrix form is formulated as follows: where V X is the correction of dynamical orbital coordinates, I is an identity matrix, and Φ t i ; ð t 0 Þ is the state transition matrix from t 0 to t i . l X i is expressed as follows:  Journal of Sensors Figure 2 shows the dynamical method for improving kinematic orbits. The details are as follows: (1) Use the Kalman filter to estimate the kinematic orbits of LEO satellites. (2) Use Lagrange interpolation to initialize the LEO satellite's initial orbital parameters. The initial dynamical parameters can be set to 0.
ð t 0 Þ and orbit coordinates X LEO;i via numerical integration of the equations of motion of the GRACE-FO satellites. (4) Use X LEO;i and X KNA;i to estimate l X . (5) Use the LS method to estimateX LEO;0 , in which X LEO;0 represents the estimated values of X LEO;0 . (6) Determine dynamical orbits by theX LEO;0 and dynamical parameters. Compare dynamical and kinematic orbits, and then delete coordinates containing outliers in the kinematic orbits. The threshold for orbital outliers is set as follows: (7) Repeat steps (2)-(4) until all dynamical orbit coordinates conform to Formula (17).

Case Study and Analysis
To verify the dynamical method used for GRACE-FO POD, a control experiment was designed based on actual data from GRACE-FO satellites, in which the experimental group achieved POD by the dynamical method and the control group used the kinematic method. The experimental data included 28 days of GRACE-FO on-board GPS data provided by GFZ (https://isdc.gfz-potsdam.de/), collected from day 32 to 59 in 2022; precise ephemeris, satellite clock offsets, and Earth rotation parameters provided by the Center for Orbit Determination in Europe (CODE), and broadcast ephemeris provided by the International GNSS Service (IGS). The data are listed in Table 1, and the dynamical models are listed in Table 2.
3.1. Reference Orbit Check. The orbits provided by Jet Propulsion Laboratory (JPL) are regarded as reference orbits. To assess orbit accuracy, GRACE-FO kinematic and dynamical orbits were compared with JPL orbits. Figures 3 and 4 show that kinematic orbit residuals in the radial, along, and cross directions are primarily distributed between −0.15 and þ0.15 m, while dynamical orbit residuals are mainly distributed between −0.05 and þ0.05 m. There are no systematic errors in kinematic and dynamical orbits. However, some kinematic orbit residuals are larger than 0.15 m, which means there are outliers. By comparison, the dynamical orbits are smoother, verifying that the dynamical method can eliminate outliers and improve orbital accuracy.
To obtain more detailed statistical information, we calculated the mean, min, max, and RMS values of residuals in the radial, along, and cross directions. Tables 3 and 4 list the statistics of the 28-day POD residuals of the GRACE-FO satellite, including the following: (1) RMS values of GRACE-FO-C and -D kinematic and dynamical orbits indicate that compared with kinematic orbits, dynamical orbit accuracy of C and D satellites improved 52.4% and 35.11% in the radial direction, and 50.01% and 33.96% in the cross direction, respectively. The statistics show kinematic and dynamical orbit accuracy of around 9 and 6 cm, respectively; the mean value of  orbital residuals in the radial, along, and cross directions was 0 m, indicating there are no biases in the kinematic and dynamical orbits and residuals conform to random distributions. However, the max and min values of kinematic orbit indicate that some residuals are larger than 0.15 m. Compared with kinematic orbit, the dynamical orbit accuracy was improved by 3 cm in the 3D direction, verifying that the dynamical method can be used to improve the accuracy of kinematic orbits. The max values of dynamical orbital residuals show that the dynamical method can weaken the influence of outliers and make the orbit smoother.
3.2. SLR Check. As a measurement method with high accuracy, satellite laser ranging is widely used to monitor GNSS and LEO satellite orbits. GRACE-FO-C and -D satellites are monitored by the International Laser Ranging Service (ILRS) SLR tracking stations, and the SLR range observations can be accessed through FTP services provided by the ILRS crustal dynamics data information system (CDDIS) and EUROLAS data center (EDC). Between days 32 and 59 in 2022, GRACE-FO-C and -D satellites were observed by 18 SLR stations, producing 3,893 and 3,929 normal points (NPs).
For SLR values, we subtracted the ranges between stations and LEO satellite orbits to get residuals. We used these SLR residuals to assess the accuracy of LEO satellite orbits and estimate the RMS of 16 SLR station residuals over 28 days. SLR residual RMS values were estimated with a 10°elevation cutoff. SLR residuals larger than 0.15 m were considered outliers. Figures 5 and 6 show the SLR residual RMS of GRACE-FO-C and D satellites. Excluding SLR stations 1824 and 7249, the SLR RMS for GRACE-FO-C dynamical orbits is between 2.0 and 6.0 cm, with most RMS values less than 5 cm. Even with stations 1824 and 7249, the RMS of SLR reached 7.7 and 8.1 cm. The SLR RMS for GRACE-FO-D dynamical orbits is mainly between 2.6 and 6.0 cm. Tables 5 and 6 show the following:

Conclusions
Poor on-board GPS data quality decreases the accuracy of LEO kinematic orbits. To solve this problem, a dynamical method was proposed. The method uses LS to estimate the initial state of the LEO satellite orbit and dynamical parameters. Coordinates containing outliers are eliminated or weakened during data processing, and dynamical orbits with higher accuracy are determined by the initial orbital and dynamical parameters. JPL orbits and SLR data were used to verify the dynamical orbit accuracy. Three conclusions can be made: (1) JPL orbit comparisons show that GRACE-FO-C and -D dynamical orbital residuals in the radial, along, and cross directions are between 2.9 and 4.7 cm, the residuals conform to random distribution, and orbital accuracy is achieved at the centimeter level.

Ethical Approval
Not applicable for studies not involving humans or animals.

Consent
Not applicable for studies not involving humans.