**1. Introduction**

Gravity Recovery and Climate Experiment Follow-on (GRACE-FO) is regarded as a new mission of gravity formation satellites launched jointly by the National Aeronautics and Space Administration (NASA) and Helmholtz-Centre Potsdam, German Research Centre for Geosciences (GFZ), with a design life of 5 years. Having been launched successfully at Vandenberg Air Force Base in California on May 22, 2018, GRACE-FO satellites attempt to replace GRACE satellites that had orbited for 15 years and retired in June, 2017 [1–3]. GRACE-FO satellites carry the same equipment as GRACE, including GPS, a K-Band Ranging (KBR) System, a satellite accelerometer, and star sensor [4,5]. Similarly, GRACE-FO satellites also adopt an orbit design similar to that of GRACE satellites, with an orbit height of about 500 ± 10 km, an orbit eccentricity of less than 0.005, and an orbit inclination of about 89◦ [5,6]. GRACE-FO satellites, mainly used to accurately measure and invert the time-varying Earth's gravitational field, are regarded as an important mission for collecting

**Citation:** Yang, Z.; Liu, X.; Guo, J.; Guo, H.; Li, G.; Kong, Q.; Chang, X. Relative Kinematic Orbit Determination for GRACE-FO Satellite by Jointing GPS and LRI. *Remote Sens.* **2022**, *14*, 993. https:// doi.org/10.3390/rs14040993

Academic Editor: Xiaogong Hu

Received: 7 December 2021 Accepted: 15 February 2022 Published: 17 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Earth's gravitational field data of high resolution [2]. The high quality satellite orbit and relative position can ensure the high quality processing of gravity data, which contributes to the estimation of gravity field models with a monthly resolution [7,8]. Therefore, the study of the precise orbit determination of GRACE-FO satellites is considered crucial. To ensure a stable relative distance in formation flight, GRACE satellites are equipped with a KBR system. In addition, GRACE-FO satellites have on-board a Laser Ranging Interferometer (LRI) [9–12]. LRI can accurately measure the intersatellite distance of GRACE-FO formation satellites and provide distance data of high quality [13].

At present, several formation satellites have been performing in orbit, including GRACE, GRACE-FO, TanDEM-X, PRISMA, SJ-9 and TechSat21 [14–21]. The existing formation satellites mainly realize relative positioning based on a GPS technique [19]. Generally, the relative kinematic orbit determination of LEO satellites is based on carrierphase differential observations (CDGPS) [14]. The differential technique promises to eliminate and weaken some common observation errors, such as receiver clock bias, satellite clock bias and ionospheric delay errors [19–30]. Gu et al. [31] combined single and double differential techniques to jointly solve the orbit of formation satellites. A comparison with satellite laser ranging (SLR) indicates that the accuracy of the orbit was improved by 25%. Van Barneveld [32] analyzed the influence of ionospheric delay on the calculation of intersatellite long baselines (>100 km) and compared existing ionospheric delay models.

The new generation GRACE-FO satellites are equipped with an LRI system, which makes it possible to use LRI observations to enhance the quality and stability of formation satellites orbit. Currently, the intersatellite distance are generally used to verify the intersatellite baseline, and to maintain the relative state of formation satellites. GRACE-FO satellites can be equipped with dual ranging system (LRI and KBR). Concerning the LRI ranging system, this manuscript describes the following activities: firstly, the inter-satellite distance observed by LRI is combined with the GPS observation data of the GRACE-FO satellites to form the GPS/LRI observations. Secondly, the estimated orbits are compared to reference orbits and KBR measurements. Thirdly, this study analyzed the relative kinematic orbit determination results of GPS/LRI and GPS only.

#### **2. Mathematical Model of GPS and LRI**

LEO satellites' orbit determination generally takes pseudo range and carrier phase of GPS as the main observed values. The pseudo range and carrier phase observation models are as follows:

$$P\_{r,j}^{s} = \rho\_r^s + c(dt\_r - dt^s) + I\_{r,j}^{s} + \varepsilon\_P^s \tag{1}$$

$$L\_{r,j}^{s} = \rho\_r^s + c(dt\_r - dt^s) - I\_{r,j}^s + \lambda\_j N\_{r,j}^s + \lambda\_j \varepsilon\_{L'}^s \tag{2}$$

where *P<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is the observed code pseudo range; *<sup>ρ</sup><sup>s</sup> <sup>r</sup>* is the geometric distance between the GRACE-FO satellite to the GPS satellite; *c* is the velocity of light; *dtr* and *dt<sup>s</sup>* are receivers and satellite clock offsets, respectively; *I<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is the ionospheric delay; *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is the observed value of carrier phase; *λjN<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is the ambiguity; and *<sup>λ</sup>jε<sup>s</sup> <sup>L</sup>* and *<sup>ε</sup><sup>s</sup> <sup>P</sup>* are observation noise. To eliminate ionospheric delay, ionospheric-free combinations are usually adopted:

$$P\_{r,IF}^{s} = \rho\_r^s + c\left(dt\_r - dt^s\right) + \varepsilon\_{P,IF}^s \tag{3}$$

$$L\_{r,IF}^{s} = \rho\_r^s + c(dt\_r - dt^s) + \lambda\_{IF} \mathcal{N}\_{r,IF}^s + \varepsilon\_{L,IF\prime}^s \tag{4}$$

where *P<sup>s</sup> <sup>r</sup>*,*IF* is the pseudo range of ionospheric-free combinations; *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*IF* is the carrier phase of ionospheric-free combinations; *λIF* is the wavelength of carrier phase *LIF*, *N<sup>s</sup> <sup>r</sup>*,*IF* is the ambiguity of ionospheric-free combinations, and *ε<sup>s</sup> <sup>r</sup>*,*IF* and *<sup>ε</sup><sup>s</sup> <sup>P</sup>*,*IF* are observation noise of pseudo range and carrier phase ionospheric-free combinations. The equations of ionospheric-free combinations are as follows:

$$L\_{IF} = \frac{f\_1^{'2}}{f\_1^{'2} - f\_2^{'2}} L\_1 - \frac{f\_2^{'2}}{f\_1^{'2} - f\_2^{'2}} L\_{2'} \tag{5}$$

where *f*<sup>1</sup> and *f*<sup>2</sup> are the frequencies of carrier phases *L*<sup>1</sup> and *L*2. Ionospheric-free combinations can eliminate most of the ionospheric delay, and loses integer characteristic.

To build double-difference observation equations, inter-satellite single-difference equations should first be constructed. The single-difference models are as follows:

$$\begin{cases} \nabla P\_{r,IF}^{s\_0,s} = \rho\_r^s - \rho\_r^{s\_0} - \mathfrak{c}(dt^s - dt^{s\_0})\\ \nabla L\_{r,IF}^{s\_0,s} = \rho\_r^s - \rho\_r^{s\_0} - \mathfrak{c}(dt^s - dt^{s\_0}) + \lambda\_{IF} (N\_{r,IF}^s - N\_{r,IF}^{s\_0})\\ \nabla N\_{r,IF}^{s\_0,s} = N\_{r,IF}^s - N\_{r,IF}^{s\_0} \end{cases} \tag{6}$$

where <sup>∇</sup>*Ls*0,*<sup>s</sup> <sup>r</sup>*,*IF* and <sup>∇</sup>*Ps*0,*<sup>s</sup> <sup>r</sup>*,*IF* are single-difference observed values of the carrier phase and pseudo range. The satellite *s*<sup>0</sup> with the largest altitude angle is selected as the reference satellite. The double-difference models are as follows:

$$\begin{cases} \Delta \nabla P^{\circ\_0, \sf s}\_{r\_0, r, IF} = \nabla \rho^{\sf s\_0, \sf s}\_r - \nabla \rho^{\sf s\_0, \sf s}\_{r\_0} \\ \Delta \nabla L^{\sf s\_0, \sf s}\_{r\_0, r, IF} = \nabla \rho^{\sf s\_0, \sf s}\_r - \nabla \rho^{\sf s\_0, \sf s}\_{r\_0} + \lambda\_{IF} (\nabla N^{\sf s\_0, \sf s}\_{r, IF} - \nabla N^{\sf s\_0, \sf s}\_{r\_0, IF}) \\ \Delta \nabla N^{\sf s\_0, \sf s}\_{r\_0, r, IF} = \nabla N^{\sf s\_0, \sf s}\_{r, IF} - \nabla N^{\sf s\_0, \sf s}\_{r\_0, IF} \end{cases} \tag{7}$$

where double-difference observed values of <sup>Δ</sup>∇*Ps*0,*<sup>s</sup> <sup>r</sup>*0,*r*,*IF*, <sup>Δ</sup>∇*Ls*0,*<sup>s</sup> <sup>r</sup>*0,*r*,*IF* and <sup>Δ</sup>∇*Ns*0,*<sup>s</sup> <sup>r</sup>*0,*r*,*IF* are doubledifference pseudo range, carrier phase and ambiguity, respectively.

For the first time, GRACE-FO formation satellites are equipped with a LRI ranging system. The model of LRI is as follows [5]:

$$\Psi = c \ast \left( -\psi\_M(t) + R\_T(t - \pi) \right) / (2f), \tag{8}$$

where Ψ is the LRI measurement (unit: m), *ψM*(*t*) is the observed phase value of LRI at *t*, *RT*(*t* − *τ*) is the time delay correction of the LRI acquisition system, *τ* is the time of signal transmission of LRI, and *f* is the LRI-frequency of the observed phase. The ground nominal value of *f* in GRACE-FO-A is 2.81616393e14 Hz, and the ground nominal value of *f* in GRACE-FO-B is 2.81615684e14 Hz. The original phase observed values (LRI1A) of LRI have glitches [33]. Glitches of LRI mainly occurs during thruster firings. JPL laboratory detects glitches of LRI by three methods. Firstly, differences in phase observed values are performed, and then glitches of LRI are detected by comparing the observed differential values with the tolerance. Secondly, small glitches of LRI are fitted by a phase filter model [33]. Thirdly, glitches are further detected by residuals of Two-Way Range [11]. In order to simplify experiments, this study adopts secondary-treated Level-1B data of JPL laboratory, and all glitches of LRI were detected and repaired.

Based on the principle of LRI and product description published by GFZ, LRI observations deliver a biased distance *NLRI*. To obtain the actual intersatellite distance, the biased distance must be gained in advance. This study proposes a simple and feasible calculation process for the biased LRI distance.

(1) The intersatellite distance is resolved by using relative kinematic orbit determination from spaceborne GPS data. Assuming that there are n epochs, the intersatellite distance between two LEO satellites *ρ<sup>i</sup> <sup>r</sup>*0,*r*(*i* = 0 ··· *n*) is obtained according to the relative kinematic orbit determination of spaceborne GPS;

(2) LRI range minus the intersatellite distance *ρ<sup>i</sup> <sup>r</sup>*0,*<sup>r</sup>* at the corresponding epoch, and the initial value *N<sup>i</sup> LRI*,0 of LRI biased distance can be introduced;

(3) The mean value *NLRI* of *N<sup>i</sup> LRI*,0(*i* = 0 ··· *n*) should be calculated;

(4) The difference between the initial value (*N<sup>i</sup> LRI*,0(*i* = 0 ··· *n*)) and the mean value (*NLRI*) is calculated, and epochs with a difference greater than a threshold are marked as outliers. This threshold is defined as follows:

$$\left| \mathcal{N}\_{LRI,0}^{i} - \overline{\mathcal{N}}\_{LRI} \right| < 3\sigma\_{rel} \tag{9}$$

where *σrel* is the precision of relative kinematic orbit determination of spaceborne GPS. For the convenience of the experiment, *σrel* was set to 0.15 m (*σrel* = 0.15 m);

(5) The marked observed values of LRI should be eliminated, and then steps from 2 to 4 are repeated until all *N<sup>i</sup> LRI*,0 pass the inspection formula;

(6) Taking into account the LRI bias, ranging values of each epoch are corrected to obtain the intersatellite distance of LEO satellites based on LRI ranging.

To evaluate the precision of LRI (actually LRI biased), the LRI biased distances for 4 months (days 121 to 242 of 2019) are calculated. The statistical results are shown in Table 1, confirming the good quality of LRI distances. More than 98.5% of LRI observations can be used to solve LRI bias.


**Table 1.** RMS of LRI bias (unit: m) and percentage of LRI outliers (unit: %).

#### **3. GPS/LRI Observation Equation and Linearization**

As observed distance values, LRI can establish the observation of the LEO satellite by combining with GPS. The GPS/LRI observation equation in matrix form reads as follows:

$$V = BX - l,\tag{10}$$

where *l* is obtained by subtracting the calculated values from actual observed values, *V* is the correction of observed values, and *B* is the linearized parameter coefficient matrix to be solved. *B* is expressed as follows:

$$B = \begin{bmatrix} a\_x^{s\_0 s\_1} & a\_y^{s\_0 s\_1} & a\_z^{s\_0 s\_1} & \lambda & 0 & \cdots & 0 \\ a\_x^{s\_0 s\_2} & a\_y^{s\_0 s\_2} & a\_z^{s\_0 s\_2} & 0 & \lambda & 0 & 0 \\ \vdots & \vdots & \vdots & 0 & 0 & \ddots & 0 \\ a\_x^{s\_0 s\_{n-1}} & a\_y^{s\_0 s\_{n-1}} & a\_z^{s\_0 s\_{n-1}} & 0 & 0 & 0 & \lambda \\ a\_x^{s\_0 s\_1} & a\_y^{s\_0 s\_1} & a\_z^{s\_0 s\_1} & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ a\_x^{s\_0 s\_{n-1}} & a\_y^{s\_0 s\_{n-1}} & a\_z^{s\_0 s\_{n-1}} & 0 & 0 & 0 & 0 \\ a\_x^{L R I} & a\_y^{L R I} & a\_z^{L R I} & 0 & 0 & 0 & 0 \end{bmatrix},\tag{11}$$

where *a <sup>s</sup>*0,*<sup>s</sup> <sup>x</sup>* <sup>=</sup> <sup>−</sup>*xs*0−*xr ρ s*0 *r* + *<sup>x</sup>s*−*xr ρs <sup>r</sup>* ; *a <sup>s</sup>*0,*<sup>s</sup> <sup>y</sup>* <sup>=</sup> <sup>−</sup>*ys*0−*yr ρ s*0 *r* + *<sup>y</sup>s*−*yr ρs <sup>r</sup>* ; *a <sup>s</sup>*0,*<sup>s</sup> <sup>z</sup>* <sup>=</sup> <sup>−</sup>*zs*0−*zr ρ s*0 *r* + *<sup>z</sup>s*−*zr ρs r* ; *aLRI <sup>x</sup>* <sup>=</sup> *xr*−*xr*<sup>0</sup> *<sup>ρ</sup>r*0,*<sup>r</sup>* ; *<sup>a</sup>LRI <sup>y</sup>* <sup>=</sup> *yr*−*yr*<sup>0</sup> *<sup>ρ</sup>r*0,*<sup>r</sup>* ; *<sup>a</sup>LRI <sup>z</sup>* <sup>=</sup> *zr*−*zr*<sup>0</sup> *<sup>ρ</sup>r*0,*<sup>r</sup>* ; *<sup>ρ</sup>r*0,*<sup>r</sup>* <sup>=</sup> - (*xr* − *xr*<sup>0</sup> ) <sup>2</sup> + (*yr* <sup>−</sup> *yr*<sup>0</sup> ) <sup>2</sup> + (*zr* <sup>−</sup> *zr*<sup>0</sup> ) 2 ; *ρs*0 *<sup>r</sup>* and *ρ<sup>s</sup> <sup>r</sup>* are the geometric distances from GRACE-FO to different GPS satellites. *xr*, *yr*, *zr*, *xs*<sup>0</sup> , *ys*<sup>0</sup> , *zs*<sup>0</sup> and *x<sup>s</sup>* , *y<sup>s</sup>* , *z<sup>s</sup>* are Cartesian Coordinates. The absolute vector *l* is given by:

$$I = \begin{bmatrix} l\_{r0,r}^{\mathbf{s}\_1,\mathbf{s}\_1} \cdot \cdots \cdot l\_{r0,r}^{\mathbf{s}\_{0,\mathbf{s}\_m-1}} \, \_\prime p\_{r0,r}^{\mathbf{s}\_0,\mathbf{s}\_1} \cdot \cdots \cdot p\_{r0,r}^{\mathbf{s}\_{0,\mathbf{s}\_m-1}} \, \_\prime \upsilon\_{LRI} \end{bmatrix} \tag{12}$$

$$\begin{cases} l\_{r\_0,r}^{s\_0,s} = \rho\_r^s - \rho\_r^{s\_0} - \rho\_{r\_0}^s + \rho\_{r\_0}^{s\_0} - \Delta \nabla L\_{r,IF}^{s\_0,s} \\ p\_{r\_0,r}^{s\_0,s} = \rho\_r^s - \rho\_r^{s\_0} - \rho\_{r\_0}^s + \rho\_{r\_0}^{s\_0} - \Delta \nabla P\_{r,IF}^{s\_0,s} \\ v\_{LRI} = \rho\_{r0,r} - \Psi \end{cases} \tag{13}$$

After the observation equation is linearized, a Kalman filter process can be used for iterative processing.

#### **4. Kalman Filter Theory Applied to Relative Kinematic Orbit Determination of LEO Satellites**

Relative Kinematic orbit determination of LEO satellites aims to obtain coordinates of satellites. Under the condition of double-difference orbit determination, the parameters to be obtained not only include the coordinates of satellites, but also the ambiguities of double-difference carrier observed values. The parameters of LEO satellites are to be calculated by setting *n* double-difference carrier observed values:

$$X = [\Delta x, \Delta y, \Delta z, \Delta \nabla N^{s\_0, s\_1}\_{r\_0, r, IF} \cdot \cdot \cdot \Delta \nabla N^{s\_0, s\_n - 1}\_{r\_0, r, IF}]\_\prime \tag{14}$$

where Δ*x*, Δ*y*, Δ*z* are components of the baseline vector between the GRACE-FO satellites (Earth-fixed reference frame) and <sup>Δ</sup>∇*Ns*0,*s*<sup>1</sup> *<sup>r</sup>*0,*r*,*IF* ··· <sup>Δ</sup>∇*Ns*0,*sn*−<sup>1</sup> *<sup>r</sup>*0,*r*,*IF* are ambiguities of the observed carrier values. The Kalman filter equations of LEO satellites are as follows:

$$\begin{cases} \begin{aligned} J &= D\_0 B^T \left( B D\_0 B^T + R \right)^{-1} \\ X &= X\_0 + Jl \\ D &= \left( D\_0^{-1} + B^T R^{-1} B \right)^{-1} \end{aligned} \end{cases} \tag{15}$$

where *D*<sup>0</sup> is the prior variance matrix of parameters; *R* is the observation variance matrix; *J* is the gain matrix; *X*<sup>0</sup> and *X* are the prior value and calculated value of parameters to be solved, respectively; *l* is the observed minus computed (O-C); and *D* is the posterior variance matrix of parameters.

Pseudo-range single point positioning is used to replace the state-transition matrix to generate prior coordinates of LEO satellites (prediction coordinates), with an accuracy of about 10 m, which requires continuous refinement of LEO satellites' coordinates according to observed carrier values and their corresponding variance information. The variance processing of parameters is regarded as an important basis for Kalman filter calculation. Therefore, this manuscript provides the detailed settings for the corresponding parameter variance matrix and observation variance matrix in the Kalman filter. The variance matrix formula of observed values *R* is as follows:

$$\begin{aligned} R \\ R\_{(2n-2)\times(2n-2)+1} \end{aligned} = \begin{bmatrix} R\_L^{s\_0} + R\_L^{s\_1} & R\_L^{s\_0} & R\_L^{s\_0} & 0 & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ R\_L^{s\_0} & R\_L^{s\_0} & R\_L^{s\_0} + R\_L^{s\_{n-1}} & 0 & 0 & 0 \\ 0 & 0 & 0 & R\_P^{s\_0} + R\_L^{s\_1} & R\_P^{s\_0} & R\_P^{s\_0} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & R\_P^{s\_0} & R\_P^{s\_0} & R\_P^{s\_0} + R\_L^{s\_{n-1}} \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \tag{16}$$

$$\begin{cases} \begin{array}{c} R\_L^s = 2 \times (a\_L^2 + \frac{b\_L^2}{\sin^2 \varepsilon l^s})\\ R\_P^s = 2 \times (a\_P^2 + \frac{b\_P^2}{\sin^2 \varepsilon l^s}) \end{array} . \end{cases} . \tag{17}$$

*Rs <sup>L</sup>* and *<sup>R</sup><sup>s</sup> <sup>P</sup>* are the phase prior variances and pseudo-range prior variances, *el<sup>s</sup>* is the elevation angle of GPS satellites, *aL* = 0.003 and *bL* = 0.003 are the error factors of phase observation (unit: m), *aP* = 0.3 and *bP* = 0.3 are the error factors of pseudo-range observation (unit: m), and *RLRI* is the variance of LRI's observed values. In the experiment of this manuscript, *RLRI* = 0.012 (unit: m2).

In this work, Kalman filter orbit determination does not take dynamic information into consideration, and pseudo-range single point positioning is used to generate a priori coordinate *X*0. The prior value of the ambiguity parameter adopts the floating-point value resolved in the previous epoch. The initial-state parameter variance matrix *D*<sup>0</sup> of LEO satellites is assigned as follows:

$$D\_0 \begin{array}{c} D\_0 \\ \cdot \\ \cdot (n+3) \times (n+3) \end{array} = \begin{bmatrix} 30^2 \\ & \cdot \\ & \cdot \\ & & 30^2 \end{bmatrix} \prime \tag{18}$$

where the variance corresponding to the coordinate parameter is reset as 302(unit: m2) in the next epoch, while the variance corresponding to the ambiguity parameter directly adopts the posterior variance of previous epochs.

#### **5. Relative Kinematic Orbit Determination Process of GPS/LRI**

The essence of relative kinematic orbit determination of spaceborne GPS/LRI is to regard LRI as a directly observed value, and build the observation equation with GPS double-difference observations for resolving the intersatellite baseline. The relative kinematic orbit determination process of GPS/LRI is shown in Figure 1.

**Figure 1.** Orbit determination process.

	- (2) Screening common-view satellites of GRACE-FO-A and GRACE-FO-B;

(3) Detecting cycle slips according to Melbourne–Wübbena (MW) and geometry-free (GF) combination [34,35]. The models of MW and GF are provided below:

$$L\_{MW} = \left(f\_1 L\_{r,1}^s - f\_2 L\_{r,2}^s\right) / \left(f\_1 - f\_2\right) - \left(f\_1 P\_{r,1}^s + f\_2 P\_{r,2}^s\right) / \left(f\_1 + f\_2\right) \tag{19}$$

$$L\_{GF} = L\_{r,1}^s - L\_{r,2}^s. \tag{20}$$

Equations (19) and (20) are usually used to jointly detect cycle slips. In this paper, the threshold of cycle slips detection between epochs (*i* and *i* + 1) is set as:

$$\begin{cases} N\_w^{i+1} - N\_w^i < 5\\ L\_{GF}^{i+1} - L\_{GF}^i < 0.05 \end{cases} (unit: m). \tag{21}$$

The detection for cycle slips is regarded as data preprocessing;

(4) Constructing double-difference observation equation according to Equations (1) to (7); (5) Updating the parameters to be resolved by Kalman filter;

(6) Correcting phase center offsets (PCO) and Sensor Offsets of LEO satellite. Taking GRACE-FO-A and GRACE-FO-B as examples, Tables 2 and 3 show the three components of PCO and Sensor Offsets in the body- fixed coordinate of GRACE-FO satellites. As for LEO satellites, their solar panel must always aim at the sun, so the satellite gesture will change according to time. Therefore, several steps should be taken in calculating PCO and Sensor Offsets correction of LEO satellites with spaceborne GPS:

(a) The coordinates of LEO satellites without being corrected by PCO and Sensor Offset are converted from Earth-fixed coordinate system to inertial coordinate system.

(b) According to the attitude of LEO satellites, the three components of PCO and Sensor offset in the body- fixed coordinate system are transferred to the inertial coordinate system.

(c) The corrections of PCO and Sensor Offset in the inertial system are added to the coordinates of LEO satellites.

(d) The coordinates of LEO satellites are transferred from the inertial coordinate system to the Earth-fixed coordinate system;


**Table 2.** PCO of GRACE-FO satellites (unit: mm).

**Table 3.** Sensor Offset of GRACE-FO satellites (unit: mm).


(7) The intersatellite baseline is recalculated according to the correction of mass center, and the relative kinematic orbit determination results are generated.

#### **6. Case Study and Analysis**

To verify the influence of LRI on the relative kinematic orbit determination of GRACE-FO formation satellites, a control experiment is designed according to actual observation data of GRACE-FO satellites. The experimental group achieves relative kinematic orbit determination by adopting GPS/LRI combinations, and the control group only uses GPS. The paper adopts reference orbit comparison and KBE check to analyze and compare the

experimental results. The observation data of GPS/LRI and KBR in GRACE-FO satellites are provided by GFZ (https://isdc.gfz-potsdam.de/, accessed on 9 June 2020), collected from day 121 to day 242 in 2019, a total of 121 days. Samples of GPS data are collected at the interval of 10 s, and those of LRI and KBR data are collected at the interval of 2 s and 5 s, respectively. Then, synchronous GPS, LRI and KBR data should be used for check. The precise ephemeris, satellite clock offsets and Earth rotation parameter files are provided by the Center for Orbit Determination in Europe (CODE), and broadcast ephemeris are provided by the International GNSS Service (IGS). More details about these data are listed in Table 4.



#### *6.1. Reference Orbit Check*

The orbit published by GFZ is regarded as the reference value. The relative orbit determination results of GPS/LRI and GPS are separately compared with the reference orbit. Figures 2–6 show the comparison results of days 121–242's reference orbits. As is shown, there are no systematic errors in the GRACE-FO satellites' relative kinematic orbit determination. The main reason is that the relative kinematic orbit determination method used in this paper essentially belongs to a geometric orbit determination without being affected by the dynamic satellite model. However, the residual errors of GPS/LRI joint relative kinematic orbit determination result in X, Y, and Z (components of the baseline vector) directions are obviously smaller than those of GPS only. Compared with relative kinematic orbit determination by GPS only, the distribution of orbit residuals calculated by GPS/LRI is closer, and the obvious spikes in the figure are weakened. The main reason is that on-board GPS is quite different from ground stations in that the number of satellites observed by the on-board GPS receiver in some arcs is smaller than five. Meanwhile, the travel speed of LEO satellites can reach 6–7 km per second, so the observed GPS satellites change frequently. Usually, a satellite can be observed for only 10 to 30 min. Therefore, the observation data quality is poor for a small part of the arc, reducing the stability of the geometric method for relative kinematic orbit determination. The Geometric Dilution Precision (GDOP) is related to the orbit determination accuracy. We analyzed the GDOP of GRACE-FO satellites for 121 to 242 days. When the value of GDOP increases, the orbit accuracy will decrease. However, GPS / LRI can limit the impact of GDOP increase on the orbit accuracy. Therefore, combined GPS/LRI observations can better eliminate GPS observations of poor quality and reduce the impact of GPS observation quality on relative kinematic orbit determination. Additionally, the results of 121 days' relative kinematic orbit determination show that introducing GPS/LRI observations is obviously better than using GPS only.

**Figure 2.** Reference orbit check for GRACE-FO in the X, Y, Z, 3D (X, Y, Z are components of the baseline vector, unit: m, (**a**–**d**)) directions (days 121–151), number of satellites for each epoch (**e**) and GDOP (**f**,**g**).

**Figure 3.** Reference orbit check for GRACE-FO in the X, Y, Z, 3D (X, Y, Z are components of the baseline vector, unit: m, (**a**–**d**)) directions (days 151–182), number of satellites for each epoch (**e**) and GDOP (**f**,**g**).

**Figure 4.** Reference orbit check for GRACE-FO in the X, Y, Z, 3D (X, Y, Z are components of the baseline vector, unit: m, (**a**–**d**)) directions (days 182–212), number of satellites for each epoch (**e**) and GDOP (**f**,**g**).

**Figure 5.** Reference orbit check for GRACE-FO in the X, Y, Z, 3D (X, Y, Z are components of the baseline vector, unit: m, (**a**–**d**)) directions (days 212–242), number of satellites for each epoch (**e**) and GDOP (**f**,**g**).

**Figure 6.** Days 121–242 orbit comparison residuals (in the X, Y and Z directions, unit: m) probability distribution for the GPS/LRI and GPS relative orbit.

To obtain more specific results of relative kinematic orbit determination, statistics over 121 days are carried out in this part. X, Y, and Z refer to the baseline vectors of the two satellites in Earth-fixed coordinates. It is shown in Table 5 that along the X direction, the 121 days accuracy of GPS/LRI relative kinematic orbit determination is improved by 12.4 mm, and the relative kinematic orbit determination accuracy along the X direction reaches 25.4 mm, with an accuracy increase of 32.8%; along the Y direction, the accuracy of GPS/LRI relative kinematic orbit determination is improved by 8.5 mm, and the accuracy along the Y direction reaches 28.6 mm, with an accuracy increase of 22.9%; along the Z direction, the accuracy of GPS/LRI relative kinematic orbit determination is improved by 15.2 mm, and the accuracy along the Z direction reaches 21.9 mm, with an accuracy increase of 41.0%; along the 3D direction, the accuracy of GPS/LRI relative kinematic orbit determination is improved by 17.5 mm, and the accuracy along the 3D direction reaches 50.1 mm, with an accuracy increase of 25.9%. LRI observed values are added to increase the number of redundant observations, enhancing the geometric strength of observed values, so the accuracy of GPS/LRI is greatly improved compared with that of only GPS for formation orbit determination. Statistical results show that the joint GPS/LRI data can effectively improve the overall accuracy of relative kinematic orbit determination of GRACE-FO formation satellites and limit orbit determination errors due to GPS observation quality, improving the stability of relative kinematic orbit determination.

**Table 5.** Statistics of reference orbit check (days 121–242) in Earth-fixed coordinate (unit: m).


To analyze the orbital residual, we set GRACE-FO-B as the reference orbit and solve the orbit of GRACE-FO-A according to the relative kinematic orbit determination results. The orbit of the GRACE-FO-A satellite and the orbit released by GFZ are compared by using the ORBCMP module of Bernese 5.2 software. Under the radial, along–track and

out–of–plane (RSW) decomposition, we analyze the residual results. Figures 7–10 show the comparison results of day 121–242's reference orbits. Figure 11 shows the probability distribution of satellite orbit residuals. In the R and W directions, the residual distributions of GPS-only and GPS/LRI are similar. However, in the S direction, GPS/LRI can significantly improve the relative kinematic orbit determination accuracy of GRACE-FO-A satellite. For further analysis, we calculated the MEAN, MEDIAN and RMS of orbital residuals in R, S and W directions. It is shown in Table 6 that along the R direction, the 121 days accuracy of GPS/LRI relative kinematic orbit determination is improved by 10.0 mm, and the relative kinematic orbit determination accuracy along the R direction reaches 29.0 mm, with an accuracy increase of 34.48%; along the S direction, the accuracy of GPS/LRI relative kinematic orbit determination is improved by 31.7 mm, and the accuracy along the S direction reaches 8.9 mm, with an accuracy increase of 78.8%; along the W direction, the accuracy of GPS/LRI relative kinematic orbit determination is not significantly improved. The RMS of the 3D vector is decreased by 18 mm which corresponds to an improvement of 26.3%.

**Figure 7.** Reference orbit check for GRACE-FO-A in the R, S, W, 3D directions (days 121–151, (**a**–**d**)).

**Figure 8.** Reference orbit check for GRACE-FO-A in the R, S, W, 3D directions (days 151–182, (**a**–**d**)).

**Figure 9.** *Cont*.

**Figure 9.** Reference orbit check for GRACE-FO-A in the R, S, W, 3D directions (days 182–212, (**a**–**d**)).

**Figure 10.** Reference orbit check for GRACE-FO-A in the R, S, W, 3D directions (days 212–242, (**a**–**d**)).

**Figure 11.** Days 121–242 orbit comparison residuals (in the R, S and W directions, unit: m) probability distribution for the GPS/LRI and GPS relative orbit.


**Table 6.** Statistics of reference orbit check (days 121–242) in RSW (unit: m).

To evaluate the GPS/LRI under sunlight and solar eclipse, we use the ORBGEN module of Bernese 5.2 software to mark the epochs in solar eclipse. GRACE-FO satellites were under solar eclipse for 25% of the total period. Figure 12 show the comparison results of day 121–242's reference orbits in the sunlight and solar eclipse. In the case of sunlight and solar eclipse, there is no significant difference in GPS/LRI orbit residuals. In both cases, the accuracy of GPS/LRI orbit determination is higher than the GPS-only. To further analyze the impact of solar eclipse on GPS/LRI, we calculated the RMS, MEAN and MEDIAN of orbital residuals. It is shown in Tables 7 and 8 that the RMS values differ by 2.2 mm and 2.1 mm in the 3D-directions. Along the X, Y, Z and 3D direction the 121-days accuracy of GPS/LRI relative kinematic orbit determination, under the conditions of sunlight and solar eclipse, is not a significant difference.

**Table 7.** Statistics of X, Y, Z and 3D residual (unit: m, days 121–242) in Sunlight and Solar eclipse.


**Table 8.** Statistics of X, Y, Z and 3D residual (unit: m, days 121–242) in Sunlight and Solar eclipse.


**Figure 12.** D (X, Y, Z are components of the baseline vector, unit: m, (**a**–**d**)) directions (days 121–242).

To evaluate the GPS/LRI in satellite orbit maneuver, we analyzed the orbit eccentricity of GRACE-FO satellites for 121 days. Landerer points out that regular orbit maintenance maneuvers will be performed throughout the GRACE-FO mission duration [36]. Figure 13 shows the orbital eccentricity of GRACE-FO satellites and their hourly variation. We marked the time of orbital maneuver in the figure with circles. Table 9 shows the date of the satellite orbit maneuver. Figure 14 shows the comparison results of reference orbits in the orbit maneuver. Table 10 shows the RMS, MEAN and MEDIAN of orbital residuals during orbital maneuver.


**Figure 14.** Reference orbit check for Orbital Maneuver in the X, Y, Z, 3D (X, Y, Z are components of the baseline vector, unit: m, (**a**–**d**)) directions (days 121–242).


**Table 10.** Statistics of X, Y, Z and 3D residual (unit: m, days 121–242) in Orbital Maneuver.

#### *6.2. KBR Validation*

The KBR system is one of the key scientific instruments carried by GRACE-FO satellites with a sampling interval of 5 s, ranging accuracy of 10 um, and range-changing rate accuracy of 1 um/s [37]. Actually, KBR and LRI are distinctly separate intersatellite ranging systems. Therefore, KBR can be used as an important external condition to verify GPS/LRI relative kinematic orbit determination results. The basic principle of KBR check is to obtain the difference between KBR-observed values and intersatellite baseline, and then evaluate the quality of the orbit difference between LEO satellites. Figures 15–19 show the residual distribution of KBR from day 121 to 242 in 2019. It shows that the relative kinematic orbit determination residuals by GPS only are mainly distributed within ±0.05 m, and the fluctuation of orbit determination residuals presents random distribution with

some residual spikes. The residual values of a few epochs exceed ±0.075 m. However, the relative kinematic orbit determination results of joint GPS/LRI are quite smooth with the residual fluctuation range within ±0.01 m, and only very few epochs have KBR residual values greater than ±0.025 m. Compared with relative kinematic orbit determination by GPS, relative kinematic orbit determination by joint GPS/LRI reduces the orbital spikes of relative kinematic orbit determination caused by GPS observation quality, and improves the orbital accuracy of relative kinematic orbit determination.

**Figure 15.** Day 121–151 KBR residuals (unit: m) for the GPS/LRI and GPS-only relative orbit and number of satellites.

**Figure 16.** Day 151–182 KBR residuals (unit: m) for the GPS/LRI and GPS-only relative orbit and number of satellites.

**Figure 17.** Day 182–212 KBR residuals (unit: m) for the GPS/LRI and GPS-only relative orbit and number of satellites.

**Figure 18.** Day 213–242 KBR residuals (unit: m) for the GPS/LRI and GPS-only relative orbit and number of satellites.

**Figure 19.** Day 121–242 KBR residuals (unit: m) probability distribution for the GPS/LRI and GPS-only relative orbit.

To obtain more detailed statistical information, we have conducted for all 121 days a comparison between KBR ranges and the relative kinematic orbit determination results. RMS, MEAN and MEDIAN values of 121 days KBR residuals are listed in Table 11. The accuracy of 121 days' relative kinematic orbit determination based on GPS only approximates 42.8 mm while that calculated from joint GPS/LRI data approximates 10.7 mm. It indicates that joint GPS/LRI improves the quality of relative kinematic orbit determination, and confirms the results of reference orbit verification. The mean value of KBR residuals in 121 days is less than 0.5 mm, indicating no obvious systematic errors in joint GPS/LRI and GPS relative kinematic orbit determination, and the results of relative kinematic orbit determination comply with the statistical law.

**Day GPS-Only GPS**\**LRI RMS MEAN MEDIAN RMS MEAN MEDIAN** 121–151 0.0415 0 0 0.0112 0 −0.0007 151–182 0.0407 0 0 0.0093 0 −0.0006 182–212 0.0458 0 0 0.0086 0 −0.0007 212–242 0.0433 0 0 0.0135 0 −0.0007 Average 0.0428 0 0 0.0107 0 −0.0007

**Table 11.** Statistics of KBR check (unit: m).

KBR can be used to evaluate the ability of GPS/LRI under sunlight and eclipse conditions. Figure 20 shows the KBR residuals of GPS/LRI and GPS-only relative kinematic orbit determination. There is no significant difference in the KBR residuals of GPS/LRI relative kinematic orbit determination under sunlight and solar eclipse. Table 12 shows the RMS, MEAN and MEDIAN of KBR residuals for GPS/LRI and GPS-only. Under the sunlight and solar eclipse, the KBR residuals RMS corresponding to GPS/LRI are 9.3 mm and 9.7 mm, and their differences are very slight. The RMS values of KBR residuals corresponding to GPS-only are 42.4 mm and 42.2 mm. Their differences are also small.

**Figure 20.** Days 121–242 KBR residuals (unit: m) for the GPS/LRI (**a**) and GPS-only (**b**) relative orbit in the Sunlight and Solar eclipse).


#### **7. Conclusions**

In this manuscript, we studied the LRI ranging system carried by GRACE-FO formation satellites, achieving high-quality relative kinematic orbit determination based on joint GPS/LRI data and the Kalman filter. Additionally, we introduced the LRI ranging systems as well as the LRI observation equation, and described the process of relative kinematic orbit determination of GPS/LRI in detail. Ultimately, relative kinematic orbit determination results of GPS/LRI were verified and analyzed by comparing to a reference orbit and verifying KBR. Accordingly, we draw two conclusions: 1. Compared with relative kinematic orbit determination by GPS only, kinematic orbit determination by GPS/LRI achieves more robust results, which can weaken some orbits' error and significantly improve the accuracy of relative kinematic orbit determination. The accuracy of relative kinematic orbit determination by GPS/LRI is improved by 17.5 mm in 3D directions, and the accuracy of relative kinematic orbit determination is improved by 25.9%, reaching 50.1 mm (compared with the reference orbit released by GFZ). 2. The results of KBR validation are ideal. The residual distribution of KBR is smooth without obvious fluctuation. The consistency between the relative kinematic orbit determination and the KBR measurements is at the +/−10.7 mm level.

**Author Contributions:** Conceptualization, writing—original draft preparation, writing—review and editing, and methodology, Z.Y., X.L. and J.G.; software, Z.Y.; validation, H.G., G.L. and Q.K.; data curation, X.C.; formal analysis, H.G.; funding acquisition, J.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by The National Natural Science Foundation of China (42174041, 41774001), the Independent and Controllable Project of National Surveying and Mapping (816-517) and the SDUST Research Fund (2014TDJH101).

**Institutional Review Board Statement:** Not applicable for studies not involving humans or animals.

**Informed Consent Statement:** Not applicable for studies not involving humans.

**Data Availability Statement:** Our sincere thanks go to the German Research Centre Geosciences(GFZ) for providing space-borne GPS data for GRACE-FO (ftp://isdcftp.gfz-potsdam.de/grace-fo, accessed on 9 June 2020); CODE and IGS for providing GPS satellite orbits, clocks and Earth rotation parameters; the GFZ for providing precise orbits, KBR and LRI data for GRACE-FO.

**Acknowledgments:** This research is supported by National Natural Science Foundation of China (grant numbers 42174041, 41774001) and SDUST Research Fund (grant number2014TDJH101).

**Conflicts of Interest:** The authors declare no conflict interest.

#### **References**

