Fusion of Laser Altimetry Data with Dems Derived from Stereo Imaging Systems

Abstract. During the last two decades surface elevation data have been gathered over the Greenland Ice Sheet (GrIS) from a variety of different sensors including spaceborne and airborne laser altimetry, such as NASA’s Ice Cloud and land Elevation Satellite (ICESat), Airborne Topographic Mapper (ATM) and Laser Vegetation Imaging Sensor (LVIS), as well as from stereo satellite imaging systems, most notably from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Worldview. The spatio-temporal resolution, the accuracy, and the spatial coverage of all these data differ widely. For example, laser altimetry systems are much more accurate than DEMs derived by correlation from imaging systems. On the other hand, DEMs usually have a superior spatial resolution and extended spatial coverage. We present in this paper an overview of the SERAC (Surface Elevation Reconstruction And Change detection) system, designed to cope with the data complexity and the computation of elevation change histories. SERAC simultaneously determines the ice sheet surface shape and the time-series of elevation changes for surface patches whose size depends on the ruggedness of the surface and the point distribution of the sensors involved. By incorporating different sensors, SERAC is a true fusion system that generates the best plausible result (time series of elevation changes) a result that is better than the sum of its individual parts. We follow this up with an example of the Helmheim gacier, involving ICESat, ATM and LVIS laser altimetry data, together with ASTER DEMs.


INTRODUCTION
During the last two decades surface elevation data have been gathered over the Greenland Ice Sheet (GrIS) from a variety of different sensors including spaceborne and airborne laser altimetry, such as NASA's Ice Cloud and land Elevation Satellite (ICE-Sat), Airborne Topographic Mapper (ATM) and Laser Vegetation Imaging Sensor (LVIS) systems, as well as from stereo satellite imaging systems, most notably from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Worldview.The spatio-temporal resolution, the accuracy, and the spatial coverage of all these data differ widely.For example, laser altimetry systems are much more accurate than DEMs derived by correlation from imaging systems.On the other hand, DEMs usually have a superior spatial resolution and extended spatial coverage.
We have developed the SERAC (Surface Elevation Reconstruction And Change detection) system to cope with the data complexity and the computation of elevation change histories, (Schenk andCsatho, 2012, Schenk et al., 2014).SERAC simultaneously determines the ice sheet surface shape and the time series of elevation changes for surface patches whose size depends on the ruggedness of the surface and the point distribution of the sensors involved.By incorporating different sensors, SERAC is a true fusion system that generates the best plausible result (time series of elevation changes) a result that is better than the sum of its individual parts.
We present detailed examples of the Zachariae Isstrøm in northeast Greenland, involving ICESat and ATM laser altimetry data, together with DEMs obtained the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and stereo aerial photogrammetry.ASTER DEMs are readily available but notorious for their accuracy behavior.The nominally stated accuracy of 15 m may occasionally reach much higher values.By embedding ASTER DEMs into the SERAC time series of elevation changes, we are able to determine plausible corrections.Thus, we can use ASTER DEMs to temporally and spatially densify the elevation change record.This is especially important on rapidly changing outlet glaciers where laser altimetry data might only sporadically be available.

OVERVIEW OF SERAC SYSTEM
Fig. 1 provides an overview of the SERAC processes to derive from laser altimetry data and DEMs complex information such as volume and mass balance change rates at individual glaciers, drainage basins, or the entire Greenland Ice Sheet (GrIS) (Schenk andCsatho, 2012, Schenk et al., 2014).
The figure is color coded; dark green boxes indicate SERAC external data sources, light green boxes contain SERAC generated data and the red boxes symbolize major SERAC processes.The chain of processes begins with obtaining laser altimetry data, such as ICESat data from NASA's satellite altimetry mission or ATM and LVIS data collected using NASA's airborne systems, from NSIDC (National Snow and Ice Data Center).After converting the NSIDC data format to the SERAC format, the first processing step determines the location of surface patches, involving as many repeat missions from ICESat, ATM and LVIS.A surface patch constitutes the areal unit used in change detection, size is about 1 km × 1 km.Fig. 2 shows a typical surface patch.
Next follows the computation of time series of surface elevation.This time series can be plotted for visualization purposes (Fig. 3).The horizontal axis denotes time, the vertical axis shows the elevation of a surface patch at any given time.Laser altimetry data essentially samples the visible surface of ice and snow.Thus, taking differences between the various time epochs gives elevation changes and because snow in its different conditions spans a The conversion from elevation to mass starts with removing the vertical crustal motion due to Glacial Isostatic Adjustment (GIA) and the elastic crustal response to recent ice loading changes.The vertical crustal motion estimates for GIA are from (A et al., 2013) provided in a 1 • × 1 • .The estimates range from 2.7 to 4.6 mm per year with errors negligible compared to the elevation change errors.The conversion from ice thickness to mass requires a firn densification model (FDM) to determine ice thickness changes related to ice dynamics, which can be converted to mass using the density of ice, typically 917 kg m −3 .The dynamic mass change is then combined with the surface mass balance (SMB) anomalies to obtain the total mass change.We refer to (Csatho et al., 2014) for details on the elevation change to mass change conversion.
For the results shown in this paper we used the firn densification model from (Kuipers Munneke et al., 2015) and RACMO2.3SMB estimates (Noël et al., 2015).
The surface elevation time series and the external data of FDM and GIA enter the module compute ice thickness changes related to ice dynamics.It produces time series of the total ice thickness, changes related to FDM, and determines changes related to ice dynamics (total change minus FDM change).Fig. 3 illustrates the results at a location near the grounding line of Zachariae Isstrøm.The upper panel shows a time series of total ice thickness; thickness change related to surficial processes (FDM) are plotted in the middle panel and those caused by ice dynamics are shown in the lower panel.These results, together with the external SMB time series enters the next stage, called compute annual rates of thickness and mass changes at time series locations, that is identical to the location of the original surface patches.We took the monthly SMB estimates from the climate model RACMO2/GR (Noël et al., 2015) that is given on 11 km × 11 km grid.This grid differs from the time series locations and the corresponding values have been found by nearest neighbor interpolation.
For convenience of computing, the final stage of SERAC takes place on regular grids.We first interpolate the annual rates computed at irregularly distributed time series locations to grids and then compute spatial integrals to get ice sheet volume and mass balance rates.The results are for the entire GrIS but they can also be analyzed on smaller regions, such as drainage basins or even individual outlet glaciers (Csatho et al., 2014).

METHOD FOR CORRECTING ASTER DEMS
This section presents a short synopsis of the method on how we correct ASTER DEMs that are introduced together with the more precise laser altimetry observations.For a more detailed description the reader is referred to (Schenk et al., 2014).
We use AST 14 because GDEM is the average of all useful ASTER DEMS at a specific location and therefore not useful for detecting surface elevation changes over time.The reported vertical accuracy of AST 14 DEMs ranges from 15 m to 60 m.It is well-known that for surfaces covered by ice or snow, these values may become even larger, because gray level correlation is unreliable if the image patches used for matching are nearly homogenous.Moreover, images with partial, undetected cloud coverage or ground fog may produce wrong elevations.In general, to render stereoscopic DEMs useful for the computation of time-series of elevation changes it is important to correct them by aligning them to the more accurate laser altimetry points.
A general method to correct stereoscopic DEMs is by way of a 3D similarity transformation.The seven unknown transformation parameters include a scale factor, 3 rotation angles, and 3 translations.Since the rotation angles and the scale factor are rather small quantities, we can simplify the 3D similarity transformation, using the differential rotation matrix, to obtain On the left side are the transformed (corrected) coordinates x , y , z , expressed as a function of the transformation parameters; ∆s, the scale; ∆ω, ∆ϕ, ∆κ, the small rotation angles about the x−, y− and z−axes of the coordinated system; and xT , yT , zT , the coordinates of the translation vector.The x, y, z are the original coordinates, that is, x and y are the grid post locations of the DEM, and z is its elevation.A traditional way to establish the seven unknown transformation parameters is using a set of known points x , y , z (Ground Control Points, GCPs).Such GCPs have known coordinates in both coordinate systems.However, there are no identical points between laser altimetry and stereoscopic DEMs and all time epochs involved are different, that is, the related observations are on different surfaces.The traditional way to use GCPs to calculate the transformation parameters does not work-a new approach is needed.
We begin by dividing the 3D transformation of Equation 1 into a horizontal and vertical transformation, assuming that the transformation parameters are small quantities and their products can be neglected.This assumption is certainly valid in our case where the transformation parameters are small corrections of DEMs.
The horizontal transformation is given as and contains the four parameters ∆s (scale), ∆κ (small rotation angle about the z−axis), and xT , yT (two horizontal shift parameters).For the vertical transformation we have where ∆ω, ∆ϕ are the two rotation angles about the x− and y−axes and zT is the vertical translation.
We apply first a 2D similarity transformation to align the ASTER DEMs horizontally, using control points and control features (distinct lines identified and measured).Using control features is very advantageous as it is much easier to identify them unambiguously compared to points (Schenk et al., 2005).Control features include ridge crests and sudden changes in slope.
After the initial horizontal alignment we now turn to the crucial question on how to introduce vertical control information to transform the stereoscopic DEMs to the laser altimetry points, bearing in mind that there are no identical points between the two sets.Moreover, the time epochs and their related surface elevations are different, too.
Fig. 4 shows the truncated time-series of elevation changes at site 5 (see Fig. 5 for site locations), for the time period 1997-2013.including laser altimetry measurements (ATM = blue dots, ICESat = red dots) and stereoscopic DEMs (ASTER = green dots).
Looking closely at the laser altimetry points we notice increasing thinning from 2001 to 2006, followed by more stable elevations until 2011.Surface elevation changes are more deterministic than random, suggesting to fit an analytical curve to the points.
Our detailed study of 130 Greenland outlet glaciers revealed that dynamic thinning events of outlet glaciers occur normally over a short time period (exceptions exist) and they can be well described by a sigmoid curve.Therefore we fit the following algebraic form of a sigmoid curve to the ATM and ICESat points.In this equation, t refers to time, h is the elevation relative to August 31, 2006 and a, b, c, d are the four parameters that have to be determined during the curve fitting process.
The upper panel of Fig. 4 shows the polynomial curve fitted through the more accurate laser altimetry points.The curve reveals a problem with the ASTER DEMs: they do not fit the curve.We consider the error between stereoscopic DEMs and the fitted polynomial as the vertical control information needed in Equation 4to calculate the 3 unknown parameters.The errors (one is labeled in Fig. 4 with residual) can easily be computed as the difference between the corresponding points on the curve and the DEM points.Examining Fig. 4 further reveals that there are seven ASTER stereoscopic DEMs involved.Each of these DEMs is assigned an error at the location of the time-series.In order to compute the three unknown correction parameters of every DEM we need to repeat the process at least at two more time-series locations, preferably many more to formulate an adjustment problem.
Suppose we have n time-series, n > 3.In this case, the three correction parameters per DEM can be found by a least-squares adjustment.Equation 4 is slightly rearranged to become a linear observation equation with ei the residual at location i, ∆ϕ, ∆ω the rotation angles, zT the vertical shift parameter, xi, yi the coordinates at location i and zi − z i the observed error.
Figure 5 illustrates the procedure.The original DEM (black) is transformed to the corrected position (red) by two small rotational angles, ∆ω, ∆ϕ, about the x and y-coordinate axes and the vertical shift, zT .The linear error model is a first approximation and future research has to address more complex, non-linear error models.
We summarize the procedure to correct stereoscopic DEMs for elevation errors, assuming that a horizontal registration has been accomplished first, for example by using Equations 2 and 3 on stationary points and lineal features.Lineal features offer the advantage of improved identification and localization but require a more sophisticated transformation algorithm than presented in Equations 2 and 3. Since we focus in this paper on correcting elevations, we do not elaborate further on the horizontal alignment.
1. Co-register all data sets involved to a common vertical reference frame, here WGS-84.
2. Select locations with repeat ICESat, ATM and LVIS laser altimetry data and stereoscopic DEMs to be corrected.
3. Compute time-series with SERAC at the selected locations and fit a sigmoid curve through the laser altimetry points.
4. Calculate the differences of DEM points to the fitted curve and enter them as observations into Equation 6 to compute the three DEM correction values (2 small angles ∆ω, ∆ϕ) and absolute elevation correction (zT ).
5. Repeat the time-series calculation with the corrected DEMs as a final check.
Finally, the panel of Fig. 4 presents the time-series of elevation changes after correcting the stereoscopic DEMs.The red line is the same polynomial curve as in the upper panel fitted through the laser altimetry points.We observe that the DEM points are now randomly distributed about the curve.However, their spread is noticeably larger than the red and pink laser altimetry points, reflecting their lower precision.

RESULTS
We present detailed examples of the Zachariae Isstrøm (ZI) in northeast Greenland, involving ICESat and ATM laser altimetry data, together with DEMs obtained the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and stereo aerial photogrammetry.Zachariae Isstrøm is one of the main outlets of the Northeast Greenland Ice Stream (NEGIS, Fig. 5).NEGIS is the only large ice stream of the GrIS, a linear body of fastmoving ice that penetrates 700 km into the interior of Greenland.Its drainage basin represents 12 % of the GrIS and contains a sealevel rise equivalent of 1.1 m (Mouginot et al., 2015).Zachariae Isstrøm is currently undergoing rapid changes in surface elevation and ice velocity, manifesting as ice front retreat, accelerating speed-up and thinning (Khan et al., 2014, Mouginot et al., 2015).
We fused the altimetry and stereo DEM record (Fig. 5) to examine the elevation changes of Zachariae Isstrøm between 1978-2013.
The fused elevation change data set allow the characterization of glacier elevation changes at different spatial and temporal scales.
The SERAC times series provided the first accurate measurements of the accelerating thinning (Csatho et al., 2014).The thinning is mostly due to ice dynamic processes, i.e., the acceleration of ice flow with a minimum contribution from increasing surface melt (Fig. 3).

CONCLUDING REMARKS
Zachariae Isstrøm, one of the three outlet glaciers draining ice from the North East Greenland Ice Stream (NEGIS), has rienced accelerating thinning and speed-up leading to increasing mass loss during the last decades (Khan et al., 2014, Mouginot et al., 2015).To investigate the mechanisms controlling ice dynamical changes, we reconstructed its elevation change history between 1978 and 2014.Our analysis shows accelerating thinning between 1974-2014 within 25 km distance from the 1996 grounding line, in a region where that glacier's bed is the deepest at 600 meter below sea level.This finding is consistent with acceleration and thinning initiated by ice loss near the grounding line, most probably caused by ocean warming (Mouginot et al.,  (2003, 2012, 2013. . in discrete steps from the calving front toward inland.Therefore, we suggest that bedrock topography or processes acting along the entire glacier, such as changes in lateral drag or ice viscosity, also influence Zachariae response to changing environmental conditions. Our analysis of ICESat time series shown that dynamic thinning extended as far as 100 km inland by 2007(Csatho et al., 2015)), indicating that an increasingly large region of the glacier is becoming vulnerable to further ocean warming and air temperature increase.We detected a pattern of "blockwise" thinning, i.e., separate large regions thinning with uniform rates (Fig. 9).This is markedly different from a simple inland propagating diffusion signal, possibly caused by rapid ice loss at the grounding line and observed for example on Kangerlussuaq Glacier in east Greenland (Schenk et al., 2014).We speculate that major bedrock topographic features act as barriers for inland propagation of thinning on Zachariae Isstrøm.A comparison of our detailed thickness changes history with numerical modeling results from the Ice Sheet System Model (Larour et al., 2012, Schlegel et al., 2014) suggests that while most of the observed rapid thinning is likely caused by grounding line retreat, a dynamic response of ice flow to surface mass balance variations might also play an important role (Csatho et al., 2015).

ACKNOWLEDGEMENT
Figure 3. series of total ice thickness (upper panel), ice thickness changes due to surficial processes (middle panel) and caused by ice dynamics (lower panel) near the grounding line of Zachariae Isstrøm.

Figure 4 .
Figure 4. Elevation time series for site 5 near the grounding line of Zachariae Isstrøm (circled in red in Fig. 5) before the DEM correction is applied (upper panel) and after the DEM correction (lower panel).Altimetry observations are from ICESat (red) and ATM (blue).Green dots mark elevations derived from ASTER DEMs.

Figure 5 .
Figure 5. Location map of data combined for elevation change reconstruction.Elevation data set includes ICESat (purple lines) and ATM (green lines) altimetry and ASTER DEMs (a total of 10; June 21, 2001 (red box); July 19, 2009 (orange box) and August 23, 2013 (blue box) are shown).Small yellow circles mark surface patches used for DEM correction.Elevation change time series for site 5 (large red circle) are presented in Fig. 4. Central flow line profile, shown Fig. 9, is marked by red line.Thick black line is 2013 calving front.Small map shows location of Zachariae Isstrøm (ZI) on ice velocity map from Rignot and Mouginot (2012).

Figure 6 .
Figure 6.Statistics and errors of DEMs used in this study.

Figure 7 .Figure 8 .
Figure 7. Statistics and errors for the 18 SERAC times series used as control information for the fusion and the correction process.Locations are shown by small yellow circles in Fig. 5.