Introduction
This paper presents an overview of the concepts and methods of differential radar interferometry applied to valley glaciers. This is done From an algorithm point of view, as the theory is well established (e.g. Reference JoughinJoughin, 1995; Reference Joughin and KwokJoughin and others, 1996a, b) with the intent of showing the progression from SAR images to a surface-velocity vector field. We also show how the constrained motion of a valley glacier can be used to resolve an inherently ambiguous component of this velocity and we discuss sources of error. In the context of this discussion of methods, this paper presents the results of an interferometric study of Bagley Icefield, showing the acceleration of the ice subsequent to the surge onset of Bering Glacier in 1993 (Fig. 1; see also Lingle and others, 1993; Reference MolniaMolnia, 1993).
The 1991 launch of the first European remote-sensing satellite (ERS-1) was followed by R. M. Goldstein's recognition that repeat-pass spaceborne radar interferometry (SRI)could be used to measure the movement of polar ice sheets (Goldstein and others, 1993). SRI makes use of the comparison of the phase of complex-valued SAR images from repeal orbits. The spatial separation between two satellite orbits functions like a stereo-vision optical base line for resolving topography on a scale of meters. A separate aspect of complex image phase records surface translations, such as ice motion, as fractions of the radar carrier wavelength on a scale of centimeters. Goldstein's original innovation made use of fortuitous satellite passes which were spatially coincident, with separations of only a few meters. These gave poor topographic resolution and so emphasized the ice-motion signal. The technique has subsequently been generalized to obtain results from more typical orbit pairs with larger base lines (hundreds of meters), giving both surface motion and topography (Reference Gabriel, Goldstein and ZebkerGabriel and others, 1989; Reference JoughinJoughin, 1995; Reference Rignot and JezekRignot and others, 1995). The generalized technique, referred to as differential SRI (here DSRI), has advanced the measurement of surface movement, deformation and topography of ice sheets and glaciers (Reference Joughin and KwokJoughin and others, 1996c; Reference Joughin and KwokKwok and Fahnestock, 1996; Reference RignotRignot, 1996; Reference RignotRignot and others, 1996; Reference Mohr, Reeh and MadsenMohr and others, 1998; Reference Joughin, Kwok and FahnestockJoughin and others, 1998).
The launch of ERS-1 was also followed by a major surge of Bering Glacier in the Chugach—St. Elias mountains of south-central Alaska. Bering Glacier together with its accumulation area, Bagley Icefield, and associated glaciers (notably Jefferies Glacier) are shown with arrows indicating flow directions in Figure 1. With an area of 5200 km−2, Bering Glacier is the largest glacier system in continental North America (Reference MolniaMolnia and Post, 1995). A full description of the ice dynamics of the 1993-95 Bering Glacier surge is incomplete but much of what is known has been derived from ERS-1 SAR images. As shown by the time line in Figure 2. the first evidence of the surge visible in a SAR image became apparent in an ERS-1 scene from April 1993 that shows surface disruptions approximately 22 km up-glacier from the Bering Glacier terminus within the lower ablation area (Reference RoushRoush, 1996). The singe front subsequently propagated down-glacier with speeds of up to 100 md−1, reaching the terminus in August 1993. Typical ice velocities varied from 10 to 20 md−1 (Reference Fatland and LingleFatland and Lingle, 1994). The surge also propagated upstream to the East/West Bagley Icefield confluence and farther up-glacier into the East Bagley Icefield (Fig. 1). The first stage of the surge ended in August 1994 with an outburst flood al the eastern edge of the terminus. The second stage of the surge began at an indeterminate time after this, with terminus advance observed in April 1995. A second outburst flood in September 1995 marked the end of the surge.
From 1992 through 1995, the Bering Glacier system was periodically imaged in the 100km wide swath of the ERS-I C-band (5.7 cm wavelength) SAR. For most of this period, the repeat-orbit interval was 35 days and the Bering Glacier system was observed from several different orbital tracks. The 35 day repeat interval is too long for useful SRI, because surface changes introduce decorrelation noise (Reference Zebker and VillasenorZebker and Villasenor, 1992), but SAR-amplilude images arc useful for tracking the progress of the surge through changes in large features such as crevasse fields and medial moraines. From January through April in both 1992 and 1994, ERS-1 was placed in an orbit (Ice I and Ice II mission phases — Figure 2) which repeated itself to generally better than 250m every 3 days, imaging a 100 km length of Bagley Icefield. Ice velocities and surface stability on Bagley Icefield were suitable for DSRI, opening the possibility of highly detailed and spatially continuous measurements of flow and topography both prior to the surge onset in the winter of 1992 and while the surge was in full progress during the winter of 1993-94. By contrast, DSRI has proven unfeasible farther downstream on the rapidly moving ice of Bering Glacier, both before and during the surge. The next section n provides a DSRI processing overview, starting with the synthesis of interferograms from SAR image pairs and continuing to the generation of topography and velocity. These are first presented in the simple cases of (i) fixed topography and a moderate satellite base line (for ERS-1, <300 m), and (ii) a moving surface with zero-length interferometric base line. DSRI is then described in the algorithm sense for the general problem of a continuously moving glacier surface imaged with a moderate interferometric base line. The algorithm is used to measure the pre-surge to surge velocity increase on West Bagley Icefield and the sources of error are discussed.
Bering Glacier, Bagley Icefield and associated glaciers, south-central Alaska.
Time line showing progression of the Bering Glacier surge (above horizontal bar) and ERS-1 observation phases (below bar). SRI data were obtained from 3day repeat-orbit periods only.
Processing Overview
This section describes the main SRI arid DSRI processing steps, summarized in Figures 3-10. An excellent resource on the derivation of the equations given here and furl her details of DSRI implementation in practice can be found in Reference JoughinJoughin (1995). This algorithm overview presents a general description of DSRI processing in order to provide a basis for evaluating the glaciological utility of the results. SAR interferometry is a relatively simple image-processing concept made complex in practice by the exigencies of the data. Figure 3 shows the interferometric imaging geometry in a plane approximately perpendicular to the instrument flight direction. Figure 4 is a flow chart for the generation of a Digital Elevation Model (DEM) from two SAR images as in the case of an unglaciated region in which the interferometric phase signal represents topography only. An example scene near Bagley Icefield, used later to determine the interferometric base line, is shown in Figure 5.
The first step in forming an interferogram is the coregis-tration of two images to a fraction of a resolution cell (Fig. 6). Resolution cells are interchangeably thought of as small blocks of pixels defined by the resolution of the SAR and the corresponding surface region in the physical scene. A resolution cell is generally several pixels in size but it is often convenient to blur the distinctions between pixels, resolution cells and corresponding patches of ground. Once two images are coregistered, the complex phases are differenced while the amplitudes are retained to form a new complex-valued image called an interferogram. The spatial separa-tion between the two imaging passes in the plane perpendicular to the spacecraft flight path is described in terms of a line-of-sight parallel component Bp and a line-of-sight perpendicular component Bn (Fig. 3). The latter is called the interferometric base line and it is analogous to an optical base line in stereo imaging. For ERS-1, both Bp and Bn are frequently less than 250 m for repeat orbits separated by 3 days, which represents good orbit control from the point of view of satellite-orbital mechanics and glacier DSRI.
Fig. 3.
SRI processing flow chart showing the generation of a digital elevation, model (DEM) under circumstances in which no glaciers are present to introduce additional translation phase into the topographic phase signal.
The intrinsic phase coherence of an interferogram is evaluated as a scalar field of autocorrelation values ρ for a small cluster of N pixels:
where Xi and yi are (complex) signal values from source images 1 and 2, respectively, and y*, i is the complex conjugate of yi . Since an interferogram pixel phase is the difference of the source-image pixel phases, interferogram pixels can be represented as Xiy* i Interferogram coherence varies throughout the image and is used as a quality guide in subsequent processing.
Having removed the random-phase variation throughout the individual SAR images by making an interferogram, the residual phase signal ψ is given by
where
and where k is the radar wave number 2π/λ, λ is the radar carrier wavelength (5.7 cm for ERS-1), R is the distance from the SAR to the scene center and ? is the image-center radar-incidence angle (Fig. 3). ΔR, the radial distance change of a resolution cell due to its spatial translation between imaging passes (Fig. 7), drives the translational phase signal ψtrans-. The topographic elevation of the resolution cell relative to a reference elevation is given by Δh, and only the topographic phase term ψtopo contains an explicit dependence on both the topography and the perpendicular base-line parameter Bn. ψnoise is related to the interferometric signal coherence ρ and is neglected for now, leaving the three terms ψgeom + ψtopo + ψtrans. These three phase signals must be separated in the general case of a moving glacier with surface topography.
Case 1: topography without surface motion
We first consider the case of a fixed surface topography with no glaciers, a moderate interferometric base line, │Bn│ <300 no and no phase noise. The total phase is then
since ψtrans = 0. Reference JoughinJoughin (1995) gave a derivation of the topographic phase signal and the conversion from SAR slant-range images to ground range. The steps for producing a ground range DEM are presented in outline in Figure 4 with some elaboration as follows:
1. Image acquisition is a matter of opportunity subject to restrictions of SAR orbital paths and other flight-agency operational issues. Suitable image pairs are reregistered to within a pixel using immobile features. These image pairs are separated by a time interval Δt which is generally short (days) to minimize signal decorrelation. Sub-pixel resampling is a “second-order” process which can significantly improve interferometric phase coherence as shown in Figure 6 (sec also Appendix).
Interferograms of a non-glacier-covered valley near Bagley Icefield. (a) Phase signal prior to removal of flat-Earth phase ψgeom.(parallel bands). (b) Phase signal after removal of flat-Earth phase, (c) Perspective rendering of figure 5b showing topographic relief with slant-range layover still present.
Coherence improvement in interferometric phase on jefferies Glacier, a tributary of West Bagley Icefield, (a) Glacier surface from two source images coregistered to within 1 pixel. (b) Signal gain after sub-pixel coregistration adjustment.
Two SAR data acquisitions, with resolution cells 1 and 2 moving during the intervening period to give radial distance changes ΔR 1 and ΔR 2.
2. An interferogram is generated by subtracting the phases of the two images on a pixel-by-pixel basis. In this case, the resulting phase signal will consist of the two terms given in Equation (6) plus a spatially varying noise contribution. Pixel phase is often represented visually using color scaled in intensity by pixel amplitude. This shows the general phase characteristics of the interferogram, as well as identifiable features.
3. Both terms of Equation (6) are dependent on the normal base-line component Bn, analogous to a stereo-vision in-terocular base line. It is thus necessary to make an initial estimate of B n by some means in order to analyze the phase signal. Image misregistration can be used for this estimate but it is generally simpler and more accurate to consult an appropriate flight agency on-line database (e.g. for ERS-1, at ESRIN: http://gds.esrin.esa.it/). It is sometimes necessary to refine iteratively the initial baseline estimate as shown in Figure 4 (see also Appendix).
4. The geometric-phase contribution to the total phase is a regular, nearly linear ramp, modulo 27π, which covers the entire interferogram (apparent in figure 5a). This is primarily due to a dependence on the base line B n which typically varies nearly linearly by only a few meters along-track from one end of an image to the other and by its dependence on the incidence angle ? which varies across-track. ψgeom is the phase signal due to the imaging geometry that would be present in the interferogram if the scene were completely flat (no surface topography; only the Earth's curvature). Consequently, ψgeorm does not contain any useful information and can be subtracted from the interferogram phase signal, leaving only residual topographic phase (fig. 5b).
Often, “flat-Earth ramp” phase subtraction is followed by low-pass filtering that reduces signal noise and produces pixels with near-unity aspect ratio. Filtering the phase signal is often necessary to some degree and may be implemented with a locally adaptive scheme which uses phase coherence and/or signal amplitude as a guide.
5. Phase unwrapping is shown in both one and two dimensions in Figure 8. It consists of the addition or subtraction of integer multiples of 2π to each pixel as necessary to eliminate phase discontinuities and is described in more detail in the next section. In noisy areas, phase unwrapping is quite problematical; such regions are often cut from the final data product or interpolated between regions of good signal (Goldstein and others, 19?Ω; Reference PrittPritt, 1996; see also Appendix).
6. An unwrapped phase signal will have a precise relationship between unwrapped phase values and corresponding surface elevations. These can be compared to a suitably selected set of tie points distributed as widely as possible throughout the scene. If the fit to the tie points is poor, then we conclude that the base line used in step 4 is in error and must be refined, for example, using the linear least-squares approximation described by Reference Joughinjoughin (1995). The base-line estimate is refined iteratively to an acceptable tolerance before proceeding further.
7. It is useful to transform SAR image products from their processed form to a “ground-range” map projection in which each pixel represents the same size surface area. Normally, single-look complex (SLC) SAR images are processed to “slant-range” format, in which each successive pixel in the cross-track direction represents a quantization of some number of meters in radial distance from the SAR antenna. For a flat-Earth surface, the increasing angle of incidence across the swath will result in successively smaller ground-range projections of each slant-range pixel, so that there is not a simple constant relationship between slant-range and ground-range pixels. The problem of transforming images to ground range is further complicated by lay-over distortion produced by scene topography. Slant-to-ground range rectification using interferometrically derived surface elevations produces a ground-range image in which all pixels have the same map-projection dimensions (Reference Joughin and KwokJoughin and others, 1996b).
8. Elevation uncertainties can be estimated directly from phase variance σψ using Equation (6) to give
The end result of the process outlined in Figure 4 is thus a DEM (here designated hjj, where h is the surface height at pixel (i,j)), consisting of an array of pixel elevations with respect to a reference datum, with associated errors. The (rectified) SAR amplitude image can be draped over the DEM to show feature correspondence, Typically, the derived DEM has an accuracy of ~5-20 m (Reference Zebker, Werner, Rosen and HensleyZebker and others, 1994), dependent on the availability of ground control.
Phase-unwrapping in one and two dimensions. Jagged nature of center-line phase plot is indicative of dala noise rather than variations in velocity.
Case 2: Surface motion with Bn
In this situation, we consider a scene containing moving resolution cells imaged twice from the same location over some time interval Δt, so that Bn = 0, giving Ψgeom=0 and Ψtopo =0. When a particular region on a glacier undergoes a "rigid" (locally non-deforming) translation Δ x relative to some other fixed region in the scene, there is generally a corresponding change in radial distance ΔR from that résolution cell to the SAR as shown in Figure 7. Spaceborne SAR has a coarse-scale two-dimensional map view) resolution of the order of tens of meters but, at the same time for each resolution cell the SAR measures radar-echo phase as a fraction of the carrier wavelength. Thus, it also acts as a range-sensitive instrument on a scale of centimeters, but with no angular spatial resolution at this scale. As a result, surface translations result in a globally ambiguous "fractional wavelength" interferogram phase signal which indicates only small radial distance changes at each resolution cell, relative to nearby resolution cells. For example, if over 3 days a resolution cell moves 9.1λ closer to the SAR, then the round-trip distance change is 18.2λ which would only be apparent in the phase signal as a fractional shift of 0.2λ or 0.2 x 2π rad in terms of pixel phase. By itself this would be indistinguishable from a radial distance change of 8.1λ, 7.1 λ, etc., so it is necessary to have a gradually varying phase signal with locally small (< π) phase increments between adjacent pixels to allow unambiguous integration of ΔR. This process of integration, starting from a known-value reference location (e.g. a fixed surface) and summing phase increments to keep track of total accumulated phase is referred to as "phase unwrapping" (Reference Goldstein, Zebker and WernerGoldstein and others, 1988). "Unwrapped" phase is no longer constrained to [—π. π] and is related, in this example, to the radar line-of-sight radial translation of surface resolution cells by the direct proportionality ψtrans = 2kΔR
Case 3: a moving surface with topographic relief and non-zero base line
Generalizing Figure 7 to a non-zero base line will add geometric and topographic phase signals ψgeom and ψtopo, both dependent on the normal base line Bn, to the translational phase ψtrans. This represents the general case for interfero-metric analysis of glaciers and ice sheets. It is of interest to distinguish the scales of translational and topographic phase, as SRI systems generally give topographic resolution on a scale of several meters (Reference Zebker, Werner, Rosen and HensleyZebker and others, 1994), whereas the fractional-wavelength translation phase measurements give translation resolution on a scale of centimeters. Furthermore, only one of the three translation-vector components is given by the translation phase; that is, only the component of ground motion in the radar line-of-sight direction is obtained. The remaining two components must be derived by other means, as discussed in the next section.
Assuming Bn has been adequately determined and the geometric phase removed from an interferogram, the remaining phase is
The separation of these two terms exploits the inherent differential scaling of the two types of information using multiple image pairs; this is shown in the flow chart in Figure 9 and schematically in Figure 10.The phase superposition and separation is shown in Figure lla-d on West Bagley Icefield. The procedure is as follows:
1. Four source images are used to generate two interfero-grams with two respective "normal" base lines B 1 and B2 . Since the random scattering phase is (ideally) eliminated in these source interferograms, there is no restriction on the time interval between pairs. However, subsequent interferogram comparisons must account for possible physical differences in the imaging conditions. For example, a glacier interferogram from winter compared with another from a different season will have phase signals reflecting seasonal differences in surface speed, violating the "constant velocity assumption" (see below).
2. It is assumed (for now) that the surface motion between images 1 and 2 is the same as between images 3 and 4, i.e. we assume that the glacier moves with constant velocity so that ΔR 1 = ΔR 2 = ΔR. To maximize the probability of this, it is best for the two image pairs to be closely spaced in time. T he result of failure of this assumption is discussed below. Differencing the phases of interferograms ψ1 and ψ2 thus cancels the translation phase component to produce a differential interferogram ψ12in which the phase represents topography only, scaled by the differential base line B 1-B 2.
where a — 2k/R sin ?. Topography-only phase ψ12 is shown in Figure l1c.
DSRI processing flow chart, showing the generation of glacier surface-velocity field and, in passing, a DEM, under the assumption of constant glacier velocity.
3. The isolated topographic phase can be unwrapped into an unconstrained scalar field and rescaled to any desired base line, particularly to B 1 and B 2. That is, ψ12 is rescaled to match the ψtopo component of interferograms ψ1 and ψ2, so that another phase difference can be used to isolate the translational phase signal:
4. Either ψ1-12, ψ2-12 or some combination (such as their average) can be taken to represent the isolated translational phase ψtrans. Because the translation phase represents relative velocity, a suitable fixed point must be chosen as the zero-velocity starting location for phase unwrapping, for example, a pixel on a mountain near the glacier margin. Because the topographic phase has been removed, all pixels on fixed features should have constant phase, regardless of pixel elevation. This dis-tinction is apparent in a comparison between figure 11a and b (which include topographic phase) with Figure l1d from which the topographic phase has been removed.
DSRI processing schematic illustrating the differentiation of topographic and translation phases. B1 and ?2 are the interfermetric base lines (normal component, Bn)for the two source image pairs. ΔR is the radial distance change to be derived from the translation phase signal ψtrans.
5. The radial distance change ΔR given by ψtranscan be converted to a surface-velocity vector by projection of the (unwrapped) ψtrans phase into an appropriately chosen flow direction (next section).
6. As with the DEM generation process, local-phase variance obtained from the data can be used (Equation (5)) to estimate uncertainties in ΔR which can be translated into velocity uncertainties.
Conversion of Radial Distance Change to Surface-Velocity Vector Field
Figure 12 shows a coordinate system with origin located at the center of a resolution cell, z axis defined vertically upwards and x axis lying in the local horizontal plane (not the glacier-surface plane) pointing in the SAR cross-track direction, which in the case of Bagley Icefield is similar but not identical to the glacier-flow direction. When the side-looking SAR is imaging this resolution cell, the SAR will have coordinates (?s,0,Zs). Between the two passes, the resolution cell moves a distance S along a flow-unit vector u The derivation of u is problematical but below it will be derived from a two-dimensional flow-unit vector f defined in the xy plane. The ice-flow vector Su results in a change in radial distance of ΔR along a line-of-sight unit vector r =(rx , 0,rz ), where rx/rz=Xs/Zs. In fact, Su can be written as the sum of the vector ΔRr and some vector A perpendicular to :
In lieu of fortuitous multiple-direction SAR data acquisitions (Reference Joughin, Kwok and FahnestockJoughin and others, 1998; Reference Mohr, Reeh and MadsenMohr and others, 1998), this method of velocity calculation necessitates some means of determining the unit vector û for surface-parallel flow in the longitudinal direction. Using the topographic gradient (u= Vni|| is inadequate, because valley-glacier topography is dynamically supported by moving ice (Reference RaymondRaymond, 1971; Reference EchelmeyerEchelmeyer, 1983). It is more feasible to calculate Û by assuming that the glacier flows parallel to the valley walls across most of its width. In this work, the two relative components of Û in the xy plane are derived first, giving the map-plane flow-unit vector f = {f 1,f 2,0}. The z component of û is then derived by taking the surface gradient of the glacier topography in the f direction. The unit vector û = {u 1,u 2,u 3} is then normalized and the translation vector ΔRr is projected in the û direction to give S. The correction from f to û is necessary to avoid introducing a 5% error in velocity. The difficult part of this process, determining f 1 and f 2 at each glacier pixel in the image, can be done "by hand" for a particular study site. Particularly problematical is the determination of f for embay-menls, tributaries and in representing small transverse velocity components. A first-order approach to deriving f at each pixel is to draw a (smoothed) center line through the image which approximates the glacier flow direction (Fig. 13). At a particular pixel along this line, fis given by the local center-line tangent. A perpendicular transect line is extended from this center-line pixel across the glacier and each pixel of this transect is assigned the same value of f. If the center-line is curved, some pixels will he visited multiple times by this technique and others will he missed entirely, necessitating some filtering to produce a continuous and smoothly varying / vector field. A second-order improvement could model the map-plane lateral convergence or divergence of this field (in accumulation and ablation areas, respectively), which empirically seems to approach 7 10° at the transition area from center-line flow to the glacier shear margin on typical temperate valley glaciers. Errors in f and thus in velocity are discussed below.
Having determined f at a particular pixel and derived û from the local surface gradient along f, the projection geometry gives
where ΔR is determined from the unwrapped translation phase ψν^ by AR = *™ans/2fc As a rule of thumb, for a glacier moving at an azimuthal angle ? with respect to the cross-track direction under ERS-1 imaging conditions, one phase fringe (a complete color cycle) from an image pair separated by 3 days corresponds to ~ 2.5/cos?cm d−1. (Here, ? is defined relative to the cross-track axis such that it satisfies the constraint |θ|<π/2.)
Interferograms corresponding to various stages of Figure 10.(a) West Bagley Icefield interfemgram, 19-22 January 1992, prior to the 1903 Bering Glacier surge onset, Ice flows from left to right, (b) Same site, data from 4-7 February 1994, with more-central phase hands indicating post-surge-onset velocity increase,(a) and (b) also have the moderate glacier-surface topography folded into the phase signal. (c) Differential interferogram in which translation phase is removed leaving only topographic phase. Phase-color boundaries are analogous to topographic contour lines, showing a typical accumulation-area profile with glacier margins higher than the center, (d) Surface translation phase only for 1994 interfergram, after removal of topographic phase. Comparison with (b) shows more central bands present.
Using the above technique, errors in flow-direction angle ? and to a lesser extent the incidence angle ? are the biggest contributors to error in the calculation of surface velocity. Figure 14 shows resultant velocity errors in per cent for given flow direction and incidence angle errors using the approximate relationship for surface speed S as a function of unwrapped phase ψ:
Errors in ?, the local incidence angle, are caused by-errors in the estimation of glacier-surface slope; these will be smaller than flow-direction errors but are also subject to greater local variability with undulations in the glacier surface. Under good conditions, DSRI will pick out such undulations from the differential topographic interferogram. Using an upper limit for the flow-direction error Δ? of 5° and a maximum acceptable velocity error of 20%, Figure 14 gives the restriction that the glacier must flow in a direction θ< 65° from the SAR cross-track image axis. That is, single-pass DSRI velocity determination works poorly or not at all on glaciers which happen to flow close lo the along-track image axis (i.e. parallel to the SAR-flight path).
Velocity-projection geometry. Observed radial translation = R, ad hoc horizontal-plane flow direction = f and derived velocity vector = Su.
West Bagley Icefield Results
figures 15 and 16 show the acceleration of West Bagley Icefield between January 1992 and February 1994. Figure 15 shows a comparison of transect velocity profiles approximately 10 km downstream from the ice divide separating West Bagley Icefield from Sieller Glacier. The 1992 velocity transect is plateau-like in the center of the glacier, whereas the 1994 transeel velocity is somewhat more rounded. In the inset graph, both velocity curves from the center tine north to the margin have been rescaled to compare curvature with a simple theoretical curve for deformational flow (flow law, n= 3, no shape factors (Reference GlenGlen, 1955)):
Here, x is distance from center line toward the margin and X1/2 is the half-width of the icefield. The 1992 profile is assumed to represent the non-surging character of the West Bagley Icefield velocity. The departure of this profile from the curvature given by Equation (17) may be indicative of the influence of the valley shape (Reference EchelmeyerEchelmeyer, 1983) and/or some degree of laterally varying basal motion.
Determination of two-dimensional flow-direction unit vector.
Velocity error as a function of flow-direct ion error. Inset plot: corresponding errors for errors in incidence-angle estimation.
The second departure in velocity-profile curvature, that from the 1992 profile to the 1994 profile, shows that the acceleration induced by the surge was stronger at the center of the icefield than at the margins. The mechanism for this acceleration is likely to be a separate lateral variation in basal motion associated with the Bering Glacier surge (as opposed to an increase in flow due lo internal deformational). Such variations have been observed previously (e.g. Reference RaymondRaymond, 1971) and in this case are presumed to be due to the manner of coupling between Bagley Icefield and Bering Glacier. One possible mechanism is a longitudinal stress impulse imparted by the removal of restraining ice downstream as Bering Glacier surged. Relative to the margins, the central part of West Bagley Icefield is less susceptible to marginal shear stress. Another potential factor in causing lateral sliding variations is a hypothetical increase in sub-glacial water pressure acting to reduce basal shear stress, e.g. Reference RobinRobin and Weertman, 1973.
Comparison of several overlaid velocity transects from 1992 and 1994. The dissimilarity in these profiles indicates variable sliding speed, probably in response to the Bering Glacier surge. The inset plot shows both velocity curves from the center line north to the margin, arbitrarily rescaled for comparison of shape with the curvature of a theoretical (flow law, n = 3) velocity profile (solid line).
West Bagley Icefield longitudinal velocity comparison, 1992-94. (a) Velocity-only phase signalfrom 1992, with, transect location and longitudinal profile shown in white, (b) Same location, 1994. Areas which are black or solid color indicate signal drop-out due to poor coherence, (c) Center-line velocity-profile comparison, (d and e) Perspective rendering of surge-related acceleration of West Bagley Icefield, looking west from the East/West Bagley Icefield confluence. Vertical relief indicates surface speed.
figure 16a and b are translation-only differential interferograms from 1992 and 1994, respectively. Considerable changes in the phase signal in the central part of the glacier show the longitudinal acceleration in the 1994 image associated with the Bering Glacier surge, also shown by the change in the center-line velocity profile in figure 16c. A remarkable aspect of this center-line comparison is that the 1994 profile is very similar to the 1992 profile multiplied by an empirical factor of 2.7 (dotted curve in figure 16c).
Near the confluence, where West and East Bagley Icefields flow together, the center-line velocity has increased from 0.36 md−1 in 1992 to 0.95 md−1 in 1994. Differential interferometry also showed appreciable acceleration over the 1 month interval from January to February 1994. At the rate of acceleration derived from this 1 month period, the influence of the surge on West Bagley Icefield can be estimated to have lasted for 4-5 months by the time the February 1994 data were acquired. With the uncertainty in surge-onset location (symbol; Fig. 1; or possibly farther up-glacier) and surge start time (apparently between late March and late April 1993, possibly earlier), these data suggest that the speed of up-glacier propagation of the Bering Glacier surge to the Bagley Icefield confluence was of the order of 200-500 md−1. Subsequently, increased longitudinal stress gradients caused the center-line velocity of West Bagley Icefield to increase by a factor of ~2.7 over a time period that may have been as little as 4 months.
Failure of The Constant Velocity Assumption
To claim that the differential interferogram shown in Figure l1c contains only information about topography, it is necessary to assume that the glacier velocity field is time-invariant so that the translation components of the source phase signals cancel. This "constant-velocity assumption" (here CVA) is generally quite valid on ice sheets (with notable exceptions; see Reference Joughin and Kwokjoughin and others, 1996c). For temperate valley glaciers with sliding speeds highly coupled to sub-glacial hydrology, it may be best to say that, if possible, the CVA assumption should be demonstrated to be valid, for example, by deriving consistent topographic phase from a time sequence of many image pairs. Most valley-glacier CVA failures arc observed to be small in magnitude, usually a fractional part of a fringe representing velocity differences of less than 1 cm d−1. A more extreme example of CVA failure is shown in the differential interferogram (nominally "topography only") in Figure 17, in which two concentric fringe patterns or bulls-eyes appear on East Bagley Icefield, each ~4km in diameter. The constituent observation intervals are 2-5 January and 11-14 January 1994 and the differential base line B 1-B 2 is 28m. Such a phase pattern, if it were to indicate topography, would represent two steep conical hills 1600 m high.
There are other possible causes for phase anomalies of this sort, most notably some manner of atmospheric or ionospheric contamination of the propagating radar signal. The effects of atmospheric water vapor on topographic phase described by Reference GoldsteinGoldstein (1995) and by Reference Zebker and RosenZebker and others (1997) are present at low latitudes but arc unlikely to have such a strong influence during winter at 60°N. Furthermore, the nature of these bulls-eyes (and others seen in similar data) argues against atmospheric causes in general as they are cleanly localized and rest symmetrically over the hypothetical deepest part of East Bagley Icefield. In fact, we suppose that the most likely explanation for these bulls-eyes is not that they represent a variation in longitudinal sliding velocity but that they represent a local rise of the surface of ~20cm over a 3 day time interval, from 2 to 5 January. This value is consistent with observations on Black Rapids Glacier, a surge-type glacier in central Alaska, where Reference Heinrichs, Mayo, Echelmeyer and HarrisonHeinrichs and others (1996) observed annual elevation cycles of ~1 m and, more to the point, elevation changes on a scale of 10-20cm have been observed over a matter of hours in association with rapid lake-drainage events (personal communication from M. Nolan). A hydrological event of this sort, normally considered unlikely in winter, may be facilitated by the dynamics of the Bering Glacier surge (for comparison see Reference KambKamb and others (1985), regarding a surge-related event in February 1983 on Variegated Glacier, Alaska). The volume represented by this hypothetical surface uplift is 7x106m3.
The hypothesis of a surface rise is only one possible type of surface translation (albeit the most obvious one) which could account for the observed bulls-eyes. This emphasizes the importance of recognizing that subtle aspects of glacier dynamics and associated CVA failures, introduces local errors into the DSRI technique for estimation of the velocity field.
Resolving Surge Draw-Down
We expect the acceleration of West Bagley Icefield to be accompanied by a draw-down of the surface but resolving this change in surface elevation using DSRI-derived DEMs is problematical. The theoretical limit to ERS-1-bascd SRI elevation is ~5 m rms (Reference Zebker, Werner, Rosen and HensleyZebker and others, 1994). However, the effects of small base lines for these data, acting in conjunction with phase noise, scattering depth uncertainly, and the potential for CVA failure (West Bagley Icefield was accelerating in 1994) require us to revise this limit upward to 20 m rms accuracy at best. At this limit of resolution, a comparison of pre-surge (1992) and during-surge (1994) elevation interferograms failed to give conclusive evidence of surface draw-down. Taking a pre-surge center-line surface velocity of 100 m a−1 and a surface slope of 0.95°, and using a conventional flow-law calculation for temperate ice with sliding speed ranging from 20% to 80% of the total speed (no valley-shape factors since the half-width is greater than three times the probable depth (Reference PatersonPaterson, 1994)), the inferred depth of West Bagley Icefield is about 500-700 m. The time of exposure to surge influence at the West Bagley Icefield equilibrium line is estimated at between 150 and 200 days. A continuity calculation for the observed speedup shown in figure 16c gives a maximum theoretical drawdown of 5 10 m, which accounts for the difficulty in resolving such a draw-down using DSRI.
Summary
The large-scale glaciological problem of characterizing the dynamic behavior of West Bagley Icefield during the 1993—95 Bering Glacier surge is addressed by differential interferometric analysis of SAR data (DSRI). DSRI is only applicable to the study of relatively slow flow, for example, on Bagley Icefield and jefferies Glacier, as it is easily rendered ineffective by decorrelation noise from rapidly moving and deforming (i.e. surging) ice. The limitations of DSRI, including its considerable computational complexities, are compensated by the capacity of this technique to produce glacier DEMs and surface-velocity vector fields at high resolution.
Due to the abundance of stable ice-free topography, interferometric base lines for valley-glacier scenes are comparatively easy to determine and refine (compared to the relatively featureless polar ice sheets) using tie-point-iteration techniques. Base-line refinement is consequently limited by tie-point errors and decorrelation noise. The determination of precise base-line data by die various flight agencies (e.g. European Space Agency) has reduced the problem of base-line refinement. Decorrelation noise also makes the analysis of temperate valley glaciers difficult relative to ice-sheet data. This decorrelation is likely due to high temperatures and proximity to coastal weather systems in the Gulf of Alaska.
Differential interferogram (nominally topography only) of part of East Bagley Icefield with two features indicating failure of constant velocity assumption, from January 1994.
Using only a single-observation geometry, an assumption of folow parallel to the axis of the valley (rather than How down the surface gradient) should be used to obtain valley-glacier surface-velocity fields, flow-direction errors arc minimized as the direction of glacier flow approaches the radar cross-track look direction but are exacerbated in regions of complex How, particularly at glacier confluences. In addition to flow-direction errors, DSRI error sources (in the ERS-1 3 day repeat context) can be summarized as follows: errors in base-line estimation will give a systematic bias to the entire scene which can easily be mistaken for real data (Reference Joughin and KwokJoughin and others, 1996a). Phase unwrapping from a non-zero-velocity starting location will introduce a small constant offset to all data within a particular scene, in general less than (2/cos?)cmd−1. A phase-unwrapping error (using the Reference Goldstein, Zebker and WernerGoldstein and others (1988) technique) may introduce discontinuities in the unwrapped phase as integer multiples of 2π. If undetected during analysis, these will lead to errors in the derived surface velocities by corresponding integer multiples of (2.5/cos?)cmd−1. Such discontinuities arc often easily noticeable and may be remedied by filtering Coherence-related phase noise introduces a high-frequency noise component in velocity or topography results. In this work, phase noise introduces a speed uncertainty of about ±(3/cos?)mmd−1 in regions of high coherence. Finally, an important case-dependent source of error is failure of the constant velocity assumption (CVA). In subtle cases, the introduced errors will be less than one fringe (< (2.5/cos?cmd−1). Drastic CVA failures can invalidate velocity and topography results but can also provide insights into glacier dynamics and subglacial hydrology.
Particular to the Bering Glacier surge, a three-fold acceleration of West Bagley Icefield is clearly observed in the comparison of DSRI results from 1992 to 1994. This observation was possible because phase coherence is maintained in 1994 scenes acquired while the surge was in progress, implying that the surface of West Bagley Icefield was fairly stable over 3 day periods during this time. Such stability. combined with moderate acceleration relative to the Bering Glacier surge, may further imply that the massive restructuring of the subglacial hydrological system associated with surging glaciers (Reference KambKamb and others, 1985) did not extend up into West Bagley Icefield. If this is the case. the 20 km longitudinal acceleration profile was probably caused by a simple reduction of downstream restraint due to lowering of ice during the surge. This downstream surface-lowering may have imparted increased longitudinal stress gradients which had greatest effect along tin· center line of West Bagley Icefield, where the resistance of marginal shear stress was minimal. The data also place a coarse limit of 200-500 md−1 on the upstream surge-propagation speed. This is considerably faster than the observed downstream propagation of the surge front at ~100 md−1(Reference RoushRoush, 1996).The pre-surge velocity was used to estimate the depth of West Bagley Icefield at 500- 700 m. The surface drawdown implied by the acceleration event was not observed, possibly because it was not within the resolution limits of topographic DSRI.
Acknowledgement
We thank the U.S. National Science Foundation for providing initial support for this work with Small Grant for Exploratory Research OPP93-19873, NASA for providing support with grants NAGW-4930 and NAGW-4068, Cray-Research, Inc., for providing additional support via their University Research and Development Grant Program, and the Arctic Region Supercomputing Center for providing computational support. We also thank the Alaska SAR Facility, University of Alaska-Fairbanks, for a graduate research assistantship to D. R. Fatland, and wc also thank I. Joughin, D. Goldstein, M. Falmestock, E. Rignot, W. Harrison, K. Echelmeycr, M. Truffer, C. Raymond and M. Nolan for valuable conversations and comments on this work.
Appendices
Image Coregistration
Images are misregistered due to scene topography, the spatial separation of the two spacecraft orbital passes (the base line) and other systemic errors. Image coregistration to within one pixel is accomplished with a straightforward linear translation of one image. Sub-pixel coregistration is commonly applied to improve phase coherence. In a system of glaciers with a variety of flow directions (Fig. 1), it is necessary to use a locally adaptive resampling technique to compensate for misregistration incurred by the ice motion on top of the other effects. For this work, 20000 local misregistration measurements are calculated on a 1/3 km grid for each 40 km x 50 km scene. This is done both from image 1 to image 2 and vice versa, with only reciprocating offset vectors retained in order to eliminate spurious results. Each measurement is made using small correlation chips over-sampled by a factor of four in both directions. The resulting misregistration grid is low-pass filtered to reduce noise, giving a seir-eonsistent smoothly varying fractional pixel-offset grid which is used to drive a local bilinear interpolation resampling of one of the two images. This technique provides extremely useful signal recovery as shown in Figure 7.
Base-Line Calculation
A linear approximation to the non-linear problem of refining the interferometric base line has been given by Reference Joughinjoughin (1995). The standard approach involves phase unwrapping an interferogram to determine phase differences between several tie-point locations in the scene. These phase differences are compared to the theoretical phase differences given an assumed base line, and the error between the two is used to refine the base-line estimate. This process is iterated to give convergence to a working base line. The difficulty for valley glaciers is that ice-covered regions with good phase-signal coherence have surface-velocity information folded into the unwrapped phase, and the areas which are not moving are generally mountains which often have poor signal coherence and are difficult lo phase-unwrap. A solution lies in performing the base-line estimation process further along the swath in an area with moderate topographic relief and no moving ice. The base line is then extrapolated back to the scene of interest, as orbital separations vary slowly and fairly linearly over distance scales of tens of kilometers.
Phase-Unwrapping and Adaptive Phase-Filtering
In general, interferograms of temperate valley glaciers contain a great deal of decorrelation noise (e.g. relative lo polar ice sheets). This is due to such influences as surface rotation in areas with curving flowlines, rapid marginal shearing between the constraining valley walls and the center of the glacier, surface deformation of rapidly moving (e.g. surging) ice, surface melting and in the case of the Bering Glacier system, wind and precipitation from mercurial coastal weather systems. Decorrelation noise necessitates adaptive phase-filtering prior to phase-unwrapping in order to recover lower-frequency phase signals that are corrupted by high-frequency noise. The trade-off in this process is a loss of resolution, so that the problem becomes one of performing the minimal amount of phase filtering.
We used Goldstein's phase-unwrapping algorithm (Reference Goldstein, Zebker and WernerGoldstein and others, 1988). Valley glaciers often include large regions of mediocre to poor coherence, so the interferometric phase is adaptively low-pass filtered, using both amplitude and coherence weighting. Adaptive filtering selectively improves coherence in noisy areas and enables the algorithm to unwrap more of the scene while keeping the full spatial resolution intact in areas with good coherence. Low-pass filtering only works in areas where some coherence remains and where the phase signal is low frequency. There is a danger in low-pass filtering areas with extremely poor coherence, as a coherent signal can "emerge" with no physical meaning. Unwrapping such a false signal using Goldstein's technique can lead lo global (multiple of 2π) phase errors elsewhere in the scene.
Manual assistance can be provided to the phase-unwrapping process. In this work, sets of points were chosen by hand such that they passed through regions of low coherence. These points were connected by lines to the edge of the image, so that the resulting pixel chains would guide the branch-cutting part of the unwrapping algorithm, acting in effect like singularity grounding cables. Manual assistance is not feasible lor large-scale SRI applications but for small-scale research projects with well-defined regions of interest it can prove useful in eliminating spurious data.
