Research ArticleOriginal Article
Open Access

A low complexity smoothing algorithm for improved GPS point solutions on board LEO spacecraft

Damon Van Buren, Penina Axelrad and Scott Palo
NAVIGATION: Journal of the Institute of Navigation March 2021, 68 (1) 185-198; DOI: https://doi.org/10.1002/navi.410

Abstract

This paper describes a novel approach for reducing the kinematic position and time errors on a GPS receiver in low Earth orbit. For dual-frequency GPS receivers, point solution errors are dominated by signal-in-space ranging errors, which introduce sharp discontinuities in the position and time estimates when GPS visibility changes. These discontinuities may be calculated and removed, greatly reducing the peaks of the errors, at the expense of introducing a slowly varying bias, which is estimated and removed. The paper gives a detailed description of the algorithm and illustrates its use on a GPS receiver using historical ephemerides. To demonstrate real-world efficacy, GRACE-FO observations are processed using the proposed approach. Real-time results show the reduction of the time error standard deviation by 16% to 1.9ns and the reduction of the 3D RMS error by 17% to 1.06m. Maximum real-time errors are reduced 26% and 32%, respectively. Delayed processing offers greater reductions of errors, near 40%.

1 INTRODUCTION

Advances in sensors and communication hardware have enabled small satellites to produce imagery and other data, which rival products from much larger spacecraft (Ali Shah & Arshad, 2013; Alminde & Laursen, 2009; Bryerton et al., 2016; Chattopadhyay, 2017; Kalemci et al., 2013; Norton et al., 2016; Peral et al., 2015; Segret et al., 2018). Cubesats have also been sent on deep space missions for the first time with the launch of the MarCo spacecraft to Mars in May 2018 (Decrossas et al., 2019; Hodges et al., 2017; Tkacenko et al., 2019). For small satellites to continue to expand into newer, more challenging roles, they must meet increasingly stringent system requirements.

Accurate position and time estimates are important for a variety of smallsat missions. Examples of missions, which require high performance timekeeping, are deep space missions using one-way ranging (Imken et al., 2017; Wallis et al., 2005; Weaver & Kantsiper, 2007), small satellite constellations functioning as sensor arrays for astronomy (Engelen et al., 2013; Rajan et al., 2015), low Earth orbit (LEO) constellations for Global Positioning System (GPS) augmentation and improved positioning (Ge et al., 2018; Reid et al., 2018), and Earth observation missions which require accurate position and time information for imagery and science data.

The GPS system provides an excellent solution for LEO position, velocity, and time (PVT) determination, eliminating the need for ground-based ranging for most LEO missions. Because the GPS signals are continuously available for LEO spacecraft equipped with an upward-looking antenna, continuous accurate geometric PVT estimates are possible. When combined with orbit determination (OD) techniques and differential GPS corrections in ground-based post-processing, orbital states may be estimated with errors as low as several centimeters (Montenbruck et al., 2005). GPS measurements have also become the de facto standard for real-time LEO OD with submeter orbit accuracies possible (Choi et al., 2009; Menzione et al., 2017; Montenbruck & Gill, 2006; Montenbruck & Ramos-Bosch, 2008; Reichert et al., 2002; Yang et al., 2016).

Kinematic OD uses only the GPS measurements to calculate instantaneous point solutions (Montenbruck, 2003) without relying on dynamic models. Because it relies entirely on the current geometry of the GPS constellation and the code and carrier phases, it is vulnerable to poor geometry and high measurement noise due to low signal-to-noise ratios.

Precise orbit determination (POD) uses force models and batch or Kalman filtering of observations to provide accurate estimates of the spacecraft’s orbital state. This approach produces a smooth estimate of the orbit, but the technique is vulnerable to errors and omissions in the dynamic model for factors such as atmospheric drag.

Reduced dynamic OD combines the advantages of the classical dynamic approach with the kinematic approach using GPS point solutions (Montenbruck et al., 2005; Montenbruck & Ramos-Bosch, 2008; Parkinson & Spilker, 1996) and balances the vulnerability of kinematic OD to weak GPS geometry and signal noise with the sensitivity of dynamic OD to model errors.

There are several challenges to reduced dynamic OD implementation on board small satellites. The process noise model must be tuned carefully to minimize sensitivity to noise in the observations, while still enabling the filter to track unmodeled forces effectively.

Differential GPS (DGPS) corrections are needed for applications requiring 3D RMS accuracies better than a decimeter as required for many Earth observation missions such as altimetry, gravimetry, SAR interferometry, and atmospheric sounding. Currently, there is no global differential service accessible by LEOs, although several regional differential services exist. NASA is planning the TDRSS Next Generation Broadcast Service (NGBS), which will broadcast differential corrections globally in S-Band from the TDRSS satellites, but eventual use of these signals will require a separate S-Band receiver (Heckler et al., 2016). For real-time applications on LEO spacecraft where DGPS corrections are not available, reduced dynamic accuracy is limited to 0.4-1.0m 3D RMS (Montenbruck & Ramos-Bosch, 2008).

Reduced dynamic algorithms also require significant processing capacity and so may exceed the processing capability of small satellites. Implementation of reduced dynamic OD requires precise numerical integration, UD factorization to improve numerical stability, and a number of other complex large matrix operations. Gravity models to achieve position accuracies of a few decimeters are between fiftieth and seventieth order (Xuewen et al., 2019; Montenbruck & Ramos-Bosch, 2008).

The proposed method is simpler to implement than the real-time orbit determination methods described above and requires less processing. It also does not require differential GPS corrections and does not use complex force models.

2 GPS POINT SOLUTIONS AND EFFECT OF GPS SYSTEM ERRORS

The GPS point solution (PS) provides an estimate of the user position and clock bias, which is optimal in a least squares sense, using a single epoch of observed pseudoranges to the GPS satellites along with GPS ephemerides and clock biases. For this research, the use of a dual-frequency receiver is assumed, allowing the generation of ionosphere-free pseudorange measurements. An example of a dual-frequency GPS receiver, which has been used on several small satellite missions, is the NovAtel OEM719 (Chang et al., 2018), which requires less than 2W of continuous power (NovAtel, n.d.). The dual-frequency pseudorange for each GPS satellite is modeled as

Embedded Image 1

where 𝑏 represents the user clock bias in meters, 𝐱(𝑘) is the Earth-centered, Earth-fixed (ECEF) position vector for the kth GPS satellite, 𝐱 is the LEO satellite position, and Embedded Image is the tracking error from the receiver combined with the GPS system errors.

Each pseudorange provides a nonlinear equation with four unknowns: the three position values and the user clock bias. A minimum of four pseudoranges are required to solve for position and time. For 𝐾 GPS satellites, the equations are solved by linearizing about the current user position estimate and iterating until the change in the user position estimate is below a predefined threshold.

The linearized pseudorange equation for the 𝑛𝑡ℎ iteration may be written compactly as

Embedded Image 2

where 𝛿𝛒𝑛 is the difference between the corrected measured pseudoranges and the pseudoranges calculated using the position and time estimates for the previous iteration, 𝛿𝐱𝑛 is the change in the estimated receiver position, 𝛿𝑏𝑛 is the change in the estimated user clock error, and Embedded Image represents the combination of the receiver measurement error, multipath error, and GPS system errors.

The geometry matrix 𝐆 is given by

Embedded Image 3

where the Embedded Image is the unit vector from the LEO spacecraft to each of the GPS satellites. The least squares solution for the overdetermined case where 𝐾 > 4 is given by

Embedded Image 4

This leads to improved estimates of the user position and clock bias:

Embedded Image 5

Embedded Image 6

After each iteration, the changes in the estimated position and time are compared to a threshold; iteration stops when the changes become sufficiently small. Typically only a few iterations are needed for the algorithm to converge on an accurate position.

2.1 Effect of GPS system errors and changes in the visible GPS constellation

The three largest sources of error in GPS-based LEO position estimates are 1) the ionosphere; 2) receiver errors caused by signal noise and multipath and carrier phase ambiguity; and 3) errors in the broadcast GPS ephemerides and clocks, also called signal-in-space ranging errors (SISRE). Dual-frequency measurements are used to eliminate the ionospheric errors. This paper describes a novel approach for reducing time and position errors caused by errors in carrier-smoothed pseudoranges and errors in broadcast GPS ephemeris and clock offsets (Misra & Enge, 2006).

Code tracking errors are caused by signal noise in the code tracking loop and result in relatively large errors (1-2m). Carrier phase observations are much more precise and have lower noise and so may be used to smooth the code measurements using a Hatch filter (Hatch, 1982). However, the initial error of the Hatch filter output may also be quite large, since initially very few pseudorange and carrier phase observations are averaged to estimate the carrier phase offset.

Real-time onboard solutions are also sensitive to errors in the broadcast GPS ephemerides and clock corrections (SISRE), which are typically on the order of 0.5m (Montenbruck et al., 2018). Differential techniques can estimate these errors using a network of ground receivers and then transmit the corrections to the user equipment via a secondary link. However, there is currently no global differential transmission service for LEO satellites.

In the absence of tracking errors and SISRE, changes in GPS visibility would not result in changes to the point solution, since any number of pseudoranges greater than three would allow calculation of the LEO satellite position and clock bias with very low errors. But since each pseudorange has different receiver tracking errors and SISRE, the point solution converges to a different value when a GPS satellite is added or removed to provide the best solution in the least-squares sense.

Figure 1 shows simulated user clock bias estimation errors caused by SISRE for the Gravity Recovery and Climate Experiment (GRACE) Follow-On (FO) (Kornfeld et al., 2019) B satellite for the last six hours of November 20, 2018. The two GRACE-FO spacecraft are separated by an average of 220km and were placed in a near circular polar orbit with an altitude of 490km. The satellites measure Earth’s gravity by precisely measuring changes in the distance between them using a microwave link, which is precise to a few microns.

FIGURE 1
FIGURE 1

Simulated SISRE time error for GRACE-FO orbit

The precise trajectory of the GRACE-FO satellite is used as the true LEO orbit, and the National Geospatial-Intelligence Agency (NGA) Office of Geomatics (NGA, n.d.) post-processed GPS ephemerides and clock biases are used as the true GPS positions and times. Ideal pseudorange measurements for GPS satellites visible to GRACE-FO are constructed from these two sources, and point solutions are calculated using the ideal pseudoranges and the actual broadcast ephemerides and clock offsets for this date and time.

The figure shows that the SISRE-induced errors are characterized by sharp discontinuities (“jumps”) connected by sloped segments, as GPS satellites are frequently added and dropped from the point solution over the LEO orbits. We note that some errors exceed 14 ns, even in the absence of measurement errors due to noise or multipath.

When the set of GPS satellites, which are included in the point solution, changes, the point solution calculation converges on a different optimum least squares solution, causing the error to jump to a different value. Updates to the broadcast ephemerides also cause jumps, but these only occur at two-hour intervals for the legacy navigation message, seen at hours 20, 22, and 24 in this plot. Between the jumps, the errors change much more gradually, as the broadcast ephemeris and clock errors slowly evolve and the GPS satellite geometry changes. Over the long term, the errors are nearly unbiased as one would expect given the nature of the GPS system (Perea et al., 2017), but over this brief 5-hour dataset short-term offsets are present, such as the negative bias after the 21st hour and the positive bias centered near the 23rd hour.

Several of the larger, more conspicuous jumps are labeled for clarity, but there are over 100 jumps in the figure, so it is not possible to label most of them. The largest jump occurs at 20.85 hours and has a magnitude of approximately 10ns. This jump occurred when the number of visible GPS satellites dropped from 8 to 7 and TDOP increased from 1.2 to 1.9. The SISRE of the GPS satellite, which is dropped, is 0.3m prior to being dropped.

3 PROPOSED APPROACH

Obviously, smoothing or otherwise addressing the SISRE-induced errors would be very beneficial. Unfortunately, the satellite-specific SISRE-induced position and time errors are never directly observable on board. The errors couple through the observation geometry into the LEO point solutions for position and user clock bias, and it is impossible to separate the SISRE from the user clock bias estimates and position estimates.

Although the SISRE cannot be measured, the values of the jumps they produce in the point solution can be directly calculated by differencing the point solution with and without a GPS signal which is newly added or dropped. Because the jumps cause the largest SISRE induced errors, we would expect the error magnitude to be dramatically reduced with the jumps removed. The ability to directly calculate the jumps in the point solution offers a unique opportunity to address the errors caused by SISRE, provided an approach which exploits this information can be developed.

The chosen approach is to subtract the accumulated jump values from the position and time estimates. This eliminates the sudden changes in the estimate and, at first glance, would appear to dramatically reduce the errors in the point solution.

For receivers using carrier smoothing, initial errors in the Hatch filter output have an effect similar to SISRE, since they cause a jump in the point solution when a GPS satellite is added or dropped. For this reason, it is assumed that a solution which addresses SISRE will also address receiver tracking errors, since biases caused by the carrier phase ambiguity have the same effect as errors in the broadcast GPS satellite clock offsets.

It is important to note that because the initial SISRE is unknown, removal of the SISRE-induced jumps will result in a fixed bias equal to the initial SISRE. This bias may be estimated by tracking the difference between the smoothed estimates and the original estimates using a phase locked loop (PLL) with a large time constant on the order of several days.

3.1 Calculation and removal of jumps

The process for calculating the jumps is shown in Figure 2. The accumulated jumps are initialized to zero. When a new sample is received, the standard point solution is calculated using the currently visible GPS satellites and the most recent broadcast ephemerides. If a GPS satellite has been added since the last sample point, or if the ephemerides have been updated, the point solution is calculated without the added satellites and/or using the previous ephemerides. The geometry matrix will have fewer rows to account for the subset of measured pseudoranges without the new one. The jump due to the GPS satellite addition or ephemeris update is found by subtracting the point solutions without the added observations or with the old ephemerides from the current point solution time and position values. The subscript “wa” indicates that the calculation is performed “without adds”, meaning without observations from GPS satellites which have just been introduced on the current sample, while the subscript “full” indicates that the full complement of visible GPS satellites is used.

Embedded Image 7

Embedded Image 8

Embedded Image 9

FIGURE 2
FIGURE 2

Jump calculation flow diagram

If one or more GPS satellites have been dropped since the previous sample time, the jump due to the dropped satellites must be computed at the previous sample time, since the dropped satellite is no longer visible. The point solution is calculated without the dropped satellites at the previous sample time, and the full point solution from the previous time is subtracted to get the jump value due to the drops. The geometry matrix will have fewer rows to account for the subset of measured pseudoranges which doesn’t include the dropped pseudoranges. The jump values are calculated for clock bias and each component of position. The subscript “wd” indicates that the calculation is performed “without drops”, meaning without observations from GPS satellites which have lost tracking on the current sample, while the subscript “full” indicates that the full complement of visible GPS satellites is used.

Embedded Image 10

Embedded Image 11

Embedded Image 12

The jumps due to the adds, ephemeris changes, and drops are then summed together to produce the total jump value for the current sample, and this total jump value is added to the accumulated jumps. Note that at most sample times, no events which would produce a jump have occurred, so that the accumulated jump value tends to remain constant for approximately 90% of the samples. Jumps for each component of the ECEF position are accumulated separately.

Embedded Image 13

Embedded Image 14

Embedded Image 15

Embedded Image 16

Figure 3 shows the effect of removing the accumulated jumps from the time estimates for the GRACE-FO orbit using broadcast ephemerides. The original error is quite choppy and has one large peak that exceeds 14ns. Subtracting the accumulated jump values results in a much smoother error and removes the 14ns peak. This simple example shows that removal of the accumulated jumps can indeed have a beneficial effect in reducing both the maximum and RMS error values. Note that this benefit has been achieved without any knowledge of the orbital dynamics of the spacecraft, or the absolute values of the user clock error. Only the relative values of the user time estimate computed at the same instant using ”before” and ”after” GPS constellations or ephemerides are needed.

FIGURE 3
FIGURE 3

Original SISRE time error, accumulated jumps, and time error with jumps subtracted

Although the removal of the accumulated jumps has resulted in the smoothing of the error and dramatic reduction in the peak errors, the remaining error after the subtraction of the jumps develops a bias of -4ns by the end of the first hour, possibly negating any benefit. Despite the fact that the jump values have a near zero mean over the full 10-days from November 16-25, 2018 (−2.6 ⋅ 10−5), the accumulated jumps sometimes develop significant offsets, up to several meters, for periods of a few hours. Subtracting the accumulated jumps imparts these offsets to the resulting smoothed time estimates.

3.2 Estimation of accumulated jump biases

To capitalize on the error reduction afforded by the jump removal, one must also derive an estimate for the bias of the accumulated jumps and add it to the estimate which has had the jumps removed (the remaining “slopes”). The offset estimate must capture the slowly varying bias of the accumulated jumps, while ignoring the rapid changes caused by each jump. Many methods exist for estimating the bias which can be broadly divided into two categories: 1) causal for missions truly requiring real-time position and time estimates, and 2) noncausal for missions which can tolerate some delay in the estimates.

For missions requiring the best real-time estimate of the position, the estimate must use only the previously accumulated jumps; this is the causal case. However, some missions do not require instantaneous knowledge of the position and time, so the estimate may be calculated after a short delay. For example, imaging spacecraft require accurate estimates of the spacecraft location when the image was captured, but do not immediately downlink the imagery. Since there is a delay between the image capture and downlink, point solution and jump values available after the image capture may also be used to improve the bias estimate; this is the noncausal or near-real-time case.

In addition to using data points following the desired estimate, the noncausal case also offers the opportunity to improve the carrier phase ambiguity estimates. If the allowable delay is long enough (>35 minutes), the full set of code and carrier phase differences from each pass of a GPS satellite may be used to estimate carrier phase ambiguity, rather than improving the estimate recursively as is typically done with Hatch filtering.

The choice of bias estimator is driven by the desire for near optimum performance in reducing the standard deviation of the errors, along with the need for simplicity. For real-time applications, the estimator should not introduce excessive latency, which would degrade performance. Clearly, one wants the estimator to ignore the high-frequency components of the accumulated jumps, while tracking the low frequency (long term) variations.

This problem was addressed by Wiener in the 1940s with resulting optimal filters for both the causal and non-causal solutions being decaying exponentials (Brown & Hwang, 1992). Casual and noncausal decaying exponential filter coefficients are shown in Figure 4 for various processing delays. The noncausal filter support is varied in length to reflect the acceptable delay to produce the time and position estimate and ranges from 15 minutes to one hour. All filters are normalized to provide 0dB gain at DC. Noncausal filters are implemented as simple finite impulse response (FIR) filters, requiring up to 360 multiply-accumulate operations for an hour of support (+/− 360 10-second samples).

FIGURE 4
FIGURE 4

Decaying exponential filters for causal and non-causal estimation of the accumulated jump bias

Interestingly, the response of a first-order PLL is equivalent to the causal decaying exponential. Here, the discipline of high stability timekeeping offers some insight into how PLLs may be used for clock stability improvement. Clocks with good short-term stability characteristics are often combined with clocks having excellent long-term stability, such as GPS, using PLLs (Green, 2012; Helsby, 2003).

The equation for the PLL bias tracker is shown below. The bias estimate is performed separately for the three elements of the ECEF position (x, y, z) and the user clock estimate 𝑏, using the corresponding accumulated jumps for each one Embedded Image. At each sample interval, the difference between the current bias estimate and the current accumulated jump value, Δ𝑎𝑐𝑐𝑢𝑚𝑡 is calculated to find the bias estimate error. The value of 𝛼 is small to reduce discontinuity. The gain 𝛼 is calculated by dividing the sample period by the desired time constant.

Embedded Image 17

Although the initial bias, 𝑏𝑖𝑎𝑠0 is set to zero, the initial SISRE is non-zero, and is unknown. This results in a fixed offset, which must be estimated by averaging the difference between the smoothed and original time and position estimates over a longer period, such as a full day. Alternatively, a second PLL with a much longer time constant may be used to estimate and remove the initial offset. For results reported in this research, the initial SISRE bias has been directly calculated and removed.

The noncausal jump bias estimation is implemented as a FIR filter, with 𝑀 filter coefficients 𝑎 equal to the exponential curves shown in Figure 4. 𝑀 is the sum of the number of noncausal coefficients 𝑔 and causal coefficients ℎ. Here, 𝜆 is the time constant and 𝜏 is the sample period. The exponential coefficients are normalized by the coefficient total prior to use to ensure that the DC gain is one.

Embedded Image 18

Embedded Image 19

Embedded Image 20

Embedded Image 21

The estimated jump mean using the first-order PLL with a 20-minute time constant is shown in Figure 5, along with the resulting time error. The bias estimate follows the mean of the accumulated jumps and mirrors the bias in the error with the jumps subtracted. The final error has near zero mean when the initial SISRE is removed and has much lower peak error.

FIGURE 5
FIGURE 5

Time error minus jumps, accumulated time jumps, PLL jump mean estimate, and improved time error

Figure 6 compares two correction approaches with the original point solution time error. The corrections use causal bias estimation via PLL with a 20-minute time constant and a noncausal bias estimation filter with a 20-minute delay and time constant. It is clear that jump removal with bias estimation has dramatically reduced the maximum errors as well as the error standard deviations.

FIGURE 6
FIGURE 6

Baseline SISRE-induced time error and reduced time error using jump removal with casual and noncausal jump bias estimation

3.3 Smoothing performance on SISRE-induced time and position errors

Various combinations of time constants and delay for bias estimation were simulated to characterize the achievable performance. Figure 7 shows the error standard deviations for a 10-day SISRE simulation covering November 16-25, 2018. As before, ideal pseudoranges with a 10-second sample period were generated using the NGA post-processed GPS ephemerides and clock error estimates along with the post-processed GRACE-FO positions from JPL, while the archived broadcast ephemerides from Nov. 16-25 of 2018 were used in the point solution algorithm. Only GPS satellites visible to GRACE-FO during this period were used in the point solution calculation. Other impairments, such as receiver noise and multipath, were not added to the simulation.

FIGURE 7
FIGURE 7

SISRE time error standard deviation vs. time constant for causal and noncausal mean estimators with delays up to 1 hour

Time constants ranged from 1 minute to 1 hour for processing delays of 1 hour, 30 minutes, 20 minutes, 15 minutes, and real-time. For the full 10-day simulation, jump removal with causal bias estimation using a time constant of 23 minutes reduced the time error standard deviation from 0.57m to 0.32m, a reduction of 43%.

It is apparent that both the causal and noncausal bias estimators enable the jump removal to provide significant improvement to the original time estimate and that they are not very sensitive to the time constant. Changing the time constant by a few minutes will provide performance within a few percent of the optimum value, as evidenced by the wide range of near minimum error at the bottoms of the curves.

As shown in Table 1, the noncausal bias estimation further reduced the error standard deviation to 0.26m with a 15-minute delay and 0.22m with a 1-hour delay, reductions of 55% and 61%. The maximum time errors were similarly reduced from 4.09m with the PLL producing a maximum time error of 1.44m, a 65% reduction, and the noncausal filters reducing the maximum error to between 1.34m and 1.01m, reductions of 67% and 75%, respectively.

View this table:
TABLE 1

SISRE time errors vs. delay

Table 2 shows the improvements in SISRE-induced position errors. Causal 3D RMS position errors were reduced by 25% with real-time processing and up to 37% with a 1-hour delay. The maximum 3D errors were also reduced by 61% with real-time processing and up to 72% with a 1-hour delay.

View this table:
TABLE 2

SISRE position errors vs. processing delay

Figure 8 shows the results of jump removal on the ECEF X error for causal bias estimation via PLL with a 20-minute time constant, noncausal bias estimation with a 15-minute time constant, and the original time error. In this case, the original estimate contains one error over 5m, which is dramatically reduced by both the causal and noncausal methods. Table 2 summarizes the reductions in the 3D RMS position errors for various time constants. It is important to note that this approach allows for the reduction of position estimation errors solely based on bias estimation of accumulated changes in position (jumps). Position changes are found by differencing point solutions calculated at the same instant in time using “before” and “after” GPS constellations. No orbital models are required.

FIGURE 8
FIGURE 8

SISRE ECEF X error reduction via jump removal and bias estimation

4 GRACE-FO PROCESSING

To demonstrate the efficacy of the proposed approach on real data, Level 1B pseudorange observations from the GRACE-FO mission were processed. GPS receivers on board the GRACE-FO spacecraft capture observations of pseudorange and carrier phase on the L1 and L2 frequencies at 10-second intervals, and these observations are made available on the NASA ftp server.

Ten days of GRACE-FO B observables, from November 16 to November 25 of 2018, were downloaded from the NASA Physical Oceanography Distributed Active Archive Center (PODAAC) (NASA, n.d.). The GRACE-FO observables captured by the 48-channel BlackJack GPS receiver include P1 and P2 code pseudoranges as well as L1 and L2 carrier phases, which were used to calculate the P12 and L12 ionosphere-free pseudoranges.

The P12 and L12 data were used along with broadcast GPS ephemerides and clock estimates to calculate point solutions. The jump removal technique was then applied to reduce position and time errors, and the results were compared to the GRACE-FO precise post-processed orbits (Wen et al., 2019).

The GRACE-FO user clock bias was calculated using the JPL GRACE-FO post-processed orbits and the NGA post-processed GPS orbits and clock errors. This clock bias was used as the “truth” for the GRACE-FO clock to evaluate the performance of the user clock estimates.

To compare the position results we calculated to those provided by JPL, the offsets from the GPS antenna to the GRACE-FO center of mass were applied to the position estimates after the jump removal algorithm was applied. A radial offset of 0.48m and an along-track offset of 0.3m were subtracted from the point solution and smoothed point solution results to reflect the GPS antenna position on the spacecraft relative to the spacecraft center of mass (Wen et al., 2019).

4.1 Preprocessing

To focus on the performance of the jump removal algorithm, we screened the GRACE-FO B observations to remove large P12 errors and L12 cycle slips. Ideal pseudoranges were calculated using NGA post-processed GPS times and positions along with the JPL precise GRACE-FO orbits. The ideal pseudoranges were then used to calculate absolute errors in P12 and to detect cycle slips in L12. P12 observations with errors greater than 12m were removed from processing, which resulted in a removal of 0.7% of the observables. L12 observations were removed if a jump larger than 1.0m occurred, which caused carrier smoothing algorithms using Hatch filtering and full pass P12-L12 averaging to reset following carrier phase slips. For the 10 days from Nov. 16 to Nov. 25 of 2018, 3,163 carrier phase cycle slips larger than 1.0m were detected, out of a total of 746,773 L12 ionosphere-free observations.

Many methods exist for detecting and correcting carrier phase slips, but most rely on double differencing data from multiple GPS receivers (Liu, 2010). For single receiver detection in real-time, the Melbourne-Wübbena linear combination provides the best option for detecting cycle slips and is often used in preprocessing to detect bad data and cycle slips (Bock et al., 2009; Liu, 2010). The choice to drop observations with relatively large slips rather than correcting the slips was made for simplicity and assumes that cycle slips greater than 1.0m are relatively easy to detect, especially when doing noncausal processing of full passes.

Both causal and noncausal carrier smoothing were applied. For the real-time calculations with and without jump removal, a real-time Hatch filter was used (Hatch, 1982). The Hatch filter continuously updates the estimated carrier phase offset based on the average difference between the P12 and L12 values for a given PRN. This method provides increasing accuracy as the number of samples with continuous phase increases, but has large errors when the GPS signal is first tracked, because very few P12 points are available for averaging. For this reason, Hatch pseudoranges were not used until at least 5 points were available for averaging. This resulted in slightly reduced pass durations for each GPS satellite but significantly improved the quality of the pseudoranges.

The Hatch filter may be implemented recursively using the following equation:

Embedded Image 22

Here, 𝑘 is the number of samples processed by the Hatch filter up to time 𝑡, Embedded Image is the current Hatch filter output, 𝜌𝑡 is the current pseudorange observation, Embedded Image is the previous Hatch output, and Φ𝑡 and Φ𝑡−1 are the current and previous carrier phase measurements, respectively.

Noncausal jump removal allows the carrier phase offset to be estimated using P12 and L12 values that occur later in time than the sample being calculated, improving the carrier smoothing. If the allowable processing delay for carrier smoothing is longer than the maximum GPS pass length, the average difference between P12 and L12 for the entire pass may be used to estimate and correct the carrier phase ambiguity. For the two noncausal scenarios evaluated, a carrier smoothing delay of 10 minutes was applied, so that P12-L12 offsets up to the current time were used to estimate the carrier phase offset for the sample 10 minutes prior. This eliminates much of the large error experienced by the Hatch filter at the beginning of each pass, when the P12 observables are still very noisy due to lower SNR.

Figure 9 shows the original L12 and P12 pseudorange errors for one pass of PRN 1, along with the real-time Hatch filter output with the first 5 samples removed, the Hatch filter output with a 10-minute delay, and the Hatch filter output with a full pass delay. Errors in pseudorange and carrier phase observations were calculated using the post-processed GPS orbits and clock errors from the NGA, along with the precise orbit estimates from JPL. The Hatch output has significantly less noise than the P12 signal, but initially has a large error due to the low SNR at the beginning of the pass, and the fact that only a few points are averaged. The Hatch output with a 10-minute latency provides a result nearly as good as the Hatch filter output with full pass delay.

FIGURE 9
FIGURE 9

Carrier phase smoothing with Hatch filter, delayed Hatch filter and full pass averaging for a single pass

4.2 Processing scenarios

Four processing scenarios, which used the GRACE-FO observations and broadcast ephemerides, were run to evaluate the jump removal performance:

  • Baseline Scenario: Hatch Filter Only (no jump removal)

  • Scenario 2: Hatch Filter and Jump Removal with Causal Jump Mean Estimation (PLL)

  • Scenario 3: Hatch Filter with 10-Minute Delay (no jump removal)

  • Scenario 4: Hatch Filter with 10-Minute Delay and Jump Removal with Noncausal Jump Mean Estimation (FIR filter) with 10-Minute Delay

The delayed processing scenarios could be used for missions which have a delay between acquiring and down-linking data, for example, imaging spacecraft which store image data after capture and only downlink it once per orbit. The Hatch filter with a 10-minute delay and no jump removal was included to differentiate the Hatch filter effect from the impact of the noncausal jump removal with bias estimation. For each scenario, the optimum time constants found in the SISRE processing were used.

Table 3 summarizes the complete results from the GRACE-FO processing. The real-time scenarios using broadcast ephemerides demonstrate that maximum time errors may be reduced by 26%, 3D RMS errors by 16%, and maximum 3D errors by 32% using the Hatch filter with a PLL for jump mean estimation. Noncausal processing with 20-minute latency reduced the standard deviations and maximum errors by roughly 50% for time and radial position and roughly 40% for 3D position.

View this table:
TABLE 3

GRACE-FO B processing results

Figure 10 shows the time and 3D RMS errors for a representative 3-hour period on November 20th, 2018. Although the real-time Hatch filter with jump removal and bias estimation shows some improvement, the non-causal processing with 20-minute latency provides significantly smaller errors. Both causal and noncausal processing dramatically reduce the maximum time and 3D RMS errors.

FIGURE 10
FIGURE 10

Time and 3D RMS errors for GRACE-FO processing with Hatch filter only, Hatch filter with jump removal and PLL for causal jump bias estimation, and delayed Hatch filter with noncausal jump bias estimation and 20-minute total delay

Figure 11 shows the radial, along-track, and cross-track position errors for the same times. The largest benefit from jump removal occurs for the radial errors with the non-causal processing providing the greatest improvement. The along-track errors show limited improvement from jump removal, particularly for real-time processing. However, even with real-time processing, the maximum errors are significantly reduced.

FIGURE 11
FIGURE 11

Radial, along-track and cross-track position errors for GRACE-FO processing with Hatch filter only, Hatch filter with jump removal and PLL for causal jump bias estimation, and delayed Hatch filter with noncausal jump bias estimation and 20-minute total delay

5 CONCLUSIONS

The proposed approach offers a method for small satellites to reduce on board GPS point solution errors in real-time and near-real-time, while minimizing the required processing. The jump removal approach allows users to implement first-order PLLs or FIR filters to dramatically reduce error standard deviations and maxima, rather than using complex orbit determination algorithms to estimate the orbital states.

Time error standard deviation is reduced by 17% for the causal case and 47% for the noncausal case, and time error maxima are reduced by 26% and 53%, respectively. Radial error standard deviations are similarly reduced by 25% and 48%, and radial maxima are reduced 31% and 54% for causal and noncausal cases. Along-track and cross-track errors also show small improvements with the exception of the along-track error standard deviation with real-time Hatch filtering and jump bias estimation, which increases from 0.46m to 0.56m. 3D RMS errors are reduced 16% and 39% for causal and noncausal processing, and 3D maximum errors are reduced by 32% and 52%, respectively.

HOW TO CITE THIS ARTICLE

Van Buren D, Axelrad P, Palo S. A low complexity smoothing algorithm for improved GPS point solutions on board LEO spacecraft. NAVIGATION. 2021;68:185–198. https://doi.org/10.1002/navi.410

  • Received May 25, 2020.
  • Revision received November 27, 2020.
  • Accepted December 13, 2020.

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. Ali Shah, S. A., &
    2. Arshad, U.
    (2013, August). General system design of Cubesat in LEO for IR imaging. 2013 International Conference on Aerospace Science & Engineering (ICASE), Islamabad, Pakistan, 1–5. https://doi.org/10.1109/ICASE.2013.6785552
  2. ↵
    1. Alminde, L., &
    2. Laursen, K. K.
    (2009, June). A strategic approach to developing space capabilities using Cubesat technology. 2009 4th International Conference on Recent Advances in Space Technologies, Istanbul, Turkey, 43–47. https://doi.org/10.1109/RAST.2009.5158239
  3. ↵
    1. Bock, H.,
    2. Dach, R.,
    3. Jäggi, A., &
    4. Beutler, G.
    (2009). High-rate GPS clock corrections from CODE: Support of 1 Hz applications. Journal of Geodesy, 83, 1083–1094. https://doi.org/10.1007/s00190-009-0326-1
  4. ↵
    1. Brown, R. G., &
    2. Hwang, P. Y. C.
    (1992). Introduction to random signals and applied Kalman filtering. John Wiley and Sons, Inc.
  5. ↵
    1. Bryerton, E. W.,
    2. Hesler, J. L.,
    3. Retzloff, S. A., &
    4. Crowe, T. W.
    (2016, July). 874-GHz heterodyne Cubesat receiver for cloud ice measurements - Flight model data. 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 5553–5555. https://doi.org/10.1109/IGARSS.2016.7730449
  6. ↵
    1. Chang, L.,
    2. Chao, C.-K.,
    3. Kuo, C.-L.,
    4. Liu, J.-Y.,
    5. Chandran, A.,
    6. Fang, T.-W.,
    7. Hari, P., &
    8. Kandi, K.
    (2018, August). IDEASSat: The ionosphere dynamics explorer and attitude subsystem satellite. 32nd Annual AIAA/USU Conference on Small Satellites, Logan, UT.
  7. ↵
    1. Chattopadhyay, G.
    (2017, December). Terahertz instruments for Cubesats. 2017 IEEE MTT-S International Microwave and RF Conference (IMaRC), Ahmedabad, India, 1–5. https://doi.org/10.1109/IMaRC.2017.8449669
  8. ↵
    1. Choi, E.,
    2. Yoon, J.,
    3. Lee, B.,
    4. Park, S., &
    5. Choi, K.
    (2009, October). Development of spaceborne GPS receiver with real-time orbit determination using unscented Kalman filter. 60th International Astronautical Congress, Daejeon, Republic of Korea, 2891–2896.
  9. ↵
    1. Decrossas, E.,
    2. Chahat, N.,
    3. Walkemeyer, P. E., &
    4. Velasco, B. S.
    (2019, July). Deployable circularly polarized UHF printed loop antenna for mars cube one (MarCO) Cubesat. 2019 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting, Atlanta, GA, 1719–1720. https://doi.org/10.1109/APUSNCURSINRSM.2019.8888027
  10. ↵
    1. Engelen, S.,
    2. Quillien, K.,
    3. Verhoeven, C.,
    4. Noroozi, A.,
    5. Sundaramoorthy, P.,
    6. Wijnholds, S.,
    7. Rajan, R. T.,
    8. Boonstra, A.-J.,
    9. Bentum, M.,
    10. Meijerink, A., &
    11. Budianu, A.
    (2013). The road to OLFAR - A roadmap to interferometric long-wavelength radio astronomy using miniaturized distributed space systems. Proceedings of the International Astronautical Congress, 6, 4357–4363.
  11. ↵
    1. Ge, H.,
    2. Li, B.,
    3. Ge, M.,
    4. Zang, N.,
    5. Nie, L.,
    6. Shen, Y., &
    7. Schuh, H.
    (2018). Initial assessment of precise point positioning with LEO enhanced global navigation satellite systems (LeGNSS). Remote Sensing, 10(7), 984. https://doi.org/10.3390/rs10070984
  12. ↵
    1. Green, C.
    (2012, April). Advantageous GPS disciplining of the OSA BVA oscillator. 2012 European Frequency and Time Forum, Gothenburg, Sweden, 263–266. https://doi.org/10.1109/EFTF.2012.6502379
  13. ↵
    1. Hatch, R. R.
    (1982, February). The synergism of GPS code and carrier measurements. Volume 2 of Proceedings of the Third International Symposium on Satellite Doppler Positioning, Las Cruces, NM, 1213–1231.
  14. ↵
    1. Heckler, G.,
    2. Gramling, C.,
    3. Valdez, J., &
    4. Baldwin, P.
    (2016, May). TDRSS augmentation service for satellites. SpaceOps 2016 Conference, Daejeon, Korea, 2016-2467. https://doi.org/10.2514/6.2016-2467
  15. ↵
    1. Helsby, N. C.
    (2003, May). GPS disciplined offset-frequency quartz oscillator. IEEE International Frequency Control Symposium and PDA Exhibition Jointly with the 17th European Frequency and Time Forum, Tampa, FL, 435–439. https://doi.org/10.1109/FREQ.2003.1275131
  16. ↵
    1. Hodges, R. E.,
    2. Chahat, N.,
    3. Hoppe, D. J., &
    4. Vacchione, J. D.
    (2017). A deployable high-gain antenna bound for Mars: Developing a new folded-panel reflectarray for the first Cubesat mission to Mars. IEEE Antennas and Propagation Magazine, 59(2), 39–49. https://doi.org/10.1109/MAP.2017.2655561
  17. ↵
    1. Imken, T.,
    2. Castillo-Rogez, J.,
    3. He, Y.,
    4. Baker, J., &
    5. Marinan, A.
    (2017, March). Cubesat flight system development for enabling deep space science. 2017 IEEE Aerospace Conference, Big Sky, MT, 1–14. https://doi.org/10.1109/AERO.2017.7943885
  18. ↵
    1. Kalemci, E.,
    2. Ümit, E., &
    3. Aslan, R.
    (2013, June). X-ray detector on 2U Cubesat BeEagleSAT of QB50. 2013 6th International Conference on Recent Advances in Space Technologies (RAST), Istanbul, Turkey, 899–902.
  19. ↵
    1. Kornfeld, R. P.,
    2. Arnold, B. W.,
    3. Gross, M. A.,
    4. Dahya, N. T.,
    5. Klipstein, W. M.,
    6. Gath, P. F., &
    7. Bettadpur, S.
    (2019). GRACE-FO: The gravity recovery and climate experiment follow-on mission. Journal of Spacecraft and Rockets, 56(3), 931–951. https://doi.org/10.2514/1.A34326
  20. ↵
    1. Liu, Z.
    (2010). A new automated cycle slip detection and repair method for a single dual-frequency GPS receiver. Journal of Geodesy, 85, 171–183. https://doi.org/10.1007/s00190-010-0426-y
  21. ↵
    1. Menzione, F.,
    2. Renga, A., &
    3. Grassi, M.
    (2017). On-board orbit determination for low thrust LEO-MEO transfer by consider Kalman filtering and multi-constellation GNSS. Acta Astronautica, 138, 242–254. https://doi.org/10.1016/j.actaastro.2017.05.029
  22. ↵
    1. Misra, P., &
    2. Enge, P.
    (2006). Global positioning system: Signals, measurements, and performance. Ganga-Jamuna Press.
  23. ↵
    1. Montenbruck, O.
    (2003). Kinematic GPS positioning of LEO satellites using ionosphere-free single frequency measurements. Aerospace Science and Technology, 7(5), 396–405. https://doi.org/10.1016/S1270-9638(03)00034-8
  24. ↵
    1. Montenbruck, O., &
    2. Gill, E.
    (2006, December). Phoenix-XNS - A miniature real-time navigation system for LEO satellites. 3rd ESA Workshop on Satellite Navigation User Equipment Technologies, Noordwijk, The Netherlands.
  25. ↵
    1. Montenbruck, O.,
    2. Helleputte, T.,
    3. Kroes, R., &
    4. Gill, E.
    (2005). Reduced dynamic orbit determination using GPS code and carrier measurements. Aerospace Science and Technology, 9(3), 261–271. https://doi.org/10.1016/j.ast.2005.01.003
  26. ↵
    1. Montenbruck, O., &
    2. Ramos-Bosch, P.
    (2008). Precision real-time navigation of LEO satellites using global positioning system measurements. GPS Solutions, 12(3), 187–198. https://doi.org/10.1007/s10291-007-0080-x
  27. ↵
    1. Montenbruck, O.,
    2. Steigenberger, P., &
    3. Hauschild, A.
    (2018). Multi-GNSS signal-in-space range error assessment - Methodology and results. Advances in Space Research, 61(12), 3020–3038. https://doi.org/10.1016/j.asr.2018.03.041
  28. ↵
    NASA. n.d. Physical oceanography distributed active archive center website. https://podaac.jpl.nasa.gov/dataaccess
  29. ↵
    NGA. n.d. Office of geomatics website. earth-info.nga.mil
  30. ↵
    1. Norton, C. D.,
    2. Millar, P. S.,
    3. Bauer, R., &
    4. Komar, G. J.
    (2016, July). New capabilities for earth science measurements with 6U Cubesat technologies. 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 5550–5552. https://doi.org/10.1109/IGARSS.2016.7730448
  31. ↵
    NovAtel. n.d. Oem719 multi-frequency GNSS receiver datasheet. www.novatel.com/assets/Documents/Papers/OEM719-Product-Sheet.pdf
  32. ↵
    1. Parkinson, B., &
    2. Spilker, J.
    (1996). Global positioning system: Theory and applications, vol. II. American Institute of Aeronautics & Astronautics.
  33. ↵
    1. Peral, E.,
    2. Tanelli, S.,
    3. Haddad, Z.,
    4. Sy, O.,
    5. Stephens, G., &
    6. Im, E.
    (2015, July). Raincube: A proposed constellation of precipitation profiling radars in Cubesat. 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Milan, Italy, 1261–1264. https://doi.org/10.1109/IGARSS.2015.7326003
  34. ↵
    1. Perea, S.,
    2. Meurer, M.,
    3. Rippl, M.,
    4. Belabbas, B., &
    5. Joerger, M.
    (2017). URA/SISA analysis for GPS and Galileo to support ARAIM. NAVIGATION, 64(2), 237–254. https://doi.org/10.1002/navi.199
  35. ↵
    1. Rajan, R. T.,
    2. Boonstra, A.,
    3. Bentum, M.,
    4. Wolt, M.,
    5. Belien, F.,
    6. Arts, M. J.,
    7. Saks, N., &
    8. Veen, A.-J.
    (2015). Space-based aperture array for ultra-long wavelength radio astronomy. Experimental Astronomy, 41, 271–306. https://doi.org/10.1007/s10686-015-9486-6
  36. ↵
    1. Reichert, A.,
    2. Meehan, T., &
    3. Munson, T.
    (2002, September). Toward decimeter-level real-time orbit determination: A demonstration using the SAC-C and CHAMP spacecraft. Proceedings of the 15th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GPS 2002), Portland, OR, 1996–2003.
  37. ↵
    1. Reid, T. G.,
    2. Neish, A. M.,
    3. Walter, T., &
    4. Enge, P. K.
    (2018). Broadband LEO constellations for navigation. NAVIGATION, 65(2), 205–220. https://doi.org/10.1002/navi.234
  38. ↵
    1. Segret, B.,
    2. Hestroffer, D.,
    3. Quinsac, G.,
    4. Agnan, M.,
    5. Vannitsen, J., &
    6. Mosser, B.
    (2018, March). In-flight orbit determination for a deep space Cubesat. 2018 IEEE Aerospace Conference, Big Sky, MT, 1–12. https://doi.org/10.1109/AERO.2018.8396756
  39. ↵
    1. Tkacenko, A.,
    2. Towfic, Z.,
    3. Stratton, M. C.,
    4. Kroll, R. T.,
    5. Abraham, D. S.,
    6. Finley, S. G.,
    7. Malhotra, S.,
    8. O’Dea, A.,
    9. Malphrus, B. K., &
    10. Conner, C. D.
    (2019, March). Insight/marco opportunistic multiple spacecraft per antenna (OMSPA) demonstration. 2019 IEEE Aerospace Conference, Big Sky, MT, 1–18. https://doi.org/10.1109/AERO.2019.8741775
  40. ↵
    1. Wallis, R. E.,
    2. Weaver, G. L.,
    3. Reinhart, M. J., &
    4. Sheng, C.
    (2005, March). An advanced synthesized ultra-stable oscillator for spacecraft applications. 2005 IEEE Aerospace Conference, Big Sky, MT, 952–960. https://doi.org/10.1109/AERO.2005.1559386
  41. ↵
    1. Weaver, G. L., &
    2. Kantsiper, B. L.
    (2007, March). Developments toward a disciplined timekeeping system for lunar and planetary navigation. 2007 IEEE Aerospace Conference, Big Sky, MT, 1–9. https://doi.org/10.1109/AERO.2007.352928
  42. ↵
    1. Wen, H.Y.,
    2. Kruizinga, G.,
    3. Paik, M.,
    4. Landerer, F.,
    5. Bertiger, W.,
    6. Sakumura, C.,
    7. Bandikova, T., &
    8. Mccullough, C.
    (2019). Gravity recovery and climate experiment follow-on (GRACE-FO) level-1 data product user handbook. https://podaac-tools.jpl.nasa.gov/drive/files/allData/gracefo/docs/GRACE-FO_L1_Handbook.pdf
  43. ↵
    1. Xuewen, G.,
    2. Guo, L.,
    3. Wang, F.,
    4. Wanwei, Z.,
    5. Sang, J.,
    6. Ge, M., &
    7. Schuh, H.
    (2019). Precise onboard real-time orbit determination with a low-cost single-frequency GPS/BDS receiver. Remote Sensing, 11, 1391. https://doi.org/10.3390/rs11111391
  44. ↵
    1. Yang, Y.,
    2. Yue, X., &
    3. Dempster, A. G.
    (2016). GPS-based onboard real-time orbit determination for LEO satellites using consider Kalman filter. IEEE Transactions on Aerospace and Electronic Systems, 52(2), 769–777. https://doi.org/10.1109/TAES.2015.140758
Back to top

In this issue

Print
Download PDF
Article Alerts
Email Article
Citation Tools
A low complexity smoothing algorithm for improved GPS point solutions on board LEO spacecraft
Damon Van Buren, Penina Axelrad, Scott Palo
NAVIGATION: Journal of the Institute of Navigation Mar 2021, 68 (1) 185-198; DOI: 10.1002/navi.410
Twitter logo Facebook logo Mendeley logo
Bookmark this article

Jump to section

Keywords