Antarctica

Remote Sensing of Environment 106 (2007) 269 – 284 www.elsevier.com/locate/rse

Refined analysis of radar altimetry data applied to the region of the subglacial Lake Vostok/Antarctica S. Roemer a,1 , B. Legrésy b , M. Horwath a,⁎, R. Dietrich a a

Technische Universität Dresden, Institut für Planetare Geodäsie, D-01062 Dresden, Germany b LEGOS, 18 av. E. Belin, Toulouse CEDEX 9 31401, France

Received 25 January 2005; received in revised form 20 September 2005; accepted 23 February 2006

Abstract Satellite altimetry is a powerful tool to map the ice sheet elevation as well as a number of other parameters linked to geometrical and geophysical properties of ice sheets. Irrespective of new instrumental developments like the laser altimeter on ICESat (Ice, Cloud and land Elevation Satellite) the well-established radar altimeter (RA) missions ERS-1,2 (European Remote Sensing satellites) and ENVISAT (ENVIronment SATellite) are unique in their temporal coverage over more than a decade and in their temporal and spatial sampling. Therefore, the full exploitation of these RA data by improved methods is imperative. Here we develop improved techniques to correct for topographically induced errors by a refined consideration of the relevant topographic conditions. Furthermore we improve the gridding procedure by adapting it to local conditions and thus preserving smaller-scale features. We apply our methods to the region of the subglacial Lake Vostok/Antarctica and derive digital elevation models (DEMs) for this region with the aim of improving/resolving smaller scales. The effect of our improvements is demonstrated in detail by comparing our DEMs and previously published DEMs to the ICESat laser measurements which are taken as a reference here. Due to our improvements, the standard deviation of the difference between RA-based DEMs of the Lake Vostok region and ICESat data decreases from more than 1.1 m to 0.5 m. This remaining error is chiefly the error inherent in the RA observations. Our RA-ICESat comparisons, supported by Fourier analyses, also reveal the presence and importance of small-scale features that can be detected by laser but not by the RA measurements. © 2006 Elsevier Inc. All rights reserved. Keywords: Antarctica; Subglacial Lake Vostok; ERS-1; Radar altimetry; Slope error correction; Ice sheet mapping; Validation of digital elevation models; ICESat

1. Introduction Satellite altimetry is a powerful tool for mapping the surface elevation and other physical parameters of the Antarctic and Greenland ice sheets with almost complete and homogeneous coverage. A number of radar altimeter (RA) missions have been operating since the 1990s (Davis, 1992), in particular the European Remote Sensing satellites ERS-1 and ERS-2 and ENVISAT (ENVIronment SATellite) which have sampled up to

⁎ Corresponding author. E-mail addresses: [email protected] (S. Roemer), [email protected] (M. Horwath). 1 Now at GeoForschungsZentrum Potsdam (GFZ), D-14473 Potsdam, Germany. 0034-4257/$ - see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2006.02.026

latitudes of 82°. Since 2003 the laser altimeter (LA) mission Ice, Cloud and land Elevation Satellite (ICESat) has observed ice sheet elevations up to 86° latitude (Schutz et al., 2005; Zwally et al., 2002). Both systems have advantages and disadvantages. Until now (partly due to the ICESat data gain being smaller than anticipated), the temporal resolution of the RA missions is higher than that of ICESat, and the cross-track resolution relevant for topographic mapping was higher for the ERS-1 geodetic mission phase than it is for ICESat. The most important advantage of the RA missions is their earlier start in observing time resulting in long and continuous time series. For these reasons, RA data should be fully exploited. An important aim is to infer the precise surface topography which, on the one hand, may serve as a reference for studying long term temporal variations and, on the other hand, is interesting on its own for rheological analyses or the study of surface processes.

270

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

Customary methods of continental or regional RA data analysis (Bamber & Bindschadler, 1997; Fricker et al., 2000; Rémy et al., 1999) can still be significantly improved. The largest improvement potential exists for the methods of correcting the slope induced error (SE) and of gridding: For an inland Antarctic region possible improvements amount to several decimeters, as will be shown in Sections 4.8 and 4.9. Errors from other sources (see Section 2) are lower, typically below 10 cm. The aim of this paper is to refine the methods of RA data analysis for precise modelling of ice elevation and additional altimetric parameters. To develop our methods we apply them to ERS-1 geodetic mission phase data in the region of Lake Vostok (Fig. 1) and derive digital elevation models (DEMs). We thoroughly evaluate and validate our results. For this purpose, we use ICESat data as independent information. We show that our DEMs derived by the refined methods are significant improvements compared to previous satellite RA-based models and also to models derived from local airborne surveys (Studinger et al., 2002, 2003; Tabacco & Bianchi, 2002). 2. ERS-1 data The ERS-1 mission (Benveniste, 1993) included the geodetic phases E/F (04/10/1994 to 03/21/1995), which by their 336 day non-repeat configuration, provided a high density of ground tracks and are therefore especially suited for surface topography mapping. At the high latitudes of Lake Vostok the track separation is about 2.0 km. Along track, the observation increment is about 350 m. These ERS-1 RA data – processed using the OSCAR2/ ICE2 algorithms (Legrésy et al., 2005; Rémy et al., 1999) and referred to precise orbits from Delft Institute for Earth-oriented Space Research (Scharroo et al., 1998) – are analysed in this study. In the Lake Vostok region, the corrections applied to the data amount to −0.1 m to 0.1 m for the earth tide, 1.4 m to 1.5 m for the dry tropospheric delay, up to 0.01 m for the wet tropospheric delay, 0.02 m to 0.05 m for the ionospheric delay and up to 0.01 m for the Doppler correction. Residual errors of these corrections are assumed to be in the order of few centimetres. The orbit error is in the order of 5 cm (Scharroo et al., 1998), and its large-scale components are removed by our analysis (see Section 3.3). The retracking of the altimeter waveform provides a mean range to the radar scatterers within a certain ground area of several square kilometres (the ‘pulse-limited footprint’, see Section 3.2.1). In case of a pure surface reflection the mean surface height over this area could be deduced. However, the 13.8 GHz Ku band microwave signal also penetrates the snow and is reflected by subsurface layers, resulting in a mean reflection horizon which is generally below the mean surface (Féménias et al., 1993; Legrésy & Rémy, 1998). Occurrence and magnitude of this effect vary in space and time dependent not only on snow characteristics (roughness, temperature, density, etc.) but, in case of a polarised signal, also on the azimuth of the signal itself (Legrésy et al., 1999). The precision of the individual retracked range observations is of the order of 0.5 m,

2

Observing continental surfaces with radar altimetry.

Vostok

Fig. 1. The Lake Vostok region in the composite Radarsat image (Jezek and RAMP Product Team, 2002). The flat lake area stands out clearly against the more undulated area of grounded ice.

as will be shown later based on the along-track variability (Section 4.6) and on the temporal variability at repeat tracks (Section 4.7). When averaging over typically 40 to 100 observations in the course of the DEM generation (see Section 3.4), the influence of this shot-to-shot precision becomes smaller than 10 cm. Besides the ellipsoidal height (h) itself, given with respect to the GRS80 (Geodetic Reference System 1980) ellipsoid, the retracking technique allows us to estimate three other parameters explained in detail in Legrésy and Rémy (1997) and Legrésy et al. (2005): namely, backscattering coefficient (Bs), leading-edge width (LeW) and trailing-edge slope (TeS). Each of these altimetric parameters contains useful information (cf. Section 4.2). They help to detect penetration effects or to interpret the results and, thus, contribute to a more complete understanding of the surface characteristics in the observed area. Beside the ERS-1 data we also use other datasets, which are introduced when needed. The DEMs derived in this study are compared to previously published DEMs (Section 4.3) and all DEMs are compared to ICESat data (Section 4.4). The temporal variability of ERS-2 data is used to assess the precision of both ERS RA systems (Section 4.7). 3. Methods of refined radar altimetry analysis 3.1. General scheme Our radar altimetry analysis comprises three steps: the correction of the slope induced error (Section 3.2), the elimination of outliers and offsets (Section 3.3), and the generation of gridded models from the irregularly distributed data (Section 3.4). The first two steps depend on a-priori topographic information, and the second step additionally depends on a-priori models of the other altimetric parameters of interest. To avoid errors in a-priori models affecting the final results, the whole analysis is iterated. The apriori models for each iteration step are obtained as follows. For the first iteration, we generate the a-priori models using the ERS-1 data without any data screening and slope error correction. In each further step, the result of the previous iteration is used as the apriori model. In what follows, the terms ‘a-priori model’ or ‘apriori DEM’ always refer to these auxiliary results within our

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

processing chain. The effect of this iteration procedure is discussed in more detail in Section 4.10. 3.2. Correction of the slope induced error 3.2.1. Review of concepts A satellite RA ‘illuminates’ an extended area, the “beamlimited footprint” (BLF). Only the echo from a small sub-area, the “pulse-limited footprint” (PLF) is used for range determination. This is the area that is closest to the satellite and, thus, provides the first radar return (Brooks et al., 1978; Fu & Cazenave, 2001). Due to surface topography the first return is not necessarily located in the satellite's nadir direction. For ERS-1 and for a flat surface, the diameters are in the order of 36 km for the BLF and 2.5 km for the PLF (Fu & Cazenave, 2001). Interpreting the observation as the range to the nadir point causes the so-called slope induced error (SE) which has to be considered over continental surfaces (Brenner et al., 1983). Fig. 2 gives an illustration of the geometry of SE correction. The SE correction already requires topographic information which has usually been introduced through simple surface parameters such as the surface slope. Such surface parameters have been provided, e.g. from an external slope data base (Bamber, 1994), or by computing slope parameters (Ihde et al., 2002) or even curvature parameters (Rémy et al., 1989) from an a-priori DEM by fitting a mathematical surface model (e.g. plane or higher order polynomial surface) over a certain area. The computed slope and curvature

271

parameters depend not only on the a-priori DEM but also on the particular mathematical surface model and on the adapted spatial scale of the surface fit (Brenner et al., 1983). Brenner et al. (1983) call the obtained surface parameters ‘effective’. They are not necessarily present in the actual topography but represent a hypothetical surface that induces the same SE as the actual topography. The ‘effective’ surface is used to correct for the SE and to estimate the location of impact (Fig. 2a): The hypothetical nearest point Ie to the satellite (in the observed distance rS) is uniquely determined by assuming an orthogonal reflection in this point on the ‘effective’ surface. Three alternative correction methods are common: the “relocation,” the “direct” and the “intermediate” method (Brenner et al., 1983; Rémy et al., 1989). Here we consider the first two of them as the ‘extreme’ concepts. A refinement of the intermediate method (which would follow mostly the improved approach of the “direct” method) may be a subject of an additional study in the future. The relocation method attempts to determine the height at the position to which the radar measurement actually refers, i.e., at the impact point It which represents the center of the PLF (Fig. 2b). This approach accounts for the geometrical principles of radar. Also the other altimetric parameters (Section 2) which constitute geophysical observations of the area around It are then properly attributed to this location. The displacement between the nadir point and the ‘upslope’ shifted impact point may amount to several kilometres (cf. end of Section 3.2.2). Consequently, the locations to which the observations refer are concentrated near local topographic maxima and are sparse near local topographic minima

Fig. 2. Illustration of slope induced error (SE) correction methods. For details not explained here refer to Sections 3.2.1 and 3.2.2. Surface slopes are largely exaggerated. (a) Overall situation with classic relocation method and direct method. (b) Refined relocation method. (c) Detail of (b) showing the effects of errors in the vertical and horizontal positions if I described by Eqs. (5) and (6).

272

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

(Bamber et al., 1998). This inhomogeneity in spatial data distribution (Fig. 4) may complicate gridding. The direct method, instead, attributes the observations to the nadir point (Fig. 2a). The avoidance of the upslope shift may, however, cause interpretation errors, since geophysical parameters are not attributed to their locations of origin. Retaining the nadir positions may facilitate gridding but does not avoid the intrinsic error (Bamber, 1994) of height estimates for positions without observations (Fig. 2). 3.2.2. Refined relocation method Here we propose a method that directly uses the topographic information from a DEM for the SE correction. It avoids the use of simple slope and curvature parameters, i.e., of intermediate products describing an ‘effective’ surface. The method is more robust and reproduces the observed range easily (see Fig. 2b for an illustration). The essential quantity affecting the location of the PLF is the range between the surface and the satellite. Hence, this location can be identified directly from a DEM. Fig. 3 shows an example of the topography and the corresponding ranges for an individual satellite position. Given an elevation model H ¼ Hðk; /Þ

ð1Þ

where λ and ϕ are ellipsoidal longitude and latitude and given a satellite position S, the range between the satellite and the surface is given as RS ¼ RS ðk; /Þ:

ð2Þ

We model the PLF as an area of fixed shape and size and search for the location of this area (i.e., of the position (λI,ϕI) of its central point I ) that minimises the mean distance R¯S to S. For each individual observation, the altimetric dataset (Section 2) contains the satellite position S (including its height hs) and the ellipsoidal surface height hN obtained under the assumption that the range was measured to nadir. Now, given the assumed impact point I, a new estimate of its height hI is computed by hI ¼ hN þ s r :

ð3Þ

The SE correction sr is (see Fig. 2b): P

sr ¼ SI−ðhS −HI Þ;

ð4Þ

where HI is the height of I according to the a-priori DEM (see P Section 3.1) and the distance SI is computed using HI. Errors of the SE correction for an assumed horizontal position of I originate from the (vertical) error of the a-priori DEM, ΔHI, at the horizontal position of I and from the error in the a-priori estimate of this horizontal position. The effect of ΔHI is negligible if the off-nadir angle α is small because ΔHI is P P approximately compensated by the resulting error D SI of SI (see Fig. 2c): P

Dsr ¼ D SI−ð−DHI Þ ¼ −DHI cosa þ DHI c0:

ð5Þ

For α = 1.3° (corresponding to the ERS altimeter beamwidth) and an error ΔHI = 10 m the effect on sr is below 3 mm. Hence,

Fig. 3. View of a pulsewidth-limited satellite RA system to the surface. (a) Topography (ellipsoidal heights w.r.t. the GRS80 ellipsoid, 2 m contour intervals) around the nadir point marked in the center. (b) Ranges (10 m contour intervals) between the satellite and the topography of Fig. 3a. The range retracked from the altimeter waveform represents a mean value over the pulselimited footprint (PLF), which is an irregular area. Our new relocation method finds the PLF position based on these ranges. The PLF is modelled as an area of fixed shape and size (marked by the square).

for the evaluation of Eq. (4), HI can be approximated by bilinear interpolation from the a-priori DEM grid. For the effect of the horizontal position error of I, Δx, on the estimated height hI we find under the assumption of a flat surface between the true impact point and I (see Fig. 2c): qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P2 P DhI c SI þ Dx2 −SI: ð6Þ To keep the effect smaller than 1 mm, Δx has to be smaller than 40 m. By estimating the footprint location with an accuracy of 10 m as described below, we meet this requirement. For a simple technical realisation, the PLF is assumed to be a rectangular area A (Fig. 3) and its side lengths Δλ,Δϕ are fixed

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

273

to multiples of the high resolution grid increments (which are 2′ in longitude and 30″ in latitude). The mean range of A is then RS ¼ P

1 A

Z

kI þDk 2

kI −Dk 2

Z

/I þD/ 2 /I −D/ 2

RS ðk; /ÞdA

ð7Þ

where dA is the surface element. In practice, the ranges RS (Eq. (2)) are given on the discrete grid of the DEM. For a rectangular grid cell j, the mean range over this cell, R¯Sj, is computed as the mean of the four corner point ranges. The mean range R¯S over an area A of n cells then arises as the mean value over those cells, weighted by their areas: P

RS ¼

n P 1X ASjRSj : A j¼1

ð8Þ

We use a simple discrete search algorithm to find the position (λI,ϕI) of A that minimises R¯S: First, we use the grid increments Δϕ,Δλ as search increments and search for an approximate minimum position over the entire BLF. If we would stop here, all observations would be relocated to grid nodes. To reach the required horizontal accuracy, to be independent of the grid resolution and to avoid stacking of observations we continue the search in the vicinity of the approximate position by iteratively decreasing the search increment to a final value of 10 m (see Eq. (6) and its discussion). The summation in Eq. (8) now includes sub-grid cells. Fig. 4a illustrates the distances of I from nadir. The maximum displacement in the Lake Vostok region is about 6 km. Fig. 4b shows the obtained pattern of relocated observations in a smaller sector south of the lake. 3.2.3. Refined direct method As discussed in Section 3.2.1, the errors of traditional SE correction methods originate both from the a-priori topographic information and from deficiencies in the way of deducing a simple ‘effective’ surface model from this a-priori information. While the adoption of such an ‘effective’ surface model can be avoided in case of the relocation method (Section 3.2.2) it is still inevitable for the direct method: Otherwise the range observation cannot be reduced to the nadir location, where no observation is available. Rémy et al. (1989) demonstrated that accounting for the curvature in addition to the slope significantly improves the result. Here we present a refined method that improves the quality of the slope and curvature parameters: Since those parameters are sensitive to the adopted spatial scale we adapt this scale to the specific local conditions by taking into account the impactpoint location as estimated in Section 3.2.2. To derive slope and curvature we fit the following 2nd order polynomial to the a-priori DEM: zðx; yÞ ¼ h þ cx x þ cy y þ cxy xy þ cxx x 2 þ cyy y 2

P

ð9Þ

where cxx, cx, cyy , cy , cxy , h are the coefficients and x, y are the horizontal coordinates centered at nadir. The area of fit should include the PLF and the nadir point because the fit should provide an ‘effective’ surface for relating the heights between

Fig. 4. Displacement of the pulse-limited footprint (PLF) from nadir due to the topography. (a) Map of the PLF displacement for the entire region of Lake Vostok. (b) Positions of PLFs in an area south of Lake Vostok. The subsatellite tracks are distorted by shifts towards local topographic maxima. This causes data gaps in surface depressions. Colours show the individual displacement between nadir points and PLF positions.

these two locations. The area should be, on the other hand, small enough to avoid smoothing of topographic features. We choose the minimal rectangular area that encloses the entire PLF (as

274

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

determined in Section 3.2.2) and that is centred in the nadir point. The last condition has proven to support a robust fit. For our region, the area size is between 6.25 km2 (the PLF's size) and about 110 km2, and in 90% of the cases it is below 23 km2. According to Rémy et al. (1989) the SE correction is then r s2 S  sd ¼  2 1−rS c− RE1þh

ð10Þ 1

where rS is the observed range, s ¼ ðc2x þ c2y Þ2 is the slope in the gradient direction, c = 2(cx2cxx + cxcycxy + cy2cyy)s− 2 the curvature in the gradient direction and RE the Earth's radius. The error of sd can then be estimated by variance propagation from the errors of the fitted parameters of Eq. (9).

3.3. Outlier and offset elimination The following procedure is repeated in every iteration step of the model generation, always starting with the whole dataset for each altimetric parameter p (h, LeW, TeS and Bs): The differences Δp between the (SE corrected) observations and the apriori model are computed. Tracks with an exceptionally high standard deviation of Δp are eliminated as a whole. We take the elimination criterion σΔp N nσ ¯Δp where σΔp is the standard deviation over the respective individual track, σ ¯Δp the standard deviation over all tracks, and n the threshold factor specified below. For each remaining track, individual outliers are identified P and eliminated according to the criterion jDp−DpjNnrDp , where P Dp is the mean difference over the track. The factor n was defined heuristically to assure that obvious outliers were eliminated: n = 10 when using the relocation method of SE correction, and n = 5 when using the direct method. About 2% of the data were excluded in both cases. The remaining data are used to analyse and correct offsets of individual tracks. First, we compute the f mean difference Dp over all observations in our region. This mean difference likely reflects an offset of the reference model (Section 3.1), e.g., due to a different SE correction applied in the f model generation. That is, Dp does not reflect an error of the observations. For every track, an offset is then computed by P f Dp− Dp, that is, by the difference between the individual track and the mean level of all observations. The data of individual tracks are then corrected for this track-related offset.

of freedom of the fit and its robustness (see Section 4.9 for a discussion of this choice). In view of the relatively smooth ice surface topography and the additional smoothing in the PLF scale which is inherent to radar altimetry this approach is appropriate. The spatial scale of the fit should be large enough to smooth observation noise but small enough to preserve the resolvable topographic variations. While previous studies have fixed this scale, we propose to keep it variable in order to adapt it pointwise to the specific local conditions. For every grid point, Eq. (9) (with x and y centered in the grid point) is fitted to the surrounding observations within a variable radius Ropt. The radius Ropt is chosen according to an optimal criterion after trying a whole sequence of radii R. For each R the leastsquares fit yields the estimated grid-point height h(R), the standard deviation σ0(R) of the residuals between the observations and the model surface, the cofactor q(R) describing the robustness of the observation and the standard pffiffiffiffiffiffiffiffiffidistribution, ffi deviation rh ðRÞ ¼ r0 ðRÞ qðRÞ of the estimated grid-point height. Here only those radii R are considered that fulfill the following three conditions: (1) R N 1.2 km, corresponding to the radius of the PLF, i.e., to the altimeter resolution limits; (2) R b 15 km as a somewhat arbitrary but sufficiently high upper border for the range of radii to be tested; (3) robustness of the fit attested by redundancy (i.e., more than 6 observations), a cofactor q b 1.5 and a moderate condition number of the normal equation matrix. Among the radii R, we finally choose the one with the minimal residual standard deviation σ0(R) of the fit. This criterion is based on considerations illustrated in Fig. 5 which shows a typical example over the flat area above Lake Vostok. As R and, hence, the number of included observations increases, the cofactor q attests an improved robustness of the fit. As long as the model surface fits the actual surface within R, σ0(R) provides an estimate of the data noise. This estimate becomes more stable for increasing radii. However, above a certain R the surface model (Eq. (9)) does not fit the actual

3.4. Gridding Various methods have been applied in ice sheet altimetry to derive a regular grid from the noisy and irregularly distributed observations: e.g., the computation of local averages (Rémy et al., 1999), the local adjustment of mathematical surface models (Ihde et al., 2002; Zwally et al., 1990) or stochastic methods (see Herzfeld (1992) for an overview) like collocation (Ekholm, 1996; Rémy et al., 1989) and kriging (Fricker et al., 2000; Herzfeld et al., 1993, 2000). Here we fit a mathematical surface model to the observations surrounding a grid point and assign the resulting model height to that grid point. We take the surface model of Eq. (9) as a compromise between the degrees

Fig. 5. Parameters of the least-squares fit of Eq. (9) to the observations surrounding a point on Lake Vostok, depending on the area radius R: cofactor q, approximation error σ0 in metres, scaled by factor 0.2, and error of the adjusted 1 height of the central point, rh ¼ q2 r0 in metres. For gridding, we choose the radius R0 for which σ0 is minimal.

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

275

Fig. 6. Illustration of the applied gridding procedure (here for the generation of the model VOSR): (a) gridding radius Ropt adapted to the local situation, (b) formal error σh(Ropt) of the derived grid-point height.

surface. Then, the residuals are the combined effect of data noise and model misfit and they increase as the misfit increases. The radius R with minimal σ0(R) is considered optimal because a smaller radius would cause a loss of robustness while a larger radius would likely cause unnecessary smoothing by fitting an inadequately simple surface model. Fig. 6 shows maps of the obtained radii Ropt and the corresponding height errors σh(Ropt). 4. Results, evaluation, interpretation The results of our investigation are the improved DEMs (Section 4.1) and the models of the additional altimetric parameters (Section 4.2) both derived from ERS-1 data. The generated DEMs are evaluated by a comparison with already existing DEMs (Section 4.3). All DEMs are validated by ICESat data and by ERS-1 observations. The validation method and the results are presented in Sections 4.4 and 4.5, respectively. In Section 4.6 the ERS-1 data and the derived DEMs as well as the ICESat data are analysed by a Fourier analysis which reveals, among other things, the error level of the used ERS RA data. In Section 4.7 this error level is alternatively derived from the analysis of ERS RA temporal variability at 35-day repeat tracks. Finally, in Sections 4.8, 4.9 and 4.10 the effects of the different improvements in radar altimetry analysis are assessed in detail by applying our validation methods. 4.1. Improved DEMs of the Lake Vostok region Two DEMs were generated using the two different methods of SE correction (Section 3.2): VOSR (Vostok region DEM from refined Relocation method) and VOSD (Vostok region DEM from refined Direct method). Fig. 7 shows the VOSR DEM.

Fig. 7. DEM for the region of Lake Vostok derived in this study. This model VOSR was generated with the refined relocation method of slope error correction and the refined gridding method. Ellipsoidal heights [m] w.r.t. the GRS80 ellipsoid. Contour intervals are 20 m, illumination is from east.

276

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

4.2. Additional altimetric parameters Maps of additional altimetric parameters, generated when applying the relocation method, are shown in Fig. 8. The backscattering coefficient Bs (Fig. 8a) varies from 11 to 16 dB from north to south which means that the reflection varies significantly across the area. This should be associated with (see Legrésy & Rémy, 1997, 1998): (1) differences in centimetre-tometre scale roughness due to different wind intensities (stronger wind leads to higher roughness which leads to a smaller Bs), (2) differences in the snow accumulation and deposition process through the intensity of volume echo controlled by stratification and snow grain size, and (3) surface undulations by a concentration of the radar echo from hollows and its dispersion over bumps. Bs is stable within 1 dB over the lake area, i.e., like the other altimetric parameters it shows the lake surface as a homogeneous flat target. The undulations as well as the flat areas are also clearly visible in the RADARSAT amplitude image (Fig. 1). The leading-edge width LeW (Fig. 8b) ranges from 1 m to more than 3 m. This parameter increases with the amplitudes of meter-scale snow dunes and surface undulations, and it also increases with radar penetration effects (Legrésy & Rémy, 1997, 1998). LeW shows a south–north gradient. As a possible explanation, snow dunes may be more frequent in the northern part due to an increase in wind intensity from south to north. This increase could, accordingly, also explain the south–north decrease of Bs as discussed above. The trailing-edge slope TeS (Fig. 8c) shows a background value about −4 μs− 1. As noticed by Legrésy et al. (1999) and

Rémy et al. (1999), this parameter provides a good measure of undulation amplitudes at the 10 to 20 km scale. Its absolute value is higher in hollows and lower at bumps. In our region, and in particular in the vicinity of the grounding line, this parameter clearly reveals this kind of undulation (cf. also Fig. 1). 4.3. Comparison to existing models For comparison with our results we use previously published DEMs of the area of investigation. All DEMs used are described in Table 1. The DEMs based on ERS altimetry are chosen to allow a comparison of the different processing approaches. Two DEMs based on local airborne measurements will be used to demonstrate the equivalent accuracy level of our DEMs. To compare these DEMs with each other their differences Δh in a certain grid are computed by using bi-linear interpolation. From Δh, the relative offsets as well as the corresponding standard deviations are derived, which are summarised for all combinations in Table 2. The mean offsets are a combined effect of systematic reference frame differences (height datum, orbits), data type (satellite/airborne, radar/laser) and retracking methods. Regarding the ERS-based DEMs, the models IHDE, REMY, VOSD and VOSR are on a similar height level, whereas RAMP is about 0.5 m below and BAMB is about 1 m above this level. The airborne-based DEMs are on separate levels near BAMB. The offset between STU1 and STU2 (about 0.7 m) may be caused by different penetration of the airborne radar signal into the surface.

Fig. 8. Additional altimetric parameters obtained by the retracking procedure and mapped using the relocation method: (a) backscattering coefficient Bs, (b) leadingedge width LeW and (c) trailing-edge slope TeS.

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

277

Table 1 DEMs covering the Lake Vostok region: Six previously published DEMs used for comparison and the two DEMs obtained in this study Model

Observation type

≈ grid-cell size [km] at 77.5S (rank number)

SE correction method

Reference

BAMB IHDE RAMP REMY STU1 STU2 VOSD VOSR

ERS-1 RA⁎ ERS-1,2 RA ERS-1 RA⁎ ERS-1 RA Airborne LA Airborne RA ERS-1 RA ERS-1 RA

5.0 × 5.0 (2) 4.8 × 5.6 (1) 1.0 × 1.0 (4) 0.8 × 3.7 (3) 0.6 × 0.6 (7) 0.6 × 0.6 (8) 0.8 × 0.9 (5) 0.8 × 0.9 (6)

Relocation (based on ‘effective’ slope) Relocation (based on ‘effective’ slope) Direct (based on ‘effective’ slope) Intermediate (based on ‘effective’ slope and curvature) – – Direct, improved (based on ‘effective’ slope and curvature) Relocation, new (based on a-priori DEM)

Bamber and Bindschadler (1997) Ihde et al. (2002) Liu et al. (1999, 2001) Rémy et al. (1999) Studinger et al. (2003) Studinger et al. (2003) This study This study

⁎In the Lake Vostok region only RA data are used. The acronym RAMP stands for Radarsat Antarctic Mapping Project, the names BAMB IHDE, REMY, STU1 and STU2 are derived from author names. The rank number corresponds to the order of the DEMs, sorted by their grid-cell sizes.

The standard deviations (σΔh) are affected by the roughness (resolution) of the DEMs and their errors. The largest values (up to more than 2 m) arise for comparisons including RAMP or STU2 despite their relatively small grid-cell sizes (Table 1). The smallest values for σΔh (0.66 m) are obtained for the differences between VOSR and VOSD. This demonstrates their consistency and similar roughness representation. Fig. 9 shows a map of the difference VOSR-VOSD. As another example, Fig. 10 shows the differences between REMY and VOSR which are of similar resolution and based on the same dataset. Apart from the effects of an outlier track which was not removed in the REMY model generation, differences at small spatial scales are visible over most parts of the region. Our validation of the DEMs (next sections) will reveal that these differences reflect topographic features which are contained in the improved DEMs (e.g. VOSR) but not in previous models like REMY.

laser beam provides an illuminated surface spot of only about 65 m in diameter spaced at 172 m alongtrack. That means that an ICESat measurement is close to a point measurement compared to an RA observation. Second, the laser signal does not penetrate the ice surface and, therefore, provides pure surface reflections. We use ICESat data as independent information to validate our radar altimetry results. We work with the GLA12 records of the Antarctic and Greenland Ice Sheet Altimetry Data release 18

4.4. Validation with ICESat observations The ICESat altimeter (Schutz et al., 2005) operates in the optical part of the electromagnetic spectrum. Hence, its characteristics differ essentially from those of the ERS RA. First, the Table 2 Intercomparison of our DEMs and previously published DEMs: Mean P difference Dh (lower triangular matrix) and standard deviation σΔh (upper triangular matrix) of the height differences between different DEMs (units: meter) BAMB IHDE RAMP REMY STU1 STU2 VOSD VOSR (2) (1) (4) (3) (7) (8) (5) (6) BAMB (2) IHDE (1) RAMP (4) REMY (3) STU1 (7) STU2 (8) VOSD (5) VOSR (6)

– 1.12 − 1.04 – 1.50 0.51 0.97 −0.08 − 0.41 −1.47 0.29 −0.79 1.10 0.05 1.20 0.16

1.77 2.42 – 0.58 −1.63 −0.93 −0.40 −0.30

0.98 0.95 1.76 1.41 1.46 2.02 1.94 1.94 2.35 – 1.14 1.84 − 1.27 – 1.51 − 0.57 0.71 – 0.12 − 1.46 − 0.77 0.22 − 1.55 − 0.86

0.73 1.92 1.85 1.21 1.18 1.88 – 0.10

0.72 1.88 1.94 1.24 1.13 1.84 0.66 –

For each pair of grids, the finer grid is interpolated to the gridpoints of the coarser grid, and the difference finer grid value minus coarser grid value is taken. The rank numbers in brackets (cf. Table 1) indicate the hierarchy of resolutions (higher number means finer grid).

Fig. 9. DEM differences: VOSR-VOSD. The standard deviation is σΔh = 0.66 m.

278

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

Fig. 11. Statistics of ICESat crossover differences. (a) offsets detected for each individual track. (b) rms of crossover differences before (grey) and after (black) correction of the offsets. See Table 3 for mean values.

Fig. 10. DEM differences REMY-VOSR (after subtraction of the mean difference). The standard deviation is σΔh ≈ 1.24 m. The undulations of the difference are mainly caused by topographic details that are contained in VOSR but not in REMY (cf. Fig. 12). In addition, outlier track 15,468 (crossing the lake) corrupts the REMY model.

(20 Feb.–18 Nov. 2003) (Zwally et al., 2003).3 Of course, the ICESat data also contain errors (Schutz et al., 2005; Shuman et al., 2006). But after the data screening and offset correction described below, our analyses and discussions (this section, Sections 4.6 and 5) give confidence that the ICESat data have a superior sensitivity in particular to smaller-scale topography and are therefore suited for validating RA data. We referenced the ICESat heights to the GRS80 ellipsoid, analysed their precision and eliminated offsets and outliers: In a crossover analysis we computed height differences at crossover points within the Lake Vostok region (95°E to 115°E and 74°S to 81°S). We determined height offsets of the individual tracks by a least squares method. Fig. 11 shows the rms crossover difference and the obtained offset for each track. The offsets are mainly of the order of a few decimeters, but in one case the offset is almost 2 m. These offsets were then removed from each track. Fig. 11b illustrates the effect of this procedure on the rms crossover differences. Irrespective of their offsets, some tracks also show a significantly higher noise in the crossover differ3

Download http://www.nsidc.org/data/docs/daac/glas_icesat_l1_12_global_ altimetry.gd.html. The material is based on work supported by the National Science Foundation under Grant OPP-9911617.

ences. Possible causes are effects of laser beam mispointing or atmospheric influences (Shuman et al., 2006). These ‘highnoise’ tracks (6 out of 214, corresponding to 3% of the data) were removed completely. Table 3 shows the effects of offset correction and elimination of ‘high-noise’ tracks. The rms crossover difference is reduced by more than a factor of 2 to 13 cm. When attributing these differences entirely to the measurement precision,pthis ffiffiffi single-shot precision estimate would be about 9 cmc 12 2d13 cm. Note, however, that the crossover differences contain interpolation errors in addition to observation errors. This could lead to an overestimation of the ICESat measurement errors. For the validation with ICESat data, surface heights for individual ICESat observation locations were generated from the DEMs by bi-linear interpolation. Then the differences between ICESat heights and DEM heights are obtained for all P positions, and finally their mean value Dh and their standard deviation σΔh are computed. Note that the ICESat dataset consists of direct observations on their original positions whereas the radar-based heights are Table 3 Crossover statistics of ICESat data: effect of offset correction and elimination of ‘high-noise’ tracks

Complete dataset Without ‘high-noise’ tracks

Before offset correction

After offset correction

rms(Δh)

max(|Δh|)

rms(Δh)

max(|Δh|)

0.30 0.22

4.12 2.07

0.22 0.13

3.30 0.86

Shown are the rms(Δh) and the maximum absolute values max(|Δh|) of crossover differences (units: meter).

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

interpolated. For an interpretation of the validation results we therefore have to distinguish between two cases: Typically, the RA observations are dense enough so that the interpolation preserves or even improves the quality of the original observations. In this case (which we term Case 1) the comparison with ICESat basically provides a comparison between both types of observations. Differences are expected to occur due to the following facts: (a) ARA measures a mean height over the PLF and hence provides the surface height smoothed at the scale of a few kilometres. ICESat, in contrast, also detects the smaller-scale variations. Therefore, differences between both datasets partly reflect the small-scale height variability; (b) Since the radar reflection contains subsurface contributions (volume scattering) while the laser reflection does not, ICESat heights can be expected to be systematically above the heights as obtained by ERS. In some regions, topographic depressions may occur in scales between the PLF and the BLF scales and with a radius of curvature smaller than the satellite height. For such depressions the bottom is never observed by the RA since the PLF is always displaced upslope. In this case (termed Case 2) the comparison to ICESat basically provides the interpolation error for an unobserved position. The statistics of the differences (ICESat-DEM) are different for the Cases 1 and 2 and should therefore be discussed separately. Here, we attempt to isolate Case 1, that is, to validate the DEM in areas where it is based on actual RA observations. We assume a Gaussian distribution for Case 1 and a higher standard deviation for Case 2 than for Case 1. Therefore we assign those P ICESat data to Case 1 for which jDh−DhjV3rDh where Δh is the P ICESat-DEM difference, Dh its mean over all ICESat data and σΔh its standard deviation. To check the effectiveness of this selection, we looked at whether those ICESat ground positions that were excluded, by applying the 3σΔh criterion with respect to the VOSR DEM (3.7% of the initial ICESat data), really are located in topographic depressions. We roughly delineate depressions by considering that differences between VOSR and VOSD (Fig. 9) should be particularly large at depressions because, at such locations, both DEMs arise from different kinds of interpolation instead of actual RA observations. We take a difference larger than 1 m (in absolute value) as the criterion for a depression. Then, 62% of the excluded ICESat data lie, indeed, in depressions, although the depressions make up only 3.6% of the total area. So, the DEM accuracies inferred from the ICESat validation are valid with the exception of the clearly identifiable areas of Case 2. Naturally, the 3σΔh criterion also eliminates data at the margin of the histogram of Case 1. According to the above numbers, about 38% × 3.7% = 1.4% extreme values of Case 1 are eliminated which, under the normal distribution assumption, reduces the standard deviation by only about 5%. In order to avoid favouring any particular DEM, the described selection of ICESat data is performed by comparison with all six spaceborne radar-based DEMs listed in Table 1. The ICESat data that passes all six selections (90.5% of the initial dataset) are then used for the validation of all DEMs.

279

Fig. 12 illustrates the validation of the model VOSR (developed in this study) and the model REMY by ICESat data in a small section south of the Lake Vostok. VOSR contains more small-scale variations than REMY. The validation shows that these variations are real topographic features. Fig. 13 summarises the validation for our DEMs and for a number of previously published DEMs. Our refined analysis significantly improves the topography recovery from RA observations. In particular, the standard deviation of differences between ICESat data and the DEM is only 0.54 m for VOSD and 0.51 m for VOSR while it is above 1 m for all other DEMs. Note, that the standard deviations are significantly smaller in the flat area over the lake. But even for this flat area, Fig. 13 clearly indicates that the RA data contain more small-scale information than obtained up to now and that we are able to partly resolve this information using our new methods.

Fig. 12. Section of the DEMs VOSR (a) and REMY (b) and their validation with ICESat. The colour code and contours (2 m interval) show the DEMs. Differences between the ICESat observations and DEMs are superimposed orthogonally to the ICESat tracks as white bars to the right and black bars to the left for positive and negative differences, respectively. For better legibility, only every second ICESat position is shown. The standard deviations σΔh of the differences in the area shown are 1.0 m for VOSR and 2.1 m for REMY.

280

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

4.6. Fourier analysis

P

Fig. 13. Validation of DEMs with ICESat data. Mean offset Dh and standard deviation σΔh of the ICESat-DEM height differences. Black: values for the whole area; grey: values for the lake area as delineated in Fig. 10.

4.5. Validation with ERS-1 observations We also performed an alternative, more internal, validation of the DEMs using the individual ERS-1 observations themselves (i.e., surface heights at particular positions). The procedure is the same as for the validation by ICESat, including the outlier and offset elimination described in Section 3.3. Concerning the SE correction we use the relocation method with the reference topography VOSR, except for the validation of VOSD. There, to be consistent with the method of model generation, we use the direct method with the reference topography VOSD. Fig. 14 summarises the validation by ERS-1 observations for our DEMs and for the set of previously published DEMs. The results (Fig. 14) support the findings of Section 4.4. The standard deviation for VOSD and VOSR is in the order of 0.5 m. As will be shown in the next two sections, this is the noise level of the RA observations. Comparable results were obtained by Bamber et al. (1998) over flat areas in Greenland. Over the lake this error level is almost achieved also by BAMB, REMY and STU1. Over the whole area, however, the standard deviations are above 1 m for all previous DEMs.

Fig. 14. Validation of DEMs using the original ERS-1 data (SE corrected). The graphs are analogous to those in Fig. 13.

To compare the spectral content of the ERS-1 and ICESat altimetric height observations and of our DEM we analyse ERS1 and ICESat tracks within our region. Along all tracks with a minimum length of 400 km we performed a Fourier analysis of (1) the altimetric observations, (2) our DEM VOSD interpolated to the same positions, and (3) of the differences between both. Prior to the Fourier analysis the ERS-1 observations were SE corrected by the direct method, maintaining a regular sampling. We chose the DEM VOSD to be consistent with the applied SE correction. All data were high-pass filtered. Fig. 15a shows the results. The mean spectrum of ERS-1 observations (Fig. 15a) exhibits a constant level for wavelengths λ b 7 km which indicates the noise level of the RA observations with a standard deviation of about 0.5 m (compare Fig. 14 and Section 4.7). The spectrum of the DEM shows that the gridding process effectively reduces this noise for wavelengths λ b 7 km. On the other hand the DEM retains virtually all RA signal at wavelengths λ N 10 km: At these wavelengths the spectrum of

Fig. 15. Power Spectral Density (PSD) of satellite altimetry observations (along their respective satellite tracks), the DEM VOSD (interpolated to the same tracks) and the difference between observations and the DEM (a) for ERS-1 and (b) for ICESat. The spectra of the DEM differ in subfigure (a) and (b) because they are computed for different sets of tracks crossing the region. The minimum wavelength is twice the sampling interval.

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

281

cycle. Since the altimeter design is identical for both satellites, the obtained results are applicable for ERS-1, too. To account for the small measurement position differences between different overflights, bicubic fits of the local topography were subtracted. Moreover, a temporal trend and a seasonal cycle were fitted and subtracted from the data. The residual variability is shown in Fig. 16. It includes measurement noise due to both instrumental noise and influences of complex topographies which induce complex waveforms. It also includes actual surface height variability, but from field observations at Vostok station it appears that no significant surface height variation can be expected for the underlying 8 year period. Therefore, the revealed variability is a good proxy of the noise level of the individual retracked range observations. This noise level varies around the value of 0.5 m. This is the same value as obtained in the previous section from the analysis of along-track variability. Spatial variations in Fig. 16 appear related to variations in topographic complexity as revealed, e.g., in the gridding error map of Fig. 6b. Our interpretation is that a more complex topography induces more complex waveforms so that their retracking contains higher noise. Fig. 16. Temporal variability of ERS-2 altimetric height measurements.

4.8. Effect of slope error correction refinement the ERS-1-DEM difference is only a small fraction of the signal. The spectrum of the ICESat observations (Fig. 15b) shows a plateau for wavelengths λ b 1 km which indicates a white-noise level in this spectral range with a standard deviation of about 0.04 m. For λ N 1 km the spectrum of the ICESat observations is likely dominated by the topographic signal. The DEM does not contain this signal at 1 km b λ b 7 km because due to the RA noise level this signal could not be extracted from the RA observations and consequently was damped by the gridding algorithm. For 7 km b λ b 15 km the ICESat observations still contain signal that is not present in the DEM. In contrast, Fig. 15a shows that the ERS-1 observations differ relatively little from the DEM at these wavelengths. We conclude that in this spectral band the reduced energy of the DEM is not mainly due to the gridding but due to the fact that ERS-1 observations themselves represent heights averaged over the scale of several kilometres and thus provide a damped topographic signal. For λ N 15 km the DEM contains the bulk of the signal present in the ICESat observations. 4.7. ERS radar altimeter precision from temporal variability analyses Here we estimate the precision of the retracked ERS RA observations from their temporal repeatability. For this aim we analyse time series of 80 35-day repeat cycles of ERS-2 in the Lake Vostok region. We use ERS-2 because it provides a more comprehensive time series: ERS-2 continuously flew in the 35day repeat configuration while ERS-1 contained the geodetic mission phase. ERS-2 always used the ‘ice mode’ of waveform sampling, which is the sampling mode of the ERS-1 geodetic phase, while ERS-1 was in ‘ocean mode’ every second 35-day

Here we want to assess the effect of our refinements in the SE correction methods. For comparison, we derived DEMs by using the ‘classic’ direct method and the ‘classic’ relocation method while maintaining all other elements of our processing. In the case of the direct method, the difference between ‘classic’ and refined version is whether to fix or to locally adapt the size of the area in which Eq. (9) is fitted for the derivation of surface slope and curvature parameters. In the case of the relocation method, the ‘classic’ version, again, utilizes the mathematical model surface of Eq. (9) fitted over a constant area. The refined method, instead, does not rely on such a polynomial model surface. As the fixed area size in the ‘classic’ case, we take 14 km × 14 km because, given the maximal impact-point displacement of 6 km in our region (see end of Section 3.2.2), the PLF is always within this area. Table 4 compares the results of the four different methods by their validation with ICESat data. Clearly, our methodological refinements improve the results of the ‘classic’ method. Presumably, the results of the ‘classic’ method would be even worse if the fixed scale was chosen larger than 14 km × 14 km.

Table 4 Effect of SE correction refinements as compared to the ‘classic’ SE correction methods: standard deviation differences between ICESat data and DEMs generated with different SE corrections σΔh [m]

Slope-error correction method Direct method, Direct method, Relocation method, Relocation method,

‘classic’ ‘new’ ‘classic’ ‘new’

0.86 (0.29) 0.54 (0.24) 0.54 (0.25) 0.51 (0.24)

Here and in the following tables values in brackets are for the lake area as delineated in Fig. 10 while otherwise the values are for the whole area.

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

Table 5 Comparison of different surface approximations (two-dimensional polynomials of order 0 to 3) tested for gridding Polynomial order

0 0 1 2 3

Number of parameters

P r 0 of the altimetric parameter:

h

Bs

LeW

TeS

[m]

[dB]

[m]

[μs− 1]

σΔh [m]

1 1⁎ 3 6 10

3.45 0.21 0.71 0.43 0.97 (0.33) 3.36 0.20 0.72 0.43 0.93 (0.32) 0.92 0.16 0.67 0.32 0.88 (0.32) 0.48 0.15 0.61 0.27 0.55 (0.24) 0.36 0.14 0.56 0.25 0.56 (0.23) r 2 Þ ⁎Observations weighted by 10ð10km where r is the distance to the grid point. P of the approximation error σ (see Section Columns 3 to 6: quadratic mean r 0 0 3.4) over all grid points. Last two columns: DEM validation by ICESat (standard deviation difference ICESat-DEM). To make the results comparable, the gridding radius was fixed to 4 km, here.

Such a choice would be necessary for continental analyses where larger impact-point displacements occur. Regarding the results from the two refined methods, their precision level is about the same (Table 4, Figs. 13 and 14). The direct comparison of VOSR and VOSD (Fig. 9) reveals that in general the pattern of height differences is correlated to the pattern of impact point displacements (Fig. 4). In flat regions (especially over the lake) both methods obtain practically the same DEM. The largest differences, in contrast, appear over concave features ‘invisible’ for a RA. With the direct method, the SE corrections are uncertain at those places. With the relocation method, the observations are shifted away from those places so that the interpolation of the gridding procedure causes large errors (Fig. 6b). This general problem can only be solved by combination of RA data with observations of higher spatial resolution (eg. ICESat or local GPS profile measurements). 4.9. Effect of gridding method refinement In Section 3.4 we explained that we compute the model height at a grid point from a polynomial fit to the surrounding observations, and a 2nd order two-dimensional polynomial fit (Eq. (9)) was chosen. This choice was made from alternative two-dimensional polynomial fits ranging from order 0 to 3. Table 5 compares the results for these fitting functions and supports the choice of the 2nd order polynomial of Eq. (9). A Table 6 Effect of adapting the DEM gridding radius to local surface characteristics Mapping radius R

σ ¯0 [m]

σΔh [m]

Ropt 3 4 5 10

0.32 0.36 0.48 0.66 1.67

0.54 (0.24) 1.95 (0.24) 0.55 (0.24) 0.59 (0.25) 1.19 (0.43)

Solutions with fixed radii are compared to the solution with the adapted radius Ropt (first line) which ranges from 2 to 15 km with a mean value of about 4 km. σ ¯0 and σΔh are the same as in Table 5. Here, the direct method of SE correction was used.

1.6 standard deviation [m]

282

1.4 1.2

X direct method,

ICESat ERS–1 direct method, relocation method, ICESat relocation method, ERS–1

X

1.0 0.8

X X

0.6

X

X

X

X

X

X

X

X

0.4 0

2

4 6 iteration step

8

10

Fig. 17. Effect of the iteration in the altimetric data processing: Standard deviations of differences with ICESat data and ERS-1 observations after each iteration step are shown for the two methods of slope error correction.

further increase of the number of parameters improves the fit to the (noisy) observations but rarely improves the resulting elevation model. We perform the fit within a radius R which is adapted to the local situation. Table 6 shows the effect of the choice of R. Since the method minimises the residual error σ0 of the local fit for every grid point, the quadratic mean of these residuals, σ ¯0 (2nd column) is minimal for an adapted radius compared to methods with fixed radius. But also the external validation with ICESat (3rd column) shows the advantage of the adapted gridding method. Although some particular fixed radii (e.g., 4 km) provide comparable results, statistical analyses would be, again, required to determine this ‘optimal fixed radius’ and it would vary from region to region. 4.10. Effect of iteration As mentioned, the processing scheme described in Section 3 is iterative. Starting with a smooth a-priori reference model derived without any corrections the DEMs are improved iteratively. Fig. 17 illustrates this improvement by the validation with ICESat and ERS-1 observations. After some iterations the results remain on a practically constant low level. These final results were given in Figs. 13 and 14. 5. Conclusions The present study reduces the error standard deviation of a satellite RA-based DEM of the Lake Vostok region from about 1.1 m to about 0.5 m. This is shown by comparison with independent ICESat observations (Fig. 13). The comparison of the DEM with individual ERS-1 observations (Fig. 14) shows that the standard deviation difference (about 0.5 m) is at the noise level of the individual ERS-1 observations as revealed in Sections 4.6 and 4.7. This indicates that our analysis is not far from extracting all available topographic information from the underlying retracked RA observations. Our refined relocation method for SE correction is very effective. It follows the principle of an RA by looking for the minimum range between the satellite and a surface area in the size of the PLF. Therefore it is more direct and robust than the common utilisation of more abstract ‘effective’ surface

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

parameters like slope and curvature. At the same time, the direct method of SE correction benefits from the impact-point displacement information obtained by our relocation method: Through introducing this displacement information, the ‘effective’ surface parameters which are still needed for the direct method can be inferred more reliably. To interpolate a regular grid from the noisy and irregular observations we used a polynomial fit (Eq. (9)). The result of this procedure was improved significantly by adapting the scale of the fit to the local conditions at each individual grid point. This proved a robust and numerically feasible refinement. Like inverse methods (Fricker et al., 2000; Herzfeld et al., 2000; Rémy et al., 1989) our approach provides some a-posteriori error estimate. Unlike inverse methods it uses only local stochastic characteristics and no regional or continental average stochastics. This is of advantage if average statistics do not adequately describe local conditions. A possible disadvantage is that spatially correlated residual errors may not be accounted for. As a future refinement of our relocation method, one could consider taking the additional waveform parameters into account. For example, the PLF size depends on surface roughness, and so does the leading-edge parameter LeW (Legrésy & Rémy, 1997). Hence, the PLF size assumed for relocating the range measurement could be locally adapted based on LeW maps. The ‘footprints’ affecting the backscattering coefficient Bs and the trailing-edge slope TeS are, strictly speaking, larger than the PLF since Bs and TeS refer to larger parts of the altimeter waveform. Hence, when relocating the observations of Bs and TeS, the assumed footprint size could be adapted accordingly. The backscattering coefficient of a surface element represents the relative weight of this element's contribution to the radar echo. That is, if Bs varies locally then the waveform and the retracked range are not simply related to surface geometry. Given the maps of Bs retrieved within the processing, one could account for the relative weighting of surface elements when interpreting the retracked range as the mean distance to the PLF and searching for the location of this PLF. In a similar way, one could account for the spatially varying influence of volume scattering for which Bs, TeS and LeW are sensitive. However, we would expect that the additional accuracy gain from such sophisticated refinements is small compared to the gain achieved in this study. The intrinsic limitations of satellite pulse-limited radar altimetry (no resolution of small-scale topographic features, mix of surface and volume echo, no observations at some positions) persist. Our comparisons with ICESat observations have illustrated these limitations. They can only be overcome by observations complementary to pulse-limited radar altimetry. Our DEM is a particularly good reference for the analysis of further RA data, and it will serve as a reference for the ENVISAT RA2 cycle by cycle ice validation, both for the height and for the other altimetric parameters. The DEM is also a good reference for any glaciological studies in the area. It is available to researchers on the website http://www.tu-dresden.de/ipg/vostok.html. The analysis presented here is a regional one. Therefore, the comparison of our DEMs with previously published continental DEMs restricted to our region may be slightly biased in favour

283

of our DEMs. However, our methods of SE correction and gridding (which are the main improvements) should work for other Antarctic regions without further adaptation. An essence of these methods is that they are self-adapting to local conditions. Only the data screening described in Section 3.3 would need adaptation to different regions or to the entire continent (e.g., concerning the threshold factor n, or the effectiveness of large-scale vertical error reduction by removing an offset). Hence, an extension of our methods to the whole Antarctic continent and to the enlarged data base of ERS-1, ERS-2, and ENVISAT observations seems feasible and worthwhile. In particular, studies related to small-scale topography would benefit from the resulting enhanced DEM. Of course, in the coastal regions with their large slopes, topographic modelling will be better served by the ICESat LA and by the CryoSat-2 RA with its interferometric synthetic aperture radar capabilities. Acknowledgements The research was supported by the German Research Foundation (DFG). Travel funding was made available by the German Academic Exchange Service (DAAD) within the PROCOPE (Programme de Coopération Scientifique) project. All four authors took part in the very productive two-way visits between LEGOS (Laboratoire d’Etudes en Géophysique et Océanographie Spatiale), Toulouse, and Technische Universität Dresden. Processed ERS RA data were provided by the OSCAR project (observing continental surfaces with radar altimetry). Comments from two anonymous reviewers were very helpful to improve the presentation of our results. References Bamber, J. L. (1994). Ice sheet altimeter processing scheme. International Journal of Remote Sensing, 15(4), 925−938. Bamber, J. L., & Bindschadler, R. (1997). An improved elevation dataset for climate and ice sheet modelling: Validation with satellite imagery. Annals of Glaciology, 25, 438−444. Bamber, J. L., Ekholm, S., & Krabill, W. B. (1998). The accuracy of satellite radar altimeter data over the Greenland ice sheet determined from airborne laser data. Geophysical Research Letters, 25, 3177−3180. Benveniste, J. (1993). Towards more efficient use of radar altimeter data. ESA Bulletin, 76. Brenner, A. C., Bindschadler, R. A., Thomas, R. H., & Zwally, H. J. (1983). Slope-induced errors in radar altimetry over continental ice sheets. Journal of Geophysical Research, C3(88), 1617−1623. Brooks, R. L., Campbell, W. J., Ramseier, R. O., Stanley, H. R., & Zwally, H. J. (1978). Ice sheet topography by satellite altimetry. Nature, 274(5671), 539−543. Davis, C. H. (1992). Satellite radar altimetry. IEEE Transactions on Microwave Theory and Techniques, 40(6), 1070−1076. Ekholm, S. (1996). A full coverage, high resolution, topographic model of Greenland, computed from a variety of digital elevation data. Journal of Geophysical Research, 101(B10), 21,961−21,972. Féménias, P., Rémy, F., Raizonville, R., & Minster, J. F. (1993). Analysis of satellite-altimeter height measurements above continental ice sheets. Journal of Glaciology, 39(133), 591−600. Fricker, H. A., Hyland, G., Coleman, R., & Young, N. W. (2000). Digital elevation models for the Lambert Glacier-Amery Ice Shelf system, East Antarctica, from ERS-1 satellite radar altimetry. Journal of Glaciology, 46 (155), 553−560.

284

S. Roemer et al. / Remote Sensing of Environment 106 (2007) 269–284

Fu, L. -L., & Cazenave, A. (2001). Satellite altimetry and Earth sciences: A handbook of techniques and applications. San Diego: Academic Press. Herzfeld, U. C. (1992). Least squares collocation, geophysical inverse theory, and geostatistics: A bird's eye view. Geophysical Journal International, 111 (2), 237−249. Herzfeld, U. C., Lingle, C. S., & Lee, L. (1993). Geostatistical evaluation of satellite radar altimetry for high-resolution mapping of Lambert Glacier, Antarctica. Annals of Glaciology, 17, 77−85. Herzfeld, U. C., Stosius, R., & Schneider, M. (2000). Geostatistical methods for mapping Antarctic ice surfaces at continental and regional scales. Annals of Glaciology, 30, 76−82. Ihde, J., Eck, J., & Schirmer, U. (2002). A digital terrain ice model of Antarctica derived by ERS radar altimeter data. Germany: Bundesamt für Kartographie und Geodäsie (BKG) Digital medium. Jezek, K., & RAMP Product Team (2002). RAMP AMM-1 SAR Image Mosaic of Antarctica. Fairbanks, AK: Alaska Satellite Facility, in association with the National Snow and Ice Data Center, Boulder, CO. Digital media. Legrésy, B., Papa, F., Rémy, F., Vinay, G., van den Bosch, M., & Zanife, O. -Z. (2005). ENVISAT radar altimeter measurements over continental surfaces and ice caps using the ICE-2 retracking algorithm. Remote Sensing of Environment, 95, 150−163. Legrésy, B., & Rémy, F. (1997). Altimetric observations of surface characteristics of the Antarctic ice sheet. Journal of Glaciology, 43(144), 265−275. Legrésy, B., & Rémy, F. (1998). Using the temporal variability of satellite radar altimetry observations to map surface properties of the Antarctic ice sheet. Journal of Glaciology, 44(147), 197−206. Legrésy, B., Rémy, F., & Schaeffer, P. (1999). Different ERS altimeter measurements between ascending and descending tracks caused by wind induced features over ice sheets. Geophysical Research Letters, 26(15), 2231−2234. Liu, H., C., J. K., & Li, B. (1999). Development of an Antarctic digital elevation model by integrating cartographic and remotely sensed data: A geographic information system based approach. Journal of Geophysical Research, 104 (B10), 23199−23213. Liu, H., Jezek, K., Li, B., & Zhao, Z. (2001). Radarsat Antarctic Mapping Project digital elevation model version 2.Boulder, CO: National Snow and Ice Data Center Digital media. Rémy, F., Mazzega, P., Houry, S., Brossier, C., & Minster, J. F. (1989). Mapping of the topography of continental ice by inversion of satellite-altimeter data. Journal of Glaciology, 35(119), 98−107. Rémy, F., Schaeffer, P., & Legrésy, B. (1999). Ice Flow Physical Processes derived from ERS-1 high-resolution map of the Antarctica and Greenland ice sheets. Geophysical Journal International, 139, 645−656. Scharroo, R., Visser, P. N. A. M., & Mets, G. J. (1998). Precise orbit determination and gravity field improvement for the ERS satellites. Journal of Geophysical Research, 103(C4), 8113−8127. Schutz, B. E., Zwally, H. J., Shuman, C. A., Hanock, D., & DiMarzio, J. P. (2005). Overview of the ICESat Mission. Geophysical Research Letters, 32, L21S01. doi:10.1029/2005GL024009 Shuman, C. A., Zwally, H. J., Schutz, B. E., Brenner, A. C., DiMarzio, J. P., Suchdeo, V. P., et al. (2006). ICESat Antarctic elevation data: Preliminary

precision and accuracy assessment. Geophysical Research Letters, 33, L07501. doi:10.1029/2005GL025227 Studinger, M., Bell, R., & Karner, G. (2002). Ice cover, landscape setting and geological framework of Lake Vostok, East Antarctica. Earth and Planetary Science Letters, 205(3–4), 195−210. Studinger, M., Bell, R. E., Karner, G. D., & Tikku, A. A. (2003). Ice flow, landscape setting, and geological framework of Lake Vostok, East Antarctica. Open-file report (pp. 107). Reston, VA, United States: U.S. Geological Survey. Tabacco, I., & Bianchi, C. (2002). Airborne radar survey above Vostok region, east central Antarctica: Ice thickness and Lake Vostok geometry. Journal of Glaciology, 48(160), 62−69. Zwally, H. J., Brenner, A. C., Major, J. A., Martin, T. V., & Bindschadler, R. A. (1990). Satellite radar altimetry over ice. Processing and corrections of Seasat data over Greenland, Vol. 1. (pp. )Washington, DC: National Aeronautics and Space Administration (NASA Ref. Pub. 1233.1). Zwally, H. J., Schutz, B., Abdalati, W., Absihre, J., Bentley, C., Brenner, A., et al. (2002). ICESat's laser measurements of polar ice, atmosphere, ocean and land. Journal of Geodynamics, 34, 405−445. Zwally, H., Schutz, R., Bentley, C., Bufton, J., Herring, T., Minster, J., et al. (2003). GLAS/ICESat L2 Antarctic and Greenland ice sheet altimetry data V018, 20. February to 18. November 2003.Boulder, CO: National Snow and Ice Data Center Digital media. Swen Roemer received the diploma in geodesy from Technische Universität Dresden, Germany, in 2000. Since 2000 he has been a research scientist at Institut für Planetare Geodäsie/Technische Universität Dresden. His fields of research are: observation of ice sheets by remote sensing (especially by altimetry) and kinematic GPS. From 2006 he is with GeoForschungsZentrum Potsdam (GFZ). Benoît Legrésy received his PhD from University of Toulouse in 1998. He received a postdoc grant from CNES from 1998 to 2000. Since then he is working as researcher for CNRS in LEGOS. He has been involved in remote sensing of ice sheets using microwave techniques. Martin Horwath (birth name: Wiehl) received the diploma in mathematics from Technische Universität Dresden, Germany, in 1998. Since 1998 he has been a research scientist at Institut für Planetare Geodäsie/Technische Universität Dresden. His fields of research are: the use of satellite gravity missions like GRACE for the study of geophysical mass variations with a focus on ice sheet mass balance, and radar remote sensing of ice sheets, including the use of reflected GPS signals. Reinhard Dietrich received his PhD in geodesy from Technische Universität Dresden, Germany, in 1977. From 1976 until 1992 he worked as a research scientist and later as a project group leader in Potsdam (Central Institute for Physics of the Earth, since 1992 at GFZ). Since November 1992 he is Professor for geodesy at the TU Dresden. His main research interest is concentrated on geodynamics and glaciology based on field observations and space techniques. He coordinated several geodetic-glaciological research projects in Greenland and Antarctica.