Research ArticleOriginal Article
Open Access

Robust Modeling of GNSS Orbit and Clock Error Dynamics

Elisa Gallon, Mathieu Joerger, and Boris Pervan
NAVIGATION: Journal of the Institute of Navigation December 2022, 69 (4) navi.539; DOI: https://doi.org/10.33012/navi.539

Abstract

In this paper, we develop new stochastic orbit and clock error models for positioning, fault detection, and integrity monitoring over time. GPS and Galileo orbit and clock data are evaluated and ranging errors are analyzed and modeled over time. This work is intended for time-sequential safety-critical navigation systems including global navigation satellite systems (GNSSs) integrated with inertial navigation systems (INSs) and Kalman filter implementations of Advanced Receiver Autonomous Integrity Monitoring (ARAIM).

Keywords

1 INTRODUCTION

Global navigation satellite systems (GNSSs) can provide continuous worldwide absolute positioning but require visibility of four or more satellites, which is not always achievable in sky-obstructed environments. Also, GNSS is vulnerable to radio-frequency interference. In contrast, inertial sensors are not directly impacted by these external factors. Inertial navigation systems (INSs) can be used as dead reckoning sensors to estimate displacements over time, but state estimation errors drift due to the temporal integration of IMU errors. Combining INS and GNSS, for example, using a Kalman filter (KF) can simultaneously limit the drift in INS positioning error and provide continuity through sky obstructions and robustness against GNSS jamming and spoofing attacks (Tanil et al., 2018). GNSS/INS integration is accomplished through measurement filtering, which requires the robust modeling of stochastic errors over time to ensure navigation integrity.

Another application in which GNSS measurements can be filtered over time is Advanced Receiver Autonomous Integrity Monitoring (ARAIM). For aircraft navigation, the baseline version of ARAIM uses carrier smoothed code (CSC) measurements at one instant in time to provide a snapshot navigation solution (Working Group C, 2012, 2014, 2016). However, in Joerger and Pervan (2020), we showed that the additional leverage of satellite motion over time provides superior positioning performance and tighter protection levels (PLs) than baseline ARAIM. Sequential ARAIM algorithms (using banks of KFs, for example) open the possibility to extend the scope of ARAIM beyond aircraft navigation to rail, harbor, or arctic operations.

To implement sequential ARAIM or inertial-GNSS integration, one must ensure that the error models implemented in the KF properly account for time correlation. In both applications, dynamic models for the three main error sources affecting ionosphere-free GNSS signals are needed: orbit and clock errors, tropospheric delay, and multipath. This paper focuses on GNSS satellite clock and orbit ephemeris errors, which are major contributors to ranging errors for dual-frequency GNSS users.

Building on the prior work of DeCleene (2000) and Rife et al. (2006), Perea Diaz (2019) employed an overbounding theory to define upper bounds on the variances of orbit and clock errors for both GPS and Galileo satellites. These error models are sufficient for snapshot positioning, but not for time-sequential implementations because they do not address the stochastic dynamics of the errors over time.

Langel (2014) derived an analytical bound on integrity risk for time-sequential linear estimators using autocorrelation function (ACF) bounding. We used experimental data to evaluate ACF bounds for GPS and Galileo orbit and clock errors in Gallon et al. (2019). But Langel’s ACF-based approach requires continuous, cumulative storage of all data and estimator coefficients over time and, except for short finite-horizon intervals, is unsuitable for KF implementations. More recently, Langel et al. (2020) introduced the concept of power spectral density (PSD) bounding, which we used for tropospheric error modeling Gallon et al. (2020). The PSD bounding method, unlike ACF bounding, is not restricted to fixed-interval implementations and is compatible with Kalman filtering.

In this work, we develop new, robust, sequential models for GNSS satellite orbit and clock errors using PSD bounding. Clock errors are analyzed per satellite clock type: Rubidium versus Cesium for GPS, and Rubidium versus passive Hydrogen masers for Galileo.

The paper is organized as follows. In Section 2, we describe the experimental data used throughout the paper, which includes three full years of orbit and clock errors for both GPS and Galileo. In Section 3, we analyze the errors’ stationarity, which is required for PSD bounding. We derive robust GPS and Galileo orbit and clock error models over time in Section 4. Concluding remarks are given in Section 5.

2 ORBIT AND CLOCK ERROR CHARACTERIZATION OVER TIME

In this section, we describe the databases and processes used to compute GPS and Galileo satellite orbit and clock errors. Orbit and clock estimates from broadcast navigation data are compared to a reference source (truth) after interpolation and coordinate transformation. In addition, we partition the data with respect to satellite clock type. 36 months of data were processed from 2018 through 2020.

2.1 Reference Orbit and Clocks

Truth data is obtained from the Multi-GNSS EXperiment (MGEX) repositories accessible at the National Aeronautics and Space Administration (NASA CDDIS, 2020). The MGEX service was initiated by the IGS to create a single GNSS data service for multiple core constellations. MGEX is comprised of several analysis centers (ACs) that independently compute their own GNSS orbit and clock products. For this work, we used precise orbit and clock data from two ACs: CODE for GPS and CNES for Galileo, both available at the IGS (2020) repository, and we consider them to be our truth reference. These reference products have an accuracy of 2.5 centimeters (specified as the one-dimensional mean root-mean-square (RMS) value over the three geocentric position components). A more detailed analysis of the consistency and accuracy of the MGEX orbit and clock products, as well as references for the GPS and Galileo reference files, can be found in Steigenberger and Montenbruck (2020).

Because we are interested in accurately characterizing orbit and clock errors over time periods of a few hours, the files’ nominal 15-minute sampling period is insufficient. The data was, therefore, interpolated at 30-second intervals. For this purpose, we used an 8th-order Lagrange polynomial following the analysis in Schenewerk (2003). Clock errors are random walk processes and should not be interpolated. Instead, we directly use the IGS reference clock products, which are provided with RMS errors lower than 70 picoseconds at 30-second sample intervals. Note that each clock product is aligned to their AC’s realization of the system time. This means that each AC product will be impacted by its own clock bias. Therefore, a common clock bias for the AC product in use must be removed from all satellites.

2.2 Broadcast Ephemerides

Broadcast ephemerides are stored in Receiver Independent Exchange (RINEX) formats that contain 24 hours of navigation messages. We used Stanford University’s sugl files for GPS satellites and brdc from CNES for Galileo, respectively, obtained from the Stanford University (2020) and CNES (2019) repositories. These institutions were chosen among several others because their data cleaning and validation algorithms ensure a minimal amount of residual file recording, storing, and labeling errors.

2.3 GPS and Galileo Orbit and Clock Errors

Satellite orbit and clock errors are obtained by taking the difference between the estimates of the broadcast ephemerides and those of the reference orbits, as shown in Figure 1. Reference orbits are provided with respect to the center of mass (CoM) of the satellite, whereas broadcast ephemerides are recorded with respect to the satellite’s antenna phase center (APC). Thus, they need to be converted to the same reference point—in this case, the reference orbits are expressed at the APC. Note that there are two conventions to define the APC: The IGS convention and the satellite metadata.

FIGURE 1
FIGURE 1

Error generation diagram

The offset for these conventions can be found at the International GNSS Service (IGS, 2020), National Oceanic and Atmospheric Administration (NOAA, 2021), and the European GNSS Service Centre (2021). More details on these offset corrections can be found in Montenbruck et al. (2014). After correcting for the offsets, orbit and clock errors were obtained by differencing reference and broadcast orbit and clocks. The final errors were then converted to the satellite-referenced local-level radial, along-track, and cross-track frame. Note that each clock product was aligned to their AC’s realization of the system time. This means that each AC product would be impacted by its own clock bias. Therefore a common clock bias must be removed from all satellites.

The user ranging error (URE) is the projection of the satellite’s orbit and clock errors for any user on Earth and is the combination of three components: Embedded Image 1

where x, w, and v are the radial+clock, along-track, and cross-track components, respectively. The weighting factors cv and cw are described below with Embedded Image.

We consider three scenarios. Given Re, the Earth radius in kilometers, and RSV, the satellite’s altitude in kilometers, the following limiting cases respectively maximize the contributions of radial+clock, along-track, and cross-track errors to URE:

  • Case 1: cv = cw = 0

  • Case 2: cv = 0 and Embedded Image

  • Case 3: Embedded Image and cw = 0

These cases are represented in Figure 2.

FIGURE 2
FIGURE 2

Projection of radial, along-track, and cross-track errors onto the user’s line of sight

We processed the orbit and clock errors of GPS (31 active satellites) and Galileo (18 active satellites) over 2018, 2019, and 2020, at a 30-sec sampling rate, which represents about 180 million data points. Figure 3 shows the orbit and clock errors of the 31 GPS (grey) and 18 Galileo (black) satellites in December 2018. Orbit and clock error time series for the entire month, for all satellites in Case 1, are plotted on top of each other.

FIGURE 3
FIGURE 3

Orbit and clock errors for all GPS and Galileo satellites in all three cases at 30-s intervals in December 2018

GPS ephemerides are nominally broadcast every 2 hours. When a new set of ephemerides is received, the GPS Interface Specifications document stipulates that the previous one is still valid (U.S. Air Force, 2020), but most users choose to use the most recently received set. The switch of ephemerides causes an abrupt change in the estimated satellite positions, which can be cumbersome to model in a KF. In response, we adopted an ephemeris switch-over method more amenable to dynamic modeling. We interpolated the broadcast ephemerides in the position domain over the 2-hour overlap between the previous and current sets to ensure continuous transitions at ephemeris switch-overs.

Figure 4 shows an example of radial error for GPS PRN07 interpolated (grey) versus not interpolated (red). The ephemeris jump can be observed for the non-interpolated case at the 2-hour time tag. The process used for Galileo ephemerides is similar, and both are further explained in Appendix A.

FIGURE 4
FIGURE 4

Example ephemerides interpolation (radial errors) of GPS PRN07 in December 2018 (upper), as well as its impact on radial+clock error ACF (lower)

Although the interpolated switch-over method does produce much smoother error transitions and is more amenable to modeling in a KF, it does not lead to error spectral characteristics significantly different from those resulting from abrupt switch-overs. The lower part of Figure 4 shows the (normalized) ACF for interpolated (grey) and not interpolated (red) orbit as well as the clock errors of GPS PRN07. The x-axis was limited here to 7 hours since it has the longest satellite pass duration (longer correlation times do not impact the user). The choice of interpolation method has minimal influence on the ACFs and, therefore, will not limit the general applicability of the final error models.

The error models developed in the rest of this paper are intended for implementations in which the ephemerides are interpolated. However, these error models also apply to applications in which the difference between interpolated versus non-interpolated ephemerides is negligible, even if the user does not interpolate ephemerides. For residual error modeling in space-based and ground-based augmentation systems (SBAS and GBAS, respectively) that rely on corrections of the un-interpolated ephemerides, this methodology could be reiterated with un-interpolated ephemerides in order to obtain appropriate models.

2.4 Impact of Satellite Clock Type on Orbit and Clock Errors

There are three main types of space-qualified atomic clocks used in GPS and Galileo: Rubidium or Rubidium Atomic Frequency Standard (Rb or RAFS), Cesium (Cs), or passive Hydrogen masers (PHMs). GPS satellites have been equipped with several combinations of clocks. GPS Block II/IIA carried two Cs and two Rb clocks, Blocks IIR and IIR-M contained three Rb clocks, and Block IIF carried two Rb and one Cs clock. Galileo satellites, on the other hand, use PHMs as their primary clocks and RAFS as secondary. Tables 1 and 2 summarize the GPS and Galileo clocks and block numbers associated with each PRN during the time period of the data we processed. Table 1 shows that most GPS satellites used a Rb clock as their main clock. Only two GPS satellites used Cs clocks. Table 2 shows that most Galileo satellites were PHM satellites and, even though all of them carried RAFS, only three used them for broadcast signals.

View this table:
TABLE 1

Clocks and Blocks of Each GPS Satellite (2018–2020)

View this table:
TABLE 2

Clocks and Blocks of Each Galileo Satellite (2018–2020)

A comparison in Teunissen and Montenbruck (2017) of atomic frequency standards was performed among the GPS and Galileo constellations over timescales ranging from 1 second to 1 day. It showed that stability could vary by as much as a factor of 10 and was generally better for the Rubidium and PHM clocks.

We can illustrate this observation with Figure 5, which shows example error time series at 30-second intervals for each satellite clock type in March 2018 for GPS (top-two charts) and Galileo (bottom-two charts). For GPS, the upper and lower plots show errors for the Rb and Cs clocks on PRN01 and PRN08, respectively. The Rb clock error oscillates within ±1 m. The Cs clock has larger error variations that reach up to ±2 m. For Galileo, both PHM and RAFS clock errors have similar behavior and remain between ±0.5 m. The variations in satellite clock errors are significantly larger in GPS than in Galileo satellites.

FIGURE 5
FIGURE 5

GPS and Galileo orbit and clock errors over December 2018

In order to model orbit and clock errors using PSD bounding, we must first ensure that they are stationary over the model’s duration.

3 ERROR STATIONARITY ANALYSIS

Consider a zero-mean random process X with autocorrelation function: Embedded Image 2

Perea Diaz (2019) established bounds on the variance, RX (0), of GNSS orbit and clock errors. By contrast, in our work, we are interested not only in RX (0), but in the entire ACF.

The following dilemma arises when dealing with sample ACFs. On the one hand, if we use too little data, the estimation of the expectation function in Equation (2) may be inaccurate. On the other hand, if we use long sets of data, the process may not be stationary over the entire data collection period.

To get some guidance on how much data is needed to get an accurate ACF estimate (and later, an accurate PSD estimate), let us analyze the example of a zero-mean first-order Gauss-Markov random process (FOGMRP) with an ACF expressed as: Embedded Image 3

where:

  • τ: the time constant of the FOGMRP

  • Embedded Image: the variance of the random process X

The variance of the sample ACF Embedded Image is given by: Embedded Image 4

(where t is the lag time, τ is the process time constant, and T is the length of data used in the estimation of the ACF). This result is derived in Appendix B.

Figure 6 presents the standard deviations of the ACF estimates for an example FOGMRP in Equation (4) with τ = 6 h and σX = 1.5 m. The curves are shown for various lengths of data T (color-coded curves) as a function of lag time (x-axis). Let us assume, for now, that this FOGMRP can be used as a rough approximation of GNSS satellite orbit and clock errors. The standard deviations of the ACF estimates (y-axis) are expressed in meters squared because the curves represent variations in orbit and clock error ACF estimation. For example, at lag time zero, the different curves capture the uncertainty in sample variance estimation error as a function of the length of data used. As T increases towards infinity, Embedded Image decreases towards zero: The longer the data set, the more accurate the ACF estimate.

FIGURE 6
FIGURE 6

Upper bound on standard deviation of ACF estimate (τ = 6 h, σ = 1.5 m)

Using 14 days of data, the standard deviation of the ACF estimate is close to 0.3 m2 at a lag time of 20 hours, whereas using 1 year of data, the standard deviation is about 0.06 m2 at the same lag. Figure 6 shows that the longer it is, the lower the uncertainty becomes. However, there appears to be diminishing benefit in using data lengths much longer than 3–6 months.

To test for the stationarity of the data, we use a combination of the Levene test and the two-sample Kolmogorov-Smirnov test. The Levene test (Levene, 1960) compares two or more sample populations to determine whether they have equal variance (homoscedastic). The two-sample Kolmogorov-Smirnov test (Massey Jr., 1951) determines whether two samples come from the same distribution. Both tests are performed with a 95% confidence level (i.e., p-value of 0.05). If both tests come back positive, the data is considered stationary.

However, both tests assume that the samples are independent. This is not the case for the actual orbit and clock error data. To approximate the effective number of independent samples, we use the properties of a FOGMRP. While the sample data processes are obviously not known to be FOGMRP a priori, the data will later verify that it is a reasonable approximation. Two samples of a FOGMRP with time constant τ can be considered independent if they are separated by a period larger than or equal to 2τ (see Gallon et al. [2021]). Therefore, to test stationarity, the data is re-sampled at regular 2τ intervals, where τ is the estimated FOGMP time constant of the data set.

Note that, for a FOGMRP, the lag time associated with an ACF estimate value of 1 / e is the estimated time constant τ used in the re-sampling mentioned above. If after final analysis of the data, the FOGMRP approximation turned out to be inappropriate, a similar approach to the one described in Gallon et al. (2021) would need to be performed to determine the number of independent samples for the specific process type.

For each PRN, the orbit and clock errors of a given satellite over 36 months of data are tested for stationarity. If the data set is deemed non-stationary, the data is subdivided into stationary data sets. Figure 7 shows the statistics of each stationary sub-data set for 2018 only (for clarity purposes). In some cases, the x-axis shows repeating PRNs. GPS PRN01, for example, is present six times. This means that stationarity tests failed until the data set was split into six different sets. Once stationarity is asserted, the partitioned data sets are treated separately.

FIGURE 7
FIGURE 7

Statistics of the 2018 GPS/GAL orbit and clock errors

4 ROBUST MODELING OF ORBIT AND CLOCK ERRORS

Unlike prior work that provided bounds on the variance of orbit and clock errors for GPS and Galileo (Perea Diaz et al., 2020), we present an approach to modeling these errors over time.

4.1 Zero-Mean Assumption

Figure 7 shows box plots of the error data for each satellite in the GPS (upper) and Galileo (lower) constellations. For clarity in exposition, the figure shows results limited to 2018; the rest of this work uses data from 2018, 2019, and 2020. The x-axis indicates the satellite’s PRN number. The color code designates the length of data used to generate the box representation as determined using the stationarity test. The red line inside each box represents the sample median and the upper and lower limits of the box represent the 75th and 25th percentiles, respectively. The dotted lines reaching away from the boxes represent the lowest and highest data points, excluding the outliers, which are represented by colored dots. A point is considered to be an outlier if it is greater than q3 + 2.7σ(q3 – q1) or smaller than q1 – 2.7σ(q3 – q1), where q1 and q3 are the 25th and 75th percentiles of the sample data. Note that ±2.7σ corresponds to 99% of the data if it is normally distributed. For both GPS and Galileo, the box plots appear to be consistent with a zero-mean assumption but will be verified below. It is worth noticing that GPS PRNs 8 and 24, the two Cesium satellites, have much larger error spreads than the rest of the GPS satellites. In addition, we can obtain the ensemble mean Embedded Image over n stationary data sets for a given constellation.

For each stationary data set xi, we can take independent samples every 2τi, and the estimates of the mean Embedded Image and variance Embedded Image for the data set can, respectively, be expressed as: Embedded Image 5 Embedded Image 6

where Ni is the number of independent samples in stationary set xi and xi,k is the k-th independent sample within xi. The variance of the error on the mean estimate Embedded Image can be written as (Bendat & Piersol, 2010): Embedded Image 7

The weighted least-squares estimate of the ensemble mean for samples within a constellation is then computed as: Embedded Image 8

Using orbit and clock error data from 2018 to 2019 (more data points for a more accurate estimation), we obtained mean estimates per constellation of Embedded Image and Embedded Image. It is important to remember that the IGS reference files were provided with an accuracy of approximately 2.5 cm. Therefore, the mean estimates obtained here are negligible. These results are consistent with those independently obtained in Perea Diaz (2019), which show that GPS and Galileo radial+clock errors for individual satellites are zero mean over one year. Therefore, we model orbit and clock error as a zero-mean process.

It is worth pointing out that we do not assert that the errors are zero-mean over the duration of a satellite pass at a given location. Often they are not, because their correlation time is typically of length similar to that of the longest satellite pass. It is the underlying parent process that is zero-mean. This is a necessary condition for evaluating PSDs in the next section.

4.2 Power Spectral Density Bounding

When it comes to estimating the PSD of stationary data, several methods exist (Bendat & Piersol, 2010). Our approach was to take the Discrete Fourier Transform (DFT) of the stationary sample ACFs. Langel et al. (2020) shows that error ACF values at time lags exceeding the duration of the Kalman filter’s operation are not relevant and can, therefore, be set to any value (zero, for example). This claim is only strictly true for zero-mean processes, which is the case in this application. To estimate the orbit and clock error PSDs, we use data collected over the entire years of 2018, 2019, and 2020, broken up into stationary segments.

We implemented the PSD estimation algorithm used in Langel et al. (2020). In order to limit spectral leakage, a tapered rectangular window (later called a tapering window) was applied to the orbit and clock error ACFs prior to performing the DFT. The red curve in Figure 8 represents the tapering window. The two other curves show an example ACF with windowing (black) and without windowing (grey). The sharp cutoffs in a simpler rectangular window would generate spectral leakage that would impact the estimated PSD (Langel et al., 2020). Instead, the tapering window smooths out the edges of an initial rectangular window to reduce spectral leakage. The window function is flat over the time interval of interest—in this case, the longest expected satellite pass duration, T1. ACF values associated with lags larger than T2 are set to zero. The farther apart T1 and T2 are from each other, the less spectral leakage is observed in the estimated PSD. Because the longest satellite pass lasts 7 hours, we used 7 and 14 hours, respectively, for T1 and T2. A sensitivity analysis on these two parameters is provided in Appendix C.

FIGURE 8
FIGURE 8

Tapering window and its impact on an ACF

The left-hand side plot in Figure 9 shows the estimated orbit and clock error PSD curves for GPS satellites over the years 2018 to 2020 and for each of the three line-of-sight limit cases selected in Section 2.3. The Cs satellites are represented with blue curves and the Rb satellites are represented with green ones. The two clock types can again be clearly distinguished with Cs curves above the others. The fact that Rb clock curves are lower means that the standard deviation of their errors is also lower, which matches our previous observations as well as those in Perea Diaz (2019).

FIGURE 9
FIGURE 9

PSDs of orbit and clock errors for (a) GPS and (b) Galileo

The right-hand side plot in Figure 9 shows the same curves for Galileo satellites. As observed before, no clear distinction can be drawn between the PHM and the RAFS Galileo satellites.

To robustly model the dynamics of orbit and clock errors over time, we upper-bound their PSD using a FOGMRP model. We chose a FOGMRP that is fully determined by two parameters, a time constant τb and a standard deviation Embedded Image, because it can easily be incorporated into a KF. Its PSD can be expressed as: Embedded Image 9

Since the two GPS clock types clearly show different trends in the frequency domain, the user may implement different bounding models for each. Figure 10(a) shows the PSD curves for the Rb satellites and its FOGMRP bound. Figure 10(c) shows the same for the Cs satellites. Figures 10(b) and 10(d), respectively, show the bounds for the Galileo satellites RAFS and PHM. If future ARAIM integrity support messages (ISMs) do not identify clock types by satellite, users could access this information via Notice Advisories to Navigation Users (NANUs). However, we only processed data for two GPS Cs satellites and three Galileo RAFS satellites. The bounds for satellites with these two clock types could benefit from validation using more data, either from more satellites or over a longer duration.

FIGURE 10
FIGURE 10

PSD bounding of (a) GPS Rb, (b) GAL RAFS, (c) GPS Cs, and (d) GAL PHM

Table 3 summarizes the parameters of the PSD-bounding Gauss-Markov random process (GMRP) models. For the GPS constellation, the PSDs of orbit and clock errors were bounded using τb = 5 hours and σb = 1.8 meters. Galileo errors were bounded using τb = 2 hours and σb = 0.65 meters.

View this table:
TABLE 3

PSD Bounding Summary: FOGMRP Parameters

These models can be readily implemented in a KF by state augmentation (Bryson, 2002).

5 CONCLUSION

In this paper, we developed new stochastic models to bound the dynamics of global navigation satellite systems (GNSSs) orbit and clock errors. The models are needed to ensure integrity in time-sequential GNSS navigation systems integrated with inertial sensors or using dual-frequency, multi-constellation, sequential Advanced Receiver Autonomous Integrity Monitoring (ARAIM).

We processed and partitioned three years of empirical GPS and Galileo orbit and clock error data. We performed stationarity analyses using two statistical tests: one on the variance of the errors, and another on their distribution.

Then, to robustly model orbit and clock errors, we used a frequency-domain bounding approach, in which the error PSDs were upper-bounded using first-order Gauss-Markov random process (GMRP) models. We derived separate high-integrity, PSD-bounding models for GPS and Galileo satellites and for their different clock types.

HOW TO CITE THIS ARTICLE

Gallon, E., Joerger, M., & Pervan, B. (2022) Robust modeling of GNSS orbit and clock error dynamics. NAVIGATION, 69(4). https://doi.org/10.33012/navi.539

ACKNOWLEDGMENTS

The authors would like to thank the Federal Aviation Administration (FAA) for their support of this research. However, the opinions in this paper are our own and do not necessarily represent those of any other person or organization.

APPENDIX A EPHEMERIS INTERPOLATION

GPS satellites typically broadcast a new set of ephemerides every 2 hours. Each set is valid for at least 4 hours. Hence, when a new set of ephemerides is received, the previous one is still valid for another 2 hours. Section 2.3 shows that ephemeris updates introduce jumps in the orbit and clock errors. If we were to account for the jumps in a batch or KF implementation, we would need to model their dynamic behavior, which would require a cumbersome hybrid continuous/discrete process model. Since smoother dynamics are more convenient to model, and since we would ideally, for ease of implementation, like to model errors using simple models (e.g., as FOGMRP), we instead use interpolated broadcast ephemerides.

A.1 GPS Ephemeris Interpolation

GPS ephemerides are valid for (at least) 4 hours—2 hours before and after the specified time of ephemeris (toe). A new set of ephemerides is normally received by the receiver every 2 hours. When the new set is received, the current one is still valid for some time (2 hours after toe [U.S. Air Force, 2020]). Users who are not concerned with stochastically modeling ephemeris error would typically decide to use the new one immediately. However, our goal was to produce such a model, so we instead interpolated the current (Xcur) and next (Xnext) sets of ephemerides in the position and velocity domains over a window Wint of 2 hours. This interpolation is represented in Figure A2, where the blue curves are the original, non-interpolated ephemerides and the green curve represents the interpolated portion over the time in which both ephemerides were valid. The equation describing this interpolation is: Embedded Image A1

where:

  • Embedded Image: the interpolated output

  • Wint: the interpolation window

  • t: the time of the jump between Xcur and Xnext

  • τ: a dummy variable within the interpolation window: τ ∈ [0: Wint ]

FIGURE A2
FIGURE A2

GPS and Galileo ephemerides interpolation diagram

Figure A1 shows the radial errors of GPS PRN07 in December 2018. The focus is on the ephemeris jump occuring at t = 2 hours. The blue dashed curve represents the set of ephemeris currently in use and whose toe was t = 0 hours. The black dashed curve represents the set of ephemeris that was just received by the user, whose toe was at t = 2 hours, but whose validity window starts at t = 0 hours. The green curve represents the result of their interpolation between t = 0 and 2 hours. At t = 0 hours, the user is currently using the blue set but has just received the black one as well; the interpolation begins at this point. At t = 1 hour, both sets are equally weighted to generate the interpolated curve (green). At t = 2 hours, the current set time of applicability has expired, and the user is now using 100% of the new set.

FIGURE A1
FIGURE A1

GPS PRN07 ephemerides interpolation

A.2 Galileo Ephemeris Interpolation

Galileo ephemerides are valid for 4 hours after the toe which, in this case, also represents the time of reception of the new ephemeris. Therefore, the interpolation process is slightly different than for GPS. It is expressed with the following formula: Embedded Image A2

where the interpolation window is now defined (in seconds) as: Embedded Image A3

Note that Galileo’s interpolation window is different to cope with the fact that the broadcast rate isn’t necessarily fixed and can be as short as 10 min. Additionally, unlike GPS, the toe is located at the beginning of validity period (see Galileo OS-SDD [2021]).

APPENDIX B VARIANCE OF AUTOCORRELATION ESTIMATES

Let us define the autocorrelation estimate of a random process X(t) as: Embedded Image B1

where RX(t) is known and δRX(t) is unknown. In the following, we derive an approximate expression for the variance of autocorrelation estimate Embedded Image.

Equation (8.103) of Bendat and Piersol (2010) show that the variance on the ACF estimate Embedded Image can be expressed as: Embedded Image B2

Approximating that the errors are derived from a zero-mean FOGMRP, the parent autocorrelation function can be expressed as: Embedded Image B3

Substituting Equation (B3) into Equation (B2) and solving the integrals, we obtain: Embedded Image B4

APPENDIX C SENSITIVITY OF PSD ESTIMATION TO PARAMETER SELECTION

The PSD estimation method used in this paper is described in Langel et al. (2020) and relies on a tapering window applied to the ACF of the errors prior to the Fourier transform computation. The tapering window, as described in Section 4.2, relies on two parameters: T1 and T2. On top of that, the length of data T plays an important role in the PSD estimation process since it drives the ACF accuracy. Hence, when estimating PSDs, we need to take into account all three time parameters: T1, T2, and T. In this appendix, we analyze the sensitivity of the PSD estimation process to our choices of T1, T2, and T.

In Section 3, we showed that the parameter T influences the overall accuracy of the ACF estimate. Since the PSD estimation approach chosen here relies on first estimating error ACFs, the same conclusion applies: The longer the data, the more accurate the ACFs and, hence, the more accurate the PSD estimate will be. Therefore, using as much data as possible would be ideal (i.e., 1 year). Note, however, that the choice of T will be entirely determined by the stationarity tests. For this analysis, we assume an average value T = 6 months.

Additionally, we know that the T1 parameter represents the lag values until which the ACF will remain unchanged. In our case, because a satellite pass lasts about 7 hours maximum, it is in our best interest not to modify correlation values during that period of time. We, therefore, chose T1 = 7 hours.

The only parameter left to analyze is T2, which will define how rectangular the tapering window will be. If T1 and T2 are close to each other, the window will be nearly rectangular, and it will result in spectral leakage in the frequency domain, decreasing the quality of our PSD estimate. Example windows are represented in Figure C1. The following section analyzes the impact of different T2 values on the PSD estimate of a FOGMRP process.

FIGURE C1
FIGURE C1

Example of tapering windows used in the PSD estimation process

C.1 Impact of the Window Shape on a FOGMRP PSD Estimate

Let us begin by looking at the impact of various tapering windows on the PSD estimate of a first-order GMRP. For that, we generated 6 months of FOGMRP data at a 30-sec sampling rate with a time constant of 6 hours and a standard deviation of 1.5 m. Because we know the theoretical expression of a FOGMRP PSD, we can compare it to the various PSD estimates. Figure C2 shows the PSD estimates using the various tapering windows represented in Figure C1.

FIGURE C2
FIGURE C2

Impact of tapering window on PSD estimate of FOGMRP in terms of spectral leakage in (a) original figure and (b) the zoomed figure

The red dashed curve represents the theoretical PSD curve of a FOGMRP. The grey curve represents the PSD estimate obtained when using a rectangular window (T1 = T2 = 7 h). Rectangular windows are known for inducing spectral leakage in the frequency domain. At high frequencies (right part of the plot), spectral leakage induces a divergence of the PSD estimate from the true PSD curve. We can see that the more T2 is increased, less spectral leakage is visible (bumpiness of the curves) and the closer the PSD estimates get to the true GMRP PSD curve.

These results suggest that PSD estimates with increasing T2 values will converge to the true GMRP PSD curve. To test this theory, Figure C3 shows PSD estimates using tapering windows whose T2 values range from 24 hours to 1 month of data. Beyond a certain length of time for T2, the curves can be seen to become noisier at all frequencies. Indeed, for values higher than 48 hours, the spectral leakage seems to be replaced by noise, most likely due to noise in the ACF estimate.

FIGURE C3
FIGURE C3

Impact of a tapering window on PSD estimate of FOGMRP in terms of data noise in (a) the original figure and (b) a zoomed-in version

The results of this section suggest that T2 = 48 h is an ideal value for the PSD estimation of a FOGMRP.

C.2 Impact of the Window Shape on Real Orbit and Clock Error Data

In this section, we verify that this theory applies to the orbit and clock errors using an example data set of PRN 01 for 2018.

The grey curve of Figure C4(b) shows the PSD estimate of this data set with T1 = 7 h and T2 = 48 h, as suggested by the analysis in Section C.1. If we were to upper bound this curve with a FOGMRP (see Section 4.2), the associated sigma would have to be inflated to upper bound the frequency bump observed at f = 2.4 × 10−5 Hz. That frequency is equivalent in the time domain to a period of 11.5 hours.

FIGURE C4
FIGURE C4

(a) Tapered ACF and (b) PSD estimates for GPS PRN 01 orbit+clock errors

Note that GPS satellites have an orbit period of half a sidereal day (23 hours, 56 minutes, and 4 seconds) and will, therefore, take about 11.9 hours to orbit the Earth. Note also that this orbital period can be observed on the grey ACF curve in Figure C4(b). The grey curve represents the tapered ACF of the data set, with T1 = 7 h and T2 = 48 h. Therefore, the ACF lobe located at the 11.9-hour time lag is barely tapered and will result, in the PSD domain, in a frequency bump that will drive the PSD bounding process.

Since a satellite pass lasts 7 hours at the maximum, tapering the ACFs after this time lag will not impact the robustness of our model. Smaller T2 values will result in an attenuated orbital period frequency bump and, therefore, a better, less conservative PSD bounding model.

Any 14h ⩽ T2 ⩽ 48h would ensure limited spectral leakage and noise. To reduce the impact of the orbital period on our PSD estimate, we chose T2 to be as small as possible within this interval.

Therefore, in this work, we set T1 = 7 h and T2 = 14 h.

APPENDIX D GPS AND GALILEO ORBIT AND CLOCK MEAN ERRORS

The tables provided in this section support the results presented in Section 4.1.

View this table:
TABLE D1

GPS Radial+Clock Means Over Stationary Data Sets From 2018

View this table:
TABLE D2

GAL Radial+Clock Means Over Stationary Data Sets from 2018

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.

REFERENCES

  1. ↵
    1. Bendat, J. S., &
    2. Piersol, A. G.
    (2010). Random data: Analysis and measurement procedures. John Wiley & Sons.
  2. ↵
    1. Bryson, A. E.
    (2002). Applied linear optimal control: Examples and algorithms (1st ed.). Cambridge University Press.
  3. ↵
    1. Centre National d’Etudes Spatiales (CNES)
    . (2019). Broadcast ephemeris repository. ftp://serenad-pub-lic.cnes.fr/SERENAD0/FROMNTMFV2/NAV/
  4. ↵
    1. DeCleene, B.
    (2000). Defining pseudorange integrity – overbounding. Proc. of the 13th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 2000), Salt Lake City, UT, 1916–1924. https://www.ion.org/publications/abstract.cfm?articleID=1603
  5. ↵
    1. European GNSS Service Centre
    . (2021). Galileo satellite metadata. European Union Agency for the Space Programme. https://www.gsc-europa.eu/support-to-developers/galileo-satellite-metadata
  6. ↵
    1. Galileo OS-SDD
    . (2021). Open service: Service definition document v1.2. (Technical Report). European GNSS Service Centre. https://www.gsc-europa.eu/sites/default/files/sites/all/files/Galileo-OS-SDD_v1.2.pdf
  7. ↵
    1. Gallon, E.,
    2. Joerger, M.,
    3. Perea, S., &
    4. Pervan, B.
    (2019). Error model development for ARAIM exploiting satellite motion. Proc. of the 32nd International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2019), Miami, FL, 3162–3174. https://doi.org/10.33012/2019.17019
  8. ↵
    1. Gallon, E.,
    2. Joerger, M., &
    3. Pervan, B.
    (2020). Robust modeling of tropospheric delay dynamics for sequential positioning. 2020 IEEE/ION Position, Location and Navigation Symposium (PLANS), Portland, OR. https://doi.org/10.1109/PLANS46316.2020.9109969
  9. ↵
    1. Gallon, E.,
    2. Joerger, M., &
    3. Pervan, B.
    (2021). Robust modeling of GNSS tropospheric delay dynamics. IEEE Transactions on Aerospace and Electronic Systems, 57(5), 2992–3003. https://doi.org/10.1109/TAES.2021.3068441
  10. ↵
    1. International GNSS Service (IGS)
    . (2020). International GNSS Service Antex files [IGS repository]. International GNSS Service File Access. https://files.igs.org/pub/station/general/pcv_archive/
  11. ↵
    1. Joerger, M., &
    2. Pervan, B.
    (2020). Multi-constellation ARAIM exploiting satellite motion. NAVIGATION, 67(2), 235–253. https://doi.org/10.1002/navi.334
  12. ↵
    1. Langel, S. E.
    (2014). Bounding estimation integrity risk for linear systems with structure stochastic modeling uncertainty [Doctoral dissertation, Illinois Institute of Technology].
  13. ↵
    1. Langel, S.,
    2. García Crespillo, O., &
    3. Joerger, M.
    (2020). A new approach for modeling correlated Gaussian errors using frequency domain overbounding. 2020 IEEE/ION Position, Location and Navigation Symposium (PLANS), Portland, OR. https://doi.org/10.1109/PLANS46316.2020.9110192
  14. ↵
    1. Levene, H.
    (1960). Robust tests for equality of variances. In I. Olkin (Ed.), Contributions to probability and statistics: Essays in honor of Harold Hotelling (pp. 278–292). Stanford University Press.
  15. ↵
    1. Massey Jr., F. J.
    (1951). The Kolmogorov-Smirnov test for goodness of fit. Journal of the American Statistical Association, 46(253), 68–78. https://doi.org/10.2307/2280095
    CrossRef
  16. ↵
    1. Montenbruck, O.,
    2. Steigenberger, P., &
    3. Hauschild, A.
    (2014). Broadcast versus precise ephemerides: A multi-GNSS perspective. GPS Solutions, 19, 321–333. https://doi.org/10.1007/s10291-014-0390-8
  17. ↵
    1. NASA CDDIS
    . (2020). Global navigation satellite system (GNSS) IGS Multi-GNSS Experiment (MGEX) products. NASA Archive of Space Geodesy Data. https://doi.org/10.5067/GNSS/GNSS_IGSMGEXP_001
  18. ↵
    1. National Oceanic and Atmospheric Administration (NOAA)
    . (2021). Antenna offset solutions [Data repository]. https://geodesy.noaa.gov/ANTCAL/LoadFile?file=ngs14.atx
  19. ↵
    1. Perea Diaz, S.
    (2019). Design of an integrity support message for offline advanced RAIM [Doctoral dissertation, Rheinisch-Westfälische Technische Hochschule Aachen]. https://elib.dlr.de/133531/1/2019_Dis_Santiago%20Perea%20Diaz.pdf
  20. ↵
    1. Perea Diaz, S.,
    2. Meurer, M., &
    3. Pervan, B.
    (2020). Impact of sample correlation on SISRE overbound for ARAIM. NAVIGATION, 67(1), 197–212. https://doi.org/10.1002/navi.346
  21. ↵
    1. Rife, J.,
    2. Pullen, S.,
    3. Enge, P., &
    4. Pervan, B.
    (2006). Paired overbounding for nonideal LAAS and WAAS error distributions. IEEE Transactions on Aerospace and Electronic Systems, 42(4), 1386–1395. https://www.doi.org/10.1109/TAES.2006.314579
  22. ↵
    1. Schenewerk, M.
    (2003). A brief review of basic GPS orbit interpolation strategies. GPS Solutions, 6(4), 265–267. https://doi.org/10.1007/s10291-002-0036-0
  23. ↵
    1. Stanford University
    . (2020). Stanford University GPS Lab ephemeris files. Stanford Department of Engineering. https://gps.stanford.edu/resources/software-tools/sugl-ephemeris-files
  24. ↵
    1. Steigenberger, P., &
    2. Montenbruck, O.
    (2020). Consistency of MGEX orbit and clock products. Engineering, 6(8), 898–903. https://doi.org/10.1016/j.eng.2019.12.005
  25. ↵
    1. Tanil, C.,
    2. Khanafseh, S.,
    3. Joerger, M., &
    4. Pervan, B.
    (2018). An INS monitor to detect GNSS spoofers capable of tracking vehicle position. IEEE Transactions on Aerospace and Electronics, 54(1), 131–143. https://www.doi.org/10.1109/TAES.2017.2739924
  26. ↵
    1. Teunissen, P. J. G., &
    2. Montenbruck, O.
    (2017). Springer handbook of Global Navigation Satellite Systems. Springer International Publishing. https://www.doi.org/10.1007/978-3-319-42928-1
  27. ↵
    1. U.S. Air Force
    . (2020). NAVSTAR GPS space segment/navigation user segment interfaces (Technical Report No. IS-GPS-200). https://www.gps.gov/technical/icwg/IS-GPS-200L.pdf
  28. ↵
    1. Working Group C
    . (2012). ARAIM Technical Subgroup Interim Report (Technical Report). EU-U.S. Cooperation on Satellite Navigation. https://www.gps.gov/policy/cooperation/europe/2013/working-group-c/ARAIM-report-1.0.pdf
  29. ↵
    1. Working Group C
    . (2014). ARAIM Technical Subgroup Milestone 2 Report (Technical Report). EU-U.S. Cooperation on Satellite Navigation. https://www.gps.gov/policy/cooperation/europe/2015/working-group-c/ARAIM-milestone-2-report.pdf
  30. ↵
    1. Working Group C
    . (2016). ARAIM Technical Subgroup Milestone 3 Report (Technical Report). EU-U.S. Cooperation on Satellite Navigation. https://www.gps.gov/policy/cooperation/europe/2016/working-group-c/ARAIM-milestone-3-report.pdf
Back to top

In this issue

Print
Download PDF
Article Alerts
Email Article
Citation Tools
Robust Modeling of GNSS Orbit and Clock Error Dynamics
Elisa Gallon, Mathieu Joerger,, Boris Pervan
NAVIGATION: Journal of the Institute of Navigation Dec 2022, 69 (4) navi.539; DOI: 10.33012/navi.539
Twitter logo Facebook logo Mendeley logo
Bookmark this article

Jump to section

Keywords