Rigorous GPS data-processing strategies for glaciological applications

Abstract

Global positioning system (GPS) data are now routinely used for many glaciological applications. In some common cases, systematic errors are unmodelled at the data-processing stage, although they are often presumed insignificant. In this paper, I investigate these assumptions for three different scenarios: (1) measurements on a moving glacier; (2) measurements on a floating ice shelf; and (3) precise height determination over large elevation ranges, such as for aircraft positioning in lidar/laser altimeter missions. In each case, systematic errors are shown to be present in the coordinate solutions that have a far greater magnitude than the formal error estimates produced by the GPS processing software, under certain conditions. If these coordinate biases go undetected, short- and long-term measurements of horizontal ice velocity or rates of ice-thickness change may be erroneous and the coordinates could not be expected to match rigorously processed data or results from different processing techniques. More rigorous processing strategies are discussed that allow for bias-free parameter estimation.

1. Introduction

The advent of an operational global positioning system (GPS) satellite constellation during the 1980s revolutionized the determination of precise coordinates on Earth’s surface. Glaciologists were amongst the first to adopt this technology (e.g. Hinze and Seeber, 1988) due to the disadvantages in using traditional survey techniques in polar regions where distances from stable survey marks with accurate coordinates were often long and conditions for manual observations typically hostile. Since that early period, GPS has been used for a variety of glaciological studies such as horizontal glacier/ice-shelf velocity measurement, ice-sheet mass-balance determination, strain-rate measurements, grounding-zone investigations and the validation of remote-sensing measurements (e.g. Hulbe and Whillans, 1994; Vaughan, 1994; Hamilton and others, 1998; Phillips and others, 1998). Achievable accuracies and precisions have been reduced from metres to millimetres during this time, with corresponding increases in possible sampling rates from once per day to once per second or higher. However, accuracies at the millimetre level (or even decimetre level) will not be achieved unless all systematic errors are accounted for. Normally it is presumed that these are each accounted for in the GPS processing software, although this may not be entirely true for applications which push the bounds of the algorithms implemented in the software, as is often the case with many glaciological applications, or where inappropriate processing techniques are used.

In this paper, I outline several possible systematic errors in typical glaciological GPS data and show that if these are not treated correctly at the processing stage they will result in biases in the derived coordinate time series. To achieve this, three case studies are presented using real glaciological data covering: (1) positioning on moving grounded ice such as a glacier; (2) positioning on floating ice such as an ice shelf; and (3) airborne elevation measurements such as in laser-altimetry/lidar projects. In each case the possible magnitude of the systematic error is demonstrated and techniques are discussed for rigorously processing such GPS data to ensure the effects of such systematic errors on derived coordinates are minimized. For background reading on GPS, its error sources and common methods of dealing with them, readers are referred to the excellent reference texts of Hofmann-Wellenhof and others (2001) and Leick (2004).

2. Case Studies

2.1. Velocity and surface height measurement of grounded ice

When attempting to monitor a moving object, a kinematic processing strategy may be adopted whereby a new set of antenna coordinate parameters is solved each and every measurement epoch. However, kinematic processing becomes increasingly problematic for baselines longer than 50–100km due to the reduction in common satellites observed simultaneously at the base and rover, and the increase in the relative ionospheric delay. This increase in relative ionospheric delay increases the noise in the solution and therefore makes ‘ambiguity fixing’ more difficult to achieve. Ambiguity fixing is the process of fixing to integers the real-value estimates of the carrier-phase wavelength ambiguities (Colombo, 1998). Importantly, it is possible to fix ambiguity parameters to the incorrect integer, especially when GPS biases or high noise levels are present. While the effects of the ionosphere may be almost entirely eliminated using the ‘ionosphere-free’ linear combination, this combination cannot be used for ambiguity fixing since the integer nature of the ambiguities is not preserved (Hofmann-Wellenhof and others, 2001; Leick, 2004). However, following the ambiguity fixing stage, the ionosphere-free combination may be used to produce the final solution.

Without ambiguity fixing, standard kinematic solutions are typically imprecise over long (>30 km) baseline distances. The only exception is a solution in which noninteger ambiguities are constant (no cycle slips) and a smoothing Kalman filter solution is used. For many glacio-logical applications in polar regions, where baseline lengths may necessarily be several hundred kilometres and ionospheric activity relatively high, the standard kinematic processing approach is not sufficient. One common approach to overcoming this limitation is to assume that for a short period of time (say, 0.5–4.0 hours) a glacier is effectively stationary and hence a static processing methodology may be adopted for this period where only one set of antenna coordinates is determined for the entire data span (e.g. Manson and others, 2000). This approach, however, violates the least-squares requirement that the functional model fully describes the relationship between the observations and the parameters. Consequently, the residuals will be biased and hence minimization of their ‘weighted squares’ will result in some parameter bias.

In order to determine the magnitude of these parameter biases and the extent to which they affect the estimates of three-dimensional coordinates, I firstly performed a sensitivity analysis of this problem (Santerre, 1991; King and others, 2003). A sensitivity analysis allows the geometry of the least-squares problem to be constructed – in this case the chosen base-and-rover combination along with the GPS satellite constellation – without the need for actual observations. It also provides a mechanism for estimating the effect of biases introduced through the ‘observed-minus-computed’ term (b) in the classical least-squares equation:

(1)

where A is the design matrix, W is the inverse of the variance–covariance matrix of the observations and x contains the corrections to the a priori parameter values. For simplicity and following Santerre (1991), we used the ‘single difference’ observable where observations from two receivers are used to difference out the satellite clock errors, while also estimating a differential receiver clock error at each epoch. The single difference is mathematically equivalent to the more common ‘double difference’ found in most GPS processing software packages, so the results shown here also apply more generally. For simplicity, the baseline was chosen to be of zero length, apart from the introduced biases.

In all, three tests were performed by introducing velocity signal with a magnitude of 1 md^–1 and, in turn, azimuths of 0°, 45° and 90°, simulating glacier flow directly north, northeast and east. The estimated parameters included the local topocentric north, east, up coordinates (n, e, u) at chosen intervals, receiver clock parameters every epoch and ambiguity parameters. To correctly simulate the effect of processing data in short sessions, one set of coordinate parameters was estimated every 1 hour and the ambiguity parameters were reset at these times. The effect of ambiguity fixing was simulated by choosing not to estimate these parameters (Santerre, 1991). Zenith tropospheric delay parameters were also estimated, but since they had no impact on the results shown here they were not estimated in the final solutions.

Figure 1 shows the biases in each parameter for each of the three velocity orientations for the ambiguity-free (real valued ambiguities estimated) and the ambiguity-fixed (ambiguities assumed fixed to their correct integer values) cases. In the case of the northward-flowing ‘ice’ (Fig. 1a, d and g), the east coordinates become biased in the ambiguity-free solution by a value approximately 5–10% of the magnitude of the horizontal velocity. There is no bias in the other parameters, and indeed once the ambiguity parameters are fixed there is no bias in any of the parameters. The superimposed dots represent the ambiguity parameter estimates mapped into the direction of each of the coordinate components. It is evident that the ambiguity parameters are absorbing the unmodelled velocity during the period of each segment, and the least-squares process together with the GPS satellite constellation results in these biases being ‘balanced’ by, in this case, a corresponding bias in the east coordinate (King and others, 2003).

Fig. 1. Biases (D) introduced in coordinate estimates due to an unmodelled 1 m d^–1 horizontal velocity being present in static 1hour GPS solutions. (a, d, g) correspond to a northward velocity; (b, e, h) correspond to a northeast velocity; (c, f, i) correspond to an eastward velocity. The solid line is the solutions where ambiguity parameters were solved as real-valued numbers, while the dashed line is the solutions where ambiguities were not estimated (fixed to integers). The dots are the projection of the solved ambiguity parameters from the satellite line-of-sight to the direction of the respective coordinate components. In all cases, the introduced velocity has been removed prior to plotting.

For eastward-flowing ice (Fig. 1c, f and i), a southward bias is evident and an additional downward bias in the height component. No bias is evident in the east component. Again, these parameter biases are balancing biases in the ambiguity parameters. Once ambiguity parameters are fixed, the biases disappear. The parameter biases for the northeast-flowing ice (Fig. 1b, e and h) are a combination of the north-flowing and east-flowing ice. For ice flowing in opposite directions to these, the signs of the biases are reversed.

In addition to the offsets evident in Figure 1, small oscillations around the mean are also evident in each of the solutions where ambiguities were not fixed. These oscillations are present even where no offset is evident (e.g. Fig. 1f). With a peak-to-peak amplitude of 0.05m with a 1md^–1 input signal, this could easily be interpreted as genuine, although confusing, periodic motion of the ice. With real GPS data containing random noise, it is doubtful, however, that this signal would be clearly observed and hence it will simply add to the apparent noise of the measurements. Fixing ambiguities will therefore reduce this ‘noise’ component significantly, over and above that normally expected from ambiguity fixing (Blewitt, 1989). However, fixing ambiguities to their correct integer values may be impossible due to the biased nature of the floating-point estimates and hence ambiguity fixing may not be a practical way to remove these spurious signals in many cases.

These results are for a baseline at a single geographical location (70° S, 0° E), although the biases at other locations were computed and it was found that these results only begin to differ at low latitudes where periodic signals begin to appear in the horizontal and vertical components with a similar amplitude to the offsets in Figure 1. The period of these signals appears linked to the mean satellite pass period, which also varies with latitude. While session lengths of 1hour were used in this study, longer session lengths did not alter the size of parameter bias. However, the parameter bias was always proportional to the magnitude of the velocity, and hence for ice with larger (smaller) velocities the bias would be larger (smaller).

Figure 2 shows GPS antenna positions determined using this short-segment static approach for data collected on Whillans Ice Stream (WIS), Antarctica, at site G2 (Bind-schadler and others, 2003b). These data are shown as circles, with the size of the circles representing the root-sum-squared (rss) error of the horizontal coordinates from the GPS software. In this case, a segment length of 4hours was adopted and the data were processed using the Precise Point Positioning (PPP) approach implemented in the GIPSY/ OASIS II software (Zumberge and others, 1997). This is equivalent to the processing performed in the AutoGIPSY service freely available on the Internet (Zumberge, 1999). The PPP approach is convenient for processing GPS data collected at remote locations since it does not require the operation of an additional base station. Instead it relies on predetermined precise satellite clock and position estimates while estimating receiver clock correction parameters along with the site coordinates. Other than the short data-segment length, the processing strategy was standard (Bar-Sever and others, 1998). These data were additionally processed using a more rigorous (for these data) kinematic PPP approach (King and Aoki, 2003), showing that WIS moves with a stick–slip behaviour (Bindschadler and others, 2003a, b). The ambiguities were not fixed to integers in this solution, but the use of a smoothed square-root information filter within GIPSY (Lichten, 1990) enabled precise kinematic results to be determined.

Fig. 2. Stick–slip motion of WIS at site G2 using 5min kinematic (black dots) and 4 hour static (grey circles) solutions (Bindschadler and others, 2003b). For the static results, the size of the circle reflects the horizontal rss error

It is apparent from Figure 2 that when the WIS site is stationary the static and kinematic methodologies produce coordinates that are in good agreement. When the WIS is moving, however, the static results are biased by several decimetres. From the sensitivity analysis above, we would expect that for ice flowing in a north-northwest direction (as occurs at this site) a bias would be present mainly towards the east and then to the south, as is evident in Figure 2. During the 4hour sessions that cover the slip events, the average velocity is approximately 2 md^–1, suggesting a bias magnitude of approximately 0.2 m which is again in agreement with the biases seen in the static estimates in Figure 2.

2.2. Velocity/tidal measurement of floating ice

The typical motion of an ice shelf outside its grounding zone is analogous to a grounded ice stream with a vertical tidal motion added. Again, static data-processing strategies are often used to determine horizontal velocity and tidal motion (Bondesan and others, 1994; King and others, 2000). Due to the increased noise introduced by the station tidal motion during each segment, ambiguities cannot normally be fixed in such a solution.

Figure 3 shows horizontal and vertical positions determined using this short-segment static approach for data collected at Halley station on the Brunt Ice Shelf, Antarctica. Processing was performed as described above for the WIS data, for both the kinematic and static solutions. The horizontal coordinates (Fig. 3a) show that both processing methodologies suggest that the ice shelf is moving in response to the tides (Doake and others, 2002), although the actual signals differ significantly. This difference is especially apparent towards the end of the time series where the tidal signal is larger. In contrast, the vertical signals are in good agreement (Fig. 3b), although a small positive height bias is evident as expected for westward-flowing ice.

Fig. 3. Horizontal (a) and vertical (b) motion at Halley station, using 5 min kinematic (black dots) and 4 hour static (grey circles) solutions (Doake and others, 2002). For the static results, the size of the circle reflects the horizontal rss error. In (a), the detrended north motion is shown for clarity since most of the motion is westerly. Instead of the east component, the non-standard westerly component is plotted to allow convenient comparison with the height component (b). The error bars on the static height solutions (b) are at the 95% confidence interval

King and others (2003) have investigated this problem for real and synthetic data for the Amery Ice Shelf (AIS), Antarctica, and show that the presence of the vertical tidal motion results in a large spurious signal, primarily in the east coordinate component, in data processed with a static methodology. Data processed using a kinematic methodology was found to be free of these signals. The source of these signals was found again to be due to the absorption of the unmodelled signal into the estimated real-valued ambiguity parameters. They disappeared when these were fixed to integers. The magnitude of the signal for the AIS was found to be approximately 40–50% of the unmodelled vertical signal. If such a spurious signal were present in the Halley station results, the biases in the east coordinate component would reach a magnitude of approximately 0.4– 0.5 m. Figure 3a, however, shows that the difference between the signal from the kinematic and static methodologies is much smaller than this, with a magnitude closer to 0.05m. This may be explained, however, by the frequency of estimation of zenith tropospheric delay (ZTD) parameters: for the AIS, these were estimated at the same frequency as the station coordinates (1 hour); for the Halley data they are estimated every 5 min. The higher frequency of estimation means that the vertical tidal motion, which is highly correlated with the ZTD parameters, is preferentially absorbed by the ZTD parameters rather than being forced into the ambiguity parameters and hence the horizontal coordinate components. Regardless of the ZTD estimation rate, the only way to obtain unbiased horizontal coordinates is to use a kinematic processing approach or, when using a static approach, fix the carrier-phase ambiguities to their correct integer values (King and others, 2003).

2.3. GPS positioning for airborne lidar/laser altimetry

Several glaciological applications require GPS heights to be computed over baselines with large elevation differences, such as for controlling airborne instrument positions for laser-altimetery or lidar measurements. Due to the centimetre-level precision of these ranging techniques, it is important that GPS heights are as accurate as possible. While most of the potential biases in GPS are not particularly sensitive to elevation differences, the ZTD is, and particular care needs to be taken in its modelling. In the absence of accurate surface meteorological data, tropospheric models poorly represent the elevation-dependent variation of the troposphere, and relative ZTD biases of 10–20mm are possible over an elevation difference of 500 m at high latitudes (Collins and others, 1996). Any mismodelling of the ZTD maps into the height component of the derived coordinates with 1 mm of ZTD bias, producing a corresponding height bias of 2-3 mm (Santerre, 1991). For base-rover combinations with similar elevations over short baselines the subsequent errors largely cancel. However, baselines with large elevation differences will be adversely impacted by any mismodelling, even over short baselines. Since accurate meteorological measurements are typically not attainable in the field, standard measurements or modelled values are normally adopted in proprietary software, sometimes with a user-selectable scale factor to attempt to account for this. Regardless, without very accurate meteorological data at both ends of a baseline, determining accurate heights over long (>30 km) baselines or large (>200 m) height differences is not possible.

The alternative to relying solely on tropospheric models is to remove the nominal tropospheric delay as estimated by one particular model and then parameterize the residual ZTD along with the site coordinates, etc. This method is employed in ‘scientific’ software (software developed in a research environment), but the most recent proprietary software does not allow for ZTD estimation. In effect, some coordinate precision is sacrificed in order to obtain greater accuracy. Figure 4a shows results from processing a baseline near the grounding zone of the AIS using a kinematic processing methodology. The ~65 km baseline consists of a base station on a rocky outcrop approximately 620 m above the grounded ice where the rover was situated, and is analogous to an airborne survey with the elevations of the base and rover reversed. The baseline was processed using two different proprietary software packages and also the Track ‘scientific’ processing software (Chen, 1998). Identical processing options were selected in the different softwares where possible, including the use of the ionosphere-free linear combination. The Saastamoinen tropospheric model was adopted in each case (Saastamoinen, 1972, 1973), but the residual ZTD was also estimated in Track. The proprietary softwares had difficulty processing these data due to the 60s sampling interval of the base station and the length of the baseline, and good solutions were not produced at every epoch. Since similar results were obtained from both of the proprietary software packages, only one solution is shown here alongside the Track solution. The smoothness of the Track height time series when compared to the proprietary solution is due to the flexibility of being able to tightly constrain the motion of the site in Track’s Kalman filter (via the process noise) during the estimation process. To emphasize the effect of the tight constraints, a more loosely constrained Track solution is also shown as dots in Figure 4a.

Fig. 4. (a) Kinematically derived heights of a glacier site near the AIS grounding zone using Track with and without ZTD estimation (thick solid line and dashed solid line respectively) and a proprietary software (thin solid line). A loosely constrained Track solution (with ZTD estimation) is shown as dots. (b) The difference between the tightly constrained Track solution (with ZTD estimated) and the proprietary solution. (c) The estimated ZTD from the tightly constrained Track solution.

A 0.05 m bias between the heights obtained from Track and the proprietary software is evident in Figure 4b. The troposphere is the only possible source of such a large bias since the antenna phase centre models are identical in the different software, multipath effects are the same and the baseline is too short to be influenced by any geophysical effects unmodelled in just one software package. Due to the more rigorous processing strategy employed in Track, we presume that this is the ‘true’ result, although no doubt some errors remain due to signal multipath effects, amongst others. The mismodelling of the troposphere has resulted in an apparent lengthening of the baseline in the proprietary solution, suggesting that too much delay has been removed from the base station data and/or too little delay removed from the rover data. This is the scenario expected when negative biases are present in the temperature value used in the hydrostatic (‘dry’) tropospheric model or when negative biases are present in the water-vapour pressure used in the wet tropospheric model (Collins and others, 1996).

The residual ZTD delay estimated in Track (Fig. 4c) does not account for the entire inter-software height bias, with the mean of 5 mm being equivalent to a mean height offset of ~10-15 mm. This may also be seen in Figure 4a where a Track solution without ZTD estimation (dashed line) is shown with the Track solution with ZTD estimation (thick solid line). The remainder of the bias may be explained, however, by different meteorological parameters being input into the Saastamoinen model in the two different software packages.

3. Discussion

In order to discuss the impacts on glaciological parameters due to the biases demonstrated above, we consider the horizontal biases separately from the vertical ones. Since horizontal coordinates are most typically used to determine ice velocities, Figure 3 alongside those in King and others (2003) suggests that, when an unmodelled tidal signal is present, velocities determined using a static technique will exhibit a time dependency, with the magnitude being a function of the distance travelled between each occupation; two 4hour occupations on subsequent days may produce ice-shelf velocities that have no connection whatsoever with reality.

GPS data collected on floating ice may not be the only situation in which the periodic horizontal signals exhibit themselves, however. Periodic vertical signals may be evident even in grounded ice due to failure to model the solid-Earth body tides and ocean tide loading (OTL; Melchior, 1966). While it is routine to model solid-Earth body tides in ‘scientific’ software, such as GIPSY, GAMIT or Bernese, the situation is less clear for proprietary software packages which are not designed for long-baseline processing. Furthermore, modelling OTL is also not routine, partly due to uncertainty about the accuracy of the numerical tide models, especially for the circum-Antarctic seas. While OTL will be largely cancelled over short baselines, this may not be true, for example, in the Antarctic Peninsula region where OTL gradients are steep (Yi and others, 2000) or when the PPP methodology is employed. Since OTL may have vertical amplitudes as large as 50–100 mm, there is potential for large, but spurious, periodic horizontal signals to appear in GPS-derived time series of grounded-ice motion.

The height-coordinate biases due to tropospheric mis-modelling or from easterly site motion are mostly of concern when used to calibrate remotely sensed measurements such as from GPS-controlled airborne laser altimetry. (Relative tropospheric errors may also be present over long baselines in ground-based surveys, with height biases as large as 0.1– 0.3 m possible.) In the case of elevation-dependent tropo-spheric mismodelling, surface height measurements made from different elevations will be biased by different amounts as the degree of tropospheric mismodelling is elevation-dependent for a given set of meteorological values. This bias becomes particularly insidious when an altimeter bias calibration is performed by comparing to ground-truth elevations (e.g. Spikes and others, 2003), since the troposphere-related elevation error will be absorbed into this altimeter bias. Consequently, the height measurements will be accurate when the aircraft is at the elevation of the calibration, but height-dependent biases will be present at any other aircraft elevation. In this case the surface height biases could be as large as those mentioned above, but with opposite sign. With even less accurate estimates of the relevant meteorological parameters, these biases could be even greater.

4. Conclusions

I have shown that processing techniques commonly applied to glaciological GPS data do not account for systematic errors, and the impact on the derived coordinates may be large enough to reduce their scientific value. GPS data collected on a moving surface must be processed in such a way that these errors are accounted for, typically by processing using a kinematic methodology or alternatively by fixing ambiguities when using short ‘static’ segments. With regard to airborne GPS surveys, residual ZTD must be correctly accounted for by applying accurate height-dependent scales to the models, or preferably, parameterizing the ZTD along with the antenna coordinates.

Fortunately the systematic errors for some surveys are time-independent, and hence velocity estimates on smoothly flowing grounded ice, for example, will remain unbiased. However, differences in the estimated quantities will occur if results that are biased are compared to results from a technique that does not contain these same systematic errors, such as from another technique, or GPS data processed taking these errors into account.

Perhaps the most important aspect of these results is that, without an independent check on the coordinates, the GPS processing software gives few indications about solution errors. Figure 2, in particular, shows that when systematic errors remain in GPS data, the formal error estimates from the processing software are wildly over-optimistic. In fact, even when state-of-the-art processing is used on a truly static site, the stochastic model in the software is still deficient, due to neglected between-epoch correlations, etc., and hence the formal error estimates will remain over-optimistic by factors of typically ~4-20 for static data; formal error estimates for kinematic data are normally more realistic due to the apparent lower number of degrees of freedom. The most convenient way to correctly scale these formal errors is to assess the repeatability of repeat measurements – something that must be factored into any rigorous observation campaign.