Registering imagery to ICESat data for measuring elevation changes on Byrd Glacier, Antarctica

GEOPHYSICAL RESEARCH LETTERS, VOL. 32, L23S05, doi:10.1029/2005GL024328, 2005 Registering imagery to ICESat data for measuring elevation changes on Byrd Glacier, Antarctica T. Schenk,1 B. Csatho,2 C. J. van der Veen,2,3 H. Brecher,2 Y. Ahn,1 and T. Yoon1 Received 7 August 2005; revised 13 October 2005; accepted 21 October 2005; published 7 December 2005. [1 ] We present a new approach to derive control information from ICESat data that enables rigorous registration of aerial and satellite imagery. The technique, based on matching terrain features identified from ICESat measurements and aerial imagery, opens the door to transform results of previous studies to a global reference frame. We demonstrate the proposed methodology with historical aerial photographs to determine surface changes between 1979 and 2004 over Byrd Glacier. This is important because there is no satellite radar altimetry coverage south of 81.5
S, which limits mass balance knowledge of outlet glaciers draining the East Antarctic ice sheet through the southern Transantarctic Mountains. Our study indicates that the grounded part of Byrd Glacier is close to being in balance. However, we observe large thinning on the floating part of the glacier, probably induced by increased basal melting. Citation: Schenk, T., B. Csatho, C. J. van der Veen, H. Brecher, Y. Ahn, and T. Yoon (2005), Registering imagery to ICESat data for measuring elevation changes on Byrd Glacier, Antarctica, Geophys. Res. Lett., 32, L23S05, doi:10.1029/2005GL024328. 1. Introduction [2] Mapping and monitoring polar regions is important for understanding current ice sheet mass balance and related sea level changes [e.g., Rignot and Thomas, 2002]. Dramatic thickness changes are occurring in ice drainages on both the Greenland and Antarctic ice sheets [e.g., Krabill et al., 2004; Thomas et al., 2004]. In order to put these changes into a broader temporal context, it is desirable to extend the change detection time line as far back as possible. [3] Since the 1940s hundreds of thousands of aerial photographs have been collected in Antarctica. Traditional photogrammetric techniques have been used successfully to determine velocities and elevations of glaciers in both the Arctic and Antarctic [e.g., Fastook et al., 1995; Brecher, 1982]. The major difficulties of this approach are (1) establishment of good Ground Control Points (GCPs) in remote areas, and (2) the need for special equipment and skilled personnel to carry out photogrammetric projects. Yet another difficulty in integrating results from earlier studies in a consistent reference frame is that the local coordinate 1 CEEGS Department, Ohio State University, Columbus, Ohio, USA. Byrd Polar Research Center, Ohio State University, Columbus, Ohio, USA. 3 Department of Geology Science, Ohio State University, Columbus, Ohio, USA. 2 Copyright 2005 by the American Geophysical Union. 0094-8276/05/2005GL024328$05.00 systems used for photogrammetric measurements are often poorly documented. As a result, only a small fraction of the available high quality photographs has been used in quantitative studies. [4] NASA’s Ice, Cloud and land Elevation Satellite (ICESat) was launched in January 2003 and since then it has collected laser altimeter data during several operational periods. ICESat measurements allow us to obtain reliable estimates of current ice-sheet mass balance and to identify regions of rapid ice sheet elevation changes, especially over the steep marginal regions. Moreover, ICESat measurements provide high-accuracy elevation data in a consistent, Earth-centered reference frame, suitable for establishing global geodetic control [Zwally et al., 2002]. [5] In this paper we describe how information from ICESat and aerial photographs can be combined to measure surface changes over outlet glaciers. First we present a new approach for registering aerial and satellite images to ICESat data, thereby circumventing the need to establish GCPs in the field. Once the data sets are registered, surface elevation changes can easily be computed. We demonstrate the feasibility of our method with aerial photographs that were acquired in 1978/79 for deriving velocities over Byrd Glacier [Brecher, 1982]. Byrd Glacier, flowing through the Transantarctic Mountains into the eastern part of the Ross Ice Shelf, is one of the largest glaciers in East Antarctica. Rignot and Thomas [2002] estimated that this glacier discharges over 20 km3/yr into the Ross Ice Shelf with a positive mass balance of 21 km3/yr, corresponding to a basin-averaged thickening rate of 2 cm/yr. A decrease in ice velocity, recently reported by Stearns and Hamilton [2005], suggests that Byrd Glacier has undergone substantial changes in velocity over the last 20 years. Establishing a rigorous relationship between the local coordinate frame used by Brecher [1982] and ICESat allows determination of surface elevation changes of the glacier between the late 1970s and 2003 – 05. Combining surface elevation changes with surface velocity observations of Stearns and Hamilton [2005] over the same time period may further our understanding of the observed deceleration. 2. Methodology [6] ICESat points have been successfully registered to known surfaces, mainly for the purpose of calibration. The transformation is established by minimizing the distance between laser points and surfaces of known elevations [e.g., Martin et al., 2005]. Another approach for assessing ICESat footprint horizontal geolocation accuracy is based on measuring the similarity between measured and simulated waveforms by shifting the ICESat laser footprint over a reference Digital Elevation Model (DEM) at regularly spaced posi- L23S05 1 of 4 L23S05 SCHENK ET AL.: REGISTERING IMAGERY TO ICESAT DATA L23S05 approximating the DEM with analytical surface functions, e.g. planar surface patches [Filin, 2003]. This would offer the advantage that all laser points can be used for determining the transformation parameters but it requires an explicit correspondence between laser points and surface patches. 3. Validation of the Methodology Figure 1. (left) ICESat measurements across a typical Antarctic nunatak with straight ridge crests, bounded by steep rectilinear slopes leading down to the ice sheet. Characteristic terrain points (triangles), corresponding to a ridge crest and the break of the slopes at the ice sheet boundary are extracted from the ICESat elevation profile. These terrain features, shown also on the oblique photograph (right), are measured stereoscopically from the aerial photographs. Characteristic ICESat points and the corresponding linear terrain features allow for a rigid body transformation. tions. The best fit is found at the location of highest similarity [Harding and Carabajal, 2005]. [7] We have developed a new method for registering images to ICESat laser points without the necessity of creating a DEM. It is based on the assumption that the orientation of overlapping aerial or satellite images with respect to a local 3D Cartesian coordinate system has been determined a priori, for example by way of aerial triangulation [e.g., Schenk, 2003]. Registering images to ICESat laser points entails establishing a transformation between the respective 3D Cartesian coordinate systems. This would be a straight-forward task if individual corresponding points could be identified from both the aerial photographs and the ICESat data sets. Since this is not possible we have to resort to terrain features that can be obtained from both data sets. Characteristic points can be extracted from elevation profiles, for example at abrupt changes of slopes [e.g., Thappa, 1987]. These characteristic points, which we call ‘‘terrain feature points’’, often correspond to landform elements, for example ridge crests or ice sheet boundaries. As computed quantities from several ICESat laser points, terrain feature points are more robust and more accurate than individual ICESat laser points. Figure 1 illustrates the principle. Terrain feature points extracted from ICESat transects should coincide with the 3D lines measured from the images. Conceptually, one can imagine shifting, rotating and scaling the local photogrammetry system until the ICESat terrain feature points agree with the measured topographic features within the expected errors. We have developed a mathematical model to incorporate this concept [Schenk, 2005]. It is based on minimizing the shortest distance between characteristic points and corresponding terrain features in a least-squares sense. [8] The proposed registration method provides the transformation parameters directly through a least-squares adjustment. In contrast, registering ICESat laser points to an image-derived DEM is an indirect method that determines the differences between laser points and DEM at arbitrary registration locations. This approach can be improved by 3.1. Registration of Aerial Imagery to ICESat Data 3.1.1. ICESat Data [9] For the registration we use Antarctic and Greenland Ice Sheet Data Product (GLA12) from Laser 2a (Release 21; October – November 2003) according to the following criteria: (1) points should be on stable ground (outcrops), (2) have small residual pointing errors and (3) small ranging errors. Points on stable ground are selected by visual inspection of the aerial photographs. As for meeting criterion (2), ICESat data from Laser 2a have the best pointing accuracy, estimated to 3 arcsec, corresponding to 9 meter planimetric accuracy on the ground [Luthcke et al., 2005]. Criterion (3) is satisfied by properly weighting ICESat observations, for example based on slope and roughness information and quality of Gaussian fitting. 3.1.2. Aerial Photographs [10] Aerial photographs over Byrd Glacier along six parallel flight lines were acquired on December 6, 1978, and January 31, 1979, with the objective of determining surface velocities and elevations [Brecher, 1982]. We selected four sites (blue boxes in Figure 3), covering outcrops along the glacier, to establish the coordinate transformation. Diapositives of aerial photographs acquired over these sites have been digitized to a resolution of
15 mm pixel size. 3.1.3. Transformation [11] The aerial triangulation, described in detail by Brecher [1982], was carried out in a local 3D Cartesian coordinate system. The sequence of transformations that Figure 2. (a) Surface elevations measured across the NE ridge of Mt. Tuatara along ICESat ground track 203(91d). Solid circles are ICESat elevations acquired on November 4, 2003 and triangles are elevations from aerial photogrammetry measured prior and after the transformation of the local system. Characteristic terrain points, extracted from the ICESat elevation profile, are shown as large circles. (b) A small area of the aerial image. White dots indicate the image positions of laser points before and the black solid circles after the transformation. Large circles mark characteristic terrain points from ICESat and white lines are terrain features measured stereoscopically on the images. 2 of 4 L23S05 SCHENK ET AL.: REGISTERING IMAGERY TO ICESAT DATA L23S05 elevation differences are now much smaller, indicating that the adjustment transforms the local system to a better location than originally obtained from prior knowledge about the local reference system. We have performed the transformation with a total of 15 characteristic lines and points in the four selected sites. The variance/covariance analysis of the adjustment revealed an accuracy (1s) of ±8 meters for the horizontal translation parameters and ±3 meters for the vertical component. These numbers agree well with the a priori error estimation of the characteristic points and the measured characteristic line in the aerial stereopair (±1 meter in planimetry, ±1.8 meter in elevation). As shown in detail by Schenk [2005], the accuracy of the characteristic points depends on the fitting error of lines or curves through ICESat laser points and also on the intersection angle of the lines. Figure 3. Average surface elevation change rates between 1979 and 2005 from photogrammetry and ICESat elevations over Byrd Glacier, depicted on AMM-1 SAR mosaic [Jezek and RAMP Product Team, 2002]. Dark thin lines show ICESat ground tracks, thick lines refer to segments used for change detection. Blue boxes indicate the 4 registration sites, red lines mark segments of ICESat tracks used for registration. The grounding zone adopted from Scofield et al. [1991] is highlighted by dark background. defines the relationship between the ICESat reference frame and the local system also includes the computation of orthometric heights because the GCPs in the photogrammetry system have been established by triangulation. [12] We now apply the registration method explained in Section 2. Figure 2a shows a section of an ICESat ground track in a distance vs. elevation representation. Also shown are photogrammetrically measured elevations at the same ICESat positions. The figure reveals that the initial approximation of the transformation, obtained from a priori knowledge about the local photogrammetric coordinate system, is off by more than one hundred meters. The ground track section shown in Figure 2a has three distinct terrain feature points that correspond to the mountain ridge and to the transition from the snow surface on the right and Byrd Glacier on the left to the steep mountain slopes. Although the terrain feature points in Figure 2a are near ICESat laser points, it is important to realize that they are quantities derived from several laser points. [13] It is quite likely that characteristic ICESat points correspond to 3D terrain features (e.g. ridge crests). These 3D features can be identified and measured in the images. Figure 2b depicts a small area of a digitized image. The ICESat laser points with the initial estimate of the transformation parameters are shown as white dots. Also shown are the three derived characteristic points (large circles), which should coincide with the measured terrain features (white lines). A rigid body transformation (3 translations, 3 rotations) is now applied for an optimal fit of characteristic ICESat points and 3D terrain features. [14] The adjustment yields improved transformation parameters. After transforming the local system we measured the elevations at the ICESat footprint locations again on the aerial images. Figure 2a shows the result. The 3.2. Ice Surface Elevation Changes 3.2.1. ICESat Data [15] For determining surface elevation changes over the Byrd Glacier we use ICESat data from the Antarctic and Greenland Ice Sheet Data Product (GLA12) for Laser 1 (Release 18; February – March 2003) Laser 2a (Release 21; October – November 2003), Laser 2b (Release 16; February – March 2004), Laser 2c (Release 17; May– June 2004) and Laser 3a (Release 22; October – November 2004). To reduce the range error induced by detector saturation we apply the correction described by Fricker et al. [2005] on the GLA12 Data. Observations with low return signal energy (apparent reflectivity less than 20% in GLA12) may have large ranging error and are therefore rejected. 3.2.2. Elevations From Aerial Photographs [16] We use elevations determined from the January 31 1979 aerial survey. The position and elevation of 601 points, located on natural features, such as crevasses, seracs, or dunes, were determined from the photogrammetry survey over the glacier by aerial triangulation to an accuracy of about 1.8 m in elevation [Brecher, 1982]. The average spacing between these points is about 2.5 km. 3.2.3. Elevation Changes [17] After transforming the photogrammetry measurements from the local reference system into the ICESat reference frame, we compare ICESat elevations (2003 – 04) with elevations derived from the 1979 photogrammetry survey. We assume that surface elevation change rates were constant during this period. First we compute elevation change rates at each ICESat point by using a linear interpolation to obtain an estimate of 1979 surface elevation from neighboring aerial survey points. To avoid large interpolation errors we only use ICESat points located within 1.1 km to the nearest aerial survey point. The map depicted in Figure 3 is obtained by kriging interpolation of these elevation change rates. [18] Since the elevation changes are derived from differences between ICESat laser points and block points established by aerial triangulation, the accuracy of elevation changes is affected by the accuracy of ICESat and aerial triangulation points, and the interpolation error. ICESat elevation errors (se,ICESat) are due to pointing and range biases, and ranging errors. Based on the global analysis of ICESat pointing and ranging biases [Luthcke et al., 2005] and our studies in the Dry Valleys where an accurate DEM 3 of 4 L23S05 SCHENK ET AL.: REGISTERING IMAGERY TO ICESAT DATA is available [Csatho et al., 2005], we estimate the combined effect of these errors over the gently sloping surface of the Byrd Glacier (2 degree maximum slope) to be 2 meters. Considering the relatively smooth surface, the interpolation error (sinterpolation) is rather small,
1 meter. Based on a variance/covariance analysis of the aerial triangulation results we estimate a remarkable point accuracy of 1 meter in planimetry and 1.8 meter in elevation (se,aerial). Another error source to consider stems from the transformation (stransformation). As mentioned in section 3.1.3, a representative accuracy is 3 meters. Combining all these errors by way of error propagation we estimate an accuracy of 4.15 m for the elevation changes: se ¼  p 2 se;ICESat þ s2interpolation þ s2e;aerial þ s2transformation : This translates into an error of 0.17 m/yr in elevation change rates for the 25 years between the observations. 4. Interpretation of Elevation Changes [19] For this study we adopted the grounding zone position shown in Figure 3 from Scofield et al. [1991], who used results from 1978 – 79 photogrammetry survey and radio-echo sounding profiles. The grounding zone coincides with the region where a transition from little change to thinning occurs, suggesting that the ice upstream is partially grounded. Over the grounded part of the glacier we estimate an average surface elevation change rate of 0.02 ± 0.17 m/yr between 1979 and 2004 with local surface elevation change rates ranging from 0.81 to +0.68 m/yr. Surface elevation changes seem to be correlated with basal drag, calculated by Whillans et al. [1989] from the 1978 – 79 velocities and elevations. Largest surface lowering is observed at areas of small basal drag, while the glacier is thickening over basal drag maxima (‘‘sticky spots’’). This supports the conclusion of Whillans et al. [1989] and Scofield et al. [1991] that some parts of the bed under the Byrd Glacier are frozen while other regions are thawed. [20] Surface lowering of 0.4 – 1.2 m/yr is detected on the floating part of the glacier. Sub-shelf circulation is generally such that the greatest basal melting rates are found in the vicinity of grounding lines well below sea level [e.g., Jacobs et al., 1992]. On Byrd Glacier, the transition from grounded to floating ice occurs at
1100 m below sea level [e.g., Reusch and Hughes, 2003]. It appears, therefore, that the surface lowering observed on the floating part of Byrd Glacier is associated with basal melting induced by ocean circulation under the shelf. It may be that the circulation has become more vigorous, or the increased basal melting may reflect a warming of the ocean water. Whatever the cause, the observed thinning should be accompanied by slowing of the glacier. On floating ice shelves, the along-flow stretching rate is strongly dependent on the ice thickness [Weertman, 1957], thus even moderate thinning will significantly reduce the velocity gradient in the flow direction. With little change observed in the grounded portion of the glacier, reduced stretching beyond the grounding line will result in a decrease in ice speed, as indeed inferred by Stearns and Hamilton [2005]. L23S05 [21] Acknowledgments. We thank W. Lee and P-L. Lai for help with data processing and H. Fricker for advice on saturation correction of the ICESat data. We thank NASA’s ICESat Science Project and the NSIDC for distribution of the ICESat data, see http://icesat.gsfc.nasa.gov and http:// nsidc.org/data/icesat/. Comments from Helen Fricker, Chris Shuman and Bob Schutz improved the manuscript. Byrd Polar Research Center contribution number C-1332. References Brecher, H. (1982), Photogrammetric determination of surface velocities and elevations on Byrd Glacier, Antarct. J. U. S., XVII, 79 – 81. Csatho, B., et al. (2005), Airborne laser scanning for high-resolution mapping of Antarctica, Eos Trans. AGU, 86, 237 – 238. Fastook, J. L., H. H. Brecher, and T. J. Hughes (1995), Derived bedrock elevations, strain rates and stresses from measured surface elevations and velocities: Jakobshavns Isbræ, Greenland, J. Glaciol., 41, 161 – 173. Filin, S. (2003), Recovery of systematic biases in laser altimetry data using natural surfaces, Photogramm. Eng. Remote Sens., 69, 1235 – 1242. Fricker, H. A., A. Borsa, B. Minster, C. Carabajal, K. Quinn, and B. Bills (2005), Assessment of ICESat performance at the salar de Uyuni, Bolivia, Geophys. Res. Lett., 32, L21S06, doi:10.1029/2005GL023423. Harding, D. J., and C. C. Carabajal (2005), ICESat waveform measurements of within-footprint topographic relief and vegetation vertical structure, Geophys. Res. Lett., 32, L21S10, doi:10.1029/2005GL023471. Jacobs, S. S., H. H. Hellmer, C. S. M. Doake, A. Jenkins, and R. M. Frolich (1992), Melting of ice shelves and the mass balance of Antarctica, J. Glaciol., 38, 375 – 387. Jezek, K., and RAMP Product Team (2002), RAMP AMM-1 SAR image mosaic of Antarctica, http://nsidc.org/data/nsidc-0103.html, Natl. Snow and Ice Data Cent., Boulder, Colo. Krabill, W., et al. (2004), Greenland Ice Sheet: Increased coastal thinning, Geophys. Res. Lett., 31, L24402, doi:10.1029/2004GL021533. Luthcke, S. B., D. D. Rowlands, T. A. Williams, and M. Sirota (2005), Reduction of ICESat systematic geolocation errors and the impact on ice sheet elevation change detection, Geophys. Res. Lett., 32, L21S05, doi:10.1029/2005GL023689. Martin, C. F., R. H. Thomas, W. B. Krabill, and S. S. Manizade (2005), ICESat range and mounting bias estimation over precisely-surveyed terrain, Geophys. Res. Lett., 32, L21S07, doi:10.1029/2005GL023800. Reusch, D., and T. Hughes (2003), Surface ‘‘waves’’ on Byrd Glacier, Antarctica, Antarct. Sci., 15, 547 – 555. Rignot, E., and R. H. Thomas (2002), Mass balance of polar ice sheets, Science, 297, 1502 – 1506. Schenk, T. (2003), Progress in automatic aerotriangulation, in Photogrammetric Week’03, edited by D. Fritsch, pp. 129 – 140, Wichmann, Heidelberg, Germany. Schenk, T. (2005), Mathematical models for registering imagery to 3D point clouds, Tech. Notes Photogramm. 23, 18 pp., Dep. of Civil and Environ. Eng. and Geod. Sci., Ohio State Univ., Columbus. Scofield, J. P., J. L. Fastook, and T. J. Hughes (1991), Evidence for a frozen bed, Byrd Glacier, Antarctica, J. Geophys. Res., 96, 11,649 – 11,655. Stearns, L., and G. Hamilton (2005), A new velocity map for Byrd Glacier, Antarctica from sequential ASTER satellite imagery, Ann. Glaciol., in press. Thappa, K. (1987), Extracting characteristic topographic points, Ph.D. diss., 118 pp., Dep. of Geod. Sci. and Surv., Ohio State Univ., Columbus. Thomas, R., et al. (2004), Accelerated sea-level rise from West Antarctica, Science, 306, 255 – 258. Weertman, J. (1957), Deformation of floating ice shelves, J. Glaciol., 3, 38 – 42. Whillans, I. M., Y. H. Chen, C. J. van der Veen, and T. J. Hughes (1989), Force budget: III. Application to three dimensional flow of Byrd Glacier, Antarctica, J. Glaciol., 35, 68 – 80. Zwally, H. J., et al. (2002), ICESat’s laser measurements of polar ice, atmosphere, ocean and land, J. Geodyn., 34, 405 – 445.  Y. Ahn, T. Schenk, and T. Yoon, CEEGS Department, Ohio State University, 2070 Neil Ave., Columbus, OH 43210-1247, USA. (schenk.2@osu.edu) H. Brecher, B. Csatho, and C. J. van der Veen, Byrd Polar Research Center, Ohio State University, 1090 Carmack Road, Columbus, OH 43210-1002, USA. C. J. van der Veen, Department of Geological Sciences, Ohio State University, 125 S Oval Mall, Columbus, OH 43210-1398, USA. 4 of 4