*Article* **GRACE-FO Antenna Phase Center Modeling and Precise Orbit Determination with Single Receiver Ambiguity Resolution**

**Biao Jin 1,2,3, Yuqiang Li 1,\*, Kecai Jiang 4, Zhulian Li <sup>1</sup> and Shanshan Chen <sup>3</sup>**


**Abstract:** Precise knowledge of the phase center location of the global navigation satellite system (GNSS) antenna is a prerequisite for precise orbit determination (POD) of the low Earth orbit (LEO) satellite. The phase center offset (PCO) and phase center variation (PCV) values for the LEO antenna obtained from ground calibration cannot reflect the error sources encountered in the actual spacecraft environment. PCV corrections are estimated by ionosphere free (IF) carrier phase post-fit residuals of reduced dynamic orbit determination. Ambiguity resolution (AR) plays a crucial role in achieving the best orbit accuracy. The single receiver AR concept is realized using wide-lane (WL) and narrow-lane (NL) bias products. Single difference (SD) observations between satellites are applied to remove the receiver dependent phase bias. SD AR and traditional double difference (DD) AR methods are applied to fix the ambiguities. The recovered SD and DD IF ambiguities are taken as pseudoobservations to constrain the undifferenced IF ambiguity parameters in the POD process. The LEO orbits based on float ambiguity (FA), SD, AR, and DD AR are investigated. One year's data collected by the Gravity Recovery And Climate Experiment Follow-On (GRACE-FO) mission and GPS precise products provided by the Center for Orbit Determination in Europe (CODE) were analyzed. Precise orbit generated by the Jet Propulsion Laboratory (JPL), independent satellite laser ranging (SLR), and K-band ranging (KBR) measurements were utilized to assess the orbit accuracy. More than 98% of SD WL and 95% of SD NL ambiguities are fixed, which confirms the good quality of the bias products and correctness of the SD AR method. With PCV corrections, the average phase residuals of DD and SD AR solutions are 0.13 and 0.41 mm, which indicates improved consistency between applied models and observations. Compared with JPL's orbit, the SD AR orbits achieve the accuracy of 6.0, 6.2, and 5.1 mm in along-track, cross-track, and radial directions. The SD AR solutions show an average improvement of 18.3% related to the FA orbits while 6.3% is gained by the DD AR approach. The root mean squares (RMSs) of SLR residuals for FA, DD AR, and SD AR solutions are 11.5, 10.2, and 9.6 mm, which validate the positive effect of AR on POD. Standard deviation (STD) of KBR residuals for SD AR orbits is 1.8 mm while 0.9 mm is achieved by the DD AR method. The explanation is that the phase bias products used for SD AR are not free of errors and the errors may degrade the KBR validation. In-flight PCV calibration and ambiguity resolution improve the LEO orbit accuracy effectively.

**Keywords:** single receiver ambiguity resolution; phase center variation (PCV) calibration; precise orbit determination; GRACE-FO satellites

#### **1. Introduction**

Low Earth orbit (LEO) satellites are considered as key technologies for space missions due to their advantages of flexibility, redundancy, efficiency, and low cost. To fulfill scientific mission requirements, precise absolute or relative positions are required. The precise orbit determination (POD) capability, based on spaceborne Global Positioning System (GPS)

**Citation:** Jin, B.; Li, Y.; Jiang, K.; Li, Z.; Chen, S. GRACE-FO Antenna Phase Center Modeling and Precise Orbit Determination with Single Receiver Ambiguity Resolution. *Remote Sens.* **2021**, *13*, 4204. https:// doi.org/10.3390/rs13214204

Academic Editor: Stefania Bonafoni

Received: 20 August 2021 Accepted: 14 October 2021 Published: 20 October 2021

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

**Copyright:** © 2021 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/).

observations, has been successfully verified on the TOPEX/Poseidon satellite [1,2]. The reduced dynamic POD technique [3,4] combined with the strength of GPS observations has shown its advantages in sampling rate and accuracy. Since then, more GPS receivers have been deployed on LEO satellites to meet the needs of various scientific missions for high precision orbit.

The quality of global navigation satellite system (GNSS) derived LEO orbits has steadily improved thanks to numerous improvements in the GNSS orbit and clock products provided by the International GNSS Service (IGS) [5], in the dynamic background models, and in modeling the carrier phase observations. LEO orbit accuracy of 1–3 cm can be fulfilled with float ambiguity (FA) [6,7]. With GPS precise orbit and clock products provided by IGS, CHAMP orbit calculated with ionosphere free (IF) observations achieves the accuracy of centimeter level [8]. GRACE orbit determined using GPS has been proven by independent satellite laser ranging (SLR) and K-band ranging (KBR) data, and 1 cm radial orbit accuracy has been achieved [9]. JASON's GPS derived orbit is evaluated based on SLR residuals, dynamic orbit produced by SLR/DORIS data, and altimeter crossover tests [10,11], and the sub-centimeter radial accuracy is verified. Similar results are also obtained in other LEO satellite missions, such as GOCE [12,13], Sentinel [14,15], TerraSAR-X [7], HY [16], and FY-3 [17,18] satellites. Ambiguity resolution (AR) and in-flight antenna phase center variation (PCV) calibration are essential to fully exploit the precision of GPS observations for POD.

In order to fix the carrier phase ambiguities, a Kalman filter modeling the relative spacecraft dynamics has been developed for the GRACE mission and the double difference (DD) ambiguities are resolved to fully exploit the inherent measurement accuracy. Finally, the resulting GRACE orbit matches the KBR measurements with an accuracy of 1 mm [19]. A GRACE POD based on undifferenced and doubly differenced GPS data is studied [6]. Different baselines, including the space baseline between two GRACE satellites, the spaceground baselines consisting of LEO satellites and ground stations, and both types of baselines together, are processed. Results show that fixing of the GPS DD ambiguities has a significant impact on the space baseline [6]. With the orbit solution constrained by resolved DD ambiguities, a GRACE baseline accuracy of 2 mm is also achieved [20].

The double difference ambiguity resolution (DD AR) approach requires two satellites to construct the baseline. Several methods have been developed to resolve the integer ambiguity for the single receiver user. The idea of a single receiver integer ambiguity resolution forms the basis of PPPRTK. Taking the ionospheric delay as the unknown parameter, the common clock model [21,22] and distinct clocks model [23,24] are proposed to perform the single receiver ambiguity resolution. For IF formulation, the integer recovery clock (IRC) model [25,26], decoupled satellite clock (DSC) model [27], and uncalibrated phase delay (UPD) or fractional cycle bias (FCB) model [28,29] are presented. The S-system theory is applied by Teunissen to establish the linkage among different PPP-RTK methods [30,31]. Their differences are shown to lie (a) in the choice of S-basis; (b) in the choice of parameterization; (c) in the choice of whether or not to eliminate the ionospheric delay.

With wide-lane (WL) and phase biases provided by the Jet Propulsion Laboratory (JPL), the single receiver AR of GRACE and JASON-2 satellites is realized [32]. GRACE baseline accuracy is improved from 6 mm of float ambiguity solution to 2 mm of ambiguity fixed solution. Based on GPS orbit, clock, and WL bias products provided by Centre National d'Études Spatiales (CNES) [26], single receiver AR is also employed to generate the orbit of Sentinel-3A satellite and the root mean square (RMS) of SLR residuals decreases from 9 to 5 mm [15]. Using the observation specific bias (OSB) products provided by the Center for Orbit Determination in Europe (CODE), AR is applied to the POD of GRACE and Sentinel-3 satellites [33]. Ambiguity fixing improves the orbit quality with validation of KBR and SLR residuals, as well as the internal consistency between the reduced dynamic and kinematic orbits. With UPDs estimated via a global distributed network, kinematic orbits of Sentinel-3A and Swarm-A satellites are determined with AR [34]. SLR residuals show an improvement of 20% for AR solution when compared with the FA solution. Different AR solutions, including DD AR, single difference (SD) AR, and integrated SD and DD AR solutions are investigated with GRACE data to access their effects on orbit accuracy [35].

Another factor affecting the orbit accuracy of the LEO satellite is the phase center position of the spaceborne GNSS antenna. Nominal antenna models obtained from ground calibration have been made available for the antennas deployed on several space missions, such as GRACE [36], GOCE [37], Swarm [38], and TerraSAR-X [39]. Such nominal antenna models, however, do not reflect the influence of error sources, which are additionally encountered in the actual spacecraft environment, e.g., the influence of near-field multipath [36]. In-flight calibration of the LEO antenna is necessary for the stringent orbit accuracy requirement, especially for the altimeter mission. Two different approaches, the residual approach and the direct approach, can be used to determine the empirical PCV correction of the LEO antenna. Using the residual approach, the PCV map of JASON-1 satellite is created. The mean post-fit phase residual is decreased from 8 to 5 mm, and 1 cm radial accuracy is demonstrated [40]. PCV correction of the GOCE satellite is determined with 154 days of data and the consistency of reduced dynamic and kinematic orbits is improved when applying the PCV map [37]. An error in the given phase center offset (PCO) of Sentinel-1A antenna has been found by comparing different PCVs and different POD approaches [14]. The influence of relative PCV on precise baseline determination (PBD) of formation flying satellites is also studied [41]. With application of the generated relative PCVs of GRACE and GRACE-FO satellites, the consistency of KBR measurements is improved by 30%. The effects of antenna PCV on orbit determination and clock estimation for CentiSpace-1 satellite, which served as a navigation satellite, were analyzed [42].

The twin GRACE-FO satellites (named GRCC and GRCD hereafter), designed as a successor of GRACE, and jointly developed by the National Aeronautics and Space Administration (NASA) and the German Research Centre for Geosciences (GFZ), were launched from California, USA, on May 22, 2018 [43]. The main goals of GRACE-FO are to continuously provide high resolution monthly solutions of the Earth's gravity field, surface mass change, and to measure the vertical temperature and humidity profiles of the Earth's atmosphere [44]. GRACE and GRACE-FO satellites adopt the same satellite appearance design, and are equipped with the KBR system, GPS receiver, SLR retroreflector, star camera assembly (SCA), and other scientific instruments [44,45]. The appearance and instrument installation position of GRACE-FO satellites are illustrated in Figure 1 [46].

**Figure 1.** Appearance and instrument installation position of GRACE-FO satellites.

Effects of PCV correction on POD of GRACE-FO satellites have been explored [35,43]. However, impacts of different ambiguity resolution strategies and in-flight PCV calibration on POD and PBD still need more investigation. In this analysis, antenna PCV models of GRACE-FO satellites are developed to further exploit the POD accuracy. Single receiver and double difference integer ambiguity resolution methods are investigated and realized. With one year of data, GRACE-FO orbits based on different AR strategies, including FA, SD AR, and DD AR, are studied and evaluated with JPL's orbit, SLR, and KBR measurements.

Section 2 provides details of the strategy and data adopted for GRACE-FO POD. The antenna PCV estimation approach and related results are presented in Section 3. Section 4 introduces the GNSS observation model and focuses on the SD and DD AR. POD and PBD results with PCV corrections and different AR methods were analyzed. The necessity of in-flight antenna calibration and ambiguity resolution is verified. Finally, discussions are made, followed by conclusions.

#### **2. POD Strategy and Data Usage**

#### *2.1. POD Strategy*

The reduced dynamic approach is employed to determine the orbits of GRACE-FO satellites. The Position And Navigation Data Analyst (PANDA) software [47], developed at the GNSS research center of Wuhan University, is modified and used for the POD process. The POD strategy is listed in Table 1. The macro model of GRACE-FO satellites [44] is applied to model the non-gravitational forces, which mainly result from atmospheric drag, solar, and earth radiation pressure. The atmospheric density values required for atmospheric drag modeling are obtained with the DTM94 model [48]. To account for deficiencies in DTM94 and the macro model, drag coefficient is estimated freely once per orbital revolution. Solar radiation pressure is calculated based on the satellite macro model and a scale factor is estimated per orbit determination arc. Additionally, one cycle per revolution (CPR) empirical accelerations in along-track and cross-track directions are estimated to compensate for deficiencies in the adopted force models. Other estimated parameters include receiver clock offsets and carrier phase ambiguities. The GRACE-FO antenna PCOs are applied according to the VGN1B product [44]. PCV corrections for IF carrier phase observations are calculated using the method described in Section 3.1.


**Table 1.** POD strategy of GRACE-FO satellites.


**Table 1.** *Cont.*

According to the ambiguity resolution methods, three kinds of orbits are calculated. The first orbit is calculated with float ambiguities (donated as FA solution hereafter), the second solution is constrained by integer DD ambiguities (DD AR solution), and the third orbit is generated with the constraint of integer SD ambiguities (SD AR solution).

#### *2.2. Data Usage*

The twin GRACE-FO satellites fly a polar orbit with an altitude of 490 km. Over the mission lifetime, the two satellites will remain in co-planar orbits and the along-track separation is about 220 km. For POD, the GRACE-FO satellites are equipped with the new generation of the GNSS space science receiver, the TriG receiver. The receiver upgrades the capabilities offered by the BlackJack receiver, which was carried on the GRACE mission. GPS C/A, P1, and P2 pseudoranges, and associated carrier phase observations, namely LA, L1, and L2 can be provided by the receiver [54]. The GPS P1 and P2 pseudoranges and L1 and L2 carrier phase observations are used in this research.

One year of data in 2019 is processed. The data can be obtained from the GRACE-FO level 1B RL04 products, which are available at ftp://isdcftp.gfz-potsdam.de/grace-fo/ (accessed on 17 October 2021). The products also include the satellite attitude (SCA1B), precise orbits provided by JPL (GNV1B), biased inter-satellite ranges measured by KBR system (KBR1B), and the positions of the GPS antenna phase center (VGN1B). The GNV1B and KBR1B data can be used to assess the orbit quality. PCOs of ionosphere free carrier phase observations and SLR reflector positions in the science reference frame (SRF) are listed in Table 2. CODE's GPS final orbits and 5 s clock products are used, and the associated OSB products are also applied to allow for SD AR [55,56].


**Table 2.** GPS antenna PCOs and SLR reflector Coordinates of GRACE-FO satellites.

#### **3. Estimation of PCV Corrections**

#### *3.1. Mathematical Models*

Modeling GNSS observation requires the computation of the geometric distance between the antenna phase center location of the GNSS satellite at signal emission time and the antenna phase center location of the receiving antenna at signal reception time. The phase center location usually differs from the mechanical antenna reference point (ARP). The difference vector is conventionally described by a set of phase center corrections. Such

a set of corrections consists of a PCO vector, which defines the position of the mean antenna phase center, with respect to the ARP and a consistent function, which models the azimuth and zenith dependent PCVs [36].

Phase center corrections have some inherent degrees of freedom. One set of corrections consisting of a PCO vector *r*<sup>0</sup> and an azimuth and zenith dependent function *φ*(*α*, *z*) can be transformed into a new set, consisting of *r* <sup>0</sup> and *φ* (*α*, *z*), which gives the same result:

$$\begin{cases} r\_0' = r\_0 + \Delta r \\ \phi'(\mathbf{a}, z) = \phi(\mathbf{a}, z) - \Delta r \cdot \mathbf{e} + \Delta \phi \end{cases} \tag{1}$$

where *α* and *z* are azimuth and zenith angles, Δ*φ* is an arbitrary offset and cannot be separated from the receiver clock. The unit vector *e* denotes the direction from the receiver to the satellite. The offset vector Δ*r* can be chosen arbitrarily. Preferably, PCVs should not induce a PCO and Δ*r* should be zero. In that case, the mean antenna phase center is explicitly defined by the PCO. This convention is particularly important if one would only apply PCO and no PCV [14].

Taking computational burden into consideration, residual approach is employed in this research. The antenna PCV is represented as piecewise linear function with respect to the zenith and azimuth angles in the antenna fixed coordinate system. The model assumes that PCV is composed of different grids, and zenith and azimuth angles are equally divided. When the observation is at the point of P within the grid ABCD, its PCV value can be linearly interpolated:

$$\begin{array}{l} \Delta\_{\text{PCV},P} = (1-\gamma)(1-\beta)\Delta\_{\text{PCV},A} + \gamma(1-\beta)\Delta\_{\text{PCV},B} + \gamma\beta\Delta\_{\text{PCV},C} + (1-\gamma)\beta\Delta\_{\text{PCV},D} \\ \gamma = (\mathfrak{a}-\mathfrak{a}\_{1})/(\mathfrak{a}\_{2}-\mathfrak{a}\_{1}), \beta = (z-z\_{1})/(z\_{2}-z\_{1}) \end{array} \tag{2}$$

where Δ*PCV*,*<sup>P</sup>* is the PCV at point P and is observation, Δ*PCV*,*i*(*i* = *A*, *B*, *C*, *D*) are the PCVs at points A, B, C, D are parameters to be estimated. *γ* and *β* are combination coefficients, *α* and *z* are azimuth and zenith angles of P, *α*<sup>1</sup> is the azimuth of point A and D, *α*<sup>2</sup> is the azimuth of B and C, *z*<sup>1</sup> is the zenith of A and B, *z*<sup>2</sup> is the zenith of C and D.

#### *3.2. Results*

The PCOs of GRACE-FO satellites provided in the VGN1B product are introduced as fixed and only the PCVs are estimated. Residual approach [36] is employed with three iterations to obtain the PCV maps. Resolutions in azimuth and zenith are 5◦. Two kinds of PCV maps are created—one is generated with the phase residuals from DD AR solution and the other is derived from the SD AR solution.

Figure 2 are PCV maps of GRACE-FO satellites produced from the DD AR solution. The azimuth of 0◦ points into the direction of flight. A similarity exists between the PCVs of GRCC and GRCD satellites, especially in the azimuth ranges of 90◦–180◦ and 270◦–360◦. The absolute maximum PCV of GRCC satellite is −14.56 mm with azimuth angle 185◦ and elevation angle 25◦. For the GRCD satellite, the absolute maximum PCV is −14.95 mm with azimuth angle 190◦ and elevation angle 30◦.

Figure 3 presents the PCV maps produced from the SD AR solution. The absolute maximum PCV of GRCC satellite is −12.57 mm with azimuth angle 185◦ and elevation angle 25◦, while the GRCD satellite is −12.15 mm with azimuth angle 185◦ and elevation angle 30◦. The PCV maps generated from DD AR and SD AR solutions have similar patterns when comparing Figure 2a with Figure 3a and comparing Figure 2b with Figure 3b. The PCVs made from the SD AR solution have less maximum values.

**Figure 2.** The 5◦ × 5◦ PCV corrections of (**a**) GRCC and (**b**) GRCD satellites based on ionosphere free carrier phase residuals from the DD AR solution.

**Figure 3.** The 5◦ × 5◦ PCV corrections of (**a**) GRCC and (**b**) GRCD satellites based on ionosphere free carrier phase residuals from SD AR solution.

PCV maps in Figures 2 and 3 are used to correct the carrier phase observations and then to perform the POD. All available carrier phase residuals are averaged and Figure 4 illustrates the results of GRCC satellite. Figure 4a presents the carrier phase residuals from the DD AR solution while Figure 4b presents the residuals from the SD AR solution. It should be noted that the limit of the color bar is 5 mm. With application of the PCV corrections, most of the carrier phase residuals of the DD AR solution are less than 1 mm and the mean value is 0.13 mm. The carrier phase residuals of the SD AR solution are a bit larger, especially in the range with azimuth of 0◦–30◦ and elevation of 30◦–60◦, the average is 0.41 mm. Similar results are obtained for the GRCD satellite and are not shown here. Carrier phase residuals can be used to measure the consistency between the applied models and GPS observations. It can be induced that the application of PCV corrections improves the orbit accuracy. The effects of PCV correction on POD and PBD will be presented in the following sections.

**Figure 4.** Carrier phase residuals of GRCC satellite derived from (**a**) DD AR solution and (**b**) SD AR solution with application of PCV corrections.

#### **4. POD with Ambiguity Resolution**

The GNSS observation model and the methods employed to resolve the SD and DD ambiguities are introduced in this part. FA, DD AR, and SD AR schemes are applied to calculate the GRACE-FO orbits. The generated orbits are assessed with JPL's orbit, SLR, and KBR data. Effects of PCV corrections on POD and PBD are also analyzed.

#### *4.1. Mathematical Models*

#### 4.1.1. GNSS Observation Model

GNSS code and carrier phase observations between a satellite and a receiver are usually described by the following equations [35,57]:

$$\begin{array}{l}P\_{r,j}^{s} = \rho\_r^s + c\left(dt\_r - dt^s\right) + I\_{r,j}^s + b\_{r,j} - b\_j^s\\L\_{r,j}^s = \rho\_r^s + c\left(dt\_r - dt^s\right) - I\_{r,j}^s + \lambda\_j\left(N\_{r,j}^s + B\_{r,j} - B\_j^s\right) + \lambda\_j\omega\_r^s\end{array} \tag{3}$$

where *P<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* and *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* are code and carrier phase observations in meters, *r* and *s* represent the receiver and GNSS satellite, *j* is the signal frequency, *ρ<sup>s</sup> <sup>r</sup>* is the geometry distance between receiver and satellite, *c* is the speed of light, *dtr* and *dt<sup>s</sup>* are receiver and satellite clock errors, *I<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is ionospheric delay [58], *br*,*<sup>j</sup>* and *<sup>b</sup><sup>s</sup> <sup>j</sup>* are receiver and satellite hardware biases of pseudorange, *Br*,*<sup>j</sup>* and *B<sup>s</sup> <sup>j</sup>* are hardware biases of the carrier phase observation, *λ<sup>j</sup>* is signal wavelength and *N<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is the integer ambiguity. *<sup>ω</sup><sup>s</sup> <sup>r</sup>* is the phase wind-up error and can be corrected by model [53], in the remaining part, this item will be omitted.

Both ionosphere free and ionosphere float models can be used for POD. The first order ionospheric delay can be eliminated by the ionosphere free model, while in the ionosphere float model, the ionospheric delay is estimated as an unknown parameter. The ionosphere free approach eliminates the ionospheric delay twice and leads to smaller redundancy. The ionosphere free method owns fewer unknown parameters and it also has the drawback of lacking flexibility for further model strengthening [59]. In this investigation, the ionosphere free model is employed.

Following the IGS convention, the IF pseudorange biases will be assimilated into the receiver and satellite clock offsets. IF observations can be formulated as:

$$\begin{array}{l} P\_{r,IF}^{s} = \rho\_r^s + c(\overline{d}t\_{I} - \overline{d}t^s) \\ L\_{r,IF}^{s} = \rho\_r^s + c(\overline{d}t\_r - \overline{d}t^s) + \lambda\_1 \overline{N}\_{r,IF}^{s} \\ \overline{N}\_{r,IF}^{s} = N\_{r,IF}^{s} + (B\_{r,IF} - b\_{r,IF}/\lambda\_1) - (B\_{IF}^{s} - b\_{IF}^{s}/\lambda\_1) \end{array} \tag{4}$$

where *P<sup>s</sup> <sup>r</sup>*,*IF* and *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*IF* are IF code and carrier phase observations, *dtr* and *dt<sup>s</sup>* are receiver and satellite clock errors, including the pseudorange hardware biases. *br*,*IF* and *b<sup>s</sup> IF* are the IF pseudorange hardware biases while *Br*,*IF* and *B<sup>s</sup> IF* are the carrier phase biases. *λ*<sup>1</sup> is the wavelength of L1. *N<sup>s</sup> <sup>r</sup>*,*IF* is the IF ambiguity, which includes the hardware biases in both pseudorange and carrier phase observations, *N<sup>s</sup> <sup>r</sup>*,*IF* is a combination of integer L1 and L2 ambiguities.

In order to fix the ambiguities, Blewitt presents a sequential integer ambiguity resolution method to resolve the DD ambiguities of long baseline. The integer DD WL ambiguities are firstly determined with a combination of pseudorange and carrier phase observations. Then the DD narrow-lane (NL) ambiguities are fixed using the estimated IF ambiguities and already fixed integer WL ambiguities [57]. For a single receiver integer ambiguity resolution, GNSS satellites SD observations are used to cancel the receiver based phase bias. Guo et al. provide a detailed description about the model to realize single receiver AR [35]. These methods are employed in this contribution to fix the SD and DD ambiguities.

#### 4.1.2. Single Receiver Ambiguity Resolution

The float IF ambiguity can be rewritten as the combination of integer WL and float NL ambiguities as the following [57]:

$$\begin{array}{l}\overline{N}^{s}\_{r,IF} = \frac{f\_1 f\_2}{f\_1^2 - f\_2^2} N^s\_{r,WL} + \frac{f\_1}{f\_1 + f\_2} \overline{N}^s\_{r,NL} \\ N^s\_{r,WL} = N^s\_{r,1} - N^s\_{r,2} \\ \overline{N}^{s}\_{r,NL} = N^s\_{r,NL} + d\_{r,NL} - d\_{NL}^s \end{array} \tag{5}$$

where *N<sup>s</sup> <sup>r</sup>*,*WL* and *<sup>N</sup><sup>s</sup> <sup>r</sup>*,*NL* are integer WL and NL ambiguities, *<sup>N</sup><sup>s</sup> <sup>r</sup>*,*NL* is float NL ambiguity, *dr*,*NL* and *d<sup>s</sup> NL* are the NL FCBs of the receiver and satellite.

The integer WL and NL ambiguities are fixed with two steps. Firstly, the WL ambiguity is resolved using the Hatch–Melbourne–Wübbena combination [60–62]:

$$\overline{N\_{r,WL}^{s}} = (\frac{f\_1 L\_{r,1}^{s} - f\_2 L\_{r,2}^{s}}{f\_1 - f\_2} - \frac{f\_1 P\_{r,1}^{s} + f\_2 P\_{r,2}^{s}}{f\_1 + f\_2}) / \lambda\_{WL} = N\_{r,WL}^{s} + d\_{r,WL} - d\_{WL}^{s} \tag{6}$$

where *dr*,*WL* and *d<sup>s</sup> WL* are receiver and satellite WL FCBs, *λWL* is WL wavelength. After correction of the WL FCBs, the integer WL ambiguity can be determined. Then integer NL ambiguity can be derived with the fixed WL ambiguity, the estimated float IF ambiguity, and the associated NL FCBs using Equation (5).

Once the WL and NL ambiguities are resolved, the IF ambiguity can be reconstructed according to Equation (5). Usually satellite dependent WL and NL FCBs can be estimated with ground network stations or provided by GNSS analysis institutions [26,55,63]. However, receiver dependent FCBs are not available and cannot be removed from Equations (5) and (6). SD operation between satellites is suggested to cancel the receiver FCBs. The SD integer WL and NL ambiguities can be calculated using Equation (7):

$$
\begin{array}{l}
\Delta N\_{r,WL}^{s\_1,s\_2} = \Delta \overline{N}\_{r,WL}^{s\_1,s\_2} + \Delta d\_{WL}^{s\_1,s\_2} \\
\Delta N\_{r,NL}^{s\_1,s\_2} = \Delta \overline{N}\_{r,NL}^{s\_1,s\_2} + \Delta d\_{NL}^{s\_1,s\_2}
\end{array}
\tag{7}
$$

where <sup>Δ</sup> is single difference operator, <sup>Δ</sup>*Ns*1,*s*<sup>2</sup> *<sup>r</sup>*,*WL* and <sup>Δ</sup>*Ns*1,*s*<sup>2</sup> *<sup>r</sup>*,*NL* are float SD WL and NL ambiguities of satellite *<sup>s</sup>*<sup>1</sup> and *<sup>s</sup>*<sup>2</sup> while <sup>Δ</sup>*Ns*1,*s*<sup>2</sup> *<sup>r</sup>*,*WL* and <sup>Δ</sup>*Ns*1,*s*<sup>2</sup> *<sup>r</sup>*,*NL* are integer ambiguities. <sup>Δ</sup>*ds*1,*s*<sup>2</sup> *WL* and Δ*ds*1,*s*<sup>2</sup> *NL* are SD hardware biases between satellites. Once SD WL and NL ambiguities are fixed, the SD IF ambiguity can be recovered as follows:

$$
\Delta \overline{N}\_{r,IF}^{s\_1,s\_2} = \frac{f\_1 f\_2}{f\_1^2 - f\_2^2} \Delta N\_{r,WL}^{s\_1,s\_2} + \frac{f\_1}{f\_1 + f\_2} (\Delta N\_{r,NL}^{s\_1,s\_2} - \Delta d\_{NL}^{s\_1,s\_2}) \tag{8}
$$

The recovered SD IF ambiguity is taken as pseudo-observation to constrain the undifferenced IF ambiguity parameters in the POD estimation using Equation (9):

$$
\Delta \overline{N}\_{r,IF}^{s\_1,s\_2} = \overline{N}\_{r,IF}^{s\_1} - \overline{N}\_{r,IF\prime}^{s\_2} \quad \mathcal{W}\_r^{s\_1,s\_2} \tag{9}
$$

where *<sup>W</sup>s*1,*s*<sup>2</sup> *<sup>r</sup>* is the weight of the pseudo-observation and will be assigned with a large value to put a strong confidence.

#### 4.1.3. Double Difference Ambiguity Resolution

With regard to the DD ambiguity, both receiver and satellite FCBs can be removed. The WL and NL DD ambiguities are theoretically integers according to Equations (5) and (6). The DD IF ambiguity can be reconstructed by:

$$
\Delta\nabla\overline{N}^{\varepsilon\_{1}\varepsilon\_{2}}\_{r\_{1}r\_{2},IF} = \frac{f\_{1}f\_{2}}{f\_{1}^{2} - f\_{2}^{2}}\Delta\nabla N^{\varepsilon\_{1}\varepsilon\_{2}}\_{r\_{1}r\_{2},NL} + \frac{f\_{1}}{f\_{1} + f\_{2}}\Delta\nabla N^{\varepsilon\_{1}\varepsilon\_{2}}\_{r\_{1}r\_{2},NL} \tag{10}
$$

where Δ∇ is double difference operator, *r*<sup>1</sup> and *r*<sup>2</sup> denote a receiver pair. The DD IF ambiguity then is also treated as a pseudo-observation to constrain the undifferenced IF ambiguities.

$$
\Delta\nabla\overline{N}^{\mathbf{s}\_1,\mathbf{s}\_2}\_{r\_1,r\_2,IF} = \overline{N}^{\mathbf{s}\_2}\_{r\_2,IF} - \overline{N}^{\mathbf{s}\_1}\_{r\_2,IF} - \overline{N}^{\mathbf{s}\_2}\_{r\_1,IF} + \overline{N}^{\mathbf{s}\_1}\_{r\_1,IF\prime} \mathcal{W}^{\mathbf{s}\_1,\mathbf{s}\_2}\_{r\_1,r\_2} \tag{11}
$$

*<sup>W</sup>s*1,*s*<sup>2</sup> *<sup>r</sup>*1,*r*<sup>2</sup> is the weight of the pseudo-observation. The precision of the DD and SD IF ambiguities is set to be 1 × <sup>10</sup>−<sup>4</sup> L1 cycles (about 0.02 mm) to impose strong constraints.

#### 4.1.4. Integer Ambiguity Validation

Several integer ambiguity resolution methods have been proposed including the integer rounding, integer bootstrapping, integer least squares, and partial ambiguity resolution methods [64]. The integer rounding approach is used in this research to resolve the WL and NL ambiguities. The ambiguity validation procedure is presented here for clearness. The reader can refer to [65] for more details.

Let *b* be the float WL or NL ambiguities and it is assumed that the probability density function for *b* is:

$$P(b|n) = \frac{1}{\left(2\pi\right)^{1/2}\sigma} \exp[-\frac{\left(b-n\right)^2}{2\sigma^2}]\tag{12}$$

where *n* is the true integer ambiguity, *σ* is the formal uncertainty from the adjustment. Let *I* be the nearest integer of *b*. In the ambiguity rounding procedure, it should be determined whether to fix the ambiguity to *I* or leave it as float. The deviation

$$
\mathfrak{x} = \mathfrak{b} - I \tag{13}
$$

must fall in the interval [−0.5, 0.5]. If the true integer ambiguity *n* is not equal to the nearest integer *I*, the probability of wrong ambiguity resolution can be calculated using Equation (14):

$$Q\_0 = \sum\_{m=1}^{\infty} \left[ \operatorname{erfc}(\frac{m-\chi}{\sqrt{2}\sigma}) - \operatorname{erfc}(\frac{m+\chi}{\sqrt{2}\sigma}) \right] \le a \tag{14}$$

$$
\operatorname{erfc}(\mathbf{x}) = \frac{2}{\sqrt{\pi t}} \int\_x^\infty e^{-t^2} dt\tag{15}
$$

where *α* is the allowable rate of wrongly fixed ambiguity, *er f c* is the complementary error function. In order to eliminate some extremes, a decision function is defined:

$$d(x, \sigma) = T / Q\_0 \tag{16}$$

$$T = \begin{cases} 0, & |\mathbf{x}| \ge T\_1 or \sigma \ge T\_2\\ (1 - \frac{|\mathbf{x}|}{T\_1})(3T\_2 - 3\sigma), & otherwise \end{cases} \tag{17}$$

where *T*<sup>1</sup> is the limit for the deviation *x* and *T*<sup>2</sup> is the limit for *σ*. In this analysis, *T*<sup>1</sup> and *T*<sup>2</sup> are both 0.25 cycles for WL ambiguity resolution, and are 0.15 cycles for NL ambiguity resolution. When *d*(*x*, *σ*) is greater than 1000, the ambiguity will be fixed to the nearest integer [28,66].

#### *4.2. Results*

One year of data in 2019 collected by GRACE-FO satellites were analyzed; 38 days were excluded due to large satellite maneuvers or data gaps. The method introduced in Section 4.1.2 is used to resolve the SD WL and NL ambiguities, and the method in Section 4.1.3 is adopted to fix the DD ambiguities. Firstly, the residuals and fixing rates of WL and NL ambiguities are presented. Then the absolute orbit accuracy of each satellite and baseline accuracy between two satellites are evaluated.

#### 4.2.1. Ambiguity Resolution Results

Ambiguity resolution performance is evaluated in terms of the residual and fixing rate. Ambiguity residual is the difference between the float ambiguity and its nearest integer. After removing the satellite and receiver FCBs, the float SD WL and NL ambiguities should be close to the integers. For DD ambiguities, the FCBs of the receiver and satellite are cancelled, which leads to an integer nature. As illustrated by Teunissen [67], the distribution of the ambiguity residuals is non-Gaussian. It is symmetric and the point of symmetry is the origin, which implies that the mean of the ambiguity residuals is zero. The distribution will become more peaked when the estimated float ambiguity is more precise. Figure 5a,b are residual distributions of SD WL and NL ambiguities while Figure 5c,d are distributions of DD ambiguities. Blue bars in the figures are percentages of the related residuals. The mean values of SD and DD WL and NL residuals are all less than 0.01 cycles and verify the symmetry of distribution. The standard deviations (STDs) of SD WL and NL residuals are 0.10 cycles and 0.07 cycles and validate the high quality of the estimated float ambiguities. For DD WL and NL ambiguities, the STDs are both 0.06 cycles. The smaller STDs of DD ambiguity residuals confirm the effectiveness of the DD operation to cancel the receiver and satellite FCBs.

**Figure 5.** Residual distributions of (**a**) SD WL, (**b**) SD NL, (**c**) DD WL, and (**d**) DD NL ambiguities of GRACE-FO satellites.

The ambiguity fixing rate is the number of fixed ambiguities divided by the total number of ambiguities. Figure 6 illustrates the daily ambiguity fixing rates of SD and DD ambiguities. Figure 6a is for SD ambiguities while Figure 6b is for DD ambiguities. The upper and lower panels of Figure 6a are results of GRCC and GRCD satellites. Blue and red dots in Figure 6 are fixing rates of WL and NL ambiguities. For the GRCC satellite, one year average fixing rates of SD WL and NL ambiguities are 98.5% and 95.8%. A similar performance is observed for the GRCD satellite; the fixing rates are 98.4% and 96.5%. The excellent fixing rates confirm the good quality of applied phase bias products. For the DD case, higher fixing rates are obtained for WL and NL ambiguities, which are 99.1% and 96.7%. Gaps in day of year (DOY) 32–51, 146–149, 205–211, and 299–302 are because of data missing.

#### 4.2.2. Single Satellite Orbit Validation

JPL's precise orbit of GRACE-FO satellites is also included in the level 1B product. The orbit is calculated with single receiver ambiguity resolution [32,44] and can be used to assess the orbits generated in this investigation. The effect of PCV correction on POD is first evaluated. Figure 7 presents the along-track orbit differences between SD AR solution and JPL's orbit. The upper panel is for GRCC satellite and lower part is for GRCD satellite. The red dots are daily RMSs of differences of the orbit determined with PCV corrections while blue dots are RMSs without PCV corrections. Using created PCV maps, the RMS of the GRCC satellite decreases from 7.2 to 6.0 mm, and from 6.4 to 5.9 mm for the GRCD satellite. The detailed results for three components of SD AR and DD AR orbits can be found in Table 3. The average improvements in along-track, cross-track, and radial direction are 12.2%, 1.6%, and 1.0% for SD AR solution, and are 10.3%, 1.5%, and 3.5% for DD AR solution. The improvement of along-track orbit accuracy is the most significant, which validates the necessity of PCV calibration, especially for PBD.

**Figure 7.** Daily RMSs of along-track orbit differences between JPL's orbit and SD AR solution with (red dots) and without (blue dots) PCV corrections. The upper panel is for the GRCC satellite and lower panel is for the GRCD satellite.


**Table 3.** Orbit improvements with application of PCV maps for SD AR and DD AR solutions.

RMSs on DOY 100, 141, and 159 of GRCC satellite, and RMSs on DOY 204 and 290 of GRCD satellite are larger than other days. It has been found that sometimes, in these days, fewer GPS satellites are available, which leads to a decrease in orbit accuracy.

Orbit accuracy with different AR methods is also evaluated. Figure 8 presents the RMSs of along-track orbit differences of FA, DD AR, and SD AR solutions with JPL's orbit. For GRCC satellite, the along-track RMS of the FA solution is 8.8 mm, the RMS is improved to 8.3 mm with DD AR and further achieves 6.0 mm by using the SD AR method. Results of three orbit components of the two satellites are listed in Table 4. Compared with the FA solution, the DD AR orbit has an average improvement of 6.3% while the SD AR result has an improvement of 18.3%. Similar to PCV correction, the along-track orbit accuracy has the most significant improvement, which is more than 30%.

**Figure 8.** Daily RMSs of along-track orbit differences of FA, DD AR, and SD AR solutions with JPL's orbit. The upper panel is for the GRCC satellite and the lower panel is for the GRCD satellite.


**Table 4.** RMSs of differences between the FA, DD AR and SD AR orbits and JPL's orbit.

GRACE-FO satellites are equipped with SLR retroreflectors. The independent SLR data can be used as external validation. SLR normal points from nine laser stations (Yarragadee, Greenbelt, Haleakala, Zimmerwald, Mount Stromlo, Wettzell, Graz, Herstmonceux, and Potsdam) are used for orbit evaluation. Residuals larger than 10 cm are rejected and an elevation cutoff of 10◦ is applied. Figure 9 shows the RMSs of SLR residuals. Red, green, and blue dots are RMSs of FA, DD AR, and SD AR orbits. The RMSs of the GRCC satellite are 11.5, 10.2, and 9.6 mm, which confirms the positive impact of AR on the orbit quality. Compared with the DD AR solution, the SD AR result has better performance, which is also validated by JPL's orbit.

**Figure 9.** RMSs of SLR residuals derived from FA (red dots), DD AR (green dots), and SD AR (blue dots) solutions. The upper panel is for the GRCC satellite and the lower panel is for the GRCD satellite.

#### 4.2.3. Baseline Validation

The GRACE-FO mission offers the possibility of validating the computed orbit, in particular, the along-track component, with the ultra-precise KBR measurements [33,68]. Figure 10 illustrates the STDs of KBR residuals (named KBR STD hereafter) calculated using FA, DD AR, and SD AR orbits. The effect of antenna PCV correction on PBD is also illustrated. Figure 10a presents the KBR STDs of the SD AR orbit with PCV corrections (blue dots), the SD AR orbit without PCV corrections (green dots), and FA orbit with PCV corrections (red dots). PCV corrections used to calculate FA orbit are generated from the SD AR solution. Figure 10b is the result of the DD AR solution.

**Figure 10.** Daily STDs of KBR residuals of (**a**) FA orbit with PCV corrections (red dots), SD AR orbit without PCV corrections (green dots), and SD AR orbit with PCV corrections (blue dots); (**b**) FA orbit with PCV corrections (red dots), DD AR orbit without PCV corrections (green dots), and DD AR orbit with PCV corrections (blue dots).

Both SD and the DD AR lead to a significant reduction in KBR residuals. Compared with the FA result, the SD AR orbit improves the KBR STD from 5.9 to 3.1 mm, and the DD AR orbit improves the STD to 1.6 mm. With PCV corrections, the KBR STD of the SD AR solution decreases from 3.1 to 1.8 mm, and the STD of the DD AR orbit decreases from 1.6 to 0.9 mm. The baseline precision of the SD AR solution gains similar performance, as presented by Arnold et al. [33]. The DD AR orbit achieves a sub-millimeter level as reported in previous research [6,19,69]. Comparisons with JPL's orbit and SLR data show that the SD AR solution possesses better performance, while in the KBR validation, the residual STD of the DD AR orbit is smaller. In other words, the SD AR orbit has better absolute orbit accuracy, while the DD AR orbit gets better relative performance. This is explained by two facts: (a) phase bias products for SD AR are not free of errors and these errors may degrade the KBR validation [35]; (b) both JPL's orbit and the SD AR orbit in this investigation adopt the single receiver ambiguity resolution method and a better consistency can be inferred.

#### **5. Discussions**

Precise knowledge of the phase center location of the GNSS antenna is a prerequisite for high precision LEO orbit determination. PCO and PCV values of the LEO antenna obtained from ground calibration cannot reflect the influence of error sources, which are additionally encountered in the actual spacecraft environment. An in-flight calibration of the LEO antenna is thus mandatory. In this analysis, antenna PCV is represented as a piecewise linear function with respect to zenith and azimuth angles in the antenna fixed coordinate system. A residual approach is employed to estimate the PCV correction. Different ambiguity resolution strategies, including SD AR and DD AR, are investigated to fully exploit the precision of GPS observations for POD. The single receiver integer ambiguity resolution concept employed here makes use of the carrier phase biases and clock products provided by CODE. One year of GRACE-FO data in 2019 are used to demonstrate the effect of integer ambiguity resolution and PCV correction on POD.

PCV maps derived from SD AR and DD AR solutions have very similar patterns. With PCV corrections, carrier phase residuals are reduced from about 10 to 1–2 mm, which implies better consistency between the observations and applied models. The average improvements of the SD AR solution in the along-track, cross-track, and radial directions are 12.2%, 1.6%, and 1.0%, and for the DD AR solution, the improvements are 10.3%, 1.5%, and 3.5%. The great enhancement in the along-track orbit accuracy validates the necessity to calibrate the PCV errors in the relative POD.

The mean values of both SD and DD ambiguity residuals are less than 0.01 cycles and confirm the symmetry of distribution. The fixing rates of SD WL and NL ambiguities are more than 98% and 95%. Compared with the FA solution, DD AR orbit accuracy has an average improvement of 6.3% while the SD AR result gains an increase of 18.3%. For high-grade SLR stations, range residuals with RMS less than 10 mm are achieved for SD AR orbit, which marks a 17% improvement compared to the FA result. Independent KBR measurements are also used as external validations. KBR residuals STD of the SD AR solution is 1.8 mm while STD of the DD AR orbit is 0.9 mm and reaches the sub-millimeter level. The better baseline accuracy for the DD AR solution is explained by the errors in phase bias products.

#### **6. Conclusions**

Integer ambiguity resolution plays a crucial role in achieving the best positioning or orbit accuracy. Single receiver and double difference integer ambiguity resolution models are introduced in this contribution. One year of GRACE-FO data were analyzed to verify the orbit improvement with different ambiguity resolution strategies. In-flight antenna PCV maps were developed to further exploit the POD accuracy.

With PCV corrections, along-track orbit accuracy of the SD AR solution is improved by 12.6%, and that of the DD AR solution is improved by 10.3%, which verifies the necessity of in-flight antenna calibration. With phase biases cancelled by the double difference operation, fixing rates of DD ambiguities are more than 96.7%. Similar performance is achieved by SD AR, which confirms the consistency of theoretical models and bias/clock products.

JPL's orbited together with independent SLR and KBR measurements are used to assess the results of different AR strategies. Both SD AR and DD AR solutions improve the orbit accuracy. SD AR solution provides the best performance of absolute orbit. Compared with JPL's orbit, the RMS of the SD AR orbit differences is better than 6 mm and the SLR residual RMS is less than 10 mm. The DD AR solution realizes the highest baseline accuracy and the STD of KBR residuals achieves 0.9 mm. In addition, sometimes there are fewer GPS satellites available, which leads to the reduction of orbit accuracy. With the application of multi-mode GNSS receivers and multi-GNSS phase bias products [70,71], orbit accuracy based on single receiver ambiguity resolution can be further improved.

**Author Contributions:** Conceptualization, B.J. and Y.L.; methodology, B.J.; software, B.J. and K.J.; validation, B.J. and S.C.; writing—original draft preparation, B.J.; writing—review and editing, S.C.; supervision, Y.L. and Z.L.; funding acquisition, Y.L. 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, grant number 12033009.

**Data Availability Statement:** GRACE-FO data in this paper can be freely accessed at ftp://isdcftp. gfz-potsdam.de/grace-fo/ (accessed on 17 October 2021). The GPS final orbit, clock, and OSB products can be found at ftp://ftp.aiub.unibe.ch/CODE/ (accessed on 17 October 2021).

**Acknowledgments:** The authors thank GFZ for providing GRACE-FO GPS data. The authors are grateful to the CODE for providing GPS products.

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

#### **References**

