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 (Reference Colombo, Kluesberg and TeunissenColombo, 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 (Reference Hofmann-Wellenhof, Lichtenegger and CollinsHofmann-Wellenhof and others, 2001; Reference LeickLeick, 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. Reference Manson, Coleman, Morgan and KingManson 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 (Reference SanterreSanterre, 1991; Reference King, Coleman and NguyenKing 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 Reference SanterreSanterre (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 (Reference SanterreSanterre, 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 (Reference King, Coleman and NguyenKing 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 (Reference BlewittBlewitt, 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-Reference Bindschadler, King, Alley, nandakrishnan and Padmanschadler 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 (Reference Zumberge, Heflin, Jefferson, Watkins and WebbZumberge and others, 1997). This is equivalent to the processing performed in the AutoGIPSY service freely available on the Internet (Reference ZumbergeZumberge, 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 (Reference Bar-Sever, Kroger and BorjessonBar-Sever and others, 1998). These data were additionally processed using a more rigorous (for these data) kinematic PPP approach (Reference King and AokiKing and Aoki, 2003), showing that WIS moves with a stick–slip behaviour (Reference Bindschadler, King, Alley, nandakrishnan and PadmanBindschadler 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 (Reference LichtenLichten, 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 (Reference Bindschadler, Vornberger, King and PadmanBindschadler 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 (Reference Bondesan, Capra, Gubellini and TisonBondesan and others, 1994; Reference King, Nguyen, Coleman and MorganKing 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 (Reference DoakeDoake 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 (Reference DoakeDoake 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

Reference King, Coleman and NguyenKing 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 (Reference King, Coleman and NguyenKing 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 (Reference Collins, Langley and LaManceCollins 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 (Reference SanterreSanterre, 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 (Reference ChenChen, 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 (Reference Saastamoinen, Henricksen, Mancini and ChovitzSaastamoinen, 1972, Reference Saastamoinen1973), 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 (Reference Collins, Langley and LaManceCollins 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.