Improved GPS-Based Single-Frequency Orbit Determination for the CYGNSS Spacecraft Using GipsyX

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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 1

CYGNSS macromodel surfaces with each surface label corresponding to the values given in Table 3

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₀, 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₀, 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2

CYGNSS code/carrier divergence (left) and detrended code/carrier (right)

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

CYGNSS code-only and phase-only clock solution comparison: DOY 191 of 2019

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)¹. 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₀ to the desired time, t_i, by initializing using b_ϕ(t₀) = 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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

Median of first difference of code minus phase (left) and estimated phase clock error (right)

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.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

CYGNSS and GRACE IRI2016 normalized election density profile

For a given LEO satellite with a height, h_s, the ionospheric pierce point, h_Ip, is defined as the altitude where 50% of the total ionosphere above the satellite is located. This can be computed using the condition:

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_Ip. The remaining fraction, α, 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.