Geometry

The basic instrumental and geometric limitations to the INSAR technique, and the need for calibrating ground truth, were established in early demonstrations of interferometry concerned mainly with measuring topography (e.g. Reference Zebker and GoldsteinZebker and Goldstein, 1986; Reference Gabriel and GoldsteinGabriel and Goldstein, 1988; Reference Gabriel and GoldsteinGoldstein and others, 1988).

Figure 3 shows a schematic imaging geometry, in which the vectors r₁, r₂, B and D lie in a plane perpendicular to the path of the sensor (or azimuth direction, denoted ζ). The range vectors r₁ and r₂ are the lines of sight to a surface element common to two images. B is the difference between the two viewpoints, known as the (interferometric) baseline, and D is the projection of any surface displacement onto the plane. The range is here defined as r = |r₁| and the range difference Δ = |r₂| - |r₁|. For a monochromatic radar with wavelength λ, the absolute phase associated with a given path-length difference 2Δ is Φ = 4π(Δ/λ).

An exact equation relating Δ to the range r, the magnitudes B and D of the baseline and surface displacement and the angles α and β they make with the line of sight is obtained by noting that r₂ = r₁ - B + D. Then

(1a)

and

(1b)

An interferogram is a map of Δ(ζ, r). When B(ζ) and D(ζ,) are known, the shape of the imaged surface may be reconstructed as θ(ζ, r) by determining α(ζ, r) from Δ(ζ, r) using Equation (lb). The function θ{ζ,r) may be discontinuous over rugged topography.

Because the range r is usually very much larger than either B or D and therefore Δ, Δ is most sensitive to topography when the component of the baseline perpendicular to the line of sight, B sin α, or simply B, is large. Range differences due to surface displacement are almost independent of the baseline, and sensitive only to displacement parallel to the line of sight. In principle, the influences on the range difference of topography and surface displacement can be separated by making use of their different dependencies on the baseline and time (e.g. Joughin and others, 1996b).

The range difference Δ cannot be measured directly, because the interferometric phase is initially known only modulo. The first part of restoring the integer wavelength part of the path-length difference is a process known as “phase unwrapping” (Reference Gabriel and GoldsteinGoldstein and others, 1988), by which phase gradients are integrated and the phase at any point related to an arbitrary datum. Calibration then establishes the arbitrary datum’s true phase.

The spatial gradients in interferometric phase that can be resolved and unwrapped depend on the size of the radar image elements, the radar wavelength and the level of phase noise. Steep slopes, large values of B or long repeat periods over surfaces undergoing rapid shear may produce phase gradients which are greater than π rad pixel⁻¹ and cannot be unwrapped.

Pairs of radar images may be collected in a single pass using two antennae (e.g. Reference Zebker and GoldsteinZebker and Goldstein, 1986; Reference ThompsonThompson and others, 1994), in which case the baseline length may be known to sub-wavelength precision (its orientation may still have to be determined by calibration). When a single antenna is used in repeat-pass interferometry, the baselines depend on aircraft paths or satellite orbits whose determinations are subject to errors equivalent to many radar wavelengths. Sufficient parameters to describe the unknown baseline must then be determined by calibration. To interpret an interferogram in terms of surface elevation and surface displacement therefore generally requires tie-points at which these quantities are known.

Insar over ice sheets using ERS-1 and -2

Those satellite-borne SAR systems flown to date have allowed only repeat-pass interferometry. Fortunately for glaciologists, the surface topography of ice sheets and typical rates of ice movement are such that ERS-1 acquisition Phases B (28 December 1991-31 March 1992) and D (1 January 1994-7 April 1994) and the ERS-1/2 Tandem Mission (April 1995-June 1996), with their respective repeat periods of 3, 3 and 1 d, were almost ideal for producing useful inter-ferograms.

An important limitation for glaciologists is the fact that any set of nominally similar SAR images can provide information about displacement in only one direction, toward or away from the radar’s path. A typical ERS-1/2 SAR look angle gives a sensitivity 2-3 times greater to vertical than to horizontal displacement, so, despite generally small surface slopes, the detail of interferograms over ice sheets may owe as much to variations in vertical motion associated with flow over undulating topography as to variations in horizontal movement. This is particularly true in the middle of fast-flowing glaciers. To extract directly any general surface motion requires three independent viewing directions. Orbit configurations which achieve this are rare, although to a very limited extent some ERS-1/2 configurations allow higher latitudes to be imaged from two or more significantly different crossing orbits. More generally, some constraint, such as the direction of the horizontal component of motion, or the assumption of surface-parallel motion (e.g. Joughin and others, 1996b) is required. Both require an elevation model of sufficient resolution, from INSAR or elsewhere, and imagery from two significantly different directions. When imagery from only one look direction is available, flow may be assumed to be “downhill” (i.e. down the greatest surface slope over an appropriate length scale) or inferred from the orientation of glacier margins, etc.

Fig. 2. (a) Interferogram. (b) SAR amplitude (backscatter) image and survey stations. Filled circles denote stations for the 1979-81 and 1984-86 surveys, crosses the 1994-96 GPS survey. Selected line features are taken from visible and SAR imagery.

Fig. 3. Orbit geometry diagram to illustrate nomenclature used in determining the interferometric baseline.

Glacier movement was first observed using INSAR by Reference Goldstein, Engelhardt, Kamb and FrolichGoldstein and others (1993) who created an interferogram of part of Rutford Ice Stream, Antarctica, from two ERS-1 SAR Single Look Complex (SAR.SLC) images. Thanks to its orientation, a 6 d interval between images and an unusually small baseline, this interferogram was almost uninfluenced by topography and contained a great deal of information about ice motion. The area covered included the Rutford Ice Stream grounding line, which Goldstein and others identified easily by interpreting the variation of interferometric phase as a combination of contributions from horizontal ice flow and vertical movement on the ocean tide. With a single interferogram and no independent elevation model, however, Goldstein and others were unable to isolate either of the two contributions from ice movement or the very small contribution from topography. Differential, or multi-baseline, interferometry techniques to achieve this have since been developed by Reference Joughin, Kwok and FahnestockKwok and Fahnestock (1996) and others. They generally require a combination of baselines (i.e. more than one pair of images) that includes at least one baseline large enough to provide the necessary elevation resolution. With suitable imagery and some simplifying assumptions, surface elevation, ice flow and tidal motion can all be isolated.

Reference Zebker and GoldsteinZebker and Goldstein (1986) identified a difference in viewing angle beyond which there can be no correlation. Expressed in terms of B at a typical look angle, the limit for the ERS-1/2 SAR systems is approximately 1100m, but good results are hard to obtain with B much greater than 600 m. Strictly, it is not B, but B - D, that is limited. If, however, the period separating coherent images over ice sheets is limited to a few days, then absolute displacements greater than a few metres will be rare in glaciological applications and B itself will be limited. Of the images acquired over ice sheets during the ERS-1 3 d repeat Phases B and D (first and second ice phases), not all pairs for which B is less than 600m (roughly 70% of the total) are correlated. This may in part be because the SAR processors at the Processing and are Archiving Facilities (PAFs) were not optimised for interferogram production, but probably the major factor is the stability of the ice-sheet surface, on which INSAR relies. Stability will depend on the amount of melting, snowfall and drift and to what extent the returned signal comes from volume as well as surface scattering. The latter is a subject of debate, but clearly the longer the time interval, the more likely are significant changes. Interferograms have been produced using images over ice-sheet surfaces taken 35 d apart (Reference Fahnestock, Joughin and KwokFahnestock and others, 1997), but success with such long periods appears uncommon.

The expense and logistical difficulties associated with operating SAR receiving stations in Antarctica allowed only a fraction of the potential Antarctic imagery from the ERS-1 acquisition Phases B and D to be collected. The ERS-1/2 SARs do not “see” that part of the ice sheet south of 80° S. North of this latitude there remain areas for which no imagery exists and areas where the combination of archived images will not allow the extraction of reliable and glaciologically useful quantities. The ERS-1/2 Tandem Mission, with its 1 d separation between instruments, has seen coverage improve dramatically, particularly within the range of the receiving station operated at McMurdo Station, but unless coherence over the 35 d repeat period proves to be common, ideal image sets may continue to be rare.

Accuracy and ground control

The factors that determine the accuracy and applicability of INSAR in a particular case include the availability of correlated image pairs from which interferograms can be formed, the degree of correlation, whether the area of interest is imaged from more than one direction, whether a separate digital elevation model (DEM) is available, and the quality of ground control available to optimise the interferometric baselines.

The applications so effective in demonstrating the glaciological potential of INSAR have sometimes lacked rigour in their use of ground truth. Joughin and others (1996a) have probably done most to address this, with their analysis of the effects of tie-point errors and tie-point distribution in simulated interferograms. Imagery over the ice sheets of Antarctica and Greenland will often contain no areas of exposed rock for use as tie-points. To what extent creating interferograms from whole swaths of images reduces the amount of ground control required in practice, and how effectively assumptions about ice flow can substitute for it, are open questions. Clearly, however, the more we can use existing surveys to calibrate interferograms in such remote areas and thus avoid revisiting them, the more effectively can the available logistical effort be used in campaigns to produce regional data sets of glaciological interest. Of equal importance, but sometimes overlooked, is the need for rational methods with which to assess the significance of inconsistencies between interferograms and between interferograms and ground surveys.

Presented here are the generation, calibration and interpretation of interferograms formed from ERS-1 SAR.SLC images of the area around Rutford Ice Stream, Antarctica. The process is described in some detail for a single interferogram using ERS-1 Phase B images and data from a series of conventional, i.e. pre-GPS, surveys of the area. Values of surface elevation and velocity obtained from these surveys contain strongly covariant errors that are taken into account along with the interferometric phase noise when optimising parameters that describe the interferometric baseline. This approach ensures that when data from recent GPS surveys and a single interferogram from ERS-1 Phase D that partly overlaps the earlier one are examined, inconsistencies between the four data sources can he interpreted reliably as variations in ice flow.

Image selection

Coherence between images is not guaranteed by the exclusion of pairs with large baselines or long repeat periods, and may be marginal where it exists. Among the potentially useful imagery in the Rutford Ice Stream area, no sets comprising three mutually coherent images or two coherent pairs of images were found. No differential interferograms from which to extract separate topographic and movement contributions could therefore be created.

The ERS-1 SAR.SLC images used (Table 1) were supplied by the German Processing and are Archiving Facility (DPAF) as quadrants of full scenes that cover a nominal 100 km by 100 km area of the Earth’s surface. Two of the images, acquired 6 d apart in February 1992, are nearly contemporary with, and overlap, those used by Reference Goldstein, Engelhardt, Kamb and FrolichGoldstein and others (1993), but cover an area further upstream on Rutford Ice Stream and include Carlson Inlet. Figure 2b shows the earlier image in conventional (backsratter or σ°) form. The Ellsworth Mountains appear along the western boundary. Note the bright (high-backscatter) ice-stream margins and the bright area in the middle of the ice stream associated with rapidly shearing flow over a prominent knoll above a step in the ice-stream bed (Reference Frolich, Vaughan and DoakeFrolich and others, 1989). Note also that the brightness in the margins continues to the upstream edge of the image, whereas visible surface crevassing, in general, does not. This indicates significant penetration of the ice-sheet surface and volume scattering within the firn. The two images from Phase D, also supplied by the DPAF, were acquired 3 d apart in January 1994.

Table 1. Scene identifications

Image registration

ERS-1 SAR.SLC images are supplied with elements arranged in the natural slant-range/azimuth coordinate system. That is, each line in the image corresponds to a particular azimuth (time or point along the orbit) and consists of pixels containing information from progressively more distant targets that share a common azimuth of closest approach. This is achieved by “zero-Doppler” processing, i.e. constructing the synthetic aperture so that each point appears on the image line that corresponds to the time at which the signal returned from the point is not Doppler-shifted. Each pixel spans roughly 4 m in azimuth and 8 m in range.

Because nominally similar images always cover slightly different areas on the ground, there are generally offsets between images in azimuth and range. Only if these offsets, which can be equivalent to 100 pixels or more, are found to within 1 or 2 pixels will the phase difference between two images have any spatial correlation. The process of finding the offsets, called registration, consists of comparing the azimuth and range and the pixel and line numbers of points common to two or more images. This is straightforward between images containing sufficient common features identifiable at a pixel level. On ice sheets, where there are usually no such features unless transponders or corner reflectors have been placed specificially for the purpose, initial estimates of the integer offsets can be obtained from a combination of the processing information contained in the SLC file header and a model for the satellite orbits.

Information contained in the European Space Agency file headers is not always entirely self-consistent, so the use of header information that cannot be verified independently is minimised here. This places particular reliance on the model used to describe the satellite orbit, in this case the “precise” orbit model provided by the Technical University of Delft (1997). The Delft model is supplied as a series of satellite positions (one per minute, roughly 100 per orbit) and software for interpolating between them. As is usual with such orbit models, the errors claimed for values of the radial component of position are smaller (on the order of 10 cm) than those claimed for the other components, which are on the order of a metre.

Relative position within an image is obtained by comparing the azimuth and range of closest approach with those of some reference. Absolute positioning in azimuth is possible to the extent that the orbit model can reliably estimate the point along the orbit of closest approach (and to the extent that the Doppler centroid frequency is indeed zero). Absolute and relative positioning in range are less straightforward, since they rely on accurate models for the refractive index between radar and target. Here the geographical position for the scene centre supplied in each header file is taken as a range reference. The offsets between images are naturally estimated at these scene centres.

An image and the orbit on which it was acquired are chosen as references. The orbit model is then used to find the time and range of closest approach on this orbit to the centre of the reference image. The time, which is usually within a few milliseconds of the scene centre time given in the header, is identified with the image’s central line, and the range is identified with its centre pixel. The process is repeated to find the range and time at which the geographical positions given for the other scene centres were imaged in the reference image. These offsets in range and time are convened to pixel and line offsets by linear interpolation assuming the scene-start and scene-end times and the pixel range length are as given in the header.

An integer search within 10 pixels in azimuth and 2 pixels in range (both equivalent to roughly 40 m on the ground) of these offsets is usually adequate to establish whether coherence is present.

Testing For correlation

Associated with each image element is a complex value C = Re^iØ , where R and Ø are the amplitude and phase of the received signal. The interferometric phase Ø ₁ — Ø ₂ is the phase of the product C ₁ C ₂* = R ₁ R ₁e ⁱ (Ø ₁-Ø ₂, where C ₁ and C ₁ are the values assigned to corresponding pixels in the two images and the superscript “*” denotes complex conjugate.

Correlation is sought by comparing small areas between images (e.g. Reference Goldstein, Engelhardt, Kamb and FrolichGoldstein and others, 1993). One common method is to search for the offsets that give the highest spectral value in a fast Fourier transform (FFT) of the interferometric phase over the test area. Another is to search for the highest value of the complex correlation Coefficient

where ... denotes expected value.

Spatial gradients in interferometric phase generally ensure that an estimate of ρ made by simple summation is biased downward. The estimates of the optimum offsets are not generally biased, but discrimination is reduced. Phase gradients may be estimated in order to optimise the estimate of , but this is usually unnecessary if the test areas are divided into sub-areas over which estimates of ρ made by simple summation before being averaged over the test area. Cases of marginal coherence (e.g. ρ ≤ 0.3) sometimes require experimentation with the size of both the test area (e.g. up to 128 pixels square) and the sub-areas (typically 16 pixels in azimuth by 4 pixels in range, roughly-square on the ground). It may also be useful in such cases to compare results between FFT and summation methods. Sometimes, of course, neither works because there is simply no correlation.

For the best results, i.e. the minimum phase noise, the search may be performed to sub-pixel precision by means of interpolation and resampling within the natural slant-range/azimuth coordinate system. Although SLC images may be resampled without significant loss of information, the process is computer-intensive, especially when offsets vary across an image because of distortions in either the processing geometry or the surface being imaged (due to ice motion, for example).

Correlation and phase noise

A map of the correlation coefficient ρ is shown in Figure 4 for the southeastern quadrant ofthe Phase B interferogram (cf. Fig. 2a). Areas of higher correlation (ρ > 0.5; phase noise of 15° or less) generally coincide with the darker areas in the SAR backseatter image (Fig. 2b). Areas of low correlation (ρ < 0.2; phase noise greater than 40°) are generally brighter in the backseatter image and coincide with the areas of noise seen in the interferogram.

Noise in the radar system, changes in the ice-sheet surface, differences in viewing angle, and inexact registration all contribute to noise in the interferometric phase. Of these, the first two are usually the most important. Phase noise hinders registration, introduces uncertainty into any “Observation” made on an interferogram and hinders the process of phase unwrapping by which each element’s phase is referred to a single, initially arbitrary, datum. Simple averaging over blocks of complex values before calculating the interferometric phase reduces the phase noise at the expense of reduced spatial resolution. When this exacerbates problems with phase unwrapping in regions of high slope or high shear, more sophisticated filtering may be appropriate.

Mosaicking

Separate interferograms were created for each quadrant of the Phase B reference image by applying those optimal inte-ger offsets in range and azimuth found by the method already outlined. The interferogram quadrants were mosaicked by a similar registration process using test areas from reference-scene quadrant overlaps. The Phase D interferogram is formed from the first quadrant of a scene overlapping the Phase B interferogram as shown in Figure 2b.

Fig. 4. SAR coherence image of the upper right quadrant of Figure 2a.

Projection

Shown in Figure 2a is a projection of the Phase B interferogram mosaic following averaging over blocks of pixels 2 in the range direction by 12 in the azimuth direction. In particularly rugged terrain, the foreshortening associated with side-looking radars such as the ERS-1 SAR can be many kilometres. Areas may even by obscured by nearby higher terrain “laid over” them. This does not occur over the gentle terrain of ice-sheet surfaces, but surface elevations must still be known throughout the image for the transformation from slant range and azimuth to geographical position to be performed accurately. To form Figure 2a in the absence of this information, the whole of the imaged surface was assumed to be 200 m above the ellipsoid surface. Since the ice-sheet surface elevation exceeds 700 m in places, projected distances from the satellite ground track may be in error by up to 1 km.

Sources of ground truth

Both interferograms used are formed from “lone pairs” of images, with no third coherent image or second coherent pair to help isolate elevation information which could then be removed to leave only motion information. Since neither of die baselines is negligibly small, both elevation and velocity must be known at the tie-points for them to be useful.

Information with which to calibrate interferograms can come from a variety of remotely sensed sources. Crevasses and even quite subtle surface features ran sometimes be tracked between pairs of visible or radar satellite images (e.g. Reference Bindschadler and ScambosBindschadler and Scambos, 1991). Satellite imagery may be used to locate ice summits (e.g. Reference Goodwin and VaughanGoodwin and Vaughan, 1995) or other places where the surface velocity is very small. Reference Goldstein, Engelhardt, Kamb and FrolichGoldstein and others (1993) calibrated their interferogram of part of Rutford Ice Stream by assuming a negligible velocity at a point off the ice stream close to the exposed rock of Flowers Hills in the Ellsworth Mountains. Some flow features, including ice-stream margins and ice divides, may tentatively be assumed to lie parallel or perpendicular to current ice flow. Where they lie respectively perpendicular or parallel to the radar’s line of sight, the contribution of ice motion to the interferometric phase will vanish. In all remotely sensed cases, some assumption must be made about the vertical component of velocity, which may be significant in areas of high accumulation rate even when the horizontal velocity is low.

Tie-point location errors and poorly known elevations will often make ground control derived from satellite imagery less precise and reliable than purpose-made surface observations. Such control is, nevertheless, perfectly usable when properly weighted using estimates of the errors involved. Although it would be instructive to attempt to use the remotely sensed control that exists in the Rutford Ice Stream area, ground truth in this work comes entirely from surface observations whose errors, while not negligible, are at least easily quantified.

The ground truth includes two major surveys of Rutford Ice Stream, from 1978-80 and 1984-86 (Reference Doake, Frolich, Mantripp, Smith and VaughanDoake and others, 1987). Both were based on angle and distance measurements using theodolite and electronic distance measurement. A number of smaller surveys, some partly intended as ground control for interferometry, were carried out using global positioning system (GPS) between 1994 and 1996. One of these crossed the Rutford Ice Stream-Carlson Inlet junction; another followed the ice-stream centre line upstream of the 1984-86 survey. A single observation of velocity on Carlson Inlet was made using GPS over the year 1995 (personal communication from A. Jenkins, 1996). In the area covered by the images used here, these surveys provide over 200 observations of surface elevation and surface velocity for potential use as tie-points (Fig. 2b).

In the Phase B interferogram, the 1978-80 survey appears in the extreme south (near range), the 1984-86 survey extends from near to far range and includes 100 km of fast-flowing ice stream (Fig. 2b), and the observation from Carlson Inlet appears at the eastern edge. The Phase D interferogram overlaps the earlier one by roughly 15 km in range at the upstream end of the 1984-86 survey. The 1994-96 GPS surveys extend from near to far range and over most of the azimuthal extent of the Phase D interferogram.

None of the surveys was contemporary with the Phase B acquisition period in early 1992. The Phase D images were acquired in February 1994, at the beginning of the 2 year period spanned by the GPS surveys. The validity of the ground control is therefore questionable, but the considerable redundancy may ensure that any problems are reflected as inconsistencies in the calibration.

The interferograms and ground survey data used are summarised in Table 2.

Observation weights and residual variance

The metric residual associated with tie-point i is defined as

An error variance matrix V is defined with elements V_ij equal to the expected values Oρ_i ρ_j V is diagonal when the residuals arise from observation errors that are independent. If, in addition, the observations carry equal weight, then estimates of the parameters describing the baseline and the phase datum, Δ_C, may be found by minimising the observed residual variance Σρ_i ² . Here, covariant errors incurred in estimating the tie-point positions and velocities, and systematic errors associated with the errors in the range and azimuth data, ensure that V is not diagonal. In such cases, the quantity to be minimised is the more general weighted residual variance, W, where

(3)

and R = [ρ_i...ρ_n is a row vector containing the residuals at the solution (Reference StrangStrang, 1986). In what follows, estimates of the contributions to V will be restricted to first-order terms of Equation (2).

Tie-point location

Tie-point location and phase unwrapping arc performed in the natural slant-range/azimuth system. The tie-points lie on snow, and none are marked with radar reflectors that appear in the images, so the slant range and azimuth of closest approach to each tie-point and the associated look angle are-found using the method described earlier in the context of image registration. Offsets in range and azimuth relative to the scene centre arc then converted to offsets in pixel and line.

Phase unwrapping and phase noise

To ease tie-point unwrapping, the raw interferogram is “flattened”. An estimate of the baseline from the Delft orbit model and a value for the look angle θ for each pixel are used in Equation (2) to synthesise the phase variation that would arise were the entire imaged surface at an ellipsoidal height of zero. Subtracting this phase removes a great deal of the phase variation (fringes) from the raw interferogram. A “spherical earth” approximation to the look angle is adequate for the flattening because the removed phase is restored to the unwrapped phases before the baseline parameters are optimised.

Here the tie-point phases are unwrapped manually by counting fringes. Tie-points for which the integer wavelength part of the phase cannot be determined are discarded. In the flattened full-resolution Phase B interferogram, only one path exists through the shear mar-gin between Rutford Ice Stream and Fletcher Promontory by which it is possible to tie the observations on the fast-flowing ice and slower ice to a common datum. Phase noise and rapid phase variation in the ice-stream margins still prevent the unwrapping of some 20 potential tie-points.

Complete rigour would require a local estimate ai each tie-point of the contribution to residual variance from phase noise. Here, for simplicity, the errors are assumed to be independent and random, with a standard error appropriate to each interferogram when the high-noise shear margins are excluded. Values used for the Phase B and the noisier Phase D interferogram were equivalent to 1.4 and 2.8 mm in Δ, respectively.

Location error

Systematic and random errors are incurred in estimating the tie-point positions, estimating the associated azimuth ζ, range r and look angle θ, and locating the points on the interferograms. Their effects depend on the distribution of phase in the interferograms, i.e. the surface topography, surface velocity and baseline, and the distribution of tie-points.

Absolute location in range and azimuth depends on the accuracy of the scene centre position given in the image file header. Uncertainties in the derived range and azimuth data, although probably no greater than 100 m, can produce systematic errors in the optimum baseline parameters that vary in rough proportion to Bδr/r. A range datum error of δr = 100 m with a typical ERS-1/2 baseline magnitude of B = 200 m could produce errors in the baseline components on the order of 1 m, depending on the interferogram and lie-point distribution. In general, the more regular the interferogram, the greater the baseline errors and the smaller the residuals generated in an optimisation based on many observations of Δ_U(r). If the imaged surface were plane and stationary, then residuals in Δ_U could be less than 1 mm. For an ellipsoidal surface they would be somewhat greater, but the range error could still be ιindétectable in the presence of phase noise and other errors. A discrepancy with precise orbit models might, however, arouse suspicion. The main effect of a systematic azimuth error is a translation in azimuth of the optimal baseline parameters. Since the variation of these is typically 10 m over the 100 km azimuthal extent of a full scene interferogram, a systematic azimuth error of 100 m would probably be undetectable by comparison with orbit models.

As well as contributing to a bias in the optimal baseline parameters, range and azimuth datum errors contribute lo residual variance through their interaction with the less regular components of the interferogram, identified here with the signatures of the surface elevation and velocity fields. This contribution is therefore estimated as the square of the typical magnitude of phase change over 100 m (maximum pixel and line errors) on the flattened interferogram (Fig. 2a). The fact that phase gradients are clearly correlated locally, particularly in regions of high shear, may be unimportant when the tie-points cover an interferogram well. Systematic orbit model error is assumed to have an effect similar to, but smaller than, the use of the file-header scene centre position as a reference. Random errors δζ, δr and δθ arise from the original ground surveys and use of the precise orbit model and contribute to residuals in Δ according to the products δζ(B/ζ) · r/r and B(ζ) · δ(r/r). Uncertainties in the lie-point positions are typically less than 5 m horizontally and vertically. Uncertainties in die point of closest approach are estimated at 1ms, or 7 m, for the azimuthal component, and 5 m for the components perpendicular to the orbit. Relative location on an interferogram should therefore be no worse than the length of a pixel (7.9m) in the range direction and two range-line widths (roughly 8 m) in the azimuth direction. For these interferograms, B is on the order of 100 m and gradients arc on the order of 10 m per 100 km. With tolal errors in r and ζ less than 10 m, the total residual contributions are unlikely to be greater than 1 mm.

Velocity contribution

Ice-sheet surface motion makes a contribution to the interferometric phase that is, to a close approximation, proportional to the component of the surface displacement parallel to the range vector. If average surface velocities over the short periods dial separated our images are the same as those obtained from the ground surveys over much longer periods, then the contribution to die range difference is

where t is the period between images, u is the three-dimensional surface velocity from the surveys and r₁ is the first range vector al closest approach. Consistency between the interferograms and the survey data will be the only evidence for or against the “steady flow” assumption.

Vertical movement

None of the tie-points has any tidal motion, but high velocities over an undulating surface ensure that vertical movement on Rutford Ice Stream contributes significantly to the interferometric phase. It is likely that such vertical movement could account for the short-length-scale discrepancies between interferogram and observed horizontal movement in Goldstein and others (1993, Fig. 5).

The ice velocity detected by phase differences between coherent images is a material one. The GPS-based surveys used aluminium poles inserted at least 1 m into the firn as markers. It is assumed here that their results, which include the vertical component of motion, represent material velocities, and any densification effect is ignored. Neither the 1978 80 nor the 1984-86 survey observed vertical movement directly, but barometric or optical levelling provided surface profiles along lines of survey markers.

A reasonable working hypothesis is that the ice sheet in the Rutford Ice Si ream area is in a “steady state”, which implies (hat the vertical component of the material surface velocity is balanced by the accumulation rate. This applies only on a smoothed-out conceptual surface, whose undulations on long length scales are generally tied through the ice dynamics to long-term accumulation patterns, basal topography and basal conditions. Climatic and dynamic changes over these length scales may produce non-steady surface topographies, with features such as kinematic waves (the detection of which could prove a powerful application for interferometry), Short-wavelength variations tend instead to be adverted with the ice motion until altered by wind or accumulation.

In areas of net accumulation, the material velocity is downward relative to the smoothed-out surface. It cities not matter whether coherence between images is taken to imply the absence of significant surface accumulation in the intervening period, since, although the surface accumulation is episodic, the ice velocity is much smoother over time. It is difficult to be precise about this component of the movement, since both the densification rate and the extent of radar penetration into the firn are poorly known and probably highly variable. Field observations on Rutford Ice Stream (Reference Doake, Frolich, Mantripp, Smith and VaughanDoake and others, 1987) suggest that it could be up to 1 m year.

With the assumption of surface-parallel flow, a minor correction to balance the long-term accumulation rate, and the neglect of any densification effect, the vertical component of motion can be estimated. This is most effective along two lines of survey stations running close to the ice-stream centre line approximately 2 km apart (Fig. 5). Though intended to lie along the flow, the station lines inevitably deviate significantly from it in places. The derivatives of cubic-spline interpolants to the two surface profiles are shown in Figure 5. Because there is some noise in the levelling data, the flow is not strictly between survey stations, and surface movement is unlikely to be parallel to surface slopes associated with length scales much shorter than the ice thickness (approximately 2000 m), residuals of up to 0.5 m are allowed in the interpolant. This has an effect similar to smoothing over a length scale comparable with the ice thickness. Slope error at the survey stations is estimated at approximately ± 0.001 rad.

The slopes along the “centre lines” are roughly normally distributed, with a mean of-0.003 rad and standard deviation of 0.005 rad. These values are used as estimates of slope and slope error for tie-points away from the centre lines but still on the ice stream itself. Slopes at tie-points on those parts of transverse lines B and C where ice on Fletcher Promontory flows towards Rutford Ice Stream (Fig. 5; Reference Frolich, Vaughan and DoakeFrolich and others, 1989) are derived from a single linear regression for each double line. Slope-related vertical movement is assumed negligible on the Carlson Inlet side of Rut ford Ice Stream margin (transverse line A).

These estimates of vertical velocity are combined with the results of the horizontal survey reduction to provide three-dimensional velocity vectors and associated error estimates in standard Earth-centred Cartesian coordinates.

Velocity error covariance

Associated with the reduction of the ground surveys of 1978-80 and 1984-86 are uncertainties whose variance and covariance are integral to the results. The uncertainties in horizontal movement over 3 or 6 d are significant--equivalent to a few radar wavelengths in extreme cases — and unequal. Although generally smaller than the uncertainties associated with vertical movement along unknown slopes, they dominate towards the ends of the transverse lines of the survey networks. This is most true along transverse line B and the eastern half of transverse line C (Fig. 2b), where surveys spanned only 6 weeks rather than a whole year and uncertainties reach 20 mm d⁻¹. Velocity errors from the GPS-based surveys are by comparison small, on the order of 0.1 mm d⁻¹, with covariances small enough to be neglected in this context.

The lowest-order contribution to Oρ_iρ_j from errors in estimating the velocity is given by

The time separation t is assumed to be known exactly. The GPS measurements are presented directly as material velocities in the Earth-centred Cartesian system. For the older surveys, only horizontal velocities are presented directly. Local vertical movement is partitioned into the product of the horizontal velocity magnitude and the surface slope in the direction of flow (over an appropriate length scale), and an accumulation-related component of motion relative to the surface. In a righthanded Cartesian system in which the horizontal velocity components of tie-point i are u _i1, and u _i2 the vertical movement is

where η_i is the surface slope, and the associated error is

Somewhat unrealistically, each surface slope is treated here as an entirely independent measurement, with an error uncorrelated with any other. The uniform accumulation rate estimate used is equivalent to a near-surface downward material velocity of u ^a = 0.5 m a⁻¹ and independent of any other measurement. In principle, however, velocity errors in all three directions inherit some covariance from the survey reduction.

Fig. 5. Elevations and slopes along two main lines of the 1984-86 survey.

The error variance matrix V

When an adequate baseline model (e.g. linear or quadratic in azimuth) is used, residuals may be expected to arise from the sources of error already outlined. Conceivable sources of significant error that arc ignored here include variations in tropospheric refractive index (Reference GoldsteinGoldstein, 1995), to which some of the comments concerning range datum error arc relevant, and non-random SAR processing errors. When potential tie-points in regions of highest shear are discarded, it turns out that location errors are generally smaller than those due to phase noise and to errors in estimating the velocity contribution. It could be argued that all of these contributions would tend to be higher in regions of moderate shear, but they are assumed to be independent of each other. To preserve the covariance in movement errors, the variance covariance matrix V is formed as the sum

with each element estimated as described above.

Optimisation

The system of non-linear observation equations whose unknowns arc the baseline parameters is generally overdetermined. Minimising W with respect to the five (or more) parameters is straightforward. The method here is to repeatedly solve the quadratic optimisation problem obtained from a pseudo-linear form of Equation (2) in which terms (B · D)/r and D ²/r are neglected. For a baseline varying linearly with azimuth ζ this is

where

Optimal values for and Δ_C are found through iteratively applying the routine E04NCF from the Numerical Algorithms Group numerical library. B ² and Δ_C ² are obtained from the previous iteration or set to zero initially. The optimisation may be left unconstrained, with the baseline parameters initially set to zero, or linear constraints may be imposed to test or reflect confidence in values obtained from a precise orbit model. The random errors in the ζ, r and θ cause a slight bias in the optimal baseline parameters (Reference Snedecor and CochranSnedecor and Cochran, 1989) which can be neglected here since the errors arc much smaller than the variation in the parameters represented in the observations. As long as the ultimate object of the calibration is to find the look angle and line-of-sight movement as functions of range and azimuth, such biases in optimal baseline parameters are of little consequence.

For a linear baseline model, the optimisation routine requires the vector of first-order partial derivatives of W with respect to the unknowns and the Hessian matrix of second-order partial derivatives of W. These are constructed, according to Equation (3), from the coefficients in each observation equation and the inverse of V, which is easily obtained from a matrix guaranteed symmetric and positive definite. The solution effectively converges within two iterations. For a linear baseline, an estimate of the baseline parameter error matrix can be obtained from the inverse of the Hessian at the solution as

If only unconstrained optimisation is required, then the MATLAB routine Iscov provides a simpler alternative that avoids inverting V.

Results

Phase B (1992) interferogram

There are 148 tie-points whose phases have been unwrapped in the Phase B interferogram: 139 from the 1978 80 and 1984-86 surveys; the Carlson Inlet GPS tic-point from 1995-96; and eight points from the 1994-96 GPS survey where it crossed the margin between Rutford Ice Stream and Carlson Inlet (Fig. 2b).

Table 3 shows results obtained with various combinations of baseline model and ground control. The baseline parameters (B _H and B _V) were assumed to vary either linearly or quadratically along the orbit, with optimal values at the scene centre given as components perpendicular and parallel to the line-of-sight. A measure of the quality of fit (though not necessarily of the accuracy of the baseline parameters themselves) is provided by the rms weighted residual. An order-of-magnitude variation in weight ensures that the rms unweighted, i.e. metric, residual is not minimised and should be interpreted with caution.

The first entry in Table 3 is the result obtained when all 148 tie-points are used and the baseline is constrained to be the best fit to the Delft precise orbit model that varies linearly with azimuth. The optimal phase datum produces a rms weighted residual greater than 2, i.e. the residuals arc on average more than twice those expected. The rms metric residual is equivalent to roughly two fringes in Figure 2a and corresponds to an elevation error of approximately 180 m, or an error in horizontal velocity towards the ground track of 9 m year.

Table 3. Optimal baselines and rms residuals {weighted and metric) for Phase B interferogram using various combinations of orbit model and ground control

If the error covariance matrix has been reliably estimated, then the Delft orbit model appears to be inadequate for the present purpose. This is confirmed, by a rms weighted residual close to unity and a halving of the rms metric residual, when the optimisation of the baseline parameters is left unconstrained. The importance of optimising the baseline parameters and including vertical movement is further illustrated in Figure 6, where (he metric residuals for the tie-points that lie along the Rutford Ice Stream centre line are shown for three cases: when the baseline parameters are free; when they are forced to be the best linear lit to the Delft precise orbits; and when vertical movement is ignored.

Although the rms weighted residual is not much greater than 1, the distribution among the tic-points shows a significant correlation with position (Fig. 7), particularly where transverse arm C (from the 1984-86 survey) crosses the shear margin and ascends Fletcher Promontory. Despite a weighting that reflects die larger errors expected at this network extremity, weighted residuals here reach values of 2-3, with a consistent sign, and metric residuals exceed 10 cm. Plausible explanations include a secular change in ice movement between 1986 and 1992 (the equivalent elevation change is large enough to be ruled out), a non-stead} element to the ice motion, and a failure of the orbit model. The required magnitude of velocity change of up to 30 m year seems unlikely in areas known to have ice-sheet-like, rather than ice-stream-like, dynamics. The small weighted residual associated with the lone lie-point on Carlson Inlet at the eastern edge of the interferogram suggests that the linear orbit model is adequate, as does the failure Γ a free optimisation of the parameters describing a quadratic variation of the baseline with azimuth to significantly reduce the rms residual. The accuracy of the original survey must therefore be questioned. What may distinguish these particular points on transverse line C is their separation by a heavily crevassed margin from the survey control along the ice-stream centre line and their survey span of only G weeks rather than 1 year. Difficulties encountered with the survey reduction where transverse line C crossed the margin did result in some observations being discarded as outliers, so it is plausible that a lack of redundancy allowed a significant systematic error to go undetected.

Residual error is also unexpectedly high on the arm of transverse line A that crosses the shear margin onto Carlson Inlet. This area, where the 1984-86 and 1994-96 surveys and the Phase B and Phase D interferograms overlap, is the subject of the next section, in which possible differences in velocity field between the various epochs arc explored.

Those tie-points for which vertical velocity errors arc small—because vertical velocity has been observed directly, the horizontal velocity is small and relatively error-free, or slope estimates are relatively accurate (i.e. along the centre tines) — form a logical subset of more accurate ground control. Although the tie-points from the 1994-96 GPS survey fall into this category, they are excluded until their possible conflict with the earlier survey is investigated in the next section. This leaves a set of 82, comprising the isolated tie-point from Carlson Inlet and those from the 1981-86 survey along the ice-stream centre lines and at the end of transverse line A on Carlson Inlet. A free optimisation with this ground control also produces a rms weighted residual close to 1. With tie-points of relatively poorly known movement now discarded, the rms metric residual is less than 1 cm (equivalent to an elevation error of approximately 30 m, or an error in horizontal velocity towards the ground track of 1.5 m year⁻¹).

Fig. 6. Residuals along the main feast) line of the 1984S6 survey. Values are for not allowing for slopes (dotted line), for an unconstrained optimisation (solid line) and for the Delft orbit (dashed line).

Fig. 7. Map of residuals. The orientation is the same as Figure 2, bill with axis units in pixels. (The x axis extends in the negative direction to accommodate the arrows marking the residuals.)

Free optimisation with either ground-control set produces baseline parameters outside the normally accepted error of the precise orbit models. The optimal baselines also differ between ground-control sets, particularly in their components parallel to the mean look angle, revealing a much shallower optimisation for B than for B. Whether these results reflect an ill-posed geometric problem or significant systematic errors such as those discussed earlier is unclear, and probably unimportant if the quality of calibration and not the baseline parameters themselves is of interest. It is interesting, however, to note how the weighted residual variance varies for the more accurate ground-control set as bounds on the baseline parameters centred on the Delft orbits are relaxed (Fig. 8). Within about 1 m of the Delft orbits.

W decreases steadily. Beyond this the optimisation is rather shallow, with a small reduction in W and little change in the final interferogram pattern, even though the optimal baseline is some 10 m away.

As an exercise, a set of live tie-points, the minimum required to determine a phase datum and a linear variation in the baseline, were chosen on the basis of minimising the associated error estimates while still covering a range of azimuths and slant ranges. These were the Carlson Inlet GPS point, one point from the Carlson Inlet side of the shear margin on transverse line A of the 1984-86 survey, and three points on fairly level surfaces al or close to the Stations along the centre line used as control in that survey. When the larger observation sets are forced to the baseline parameters obtained using these five tie-points, the results in terms of the rms residuals are marginally worse than those using free optimisation, but far better than those obtained from the Delft orbits. This is encouraging for those who would minimise fieldwork, but the tie-point selection, while not arbitrary, could clearly have been different, with potentially different results. To conclude that tie-point redundancy is unnecessary would be unwise.

Fig. 8. Variation of weighted residual variance as a function of distance from Delft orbit.

With the Phase B (1992) interferogram and more accurate ground-control set, a fit generally in line with the estimated errors and without obvious problem areas suggests no evidence for significant change in ice movement along the Rutford Ice Stream centre line or over Fletcher Promontory between 1978 and 1992. The area where the Phase B interferogram and the 1984-86 survey are overlapped by the Phase D interferogram and a roughly contemporary ground survey may be an exception that is amenable to closer examination.

The Rutford Ice Stream-Carlson Inlet shear margin

A suspicion of velocity changes in the ice-stream margin was aroused when the 1994-96 GPS survey results were included in a calibration, described above, of the 1992 Phase B interferogram. Inconsistencies were revealed in the form of residuals greater than the known observation errors. To investigate this, the natural choice of reference is the calibration of the Phase B interferogram with the smaller, more accurate, set of tie-points, from which the 1994-96 GPS tie-points were excluded. Additional data sources are the Phase D (1994) interferogram and that part of the 1994-96 GPS survey that continued upstream along the ice-stream “centre line” from the end of the 1984-86 survey. The Phase D interferogram partly overlaps (he Phase B interferogram and includes observation sites from both the 1984-86 survey and the 1994-96 GPS survey. The 1994 Phase D interferogram was calibrated with a total of 35 tie-points: 27 from the 1984-86 survey and eight from the 1994-96 GPS survey. With a linear baseline model, free optimisation produces a rms weighted residual of 1.44 for the Phase D interferogram. Some other results are shown in Table 4. Excluded from the calibrations were those tic-points sited in or near the 5 km wide zone of greatest shear. Of these, 11 tie-points from the 1994-96 GPS survey and four from transverse arm A of the 1984-86 survey were unwrapped in both interferograms.

The excluded tie-point observations and the calibrated interferograms arc compared in Figure 9a. Residuals are expressed in terms of the weight the observations would have been given if included in the calibrations. The observations from the 1994-96 GPS survey produce residual weights as high as 8 in the Phase B interferogram and up to 3 in the Phase D interferogram. On the Carlson Inlet (out-er) side of the shear margin, however, residual weights have magnitudes less than unity in both interferograms, i.e. the ground truth and interferograms arc consistent. The 1984-86 survey observations produce residuals with a similar pattern: low on Carlson Inlet (these were included in the cali-bration), high near the middle of the margin and decreasing towards the centre of the ice stream.

Table 4. Optimal baselines and rms residuals for Phase D interferogram using various orbit model

Fig. 9. (a) Weighted residuals across shear margin: solid line is for 1992 interferogram, dashed line for 1994 interferogram; crosses are 1994-96 GPS stations, circles 1984-86 survey stations. (b) Relative velocity (with respect to 1992): 1994-96: solid line with crosses, error ±1 m year⁻¹; 1984-86: solid line with circles, error ±3 m year⁻¹; 1994: dashed line with crosses, error ±2 m year⁻¹; 1994: dashed line with triangles, error ±6 m year⁻¹. (c) Velocity profiles: solid line with crosses marks GPS stations occupied in 1994-96; dashed line is inferred 1992 velocity.

The covariant nature of the 1984-86 survey results means that no subset of them is strictly independent of any other, but any bias reducing the residuals from the excluded tie-points is expected to be small. The 1984 86 and 1994-96 surveys arc linked only through a common (rock) station. whose velocity is known exactly, so their velocity errors are independent. The patterns of residuals in Figure 9a are therefore taken to reflect real changes in velocity between the various epochs. The possibility of equivalent, hundreds of metres, changes in elevation is discounted.

Metric residuals may be expressed as changes in the components of horizontal velocity perpendicular to the satellite ground tracks if vertical velocities are assumed not to change. In Figure 9b the horizontal velocity for the shear margin lie-points is shown relative to that in 1992 for the nominal epochs 1984 86, 1994 and 1994-96. For the 1984-86 and 1994-96 epochs the relative velocity is obtained directly from the respective sets of residuals in the 1992 Phase B interferogram (curves c and a in Figure 9a). For 1994 it is found by differencing the residuals for each tie-point in the Phase B and Phase B interferograms (curve c minus curve d, and curve a minus curve b in Figure 9a). In allowing for the different orientation of the two interferograms, it is assumed that velocity changes are parallel to the velocities themselves. For comparison, the horizontal velocity profile through the shear margin measured by the 1994-96 GPS survey is shown with parts of the inferred 1992 profile in Figure 9c.

Concluding Remarks

Realising INSAR’s enormous potential requires various degrees of sophistication, Interferograms can be created without ground control, using a precise orbit model to estimate integer offsets between coherent images and the same precise orbit model to estimate the interferometric baselines. Such interferograms generally contain a mixture of topographic and kinematic information, but can still be extremely effective in identifying and roughly locating features such as grounding lines and boundaries between fast and slow flow. When these interferograms are accurately calibrated they can also be useful for monitoring ice sheets for change. To assess the significance of putative changes between interferograms and between interferograms and terrestrial observations requires a rigorous approach to error analysis such as that attempted here. Redundancy in ground control is essential, even when, as in this work, the variation of baseline with azimuth can be assumed linear.

At the other extreme of sophistication are accurately calibrated differential interferograms from which elevation models and three-dimensional surface velocity fields can be obtained. These require extensive image sets that will be far from ubiquitous. Less complete results, provided they are properly calibrated, may still aid quantitative work such as ice-sheet modelling. For example, interferograms that provide only line-of-sight velocity or a single horizontal velocity component could prove extremely useful as constraints in applications using optimal control methods (e.g. Reference MacAyealMacAyeal, 1993).

There remains the question of ground truth; how much is required and how accurate it needs to be. The answers depend on the objects of a particular application. It is unlikely that changes in the margin of Rutford Ice Stream could have been identified without ground control in the area to optimise the baseline parameters. No attempt was made, however, to use natural indicators of ice motion, such as ice divides or longitudinal flow features, as control.

The significance of vertical motion distinguishes interferometry from most other remotely sensed sources of velocity information. This is sometimes an advantage, but not usually on grounded ice, where knowledge of surface slopes in the direction of flow, and possibly local accumulation rates, may be essential in order to obtain accurate horizontal velocities. The appropriate length scale over which to measure surface slopes is locally determined by the ice dynamics. For this reason, a powerful application of interferometry may be in identifying areas of glaciological significance where current flow is not parallel to current surface slopes.