Precise Orbit Determination of LEO Satellites Based on Undifferenced GNSS Observations

The kinematic POD method estimates the orbit by processing only the GNSS measurements with a precise positioning approach without considering the dynamic models. Hence, it is susceptible to satellite geometry, measurement errors, noise, and data outages, and orbit estimation will not be accurate or even available when the GNSS data have large biases or gaps. This something could happen in small LEO satellites with limited power capability, such as CubeSats (Wang et al. 2020). Kinematic POD is a popular method in recovering the gravity field information because a priori knowledge of the involved physical forces is not required (Jäggi et al. 2016). The use of this method for the gravity field determination was demonstrated by Švehla and Rothacher (2003b), and the first gravity field model based on kinematic orbits was calculated by Gerlach et al. (2003).

The kinematic POD in postmission mode starts with the quality control of the observations, followed by orbit determination based on the precise point positioning (PPP) approach. These steps are discussed subsequently.

As discussed previously, the dynamic and kinematic POD have some limitations. The former requires comprehensive dynamic models that are not precisely available for LEO satellites (especially the nongravitational-force models), and the latter is sensitive to data outages and weak GNSS satellite geometry. The reduced-dynamic POD, however, can benefit from both sides, i.e., the availability of dynamic models, as well as GNSS observations. In this method, the initial conditions, force model parameters, and possible stochastic accelerations are estimated instead of the epochwise LEO positions in the kinematic POD.

Because this POD method delivers more accurate orbital positions than the other two approaches, it has been used as the primary method to determine the orbit of different LEO missions, such as TOPEX/Poseidon (Yunck et al. 1994), CHAMP (Švehla and Rothacher 2003b) Jason-1 (Haines et al. 2004), GRACE (Kang et al. 2003), GOCE (Bock et al. 2014), Swarm (Van Den Ijssel et al. 2015), HY-2A (Guo et al. 2015), FY-3C (Zhang et al. 2018), Sentinel (Fernández 2019), and GRACE-FO (Kang et al. 2020). In the following, different steps of this approach are discussed.

The kinematic and the reduced-dynamic POD using single- and dual-frequency observations are used in different LEO satellite–related studies. Table 2 presents a selection of LEO POD accuracies achieved with the different validation methods based on undifferenced dual-frequency measurements with both the kinematic and reduced-dynamic methods used in different missions. These results were obtained from different studies over a rather long timespan (2005–2020). The presented POD accuracies are therefore affected by the progress in the gravitational force modeling achieved by different gravity field missions over these years. Satellites in higher orbits often have lower root mean square (RMS) errors, mainly due to the harsh dynamical environment of satellites in the lower altitudes (e.g., Swarm-B satellite in Table 2 that has a higher altitude than Swarm-A and Swarm-C).

The orbital accuracy of the kinematic POD method is correlated to the epochwise receiver clock offset determination. There is a high correlation between the radial component of the orbit and the receiver clock parameters (Weinbach and Schön 2012). Stable oscillators thus enable a robust clock model and can stabilize the radial orbital component. This concept is simulated for the GRACE mission to model the receiver clock offset, and a 40% improvement in the radial direction is shown (Yang et al. 2014).

In all the aforementioned studies, the first-order ionospheric delay is removed by forming the IF linear combination. However, as discussed in the “Measurement and Error Models of Kinematic POD” section, the uncombined observations can also be applied in the POD. The kinematic POD with raw measurements, considering the higher-orders of the ionospheric delays, provided similar orbital accuracy compared with using the IF combination (Zehentner and Mayer-Gürr 2016).

Although this paper reviews LEO POD methods based on using undifferenced GNSS measurements, it is worth mentioning that differencing of the GNSS observables were also used for LEO POD, where GNSS observations collected on a LEO satellite are differenced with GNSS observations at a known ground reference station. The double differencing of GPS phase measurements with ambiguity resolution for the kinematic POD of CHAMP was carried out by Švehla and Rothacher (2003a). The double differences (DD) were formed between the IGS permanent stations and the LEO satellite to reduce the orbit- and clock-related errors from the observation model. However, with such long baselines, the remaining orbital and clock errors still need to be considered. Besides, tropospheric errors need to be considered for the observations at the ground stations. Furthermore, the concept of DD between two GRACE satellites formed the first GPS baseline in space with the millimeter-accuracy after performing ambiguity resolution (Švehla and Rothacher 2004). Such a baseline can be used for the orbit validation and gravity recovery in the GRACE-FO mission (Kang et al. 2020).

Single receiver integer ambiguity resolution (IAR) can effectively improve the positioning accuracy; however, it requires corrections for the hardware biases to recover the integer nature of the ambiguities. Different IAR methods were reviewed comprehensively by Teunissen and Khodabandeh (2015). Among the available single-receiver-based IAR methods, the methods proposed by Laurichesse et al. (2009) and Bertiger et al. (2010) are usually used for the LEO POD. The former method is based on identifying the satellite wide-lane biases, using them to estimate the undifferenced wide-lane ambiguities for a network, and finally fixing L1 ambiguities for these receivers.

A similar approach has been tested for GRACE and Jason-1, i.e., their wide-lane measurements were corrected by the GPS fractional wide-lane corrections provided by the network. Next, the GPS satellite integer clocks derived from the network were used to estimate the IF phase residuals. Finally, by forming between-satellite single differences of these residuals, the receiver clock offset was removed and the integer ambiguities were estimated. The orbits of these satellites determined by the unambiguous measurements were shown to be more accurate than the orbits from float solutions, i.e., those estimated with float ambiguities.

In the method proposed by Bertiger et al. (2010), GPS orbit errors and clock offsets, dual-frequency phase biases, and corresponding wide-lane values were estimated in a network solution and were used to form integer constraints in the LEO data processing. The method has been tested for the Jason-2 mission, and it was shown that the radial orbital accuracy increased significantly, i.e., reaching below 1 cm, fulfilling the altimetry requirements.

The wide-lane ambiguities have a longer wavelength than the original ambiguities, which eases the ambiguity resolution by reducing the number of potential ambiguity candidates. In addition to theaforementioned methods, the IAR for Sentinel-3A data has been performed by Montenbruck et al. (2018) using the wide-lane bias products and the integer phase clock offset corrections from the Centre National d’Etudes Spatials (CNES). The orbit generated with fixed ambiguities outperformed the orbit of float solutions. The RMS value (validated by the SLR) reached 5 mm, which has a 30% improvement. The systematic biases in the cross-track direction were also revealed through the IAR due to the errors in the antenna PCOs, which can be absorbed by the ambiguities in the float solution.

In addition to the kinematic and the reduced-dynamic methods, a fourth POD approach was introduced by Švehla and Rothacher (2005a), called the reduced-kinematic POD. However, it has not been developed further because it was used inside the kinematic POD only for the gravity field determination (Švehla 2018).

In summary, both the kinematic and the reduced-dynamic POD in the postmission processing can deliver orbital accuracy at the centimeter-level, which can be improved by the IAR. This accuracy fulfills the requirements of most postprocessed LEO applications (Montenbruck 2017). The reduced-dynamic POD makes use of the strong dynamic models and can normally deliver higher orbital accuracy than the kinematic mode. However, the extensive dynamic modeling requires high computational load, which is not a big issue in postmission by the ground-based processing centers, but could be an issue for onboard processing, especially for small satellites such as CubeSats.

The state vector for the kinematic POD at each epoch $t$ is expressed as $\mathbf{x}={(\begin{array}{c}{\mathbf{r}}^T,dt,{\mathbf{b}}^T\end{array})}^T$. The processing is initialized with the initial values of the filter state (${\mathbf{x}}_0$) and its covariance matrix (${\mathbf{Q}}_0$), which are typically estimated using the least-squares adjustment. Because no dynamic model was used in the kinematic POD, an identity matrix is considered as the transition matrix (${\mathbf{\varphi }}_t^{\mathrm{KIN}}=\mathbf{I}$, where KIN refers to kinematic) to update the clock offset and ambiguities combined in ${\hat{\mathbf{x}}}^{\prime }$ as follows:where the process noise matrix is diagonal, i.e., ${\mathbf{Q}}_t^{\mathrm{KIN}}=\mathrm{diag}({\sigma }_{dt}^2,{\sigma }_{\mathbf{b}}^2)$, which represents the process noise of the clock offset modeled as random-walk, with the term $\mathbf{b}$ due to the hardware biases included in it. If cycle slips take place, the temporal link of the ambiguities is interrupted.

Before applying the measurement update, data screening process discussed in the “Data Screening for Kinematic POD” section should be performed. During the measurement update, the predicted states can be considered as pseudo-observations that aid GNSS observations of the current epoch, and the unknowns can be estimated using a sequential least-squares adjustment, with the measurement and stochastic models expressed as follows:where ${\mathbf{J}}_{\mathbf{r}}$ and ${\mathbf{J}}_{{\mathbf{x}}^{\prime }}$ = partial derivatives of the observations with respect to the kinematic orbits and the combined clocks and ambiguities, respectively. Although expressed in different forms, the aforementioned sequential least-squares adjustment is equivalent to EKF (Humpherys and West 2010), provided that the same temporal link is considered.

The vector $\mathbf{x}={(\begin{array}{c}{\mathbf{r}}_0^T,{\mathbf{v}}_0^T,{\mathbf{p}}^T,{\mathbf{\alpha }}_{\mathbf{E}\mathbf{M}\mathbf{P}}^T,dt,{\mathbf{b}}^T\end{array})}^T$ is the EKF state vector in the reduced-dynamic POD, where $\mathbf{p}={(\begin{array}{c}{c}_r,c_d\end{array})}^T$ contains the SRP and drag coefficients. For $m$ number of satellites, there are $12+m$ unknowns at each epoch (six satellite state components, two coefficients for the SRP and air drag, three empirical accelerations, one clock offset, and $m$ phase ambiguities). The propagation in time update involves a simple integration method to compute both the satellite trajectory and the partial derivatives concerning the state vector and the term $\mathbf{p}$. The empirical acceleration parameters ${\mathbf{\alpha }}_{\mathbf{EMP}}$ are updated using a Gauss-Markov process, expressedwhere ${\mathrm{\Delta }}_{t,t-1}$ = time difference between the current and previous epochs; $\tau$ = time constant; and $l$ = exponential damping factor. Other state parameters of the filter remain constant during this step with a white process noise model. Considering ${\partial }_{(x,y,z,v_x,v_y,v_z)}{\mathbf{\alpha }}_{{\mathbf{EMP}}_{t|t-1}}=0$ and ${\partial }_{{\mathbf{\alpha }}_{{\mathbf{EMP}}_{t-1|t-1}}}{\mathbf{\alpha }}_{\mathbf{E}\mathbf{M}{\mathbf{P}}_{t|t-1}}=l$, the full structure of the reduced-dynamic EKF transitions matrix ${\mathbf{\varphi }}_t^{\mathrm{RD}}$, where the superscript RD denotes reduced dynamic, is expressedwhere $\tilde{\mathbf{\Phi }}$ and $\tilde{\mathbf{S}}$ are obtained from the numerical integration of the variational equation. The covariance update is expressedwhere the process noise matrix ${\mathbf{Q}}_t^{\mathrm{RD}}$ is determined based on assumptions and simulations that are designed to prevent filter divergence. An example of such a matrix has been given by Montenbruck and Ramos-Bosch (2008), where the empirical accelerations and the receiver clock offsets were modeled using Gauss-Markov and random walk processes, respectively.

Next, the measurement update is performed to improve the filter state vector with the newly available observations, after data screening and cycle slip detection process, as follows:

The gain matrix ${\mathbf{K}}_t$ depends on the covariance matrix of the measurements ${\mathbf{Q}}_{{\mathbf{y}}_t}$ and the propagated covariance from the time update step ${\mathbf{Q}}_{t|t-1}$. The term ${\mathbf{J}}_t$ is the partial derivative for the filter state vector, and ${\mathbf{j}}_t({\hat{\mathbf{x}}}_{t|t-1})$ contains the computed values of the observations.

Finally, an integration method, such as the fourth-order Runge-Kutta (Butcher 2016), can be used to propagate the trajectory of the LEO satellite over a fixed time step to the desired output sampling interval (Montenbruck and Ramos-Bosch 2008). A flowchart of these steps is given in Fig. 4. It is possible to simplify the aforementioned algorithm to meet the computational budget of the small satellites such as CubeSats by considering drag and SRP coefficients as modeled values rather than estimating them, and not estimating empirical accelerations, as explained by Yang et al. (2016). However, it reduces the orbital accuracy to a meter level.

Different studies were carried out to overcome the limitations of real-time POD. One of these limitations is the poor accuracy of GPS orbits and clock offsets provided by the broadcast ephemeris. It was reported by Montenbruck and Ramos-Bosch (2008) that depending on the satellite altitude, an accuracy [defined by the daily three-dimensional (3D) RMS] of 40–60 cm was achieved using GPS dual-frequency observations and broadcast ephemeris. The simulated multi-GNSS measurements from GPS, Galileo, and Beidou-3 satellites were tested by Hauschild and Montenbruck (2020) with a flight-proven algorithm developed by Montenbruck et al. (2008). It was reported that combining all measurements deliver LEO satellite orbits with 10.4-cm 3D RMS error when using broadcast ephemeris.