## Summary

This paper presents methods for the precise orbit determination (POD) of a satellite in the CYGNSS constellation based on available single-frequency GPS code and carrier measurements. The contributions include the development and evaluation of procedures for single-frequency POD with GipsyX, improvement of CYGNSS orbit knowledge, and an assessment of its final accuracy. Ionospheric effects are mitigated using the GRAPHIC processing method, and spacecraft multipath effects are calibrated with an azimuth/elevation-dependent antenna calibration map. The method is demonstrated using comparable data from the GRACE mission, from which we infer the expected accuracy of the CYGNSS results. Processing more than 170 days of data from each mission, a 1*σ* CYGNSS orbit accuracy of 2.8 cm radial, 2.4 cm cross-track, and 6 cm in-track is demonstrated. We expect that achieving this level of performance will expand the set of future scientific investigations that can be undertaken using satellites equipped with single-frequency GNSS.

## 1 INTRODUCTION

Global navigation satellite system (GNSS) tracking for spacecraft in low Earth orbit (LEO) has been in use for more than 30 years. The TOPEX/POSEIDON ocean surface altimetry mission, launched in 1992, achieved orbit accuracies better than 3 cm (radial RMS; Bertiger et al., 1994) when estimated with a reduced dynamic technique (Wu et al., 1991) using GNSS tracking alone. Since then, follow-on missions have achieved radial accuracies of 1-cm accuracy for Jason-1 (Haines et al., 2004) and sub-centimeter for Jason-2 (Bertiger et al., 2010a). Precise orbit determination (POD) accuracy improvement in the most recent missions is due, in part, to improved GNSS orbit and clock products and calibrated transmitter antenna patterns with improved phase modeling (Schmid et al., 2005, 2007).

POD accuracy is a critical component of satellite geodetic techniques such as ocean surface altimetry, gravity recovery, and environmental sensing with GNSS radio occultation. These missions are typically equipped with scientific-grade (dual-frequency) GNSS receivers that allow for the direct removal of the first-order ionosphere delay. Even some small satellites like the 3U CubeSats flown by Spire (Nguyen et al., 2020) have this capability, which is used to produce a variety of science data products. However, there are also missions where precise orbits are not prioritized, and satellites are equipped with only single-frequency GNSS capability to support basic navigation and geolocation functions. Space-capable, single-frequency receivers with accuracies of 10 m are readily available (Montenbruck et al., 2012) and function well, meeting the lower cost, weight, power, and telemetry constraints for small spacecraft.

CYGNSS is a constellation of eight small satellites, each equipped with a GNSS reflectometry (GNSS-R) instrument consisting of a receiver, a zenith antenna for traditional direct signal tracking and formation of navigation solutions, and two downward-looking antennas capable of measuring GPS reflections from the surface of the Earth (Ruf et al., 2012). While CYGNSS was designed to use GNSS-R to measure ocean surface wind speed, it also provides an opportunity to study the feasibility of using GNSS-R for ocean altimetry. For ocean altimetry application, the delay of the reflected GPS signal must be precisely and accurately measured. Previous work (Li et al., 2018, 2019; Mashburn et al., 2020) identified uncertainty in the knowledge of CYGNSS satellite orbits as a significant source of error for GNSS-R altimetry. Consequently, a focused effort to analyze and improve the CYGNSS orbit estimates was undertaken.

LEO satellite missions equipped with dual-frequency receivers can achieve 3D orbit accuracies of 3–5 cm (Hackel et al., 2017). This level of accuracy is achieved by processing dual-frequency pseudorange and carrier-phase observations in combination with state-of-the-art gravity fields, satellite surface macromodels, antenna calibration maps, and reduced dynamic processing. Antenna calibration maps can be generated either by residual stacking or direct estimation within a POD filter. Residual stacking generates the map by averaging post-fit residuals in azimuth and elevation bins based on the transmitter line of sight (Haines et al., 2004). Jäggi et al. (2009) showed that the direct approach is superior to residual stacking, due to the simultaneous estimation of the map along with the receiver clock and phase ambiguities. The implementation of single-receiver integer ambiguity resolution has also been shown to improve orbit estimates (Bertiger et al., 2010b; Montenbruck et al., 2018) with recent results by Mao et al. (2021) demonstrating 1-cm precision.

Many of these techniques have also been applied to POD with single-frequency observations; however, it requires an alternative approach for the mitigation of ionospheric effects and cannot take advantage of carrier-phase integer ambiguity resolution. The latter limitation reduces the effectiveness of purely kinematic solutions, with prior studies (Montenbruck et al., 2005; Wu et al., 1991) showing reduced dynamic filtering to be the most effective for single-frequency POD. To address the ionospheric effects, three basic approaches are well-established—modeling based on ionospheric density databases, estimation based on the difference in code and carrier observations, and direct removal of ionospheric effects by combining code and carrier observations.

Model-based methods approximate ionospheric effects using measurement-based global ionospheric maps (GIMs; Hernández-Pajares et al., 2009) or empirical models such as the International Reference Ionosphere (IRI; Bilitza et al., 2017; Bock et al., 2009; Montenbruck & Gill, 2002). The Group and Phase Ionosphere Correction (GRAPHIC) method (Montenbruck, 2003; Yunck, 1993) eliminates first-order ionospheric effects directly by adding code-based pseudorange and carrier-based measurements to produce an ionosphere-free observable with a fixed ambiguity per pass and error dominated by half the noise in the pseudorange measurement. This approach is effective and widely used.

Early single-frequency POD of the TOPEX/POSEIDON spacecraft achieved radial accuracies of 4 to 5 cm (Muellerschoen et al., 1994) with GRAPHIC. Subsequent studies (Bock et al., 2009; Montenbruck et al., 2005) achieved full 3D accuracies below 10 cm. Because GRAPHIC incorporates a code-base pseudorange as part of the primary observable, the noise and group delay in the code has a more substantial influence here than in an integer ambiguity resolved dual-frequency POD solution. Shao et al. (2019) demonstrated further improvements by reducing systematic errors with an antenna calibration map for the GRAPHIC observable generated with a residual stacking approach. They achieved 3D accuracies better than 5 cm for the GRACE-B satellite.

The goal of our work is to improve the CYGNSS orbits and quantify the precision and accuracy of the results in support of GNSS-R altimetry research. As part of this effort, we introduce a method for correcting a clock-like signal in the CYGNSS raw carrier measurements that is not present in the pseudorange. This correction allows for the implementation of the GRAPHIC orbit solutions. We estimate code-only and GRAPHIC observable antenna calibrations directly within a dynamic POD filter, and demonstrate performance of single frequency, model-corrected, code-only, and GRAPHIC orbit solutions in the GipsyX software environment.

The GipsyX software suite was developed by NASA’s Jet Propulsion Laboratory (JPL) for positioning, navigation, and timing using space geodetic data (Bertiger et al., 2020). While GipsyX is primarily used with dual-frequency GNSS measurements, we show that it can also be configured to produce accurate single-frequency POD. To demonstrate the methods for CYGNSS, we have selected GRACE-B as the primary reference satellite. The GRACE mission, at a similar altitude to CYGNSS, provides extensive sets of high-quality dual-frequency observations that are readily available, as well as previously established high-precision orbits. The GRACE satellite orbit and GPS measurements are used to directly evaluate errors associated with L1 single-frequency POD strategies that we apply to CYGNSS. In addition to GRACE, we also consider the higher-altitude Sentinel-6 Michael Freilich (MF) spacecraft to provide an additional dual-frequency reference. Like GRACE, the Sentinel-6 MF altimetry mission is highly dependent on accurate orbits and provides a reference for the L1 single-frequency POD errors. While our results for this mission are less extensive, it provides a useful reference point because it flies at a 66-degree inclination, more similar to CYGNSS than the GRACE polar orbit.

The remainder of this paper is organized into three main sections. Section 2 describes the background for LEO POD within GipsyX, including the specifics of the GRACE and CYGNSS missions, estimation strategies, antenna calibration, and issues encountered with the CYGNSS phase measurements. Section 3 presents the results, including estimated antenna calibration maps, followed by an orbit precision and accuracy assessment using internal metrics of residual fit root-mean-squre (RMS) and orbit overlaps for GRACE and CYGNSS, and external metrics for GRACE by comparing the dual-frequency solutions. We also use a small set of Sentinel-6 MF orbits to support the assessment of CYGNSS GRAPHIC solutions. Finally, the conclusions in Section 4 summarize the findings of this work.

## 2 METHODS

### 2.1 CYGNSS, GRACE, and Sentinel-6 MF Comparison

Table 1 shows a comparison of the key features of the GRACE-B, CYGNSS FM05, and Sentinel-6 MF satellites, orbits, and onboard GPS receivers. Both the CYGNSS eight-satellite constellation and GRACE A and B pair of satellites operate in the LEO regime in near circular orbits, with CYGNSS about 50-km higher in altitude. They have similar attitude profiles with a nominal yaw-fixed orientation relative to the flight direction. The biggest difference between the GRACE and CYGNSS orbits is that GRACE is in a polar orbit, while CYGNSS is in a lower (35-deg) inclination orbit driven by its mission to observe ocean surface winds. Sentinel-6 MF is much higher than both GRACE and CYGNSS but, like CYGNSS, is not in a polar orbit (inclination of 66-degrees).

Each CYGNSS satellite is equipped with a Surrey Satellite Technology Ltd. (SSTL) Space GPS Receiver Remote Sensing Instrument (SGR-ReSI) receiver that determines the satellite position in real time using civil single-frequency GPS (L1 C/A) signals (Unwin et al., 2013). Since the primary mission is not reliant on precise knowledge of the CYGNSS orbit, only this onboard GPS instantaneous navigation (Nav) solution was originally downlinked. The Nav solution is noisy with 1*σ* 3D position errors of approximately 3 m. To enable an evaluation of altimetry applications, in July 2019, the satellite flight software was updated to also downlink the individual GPS satellite pseudorange and carrier-phase measurements at 15-second intervals. At the time of writing, these lower-level GPS data are only available in the CYGNSS Level-0 data product, accessed by special request for this investigation. The availability of these raw tracking measurements provides the opportunity to estimate the CYGNSS orbits on the ground using advanced postprocessing techniques in GipsyX or other POD software.

The Blackjack GPS receiver flown on GRACE was designed specifically for geodetic-quality POD. It uses advanced codeless techniques to track L1 C/A and P(Y) pseudorange and carrier as well as L2 P(Y) pseudorange and carrier (Dunn et al., 2002; Montenbruck & Kroes, 2003). The receiver is coupled with a flush-mounted choke-ring antenna to reduce multipath reflections. The overall measurement approach promotes more precise measurements than typical receiver implementations geared to support real-time onboard navigation and timing. The Sentinel-6 MF spacecraft is equipped with a TriG receiver, also capable of tracking L1 C/A, P(Y), and L2 P(Y) measurements (Donlon et al., 2021). With both CYGNSS, GRACE, and Sentinel-6 MF tracking the L1 C/A code and carrier, we can apply the same data processing and orbit estimation methods to the three missions. The observed C/A code noise for both GRACE and Sentinel-6 MF is below 30 cm based on post-fit residual RMS. For CYGNSS, it is significantly larger with a value around 1 m. Finally, to ensure comparable ionospheric effects for both GRACE and CYGNSS, we selected data sets for each mission near the solar minimum, using data from early 2009 for GRACE and the last half of 2019 for CYGNSS.

### 2.2 GRAPHIC Observable

The GRAPHIC method is well established as an effective approach for creating an ionosphere-free observable from single-frequency GNSS measurements (Bock et al., 2009). The GRAPHIC observable is formed from the pseudorange and carrier by computing the mean of the two measurements as follows:

1

2

3

This results in an ionosphere-free phase-like observable that contains half the code noise, half the phase noise, and a half-wavelength integer ambiguity (*N*/2). While the GRAPHIC observable only contains half the code noise, this is substantially larger than what would be observed for typical carrier measurements. Bock et al. (2009) showed that a higher data sampling rate can mitigate the effect of the noise. Because GRAPHIC clock and ambiguities cannot be separated, the raw (undifferenced) pseudorange is included to anchor the clock solution; but it is de-weighted relative to the GRAPHIC observable such that it does not significantly impact the final orbit solution. The half-integer ambiguity that remains in the GRAPHIC observable can be estimated as a constant over the duration of the satellite transmitter pass.

### 2.3 GipsyX Orbit Modeling

This section describes the relevant aspects of POD within the GipsyX 1.1 software environment. GipsyX is used for nearly all orbit determination operations and research activities using GNSS at JPL. It derives its heritage from the legacy GIPSY/Oasis software (dating back to the 1980s), and is thus well-established for applications of GNSS to positioning, timing, navigation, and geodetic science (Bertiger et al., 2020). Our single-frequency POD employs the high-fidelity EIGEN-GRGS.RL03-v2. MEAN-FIELD gravity model, solid Earth tide, pole tide, ocean tides, and third-body effects as listed in Table 2. Non-gravitational effects are computed based on a custom macromodel constructed for each spacecraft, comprised of individual plates defined by their body frame normal vectors, areas and reflection properties (specularity and diffusivity). This model, along with knowledge of satellite attitude, provides input for the computation of drag and radiation pressure forces on the spacecraft. Drag force is computed based on the DTM2000 model, driven by the F10.7-cm solar flux and Kp geomagnetic indices to compute the atmospheric density (Bruinsma et al., 2003). Radiation pressure force is modeled for both solar (Milani et al., 1987) and Earth albedo (Knocke et al., 1988) radiation sources. GPS satellite ephemerides and clock solutions, including Earth orientation parameters, come from the JPL International GNSS Service (IGS) Analysis Center final products (Johnston et al., 2017).

The GRACE macromodel has already been well-established in the GipsyX environment (Bertiger et al., 2020). For CYGNSS, we created a new macromodel with eight surface areas based on engineering drawings. Coefficients for their specular and diffuse reflectivity were estimated in GipsyX using dynamic orbits over approximately 170 days. We then computed a weighted mean of the daily results, using the inverse of the filter formal errors to arrive at the fixed values for CYGNSS specularity and diffusivity per surface element, as listed in Table 3 and illustrated in Figure 1.

### 2.4 POD Processing Strategy

GipsyX offers a variety of options for orbit determination such as the use of empirical acceleration estimation, configurable reduced-dynamic models, and antenna calibration maps (Bertiger et al., 2020). Reduced-dynamic approaches have been shown to improve orbit solution accuracy by allowing the filter to estimate accelerations that are not captured in the dynamical models (Wu et al., 1991). They also provide the ability to estimate these accelerations in a stochastic manner by setting the correlation time and update frequency along with the a-priori uncertainty and the amount of process noise injected at each update. This can be applied to any estimated parameter such as custom empirical accelerations, providing a flexible method to account for the mismodeling of forces due to drag, radiation pressure, and higher-order gravity.

We apply the same GipsyX POD strategy to both the code-only and GRAPHIC methods. Daily POD solutions are generated using 30-hour data sets centered at noon GPS time. A 30-hour measurement set allows for orbit evaluation using orbit overlaps as well as avoiding edge effects in the central 24 hours of the solution (Bertiger et al., 2010a; Haines et al., 2004). GipsyX is initialized by first estimating a dynamic orbit fit to the Nav solution. This is followed up with a two-step iteration process using pseudorange and phase measurements. The first iteration is comprised of three dynamic passes where the initial state and dynamic properties are updated after each iteration. This is followed by a final reduced-dynamic pass. The parameters for each estimation strategy are shown in Table 4. The a-priori *σ* sets the initial uncertainty of the filter. The predominant orbit errors, including those due to mismodeled dynamical effects such as drag and solar radiation pressure, manifest at a once-per-orbit frequency. To account for these effects, empirical accelerations in the cross-track and in-track directions are estimated in terms of sine and cosine components using:

4

where *C* and *S* are the estimated sine and cosine amplitudes, *θ* is the argument of latitude, and is in the direction of the desired component (in-track and cross-track). These parameters, along with the drag coefficient, are estimated as constant biases for each 30-hour data set in the dynamic passes. This is followed by a reduced-dynamic pass where stochastic radial, cross-track, and in-track accelerations are estimated with a correlation time of 6 hours and held constant across the update interval of 30 minutes. The constant amplitude once-per-orbit cross-track and in-track estimates from the dynamic pass are augmented by stochastic states with correlation times of 6 hours and update intervals once-per-orbit. The reduced-dynamic pass enables the final orbit solution to follow the measurements more closely by allowing the stochastic estimates to slowly change during the 30-hour data set.

### 2.5 Antenna Calibration Maps

Calibration maps are used to compensate for angle-dependent measurement biases (e.g., multipath) caused by the receiver antenna and its electromagnetic interaction with the spacecraft (Haines et al., 2004). The maps apply corrections to GPS observations using a popular model that decomposes the biases into a phase center offset (PCO) and an azimuth- and elevation-dependent phase center variation (PCV; Jäggi et al., 2009). The PCO can be described as a vector, **r**_{0}, that represents a fixed adjustment to the spacecraft antenna reference point. The PCV is a function, *ϕ*(*θ, φ*), that applies corrections based on the line-of-sight azimuth, *θ*, and elevation, *φ*. Rothacher et al. (1995) showed that there are inherent degrees of freedom with this PCO-PCV model, in that one can shift the phase center offset and compensate for it with a corresponding adjustment in the phase center variation function, resulting in the same antenna calibration map. That is, a PCO vector, **r**_{0}, and PCV azimuth- and elevation-dependent function, *ϕ*(*θ, φ*), can be transformed into a new set of and *ϕ*′(*θ, φ*) that will give identical results when using GPS data as given by:

5

6

Here, the unit vector, **e**(*θ, φ*), points from the receiver to the transmitter, and Δ*ϕ* represents an angle-independent delay offset. This delay offset cannot be disambiguated from the receiver clock bias; instead, it can be uniquely defined by artificially constraining either *ϕ*′(*θ, φ*) in the boresight direction to zero or by constraining the mean value of *ϕ*(*θ, φ*) to zero. Estimation of the antenna calibration map is done using dynamic orbits without constant, body-fixed empirical accelerations as these could induce undesirable offsets in the map. For this analysis, it is assumed that there is no a-priori information regarding either the PCOs or the PCVs. This is potentially problematic for yaw-fixed spacecraft, such as GRACE and CYGNSS, as the in-track direction of the map experiences reduced observability. In our approach, we rely on the orbits to inform the use of additional constraints on the PCOs.

The PCO-PCV model can be applied to both the phase and code data types with the code being highly susceptible to multipath errors, resulting in much larger PCVs (Haines et al., 2004). The antenna calibration map within GipsyX is applied to each data type using a grid of azimuth and elevation vertices with a corresponding calibration value. Application within the filter, then, uses a bilinear interpolation for individual measurements. The estimation of the antenna calibration map is done using the direct approach recommended by Jäggi et al. (2009), where corrections to the a-priori PCV values are estimated within the orbit determination filter. The advantage of this method is the ability to estimate the orbit, receiver clock, carrier ambiguities, and PCVs simultaneously. The GipsyX estimation strategy uses 24-hour dynamic orbits. Orbit estimates are first processed with estimated drag coefficient and once-per-rev empirical accelerations for the cross-track and in-track components. A final pass is performed where only the satellite orbit, clock, and calibration values are estimated, holding the previous dynamical estimates as fixed parameters. To account for the reduced measurement density in the zenith direction, the estimated map vertices use a fixed spacing of 2 degrees for elevation and a variable azimuth spacing as shown in Table 5. This results in a total of 4,242 vertices that are estimated from each 24-hour solution. The output of each daily solution is saved, and then accumulated into a final solution over the desired time span. It is during accumulation that constraints are applied to the PCVs. Here, we constrain the sum of *ϕ*(*θ, φ*) above 30 degrees of elevation to be zero.

### 2.6 CYGNSS Clock Inconsistency

The first attempts at processing the CYGNSS measurements revealed an unexpected inconsistency between the pseudorange and carrier-phase measurements. Typically, GPS pseudorange and carrier-phase measurements made on the same receiver channel have the same clock offset, and diverge from each other systematically only due to the difference in the sign of the ionospheric effect as follows:

7

With CYGNSS however, we observed different apparent clock offsets in the pseudorange and carrier measurements that were common across all simultaneously tracked GPS satellites. An example of the resulting pseudorange/carrier divergence is shown in Figure 2 (left) for a single satellite pass. This rate of divergence is much larger than would be expected due to ionosphere or noise alone, and it contains significant ramp/jump structure. When the code and phase are separately detrended, as shown in Figure 2 (right), the ramp/jump behavior is shown to be present in the phase measurements and not the code. To further isolate the effect, it is possible to process the pseudorange and carrier separately in GipsyX, solving for independent clock solutions, as shown in Figure 3. The left panel shows solutions for receiver clock offset over a full day using only the pseudorange (blue) and only the phase (red). The clock solution behavior revealed by the pseudorange data shows small offsets from zero of only a few meters with little variability. This is quite typical for a GPS receiver and suggests that the receiver’s internal oscillator is being steered to GPS time or that the measurements are compensated onboard for the estimated offset of the oscillator. In contrast, the clock behavior revealed by the phase data features a large linear drift of tens of km per day with jumps of tens of meters about the linear trend. A zoomed-in view of Figure 3 (right) shows similar behavior as those observed in the detrended phase seen earlier in Figure 2.

These features in the phase-only clock solution suggest that the phase measurements are either uncorrected for the clock or corrupted in a clock-like way (or both)^{1}. This precludes the use of conventional pseudorange/carrier divergence methods to estimate and remove the ionospheric effect and limits our analysis to just two approaches for compensating for ionospheric effects in POD: single-frequency code-only observation processing with ionospheric model corrections or the GRAPHIC method. The GRAPHIC data combination will result in an apparent clock offset that is the average of the pseudorange and phase. With the SGR-ReSI phase vs pseudorange discrepancy, the issue that arises is that the phase measurements appear inconsistent with the time tags. Timing errors are correlated with in-track orbit errors, and in-track differences greater than 1 m were initially observed in GRAPHIC overlaps of the 30-hour POD solution arcs. Because of this, an effort was undertaken to correct the most significant portion of the clock inconsistency in the phase. Starting with Equation (7), the difference, *d*, is now modified to include the non-common clock-like phase error, *b _{ϕ}*:

8

To remove the integer ambiguity, we take the time difference for all tracked satellites, *k*, at time *t _{i}* and

*t*

_{i+1}as follows:

9

Because the Δ*b _{ϕ}* term is common to all satellites, we can use the assumption that the median values at each measurement interval of Δ

*I*and Δ

*ε*are much smaller than Δ

*b*and are near zero. We use the median because it is less susceptible to large outliers that can occur for low-elevation satellites that experience significant changes in ionosphere over short periods of time. The phase clock term at any given time can then be constructed recursively using the median from

_{ϕ}*t*

_{0}to the desired time,

*t*, by initializing using

_{i}*b*(

_{ϕ}*t*

_{0}) = 0 and using:

10

The left plot in Figure 4 shows a 10-minute span with the computed Δ*d ^{k}* (orange) and the median values (blue), along with the reconstructed phase clock term (right). The reconstructed clock term is, then, removed from the phase measurements. This creates a clock bias time history for the phase measurements that is much more consistent with the pseudorange measurements and, more importantly, consistent with the measurement time tags. The assumption of Δ

*I*and Δ

*ε*being near zero does not need to hold true in the short term; rather, the sum only needs to stay near zero in the long term. It will still contain some residual ionosphere and measurement noise at each epoch, and so this clock effect prevents us from determining the ionospheric effect based on code/carrier divergence. Fortunately, with this adjustment to the clock bias, we can now reliably employ the GRAPHIC approach to remove the ionospheric effects.

### 2.7 Code-Only Corrections

The processing of the code-only measurements requires additional corrections for the timing group delay (TGD) bias and CA-P1 differential code biases (DCBs). These corrections are necessary to make the pseudorange measurements consistent with the standard clock products used in GipsyX that are referenced to P(Y) dual-frequency ionosphere-free measurements. An ionosphere correction is made using the IGS GIMs to model the total electron content (TEC; Hernández-Pajares et al., 2009). We follow the method described by Montenbruck and Gill (2002) with a few modifications. Where Montenbruck and Gill integrated a Chapman profile fit to the IRI model, we instead numerically integrate a single representative IRI2016 (Bilitza et al., 2017) distribution directly, with the purpose of defining the ionosphere pierce point. From this, the portion of the vertical TEC (vTEC) based on the GIM is then mapped to correct the ionospheric delay in the LEO pseudorange. For both GRACE and CYGNSS, we select a representative profile near the daily global TEC maximum. For GRACE, the reference IRI2016 profile is from April 9, 2009, at 12:00 UTC located at N 0 deg, E 20 deg, and for CYGNSS, October 2, 2019, at 12:00 UTC, located at S 5 deg, E 25 deg. These two profiles are shown in Figure 5. One limitation of the IRI2016 model is that it does not converge to zero at the upper altitude limit of the model (2,000 km). The remaining portion above 2,000 km is small, however, and we approximate it as zero above this altitude.

For a given LEO satellite with a height, *h _{s}*, the ionospheric pierce point,

*h*, is defined as the altitude where 50% of the total ionosphere above the satellite is located. This can be computed using the condition:

_{Ip}11

where *vTEC _{i}* is the normalized electron density at altitude

*i*from the IRI2016 profile in Figure 5. This relationship, along with the satellite altitude, allows for computation of the height of the ionospheric pierce point,

*h*. The remaining fraction,

_{Ip}*α*, of the

*vTEC*from the IGS GIM above the satellite can then be computed using:

12

Finally, the L1 C/A delay can be computed using the mapping function as follows:

13

where *E _{IP}* is the calculated elevation of the signal ray path at the intersection of the ionospheric pierce point. This produces an estimate of the ionosphere delay that can be removed from each measurement in a pure single-frequency code solution.

## 3 RESULTS AND ANALYSIS

To demonstrate the single-frequency orbit solution methods proposed for CYGNSS, we processed 200 days of data for GRACE-B and 175 days of data for one of the eight CYGNSS satellites, FM05. We limit the analysis here to a single satellite from the CYGNSS constellation to demonstrate the POD evaluation, but with the understanding that it can readily be applied to the rest of the CYGNSS constellation. For GRACE, 30-second C/A L1 measurements from DOY 1 – 200, 2009, were processed using the code-only and GRAPHIC methods. For CYGNSS FM05, 15-second C/A L1 measurements from DOY 190 – 365, 2019, were used. The CYGNSS data sets have several large gaps of approximately a day in length on DOY 212, 252, and 353. Large data gaps such as these prevent the formation of 30-hour arcs of continuous data on either side of the missing data, resulting in 3 missing days for each gap. Additionally, we use a single day of Sentinel-6 MF data to investigate a cross-track anomaly observed in CYGNSS but not in GRACE. First, we present the estimated code-only and GRAPHIC antenna calibration maps for both GRACE and CYGNSS. Next, we examine the internal metrics consisting of daily post-fit measurement residual RMS values along with the orbit overlap statistics. The GRACE code-only and GRAPHIC orbits are compared to the dual-frequency orbits to evaluate the associated errors. Finally, we use the GRACE and Sentinel-6 MF errors to understand the potential errors associated with the CYGNSS POD.

### 3.1 Antenna Calibration Maps

Figures 6 and 7 show the calibration maps estimated in GipsyX using the code-only and GRAPHIC methods for GRACE and CYGNSS, respectively. Both sets of antenna calibration maps exhibit significant systematic measurement biases. For GRACE-B, the code multipath effects can cause errors on the order of 30 to 40 cm and, for CYGNSS, the variability is as high as several meters. Given that the GRAPHIC observable contains half of the code multipath errors, it is not surprising to see that the GRAPHIC antenna calibration map is dominated by the code multipath for both GRACE and CYGNSS, but with half the magnitude (notice map scale). For GRACE, this deviates to some degree at the lower elevations where it is likely that deficiencies in the code-only ionosphere correction result in ionospheric errors being included in the map. The effect is more pronounced for GRACE due to its lower orbit altitude (more ionosphere) and smaller multipath effects.

In Figure 7, for comparison, we also provide a group delay map produced from anechoic chamber measurements using a CYGNSS zenith antenna mounted on an approximate CYGNSS satellite mockup. For this map, the group delay is computed directly as the derivative of the measured right-hand circular polarized (RHCP) phase with respect to frequency at the L1 center frequency. We observe that the CYGNSS on-orbit and anechoic chamber antenna maps are qualitatively very similar, which gives us confidence in the results produced by GipsyX. Since both the specific antenna and satellite differ from the CYGNSS FM05 on-orbit case, we should expect some differences in the two antenna maps. For POD, we will utilize the GipsyX maps based on actual on-orbit data since these should be more accurate.

Initial evaluation of the GRACE orbits showed a smaller than 2 cm in-track bias (described in more detail below) for the GRAPHIC orbits. However, the code-only antenna calibration map produced a consistent bias larger than 6 cm. A similar bias was observed between the CYGNSS code-only and GRAPHIC orbits. Given these results, we chose to constrain the spacecraft flight direction PCOs of the code-only map for both GRACE and CYGNSS to zero.

### 3.2 Orbit Quality Assessment

The internal metrics presented here provide a means to evaluate the overall consistency of the daily POD solutions. Two different metrics are shown here: the daily residual RMS and the orbit overlap differences. The daily residual RMS is a measure of the quality of fit of the daily reduced dynamic POD strategy. Figure 8 shows the statistical distribution of the daily GRACE-B residual RMS values for the time frame from January 2009 – July 2009. The daily GRAPHIC RMS values are very consistent with a mean and standard deviation of 3.5 ± 0.04 cm. This would appear to be consistent with code noise of about 7 cm. However, the code-only residual RMS values are significantly larger at 27.6 ± 3.2 cm. This is due primarily to errors in the ionosphere correction and remaining DCBs. This gives GRAPHIC a distinct advantage, given it is ionosphere-free and able to absorb the code DCBs in the carrier bias estimate. Turning now to CYGNSS, Figure 9 shows the statistical distribution of the daily residual RMS from July 2019 – December 2019. For CYGNSS, the daily GRAPHIC residual RMS values are much larger than for GRACE but still consistent with a mean and standard deviation of 54 ± 0.8 cm. This would correspond to an expected code noise for CYGNSS of 108 cm, which is close to the code-only mean and standard deviation of 116 ± 2.0 cm. When compared to GRACE, the CYGNSS GRAPHIC residuals are more than an order of magnitude larger, representing a substantial increase in measurement noise. This is understood to be a combination of the quality of the receiver measurements and remaining multipath not adequately captured in the antenna map.

Orbit overlaps are a way to evaluate the precision or consistency of the POD solutions over time and are computed from the six hours of overlap that occur for 30-hour solution arcs centered at noon GPS time. Daily RMS statistics are computed from the component differences of the central four hours of the six hours of overlap to avoid edge effects from the reduced dynamics. Figure 10 shows the distribution of the GRACE code-only and GRAPHIC statistics. Here, we can see that the GRAPHIC POD solutions are much more consistent than the code-only method for all components. Table 6 shows the average and standard deviation of all the component statistics. The mean and standard deviation of the GRAPHIC overlap RMS statistics show very consistent results of 0.5±0.2 cm for radial and cross-track and 1.2±0.5 cm for the in-track. This is a significant improvement over the code-only method with a more than six-fold improvement in all three components.

We compute the same orbit overlap RMS statistics for CYGNSS, shown in Figure 11 and Table 7. Like GRACE, the GRAPHIC POD solutions are more consistent than the code-only method for all components, but in this case, they show reductions by a factor of 1.5 for radial, 5 for cross-track, and 1.8 for in-track. The improvement is not as substantial for the CYGNSS GRAPHIC results as for GRACE because of the higher noise level of the code measurements that do remain, reduced by a factor of two, in the GRAPHIC observable. These results show a substantial improvement in orbit consistency with GRAPHIC, but cannot confirm the improvement in accuracy. For that, we turn to an external reference comparison using both GRACE and Sentinel-6 MF missions.

### 3.3 GRACE Single-Frequency Errors

Because the orbit knowledge of the GRACE satellites is critical to its primary mission of mapping the gravity field, reduced-dynamic dual-frequency orbits have been described by Bertiger et al. (2010b). These orbits have an overlap precision of 1.4 mm radial, 2.2 mm cross-track, and 2.3 mm in-track, as well as an accuracy of 2 mm for the distance between the satellites, as confirmed by comparison with independent K-band inter-satellite range measurements. These validated orbits can be treated as truth orbits when compared to the single-frequency solutions. Here, we first examine the biases associated with code-only and GRAPHIC POD solutions relative to the dual-frequency solution by computing component differences of the central 24 hours of the 30-hour arcs. The bias value is useful for evaluating the antenna calibration map where it is possible to both introduce a bias or correct for one if it is known. The daily bias statistics for the radial, cross-track, and in-track components are shown in Figure 12. Both the code-only and GRAPHIC solutions show sub-cm bias statistics for the radial and cross-track components, indicating very good performance here. The in-track bias is the largest for both the code-only and GRAPHIC solutions.

The solution RMS errors are computed from the same component differences. Figure 13 shows the daily RMS statistics for the radial (left), cross-track (middle), and in-track (right). The GRAPHIC solution error RMS is markedly smaller than the code-only solution by more than a factor of two for all components. Table 8 contains the overall mean and standard deviation of the daily statistics of the code-only and GRAPHIC errors. The GRAPHIC method performs significantly better than the code-only method with an average 3D RMS error of 3.1 ± 0.4 cm with an overall factor of three improvement when compared to the code-only solution errors at 9.3 ± 3.0 cm. Relating these errors back to the overlap statistics for each method, we can see that the average code-only overlap RMS values are consistent with the overall error RMS, but that for GRAPHIC, the overlap RMS statistics are optimistic by about a factor of two. Overall, the GRAPHIC method shows much better performance than the code-only method with a smaller in-track bias and reduced errors in all three components.

### 3.4 Code-Only and GRAPHIC Comparison

In the absence of dual-frequency measurements or other independent metrics to evaluate the CYGNSS single-frequency orbit accuracy, we leverage the GRACE results to understand the errors in the CYGNSS POD. We can conclude from the internal metrics that the GRAPHIC solution is more precise, and, from the GRACE external orbit comparison, that errors are smaller in all three components when compared to the code-only solution. While these conclusions are useful, they do not provide a way to directly quantify the size of the CYGNSS GRAPHIC orbit errors except to say that the GRAPHIC overlap metrics are optimistic and place a lower bound on the errors.

An additional way to quantify the uncertainty in the GRAPHIC solution is to compare it directly to the code-only solution. Figure 14 shows a comparison of the statistics for the component differences in the GRACE code-only and GRAPHIC solutions along with the associated errors. The solution differences are dominated by the code-only errors and can be used to place an upper bound on the GRAPHIC errors. For GRACE, the optimistic lower bound derived from GRAPHIC overlaps, combined with the pessimistic upper bound from the orbit solution differences between code-only and GRAPHIC solutions, produces a range of possible GRAPHIC orbit errors between 0.5–2.3 cm for radial, 0.5–3.0 cm for cross-track, and 1.2–8.4 cm for in-track. This is in good agreement with the actual GRACE GRAPHIC errors of 0.9 ± 0.13 cm for radial, 1.1 ± 0.14 for cross-track, and 2.5 ± 0.35 for in-track, derived from comparisons with dual-frequency reference orbits. The same comparison between the code-only and GRAPHIC solutions, but for CYGNSS, is shown in Figure 15. Here, the component differences for CYGNSS are much larger than those seen for GRACE. It is also notable that the component with the largest mean RMS is the cross-track, resulting in a more than 30-cm difference. A summary of the mean and standard deviation of the daily bias and RMS statistics is shown in Table 9. Once again, we can use the GRAPHIC overlap statistics from Table 7 to approximate the CYGNSS GRAPHIC solution lower error bound and the code-only and GRAPHIC solution differences from Table 9 for the upper error bound. This produces a potential error range for the CYGNSS GRAPHIC orbits of 1.4–5.3 cm for radial, 1.2–31 cm for cross-track, and 3.0-15 cm for in-track.

It is interesting to observe that the cross-track differences between the CYGNSS GRAPHIC and code-only solutions are consistent with a very small inclination change. The cross-track difference has a frequency of once-per-orbit with a relatively constant amplitude, and the largest difference occurs at the highest/lowest latitude of the CYGNSS orbit. These differences persist in dynamic-only orbit solutions with and without any cross-track empirical accelerations. This implies a persistent signal difference between the pseudorange and GRAPHIC observable which would indicate either an ionosphere-related signal, or some behavior related to the GPS transmitter DCBs, transmitter antenna calibration maps, or clock solutions. To test the theory that the code-only solution is the one that contains the error, a single day of Sentinel-6 MF is examined. It is in a roughly circular orbit at an altitude of 1,336 kilometers with an inclination of 66 degrees. Figure 16 shows a time history comparison of the CYGNSS cross-track GRAPHIC and code-only differences (left) and Sentinel-6 MF cross-track differences for a GRAPHIC and code-only solution relative to the dual-frequency solution. Here, we see a similar pattern in the Sentinel-6 MF code-only cross-track differences with a once-per-orbit variability indicating that large differences observed for CYGNSS are a result of errors in the code-only solutions.

The results in Figure 16 are consistent with the conclusion that the code-only solution cross-track errors are much larger than the GRAPHIC cross-track errors. Here, the CYGNSS GRAPHIC/code-only difference has an RMS of 28 cm compared to the Sentinel-6 MF code-only error RMS of 11 cm. This is in contrast to the Sentinel-6 MF GRAPHIC error RMS which is an order of magnitude smaller at 1 cm. This observation shows that the cross-track signal in the code-only solution is not present in the GRAPHIC solution. If it were, we would expect to see a cross-track error in the GRAPHIC solution with half the amplitude, but this is not the case. From this, we conclude that the CYGNSS GRAPHIC/code-only differences, like GRACE, are dominated by the code-only errors. Further investigation of the fundamental causes for this behavior would need to be done to fully explain the observed cross-track differences of the code-only solutions. Such an investigation is beyond the scope of this paper, but this suggests that potentially large errors (>25 cm) in a code-only solution can exist that would be hidden when assessing the orbits through overlaps alone. Without any method of independent evaluation of the CYGNSS GRAPHIC orbits, we rely on the lessons learned from GRACE and Sentinel-6 MF. Applying the fact that the GRACE GRAPHIC overlaps were optimistic by a factor of two places the CYGNSS orbit errors at 2.8 cm radial, 2.4 cm cross-track, and 6 cm in-track. Overall, the GRAPHIC method is far more precise and accurate than the code-only solution in the GipsyX environment and is the preferred single-frequency method for producing orbits.

## 4 CONCLUSION

Effective strategies for POD with single-frequency CYGNSS data in GipsyX have been introduced and demonstrated on more than 5 months of data from the CYGNSS FM05 satellite in the context of assessing its applicability for altimetry. Parallel analysis of data from the GRACE mission provides a basis for quantifying the accuracy of the CYGNSS solutions. With GRACE, the GRAPHIC method, when compared to the code-only method, is shown to be both more accurate and precise. A bound on the 1*σ* CYGNSS GRAPHIC orbit errors can be determined from the GRAPHIC overlap statistics and the differences between the two solutions. This approach places a potential 1*σ* range of 1.4 – 5.3 cm for radial, 1.2 – 31 cm for cross-track, and 3.0 – 15 cm for in-track with a 3D position error of 3.5 – 40 cm. Because the solution differences are dominated by the code-only errors as shown in the Sentinel-6 MF example, GRAPHIC solutions will be closer to the lower bound, particularly for the cross-track component. Based on analysis of GRACE POD solutions using the same methodology, and in particular the GRAPHIC overlap statistics, we estimate that the 1*σ* CYGNSS GRAPHIC orbit accuracy is 2.8 cm radial, 2.4 cm cross-track, and 6 cm in-track, placing the 3D accuracy at less than 10 cm. With this method and the assessment of the improved orbit accuracy in place, researchers will be able to better determine the potential utility of the CYGNSS constellation for GNSS-R-based altimetry.

## HOW TO CITE THIS ARTICLE

Conrad, A. V., Axelrad, P., Haines, B., Zuffada, C., & O’Brien, A. (2023). Improved GPS-based single-frequency orbit determination for the CYGNSS spacecraft using GipsyX. *NAVIGATION, 70*(1). https://doi.org/10.33012/navi.565

## CONFLICT OF INTEREST

The authors declare no potential conflict of interests.

## ACKNOWLEDGMENTS

The research reported here is supported by the NASA ROSES PO award NRA NNH16ZDA001N, 2017-2020. Part of the research was carried out at the Jet Propulsion Laboratory, Caltech, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). The authors would like to thank the CYGNSS project for providing the ability to downlink the raw GPS measurements that have made this work possible.

## Footnotes

↵1 We note that the phase data from the SGR-ReSI receiver are not needed to satisfy any of the mission requirements for the primary CYGNSS mission. The onboard GPS navigation solutions underpinning the primary science mission are conventional and use only the pseudorange measurements. That the phase measurements can be recorded and downlinked, however, is a useful feature for POD studies focused on value-added science (e.g., ocean altimetry).

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.