APPLICATION OF GLAS LASER ALTIMETRY TO DETECT ELEVA TION CHANGES IN EAST ANTARCTICA

In this paper the use of ICESat/GLAS laser altimeter for estimating multi-temporal elevation changes on p lar ice sheets is afforded. Due to non-overlapping laser spots during repeat pa sses, interpolation methods are required to make co mparisons. After reviewing the main methods described in the literature ( c ossover point analysis, cross-track DEM projection, space-temporal regressions), the last one has been chosen for its capability of prov iding more elevation change rate measurements. The standard implementation of the space-temporal linear regression technique has been revisited and improved to better cope with outliers and to check the estimability of model’s parameters. GLAS data over the PANDA route in East Antarctica h ave been used for testing. Obtained results have been quite meaningful from a physical point of view, confirming the trend reported by the lit rature of a constant snow accumulation in the area during the two past decade s, unlike the most part of the continent that has b een losing mass. 1 Corresponding author


INTRODUCTION 1.1 GLAS laser altimeter on-board of ICESat satellite
The GLAS (Geoscience Laser Altimeter System) sensor onboard ICESat satellite offered unprecedented elevation data of Earth surface topography from [2003][2004][2005][2006][2007][2008][2009], when it stopped working due to the failure of its primary instrument (NSIDC, 2013).The main aim of GLAS sensor was the measurement of elevation changes over ice sheets, especially for the two major polar regions of Antarctica and Greenland.Considering the difficulty of operating airborne laser scanning flights over such areas (Young et al., 2008), GLAS provided the most precise laser altimetry data so far.These have been exploited for many purposes (see Moholdt et al., 2010) but the production of digital elevation models (DEM) and the measurement of elevation changes are undoubtedly the most important achievements (Shepherd et al., 2012).DEMs derived from ICESat data can be considered as the highest-precision surface models of Antarctica today available (DiMarzio et al., 2008).Considering the global warming is melting the ice sheets (Kerr, 2012), the accurate evaluation of elevation change is the first step for computing ice mass balance.This can be successively adopted for estimating the contribution of ice-melting to sealevel rise in the next decades.This paper would like to focus on this aspect and in particular on the change detection techniques for comparing multi-temporal GLAS data.As well known, comparing laser scanning data is in general a non-trivial task, due to the imprecise relocation of laser beams and to the contribute of several error sources which very often result in a measurement uncertainty larger than elevation variations to evaluate (Scaioni et al., 2013).In the case of GLAS data, this problem is emphasized due to atmospheric disturbances, signal saturation, large spot size and the large sampling distance, especially across tracks.Consequently, different portion of grounds are measured during repeat passes and a pointwise comparison is not meaningful.The worse situation happens in significantly sloped terrains.The low resolution of GLAS data also prevents the use of DEM interpolation for comparison, being this solution often adopted in laser scanning practice to overcome the problem of point relocation.In Figure 1an example of data acquisition during the same ICESat GLAS campaign is reported, highlighting the Amery Ice Shelf region where the test reported in this paper is focused.Here the presence of large areas without any data is clear, especially in the coastal zones which are mostly interested by the larger elevation changes along time.Comparing DEMs obtained from ICESat data may be worth only in the case of very large study areas, like when working at the continental level (see e.g., Bamber and Gomez-Dans, 2005).The aim of this paper is to define a methodology for computing elevation changes from GLAS data over Antarctica.At the moment, a smaller study area has been selected, which is located in Mac Robertson Land (East Antarctica) on the west side with respect to Amery Ice Shelf.This region spreads around the so called 'PANDA' route from the Chinese polar station of Zhongshan to Dome Argus.The route has been investigated since 1997 by several scientific expeditions (Ding et al., 2011), which provided several datasets concerning surface and sub-surface data.The integration of elevation changes into this dataset is expected to improve the knowledge of the dynamics of ice flows and other glaciological processes in the area.Starting by reviewing some state-of-the-art methods for computing elevation changes (Subsect.0), the paper proposes some improvement to a method based on space-temporal regression, which is applied to the study area (Sect.2).

Methods for comparing GLAS data
The main methods that have been adopted to compute elevation changes in satellite laser altimetry relied on local spatial interpolations (crossover point analysis and cross-track DEM projections), or on the estimate of space-temporal regression.
In Moholdt et al. (2010) a review and some comparisons between different methods over the same area are reported.
Other approaches can be found in the literature, which are based on more involved techniques for time series analysis (Ferguson et al., 2004;Nguyen and Herring, 2005;Smith et al, 2005), but basically are similar to one of the three main methods.
The basic concept of the first two approaches is to interpolate the original laser points to compute comparable elevations in the same positions.With reference to Figure 2-a, in the crossover point analysis, one ascending and one descending tracks are compared in the points where the track lines overlap.
Here elevations (H A , H B ) are linearly interpolated on the basis of the two closest points in each track, then the difference dH AB = H B -H A is computed.This technique is the most precise (Brenner et al., 2007), because the short along-track spacing of two GLAS points allows comparing points which are very close in space.Consequently the interpolation error is quite small, being average slopes on ice-sheets quite smooth (less than 2°).For this reason, crossover point analysis is also used to evaluate the intrinsic data quality by using data collected in a close time (see Subsect. 2.2).On the other hand, the number of crossover points is quite limited, as can be seen in Figure 1, and this method is not suitable for detailed analysis of small regions.The cross-track DEM projection method is based on the use of an independent DEM for compensating the unknown topographic difference (dH DEM in Fig. 2-b) between different points (Slobbe et al., 2008).To limit the slope-induced error, which depends on the quality of the adopted DEM and the distance between points to compare, the analysis is restricted to repeat-tracks featuring a small separation between them (usually less than a few hundred metres, depending on the slope).
The third approach (space-temporal linear regression) tries to overcome the limitation of the previous techniques, extending the area where information can be retrieved and avoiding the use of an external DEM.Indeed, over large polar regions the most existing DEMs are based on data from radar altimetry, which suffers for large non-modelled errors in areas with steeper slopes (Bamber and Gomez-Dans, 2005).Being the GLAS data too sparse to provide a sufficient coverage for continental-scale DEMs, ERS 1-2 radar data have been also integrated to laser altimetry to produce a topographic model of Antarctica that exploits the higher precision of the latter and the better spatial resolution of the former (Bamber et al., 2009).On the other hand, this improved model is still not fully reliable in sloped areas.Howat et al. (2008) proposed a method based on the contemporary estimate of local topography and temporal elevation changes on the basis of a laser point dataset collected over time into a spatial window (see Fig. 2-c).A linear regression model is used for the Least Squares (LS) estimate of a local average plane (parameters α E and α N that correspond to slopes in East and North directions, respectively).Under the assumption of constant elevation change rate (dH/dt) in each window, model parameters are related through the model: In Eq. ( 1), dE i , DN i , and dH i are East, North mapping coordinates and ellipsoidal elevations referred to the average values in the selection window.The advantage of such method is to enlarge and to densify the number measurements, although it cannot output elevation changes but the corresponding velocities computed across more observation epochs.On the other hand, the interest of research is more oriented on detecting medium and long term trends than having precise changes on limited time span (Shepherd et al., 2012).This approach accomplishes independent analyses over windows.The position of each window is selected along a regular orthogonal grid, which may be aligned to the main track direction.In Moholdt et al. (2010), some details on the implementation of this technique are given, in particular on the need of removing outliers from the dataset by establishing an empirical threshold on the estimated residuals v i in Eq. ( 1) after LS regression.

GLAS/ICESat data characteristics
The ICESat/GLAS instrument operated 17 observation campaigns of approx.35 days each from 2003 to 2009.Orbits were planned to obtain repeated ground tracks (5-10 profiles are usually available within a ground swath of a few hundred metres).Unfortunately, some tracks are not complete due to thick cloud absorption.Three annual campaigns were accomplished from February-April, May-June, and October-November (see Table 1).The acquisition of altimetry data was accomplished by using 1064 nm laser pulse transmission and measuring the time delay of surface echo returns (Zwally et al., 2002).The foot-print size was 52x95 m until summer 2004 (Laser 1/2), and 47x61 m since fall 2004 (Laser 3).The spacing of laser points along ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey tracks is 172 m, while it is much larger across tracks (e.g., at 70°S Lat is about 20 km).The accuracy design of elevation measurements was about 15 cm over gently sloped terrains.Fricker et al. (2005) proved that accuracy about 5 cm could be reached, but this performance sharply dropped down over sloped terrains and in the case of unfavourable atmospheric conditions resulting in forward scattering and detector saturation (Yang et al., 2010).Different datasets can be downloaded from NSIDC ( 2013), whose selection depends on the application environment: ice, cloud, water land topography (Duong et al., 2007;Shtain and Filin, 2011).
In this study we have adopted the dataset GLA12-Level 2, release 32.No pre-processing and corrections of downloaded data have been carried out, being the aim of our investigation to apply existing satellite altimetry products.

APPLICATION OF SPACE-TEMPORAL REGRESSION
In some cases a high density of elevation change observations is required.This occurs for instance when small regions (e.g., single glaciers or coastal areas) have to be investigated.In the case of the PANDA route depicted in Figure 2, precise elevation changes are not mainly motivated by the extension (ca.263⋅10 6 km), but by the need of in-depth analyses focused on looking at correlations between diverse observations in the area (bedrock map, ice-flow velocity, slope, elevation changes, snow density, wind speed and others).
For this reason the attention was put on extracting denser, more precise and reliable information on elevation changes from GLAS data.

Method description
The method applied here was based on the space-temporal linear regression described in Moholdt et al. (2010) and reviewed in subsection 1.2.Some improvements have been introduced to better cope with errors in GLAS data and with approximations related to the adopted method.Moreover, the functional model described in Eq. ( 1) is based on the assumption that the same elevation change occurs within each sampling window.This hypothesis is quite acceptable in the study area, but some exceptions are expected in regions with deeper slope or rougher surface (e.g., due to snow dunes of crevasses).
The reasons for preferring this method to others were: (i) crossover point analysis gave too sparse results; (ii) no precise DEMs were available in the region of PANDA route, especially in the coastal region; and (iii) the chance to have a larger redundancy to better locate outliers and understand nonmodelled systematic errors.
The workflow of the adopted procedure is shown in Figure 4.In the first stage, GLAS data are prepared to be input to the regression estimate, including some basic tasks like coordinate transformation (from geographic to UPS mapping grid system) and point selection into the Region-Of-Interest (ROI).
Ellipsoidal elevations with respect to WGS84 ellipsoid are used.
In addition, the input data do not entail the track from which each point comes from.This information is recovered from the analysis of geometric and temporal point proximity.Track identifiers will be used in the following processing stages.
The ROI is subdivided into regular window sizing ∆E and ∆W.
The spacing between different windows (δE and δW) can be as large as the window size, can be smaller with partial overlap between windows, or can be larger with consequent decimation of sampled data.General considerations on methods for 'areabased' deformation measurement hold here and are not discussed (see e.g., Scaioni et al., 2013).ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey Data are selected into each sampling window.Some tests are applied to check the consistency of the dataset, according to the number of points, the number of involved tracks, and the time span of selected data.This check is aimed to avoid instability in the estimate of the LS solution, which is carry out in standard way after normalizing point coordinates with respect to their average values inside each selection window.Two different tests are adopted to cope with outliers.The first test is based on the use of an empirical threshold for removing blunders, as also reported in Moholdt et al. (2010).The expected number of such errors is quite small and they are removed to speed up the process.A second step in the outlier rejection has been included to iteratively reject smaller errors on the basis of standard data snooping (Baarda, 1968).This solution resulted in increasing the number of removed points.
A third group of tests concerned the statistical significance of the estimated parameters, and the analysis of correlations between them.In the case one of the slope coefficient α E and α N cannot be estimated with statistical significance, its value is set to zero and the solution constrained by introducing a weight matrix W c .The solution can be worked out as follows: ( ) where x ˆ is the solution vector, A is the design matrix computed on the basis of functional model in Eq. ( 1), W is the weight matrix of observations in vector dH.Here W has been set equal to unary matrix, because there was not any evidence of assigning different weights.The slope parameters have been also constrained when they showed high correlations with elevation change rate dH/dt (larger then 0.9).

Experimental results
The approach described in the previous subsection has been applied to the region around the PANDA route shown in Figure 3. Two series of 7 campaigns of GLAS data have been adopted.The first one included data gathered in the period between February-April, i.e., at the beginning of the austral autumn (hereafter referred to as 'Aut' dataset).The second period spanned between September-November (austral spring -'Spr').
A preliminary evaluation of the data quality had been carried out by computing elevation differences at crossover points by using all data belonging to the same campaign.Results are reported in Table 1, which also gives a description of the properties of different GLAS datasets.As can be seen, the data quality seems to improve from the first to last missions, as probable consequence of the adjustment of satellite and sensor parameters.The same processing method had been applied by the NSDIC to all datasets.Indeed, in any new data release, all existing campaigns are re-processed to maintain a homogenous data quality.Another important fact to note is the presence of very large maximum absolute values (several metres), which means some local disturbances exist.
The processing scheme was independently applied to datasets collected during different seasons of year.This choice was motivated by the presence of seasonal changes due to accumulation during winter and ablation in summer.The analysis has been first applied to the entire region around 'PANDA' route.In a following step, a smaller area around Dome Argus has been investigated in greater detail.
The selection of processing parameters has revealed to be a crucial point in the application of this technique.In Table 2 an overview of the adopted parameter settings per each experiment run is shown.
The first group of critical parameters incorporates the window size (∆E, ∆W) and the spacing between different windows (δE, δW).Also the orientation (α) of the grid adopted for window selection can be aligned to the axes of mapping coordinate system, or can be rotated to follow the direction of ICESat's tracks.Indeed, in Moholdt et al. (2010) the application of this method was based on quite small windows (700 m along track and a few hundred metres across track, depending on the separation distance between different tracks) whose direction was aligned to the track direction.A 50% overlap was used.No mention to the joint use of both ascending and descending tracks is done.On the other hand, the processed datasets concerned small glaciers only.Here a different strategy has been preferred to cope with larger areas over ice sheets.When processing the whole ROI along the PANDA route, a larger selection window has been setup.After different tests, a window size ∆E=∆W=2000 m with a partial overlap of 500 m has been used.This means that 25% of points is shared between two adjacent windows in both row and column directions.The direction of the grid driving the window selection has been established along the orthogonal axes of the mapping coordinate system.Both ascending and descending tracks have been included in the selection.This has resulted in some problems, as shown later.The adopted window size has been a trade-off between the assumption of constant elevation change inside each window, which established a superior limitation, and the need of a consistent set of point to allow a stable solution of the LS regression (lower limitation).A smaller window size was used to analyse the area around Dome Argus (see Table 2).This was motivated by some bias on the average estimated residuals that has been found after processing the whole ROI.
The second group of parameters (of selected points n min , min. no. of tracks n trmin , and min.time span ∆t min in a window) also has resulted to be quite influencing on the final outcomes.A too strict selection of these parameters has turned out to include a limited number of regions in the analysis of elevation change.
On the other hand, too loose values resulted in instability problems.The whole datasets 'Spr' and 'Aut' have been also processed by considering half the whole time span of GLAS data (2003-2006 and 2006-2009) As can be seen from Table 2, the space-temporal interpolation adopted here provided a large number of velocities dH/dt w.r.t.crossover point analysis reported in Table 1.

Analysis of the whole 'PANDA' route area.
The average computed velocity is slightly different (3.7 cm/y for 'Spr' and 6.4 cm/y for 'Aut', respectively).Although this parameter has to be carefully considered because of the large variability of estimated velocities in the ROI, the obtained values are meaningful.First of all, the higher value for the 'Aut' accounts for the major accumulation typical of the colder seasons.
Secondly, Zwally et al. (2005) reported results on the evaluation of mass-changes in Antarctica for the previous decade 1992-2002.Their outcomes highlight that in contrast with the general trend of the continent which is melting, the region East of Amery show a tendency to accumulate ice with a rate of ca. 5 cm/y.On the other hand, the variability of these results is quite evident from plots in Figure 5, including areas of accumulation and loss.In particular, two aspects should be noted.
The first concerned the presence of large value of estimated dH/dt on areas where ascending and descending tracks are overlapped.The problem is more evident in the northern area closer to the see, where the spatial density of observation is lower.This suggested that systematic errors largely affected some tracks, probably due to satellite orbit determination.In these areas an in-depth analysis of theoretical accuracy of the estimated velocity showed some values larger than the average one.While the average residual over all windows is quite small, some large biases (tens of cm) have been found in the region of Dome Argus.For this motivation, a specific analysis of that region in described in next paragraph 2.2.2.In a second stage, the analysis of the whole ROI by considering two distinct periods (2003-2006 and 2006-2009)  Similar results in terms of average velocities have been obtained in the two periods for 'Aut' and 'Spr' datasets.In the case of 'Spr,' estimated velocities are very close between them (4.1 and 3.9 cm/y) and also near to the value (3.

2.2.2
Analysis of the Dome Argus region.The area around Dome Argus is the one presenting the steepest slopes (1°-2°) in the region of 'PANDA' route.Moreover, this is also the area with major snow accumulation, being nearer to the South Pole.To better investigate this region, an independent analysis has been performed, after tuning the parameter settings to deal with the smaller dataset.Data from 'Spr' campaigns spanning the whole period 2003-2009 have been exploited.The results have been quite satisfying, either in term of number of elevation change rates computed (623), and in term of precision.An average dH/dt=5.6cm/y was found, which resulted to be similar to values found in other datasets when using the same 'Spr' campaigns.On the other hand, the slightly larger value obtained in Dome Argus is motivated by the major snow accumulation here.Dispersion of velocities is 2-3 times smaller than found from other dataset, as expected considering the smaller size of the study area.The most significant achievement consisted in the improvement of the average residuals, with a significant reduction of large biased points.ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey

FINAL DISCUSSION AND CONCLUSIONS
This paper has discussed the application of the space-temporal linear regression method proposed by Howat et al. (2008) for the analysis of elevation change on Antarctic ice-sheets.Some improvements have been introduced into this technique to increase the capability of coping with outliers.Application of ICESat/GLAS laser altimeter data over a region in East Antarctica demonstrated the method can provide a denser map of elevation change rate estimate than other existing techniques.On the other hand, results are still locally influenced by anomalies in laser data.Improvements are needed to incorporate in the method some empirical models to compensate for errors (e.g., orbit errors) and to identify and discard noisy tracks.The use of the latest data release (33) of GLA12 products and the correction of inter-campaign biases recently investigated by Borsa et al. (2013) are both expected to improve the quality of results.
At the current state-of-art, regressions are independently computed in running sampling windows.The cross-comparison of different results in a kind of global adjustment could be used for estimating correction models.Improving the robustness of this technique would open the chance to be applied also to radar altimeter data, which are much noisy and affected by relevant systematic errors.On the other hand, the better spatial resolution of radar altimetry and the current operation of Cryosat satellite draw the attention on these kinds of data and on methods to use them for estimating accurate elevation changes.
Eventually, the analysis of Dome Argus suggested that the adopted space-temporal linear regression method requires a precise tuning of parameter settings to cope effectively with specific characteristics of different regions.
The paper also showed the analysis of GLAS data over the PANDA route in East Antarctica.Although the quality of results must be improved, they confirmed the trend of accumulating mass typical of this area in East Antarctica during the last decades.

Figure 2 .
Figure 2. Three main methods adopted for calculating elevation changes (figure from Moholdt et al., 2010); a) crossover point analysis; b) cross-track DEM projection; and c) space-temporal linear regression over regularly spaced windows.

Figure 3 .
Figure 3.The location of the study area (with red borders) along the PANDA route from Zhongshan Station to Dome Argus in East Antarctica.

Figure 4 .
Figure 4. Workflow of the space-temporal linear regression method for estimating elevation change rate over icesheets.

Figure 5 .
Figure 5. Plots of estimated elevation change rates (in cm/y) over the whole ROI.In the upper and lower plots the results obtained from 'Spr' and 'Aut' campaigns are shown, respectively.Colours of different points represent the estimated average velocity, according to the scale in colorbars on the right.

Figure 6 .
Figure 6.Plots of estimated elevation change rates (in cm/y) in the area of Dome Argus.

Table 1
. Results are separately discussed in following paragraphs.ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey . Characteristics of different campaigns of GLAS data adopted in experiments and results of crossover point analysis to check data quality ('tr': tracks; 'max abs': max absolute value).

Table 2
. Dataset properties, input parameter settings and results of experiments described in Subsect.2.2.Statistics of results have been averaged over all windows (PAN: 'PANDA' route area; 'D-A': Dome Argus area; 'p.w.': per window).