Open AccessArticle

Using Deep Learning to Model Elevation Differences between Radar and Laser Altimetry

Alex Horton

^1,*,

Martin Ewart

¹,

Noel Gourmelen

Xavier Fettweis

and

Amos Storkey

⁴

Earthwave Ltd., Edinburgh EH3 9DR, UK

School of Geosciences, University of Edinburgh, Edinburgh EH8 9YL, UK

SPHERES Research Unit, Department of Geography, University of Liège, 4000 Liège, Belgium

⁴

School of Informatics, University of Edinburgh, Edinburgh EH8 9YL, UK

Author to whom correspondence should be addressed.

Remote Sens. 2022, 14(24), 6210; https://doi.org/10.3390/rs14246210

Submission received: 13 September 2022 / Revised: 14 November 2022 / Accepted: 17 November 2022 / Published: 8 December 2022

(This article belongs to the Special Issue Multi-Source Data with Remote Sensing Techniques)

Download

Browse Figures

Figure 1
(a) The geographical distribution of the ~29 million collocated SARIn vs. LIDAR observations in this study and (b) the statistical distribution of the same observations with a fitted normal inverse Gaussian distribution trend line. (c.f. <a href="#sec3dot2-remotesensing-14-06210" class="html-sec">Section 3.2</a> for details on the standard error calculation). "> Figure 2
Summary of the observed and modelled LIDAR-SARIn differences as a function of a selection of input parameters for all 3.6 million unseen test data points across GrIS. The distribution of observation counts per parameter is shown in the grey histogram for each plot. "> Figure 3
Gridded maps of observed LIDAR-SARIn differences for February 2020 (a) and August 2020 (d) and their difference (g), and equivalent model predictions from the ordinary least squares model (b,e,h), and from the neural network model (c,f,i). Note that the scale bars have been adjusted between observations and model predictions to allow for a clearer visual comparison of broad trends. "> Figure 4
LIDAR-SARIn adjustment as a function of distance to POCA broken down by 6 quantile bands of coherence, power (dB), and slope across. "> Figure 5
Spatial plots of MAR model outputs across GrIS for (a) cumulative snowfall in February 2020 and the prior five months; (b) cumulative snowfall in August 2020 and the prior five months; (c) the change in accumulative snowfall between August 2020 and February 2020: i.e., the difference between (b) and (a); (d) the approximate average snow density in the top 2 m of the snowpack for February 2020; (e) the average surface mass balance over 12 months prior to and including February 2020; (f) the total summer melt for June, July and August 2020. "> Figure 5 Cont.
Spatial plots of MAR model outputs across GrIS for (a) cumulative snowfall in February 2020 and the prior five months; (b) cumulative snowfall in August 2020 and the prior five months; (c) the change in accumulative snowfall between August 2020 and February 2020: i.e., the difference between (b) and (a); (d) the approximate average snow density in the top 2 m of the snowpack for February 2020; (e) the average surface mass balance over 12 months prior to and including February 2020; (f) the total summer melt for June, July and August 2020. ">

Review Reports Versions Notes

Abstract

Satellite and airborne observations of surface elevation are critical in understanding climatic and glaciological processes and quantifying their impact on changes in ice masses and sea level contribution. With the growing number of dedicated airborne campaigns and experimental and operational satellite missions, the science community has access to unprecedented and ever-increasing data. Combining elevation datasets allows potentially greater spatial-temporal coverage and improved accuracy; however, combining data from different sensor types and acquisition modes is difficult by differences in intrinsic sensor properties and processing methods. This study focuses on the combination of elevation measurements derived from ICESat-2 and Operation IceBridge LIDAR instruments and from CryoSat-2’s novel interferometric radar altimeter over Greenland. We develop a deep neural network based on sub-waveform information from CryoSat-2, elevation differences between radar and LIDAR, and additional inputs representing local geophysical information. A time series of maps are created showing observed LIDAR-radar differences and neural network model predictions. Mean LIDAR vs. interferometric radar adjustments and the broad spatial and temporal trends thereof are recreated by the neural network. The neural network also predicts radar-LIDAR differences with respect to waveform parameters better than a simple linear model; however, point level adjustments and the magnitudes of the spatial and temporal trends are underestimated.

Keywords:

SARIn; interferometry; CryoSat; swath; ICESat-2; IceBridge; artificial intelligence (AI); Greenland; cryosphere; altimetry

1. Introduction

Over the past 30 years, altimetry has revolutionised our ability to monitor changes across the world’s cryosphere. Radar altimetry is now routinely used to quantify changes in the world’s ice masses and their impact on sea level and water availability and to provide key observations to better understand climatic and glaciological processes behind those changes [1,2,3,4,5,6,7,8,9,10,11,12]. With two high-resolution, satellite-based altimeters currently active—the interferometric radar altimeter on board CryoSat-2 (CS2) and the laser altimeter on ICESat-2 (IS2)—the present period offers a unique opportunity to co-exploit the observations made by the two sensors and improve the monitoring of ice height and trends, as well as reveal insights into the properties of the firn layer and hence surface conditions. The ever-evolving roadmap for future generations of satellites to be made operational, such as the Copernicus Polar Ice and Snow Topography Altimeter, CRISTAL, will further compound the need to combine measurements from different sensor types meaningfully.

Recent advances in swath altimetry processing, using the full waveform of the interferometric synthetic aperture radar (SARIn) mode of CS2, have enabled improved spatial resolution of surface elevation [13]. Meanwhile, IS2 provides enhanced resolution compared to the previous generation of laser altimeters thanks to its six laser beams and sharper footprint. Due to the fact that radar and laser altimeters have different intrinsic properties and behaviours, joining and interpreting their combined measurements requires careful consideration of factors such as differences in electromagnetic interaction with the surface, footprint size, weather conditions, and topography. Several studies have focused, in particular, on the variability in the scattering depth of radar altimeters and on ways to mitigate its impact [14,15,16] and exploit it to retrieve geophysical information about subsurface and climatic conditions [17,18]. However, equivalent studies have not been conducted for SARIn altimetry to enable the calculation of adjustments at the sub-waveform level for large acquisition regions such as the Greenland Ice Sheet (GrIS). Such an advancement would enable the meaningful combination of Ku-band SARIn and LIDAR altimetry allowing higher-resolution reference digital elevation models (DEMs) and more accurate changes in elevation over time. Such a method would also be directly applicable to future Ku-band SARIn altimeters, as proposed for the upcoming CRISTAL mission.

Slater et al.’s study into compensating changes in the penetration depth of pulse-limited non-interferometric radar altimetry [16] demonstrates the expected depth distribution of Ku-band radar penetrations over the simple topographies of the central GrIS. A key finding was the change in signal penetration depth of over 2 m on either side of the unique pan-GrIS melt event in mid-2012, which highlights the radar measurement’s dependency on the firn state. The numerical solution in the study focuses on the best fit parameterisation of a complex analytical waveform function [19]. Still, it does not easily extend to sub-waveform level SARIn adjustments with the addition of coherence, nor does it take additional factors such as complex topography or firn layers into account, such as are more prevalent around the GrIS periphery. Gray et al. studied the complications in using SARIn measurements for monitoring time series elevation changes across GrIS due to changing surface conditions [20] with a suggestion for using data from late summer and fall to minimise the possibility of varying bias between real surface height change and that reported by CS2’s SARIn altimeter. Additionally, non-trivial sub-waveform corrections are known to be unaccounted for in current CS2 SARIn products from the off-nadir ‘mis-pointing’ of the altimeter [21]. In summary, a complete solution for the combination of LIDAR and SARIn altimetry does not currently exist in numerical or analytical form, and if it did, such a solution would require a multivariate, higher order, numerical, or semi-analytical construction.

While neural networks (NNs) have been applied to a multitude of scientific disciplines, their use in the cryosphere-related domain is comparatively new. Successful applications range from improved snow water equivalent products over British Columbia, Canada [22] to the identification of potential Antarctic meteorite sites [23]. Further studies of particular relevance are the regression-based neural network solution for the estimation of snow depth on Arctic sea ice using multi-band observations from the Copernicus Imaging Microwave Radiometer [24] and the enhanced waveform retracking of lake surface elevation using deep learning [25].

Here we use a deep NN model to combine elevation measurements acquired by CS2 with those from Operation IceBridge (OIB) and IS2 between 2011 and 2021, utilising SARIn waveform parameters with additional local topography products over the GrIS. We explore the difference between interferometric radar and laser altimetry and its relationship with the surface condition, the impact of penetration of radar waves into the snow and firn, and the respective measurement uncertainties. The modelled elevation corrections are used to generate time-dependent DEMs for the required LIDAR-SARIn adjustment.

2. Data

Point data were collected for CryoSat-2 using its Level 2 (L2) SARIn product, Cryo-TEMPO-EOLIS [13], acquired between 2010 and 2021; the study focused on the full swath data provided by the SARIn mode in order to utilise the higher resolution offered when compared to the traditional point of closest approach (POCA) product [13,26,27].

Operation IceBridge ATM Level-2 Icessn Elevation, Slope, and Roughness (ILATM2) dataset [28] acquired using the NASA Airborne Topographic Mapper (ATM) instrumentation was collected between 2010 and 2016 from which the WGS84_Ellipsoid_Height measurement was taken for elevation measurements.

The ICESat-2 ATLAS L3A Land Ice Height, Version 2, ATL06 was gathered from the satellite’s launch in October 2018 up to February 2021 [29]. From this dataset, elevation measurements and expected RMS segment misfit were extracted (c.f. h_li and h_li_sigma, respectively, from the product handbook).

Lastly, the ArcticDEM digital surface model dataset [30] was gathered.

3. Method

3.1. Data Filtering and Preparation

3.1.1. CryoSat-2 Data

To ensure that data of sufficient quality is used, a filter is applied that consists of two physical parameters; the power in decibels must be greater than −160 dB, and the interferometric coherence must be greater than 0.8. In addition, a comparison of the elevation to the Arctic DEM with the absolute difference (DEM-diff) being less than 100 m and the median absolute deviation of the difference (DEM-diff-MAD) across the whole waveform needs to be no more than 10 m. The DEM-diff filters are used to identify errors in the derived phase ambiguity solutions and remaining outliers after power and coherence screening, consistent with Gourmelen et al. [13].

The following parameters were selected as inputs to the deep NN model because of their potential sensitivity to scattering processes: Coherence which is expected to generally decrease with increasing volume scattering; power with similar properties with respect to scattering and which has been used in previous scattering depth analysis [16] and also yields information to the position of the sample in the waveform; the horizontal scalar distance of each point measurement to its corresponding POCA measurement was calculated (distance to POCA) which was selected to provide information about incidence angle as well as waveform sample dependent characteristics such as recorded power; the waveform leading edge width which is expected to be larger with increased volume scattering [15,31]. The following two derived waveform quantities were also introduced to the model inputs to provide information about localised gradient and statistical uncertainties around a neighbouring sample point: relative-elevationⁱ and relative-elevation-mean^j. For a given waveform sample number n, the relative-elevationⁱ is defined as elevationⁿ⁺ⁱ–elevationⁿ where i = {−3, −2, −1, 1, 2, 3}. For the same sampled waveform point, n, the relative-elevation-mean^j is the arithmetic mean of relative-elevationⁱ selecting the neighbouring j points where −

\frac{j}{2}

≤ i <

\frac{j}{2}

and j = {10, 20}. The direct elevation measurement was omitted from model inputs to avoid spatial overfitting with the aim of improved spatial transferability.

3.1.2. ICESat-2 Data

The IS2 data were filtered to exclude expected RMS segment misfits greater than 0.03. Acquisition data from all six beams were used in the study, but the ATLAS spot (beam) number was not retained or considered for specific analysis.

3.1.3. ArcticDEM Data

ArcticDEM was used to derive the slope parallel to the sensor heading (slope along), slope perpendicular to the sensor heading (slope across), and aspect represented as a south/north vector component between −1 and 1 respectively, and east/west aspect vector component between −1 and 1 respectively for each CS2 point observation. The slope along and slope across were selected as model inputs to provide a reliable reference estimate of localised topography. They were incorporated into the model to account for physical effects from electromagnetic scattering as a function of incident angle as well as providing information to help correct for additional, non-trivial relationships observed between SARIn and LIDAR measurements for specific topographies. These values were derived using a cubic bivariate spline approximation (c.f. scipy.interpolate.RectBivariateSpline [32]) given the coordinates for each CS2 observation.

A summary of all model inputs is given in Table 1.

3.1.4. Spatio-Temporal Join

Radar and laser measurements used in the analysis are selected based on a spatio-temporal join. For each CS2 elevation measurement, the nearest neighbour is found in both the OIB and IS2 datasets such that the joined dataset retains point data that is coincident within both a 10-day time window and a 50-metre radius. These criteria were selected to maximise dataset volume while minimising the impact of changes in ground conditions and local topography in space and time. To minimise the impact of local gradients, within the 50 m join distance, on the difference between almost-coincident SARIn and LIDAR elevation measurements, the ArcticDEM was used again to derive an expected elevation difference due to slope using a cubic bivariate spline approximation [32]. This slope correction was then removed from the observed elevation difference between the radar and LIDAR sensors.

After filtering and joining, the full observation dataset was comprised of 1.7 million points from OIB and 27 million from IS2. The dominant contribution of points from IS2 explains the higher density of observations seen in the north of the GrIS in Figure 1a. When looking at a histogram of the observed differences between SARIn and LIDAR elevation measurements, per Figure 1b, a wide-tailed, non-symmetric distribution is observed for which a normal inverse gaussian distribution was found to fit well.

3.2. Statistical Assumptions

It has been shown that there exists a non-zero spatial correlation in the underlying dataset that drops to zero after 3.2 km (c.f. CryoTEMPO-EOLIS ATDB [33]). A conservative approach was taken here, assuming an autocorrelation of 100% below 3.2 km. With this assumption, and for the purpose of statistical analysis only, the data were clustered using 100% correlation up to 3.2 km. Then, an uncorrelated standard error was calculated across the clustered data. This is consistent with McNabb et al. [33]; however, in this case,

r = L = 3.2 km

, simplifying the calculation to:

n_{e f f} = \frac{n_{t o t} r^{2}}{L^{2}}

where

n_{t o t}

is the total number of clusters. This approach gives a conservative estimate of the standard error given that, in practice, the data is not 100% correlated below 3.2 km.

3.3. CryoSat-2 SARIn vs. ICESat-2 LIDAR Maps

A gridded monthly raster product of differences between CS2 and IS2 was generated between October 2018 and February 2021. For each month, centred around the 15th of each month, 90 days of data were used to create a 2 km resolution product for CryoSat-2 and ICESat-2, respectively. The difference was then taken, excluding pixels where no data were present for either dataset. A median filter was subsequently applied using a square kernel of side length 5 pixels (10 km) which excludes no-data values when considering the median calculation.

3.4. Models

A series of neural networks (NNs) were created and compared to a simple, multivariate, ordinary least squares regression model (OLS). The OLS model was selected as a simple benchmark similar to linear models previously employed in investigating the impact of backscatter on fluctuations in the range [11]. To explore the impact of neural network design on LIDAR-SARIn prediction effectiveness, a neural network framework was constructed using pytorch [34] that allowed the flexible creation of a given neural network by specifying: the number of inputs, number of hidden layers, number of nodes per hidden layer as a fraction of number of inputs, and the activation function applied for hidden layers (None, Relu [35], LeakyRelu [35], Selu [36]), the dropout ratio (0->1) [37], and the hidden layer normalisation method (None, BatchNorm [38], LayerNorm [39]). There was only 1 output node for the normalised, regression-based prediction of the LIDAR-SARIn adjustment. As well as dynamic model creation, the SGD [40], Adam, and AdaMax [41] optimisers were investigated alongside the choice of L1, MSE, and Huber [42] loss functions. Each training epoch was subdivided and optimised via mini-batch gradient descent. Input and output data were normalised using a standard scalar approach [43]. Performance across each of these degrees of freedom in the model configuration space was explored using a standard grid search approach.

The data were separated into training and test datasets to focus on temporal transferability across the GrIS: OIB data prior to 2016, and IS2 data prior to September 2020, were used as training data. OIB data for 2016 and IS2 data from September 2020 to December 2021 were used as unseen test data for performance analysis. For the neural network training, 10% of the training data were selected randomly to be used as validation data and were thus excluded from the loss function optimisation. The optimised model weightings were taken from the model with the lowest loss against the validation data to avoid overfitting.

The output from each model configuration is a collection of predictions for the elevation adjustments that would need to be added to a LIDAR measurement in order to retrieve the CS2 equivalent. These adjustments were compared to the observations in the training dataset and the unseen test dataset. The following simple statistical assessments were calculated to assess model accuracy: mean, abs(mean), median, RMSE, median absolute deviation and Pearson’s correlation coefficient (r). The final model selection was focused on minimising the RMSE, mean, and median while maximising the r.

The non-exhaustive grid search approach highlighted ranges of model configuration and hyperparameters that produced similar levels of model performance. The model consequently selected and discussed here was constructed as follows: 16 input nodes, 6 hidden layers, with each layer reducing in node count by 85% and rounded up to the nearest integer, 1 output node; a mini-batch size of 32,000; the AdaMax optimiser; a smooth L1 loss function; and dropout was not used. This configuration will be referred to simply as the neural network or NN hereafter.

The OLS and NN models were used to predict elevation adjustments using all CryoSat-2 from 2010 to 2020 over the GrIS. The input CS2 SARIn data had the same quality criteria applied per the training data. A prediction is calculated for each individual SARIn point elevation measurement and then converted into a monthly raster product at a 2 km resolution. Each pixel in the rasterization is generated using a weighted average where the weights are defined by the inverse of the square distance of each point’s location to the centre of each pixel, as used by the CryoTEMPO-EOLIS v1 product [44]. The monthly raster products were used to build maps of predicted elevation adjustments and changes in predicted adjustments between the winter and summer months, represented by February and August, respectively.

4. Results

The mean difference of −1.01 m and standard deviation of 3.45 m (c.f. Figure 1b) between the radar and laser height measurement across the entire SARIn mode sector of the Greenland Ice Sheet, for the coincident data included in this study, is in line with previous more localised validation of swath data in Greenland [13,45]. This reflects the tendency of the Ku-band radar to scatter further within the snowpack than LIDAR, stochastic errors in the radar and LIDAR measurements, errors in the reference DEM used to slope-correct the radar-LIDAR join, and real height differences due to the temporal and spatial mismatch of the radar-LIDAR join.

Per Table 2, the unseen test dataset has a mean LIDAR-SARIn difference of −1.05 m. The neural network (NN) predicted this adjustment with a mean error of −0.048 m; the Ordinary least squares (OLS) method predicted this adjustment with a mean error of −0.175 m. To explore the statistical robustness of the improvement in the NN vs. the OLS method, conservative standard errors were taken (c.f. Section 3.2). Additionally, a two-sample t-test was conducted on the 3.2 km clustered data giving a near-zero p-value of 5.57 × 10⁻¹³—i.e., significantly less than 0.05. Lastly, a 34% improvement in the Pearson correlation coefficient (r) demonstrates the NN tends to better predict the observations.

Complex, non-linear trends are seen when comparing the LIDAR elevation measurements to sub-waveform SARIn point elevation measurements as a function of each input in the study per Figure 2. This figure focuses on the 3.6 million test dataset points to demonstrate model performance against unseen data. While power and coherence display a generally positive correlation with height differences, some behaviours are more complex. For example, the observed LIDAR-SARIn adjustment as a function of distance to POCA demonstrates multiple local minima and maxima and a general trend in decreasing magnitude with increasing distance to POCA.

When assessing the median LIDAR-SARIn adjustment prediction as a function of key quantities, it is clear that the NN model outperforms the OLS at replicating complex, non-linear trends per Figure 2. Model predictions vs. power and distance to POCA are good examples in this regard. However, it must be noted that the NN method does not outperform the OLS method for all quantities assessed or all ranges of a given quantity. For example, the OLS yields a close median agreement to the observations for lower coherence and more negative slope-across ranges as well as at a range of lower elevations. Such behaviours typically occur when the training and test datasets have smaller sample sizes.

When the observation data is aggregated into a 2 km raster product, LIDAR-SARIn adjustments of up to circa −4 m are seen per Figure 3a. This figure also shows the correlation between adjustment magnitude and the south/north aspect, as seen in Figure 2. However, localised regions exhibiting a larger adjustment are also seen for example, the western Rink Isbræ glacier and, to a larger extent, the northeastern region of King Frederick VIII Land. Figure 3d shows that LIDAR-SARIn adjustments typically decreased in magnitude in the summer of 2020, with the exception of the most northerly region, localised grounding line zones, and patches of the eastern GrIS. This is seen more clearly in subplot Figure 3g, which contains the difference in the observed adjustment between August 2020 and February 2020.

Figure 3b,e,h show the equivalent OLS model predictions, whereas plots Figure 3c,f,i show the equivalent NN model predictions. It can be seen that broad spatial trends are reproduced by both OLS and NN models. Regions of distinctly higher observed adjustments, such as Rink Isbræ and King Frederick VIII Land, are more accurately predicted by the NN model.

Both models underpredict the LIDAR-SARIn adjustment by up to ~3.5 m, for example, in the southerly aspects in the northwest of the GrIS and overpredict the adjustment for the most northerly facing aspects by up to ~1 m. When comparing observations to the model predictions of the change in LIDAR-SARIn adjustments over 6 months—as shown in plots Figure 3g–i—the OLS and NN models show a broad directional agreement with the observations. Still, the magnitude of this change is again underpredicted in the models by up to a factor of ~2.5 in the case of the NN model and over 3 for the OLS model.

5. Discussion

5.1. Observations

This first assessment of differences between the swath processing of CS2 and LIDAR onboard OIB and ICESat-2 at the ice sheet scale highlights the impact of both changing surface conditions and waveform parameters on height differences. An approach for predicting and adjusting for such LIDAR-SARIn altimetry differences is therefore critical in combining these datasets to utilise both spatial and temporal coverage fully. It also highlights the potential of such analysis for mapping and quantifying processes related to accumulation and melt at the ice sheet’s surface. The CS2 altimeter in SARIn mode uses 13.575 Ghz Ku-band microwave pulses close to the nadir and records backscattered power, coherence, and phase as a function of time. A single Ku-band pulse can be reflected from the topmost layer of ice, firn, or snow in surface scattering or penetrate several metres beyond the snow surface in volume scattering [46,47,48]. Remy et al.’s review of radar altimetry [14] also provides a good overview of related literature detailing surface backscattering as a function of incidence angle, snow temperature, grain size, roughness and density. However, while simplified electromagnetic penetration into snowpack has been well studied, a comprehensive understanding of the nature of backscattering in a non-homogenous snowpack across a complex elevation topography in the context of interferometric synthetic aperture radar altimetry and full waveform swath processing does not exist.

In a prior study [16], the scattering horizon for CS2, in non-interferometric low-resolution mode, was predicted to be between 2.34 ± 1.41 m and 3.28 ± 1.13 m beneath the surface on average in sectors of the ice sheet above 2000 m in altitude. Due to lower altitude and more frequent snowfall and melt in the GrIS periphery, and the introduction of interferometric processing, values on either side of this range were expected in this study and have been observed per Figure 1 and Figure 2. Conversely, the shorter wavelengths used by the LIDAR acquisition result in significantly less penetration into the snowpack with a measurement accuracy of ±10 cm [49,50]. A residual inter-beam range bias has been documented between ICESat-2’s six beams of circa from −2.8 cm ± 8.9 cm to +2.7 cm ± 8.1 cm [51,52]. However, given that these differences are more than an order of magnitude less than the observed LIDAR-SARIn differences, the LIDAR-specific contributions are assumed to be insignificant in the context of the analysis in this study.

5.1.1. LIDAR-SARIn Differences-Observed Contribution Drivers

A given sub-waveform, point-level, LIDAR-SARIn adjustment could theoretically be decomposed into a series of contributing components—we analyse the observed adjustments (c.f. Figure 1 and Figure 2) with respect to the input quantities and three high-level contributors: physical processes, systematic processes, and stochastic noise. In the following sections, the terms increasing and decreasing are assumed to refer to the magnitude of the observed and modelled adjustments unless otherwise stated. That is, as the typical median adjustment in the results plots is negative, an increase is considered to imply more negative, and a decrease is to imply less negative. The terms adjustment and observation imply the observed LIDAR-SARIn adjustment.

Figure 2 highlights several high-level trends tied to physical processes and interpretation. The adjustment is seen to typically increase with power. The increase in adjustment implies SARIn signal penetration, and thus volume scattering, higher power absorption, and lower power return. Similarly, volume scattering is expected to result in a decrease in coherence which agrees with our observations. However, for coherences below ~0.925, this trend weakens. This may be related to the coherence drop seen for higher sample numbers and, thus, increasing the distance to POCAs, which are seen to lead to a decrease in adjustment. This distance to the POCA trend, with adjustment increasing with distance to POCA while exhibiting multiple maxima and minima, is counterintuitive. A naïve argument would say that power and coherence decrease as sample number and thus distance to POCA increase. Therefore, an increase in adjustment should be seen for increases in distance to POCA, which is not the case. It is feasible that the complex interactions of the adjustment as a function of slope across, i.e., incidence angle, contribute more significantly than that of power and coherence to the overall behaviour of adjustment as a function of distance to POCA. Figure 4 demonstrates the sensitivity of the adjustment to different ranges of slope-across values and suggests that this might explain the local minima/maxima seen in the overall trend in Figure 2. The trend of adjustment as a function of slope across is not well understood but can be related to the impact of the cross-track angle on the accuracy of the swath processing [26] or to the fact that CryoSat is slightly ‘mispointing’ [21]. The correlation of adjustment with respect to the south/north aspect component can be attributed to increased accumulation on the southerly aspects. At the same time, no significance is drawn from the lack of correlation in the east/west component.

The relative elevation quantities were added for three reasons: to allow the models to capture information about localised slope across gradients according to the CS2 sensor, to grant information to the model about how this may differ from the reference DEM slope across, and to grant the models information about the stochastic nature for a localised section of the waveform. To aid the interpretation of these quantities, consider seven sequential waveform samples over a flat terrain where all measure the same elevation except the central point, which is 10 m lower: RelativeElevation¹ would be 10 m in this case. This would require a LIDAR-SARIn correction of −10 m to be accounted for as well as additional physical and systematic contributors. If relative-elevationⁱ captured only stochastic noise, one would expect the gradient of adjustment vs. relative elevation to be −1. A gradient of circa −0.55 is observed for relative elevations beneath its modal value and circa −0.90 otherwise. One interpretation is that relative elevations above an expected amount are more likely to be the result of stochastic drivers. In contrast, those less than the expected amount are made up less of stochastic noise and more of other physical or systematic contributions. The behaviour of adjustment as a function of relative-elevation⁻¹ is an x-shifted equivalent of that of relative-elevation¹, which is to be expected on an aggregated basis. The relative mean elevation quantities can also be interpreted in a similar way and exhibit similar behaviour. However, the pre-modal gradients for the adjustment as a function of relative-mean-elevation¹⁰ and relative-mean-elevation²⁰ are circa −0.63 and −0.70, respectively, suggesting they are even stronger stochastic correction indicators.

Figure 1b shows that the study data has a wide-tailed, asymmetric distribution of observations. In total, 34% of the observations relate to a LIDAR-SARIn adjustment greater than zero, implying either significant stochastic and/or unaccounted-for systematic processing contributions in the data. To better understand at least part of these systematic processing contributions, further study into LIDAR-SARIn differences over a suitable solid land body, i.e., land not covered in snow/ice/firn, could be of interest. To best mimic an ice sheet, the land body would need to be carefully selected: for example, a large region of non-flat but low roughness topography would be worth considering.

5.1.2. LIDAR-SARIn Differences–Spatial and Temporal Trends

Once aggregated into a spatiotemporally gridded product, per Figure 3a,d, the sub-waveform level and stochastic contributors are largely eliminated. Therefore, it is assumed that the aggregate signal reflects the LIDAR-SARIn adjustments related to physical conditions. The possible introduction of a net bias into this aggregate signal from systematic and/or stochastic components is not explored further here.

It has previously been observed that CS2 SARIn measurements can require adjustments of up 2 m in the southeastern region of GrIS in relation to the accumulation and snowpack conditions. In contrast, the northerly aspects and the southwestern region GrIS typically see substantially less accumulation and expect less CS2 SARIn penetration [20]. This study is in general agreement. However, our data suggest adjustments up to circa 4 m due to the use of the full waveform records (c.f. Figure 3a). Typical seasonal accumulation and melt across GrIS have also been reported previously [18] and are shown again in this study.

A selection of outputs from the regional climate MAR model forced by the ERA5 reanalysis for GrIS [53] is shown in Figure 5 as indicators of the snowpack state. Plots (a), (b), and (c) in Figure 5 offer a crude estimate of the snowpack state for February 2020 and August 2020 and snowpack change over 6 months. Broad qualitative spatial correlations are observed between accumulated snowfall, observed LIDAR-SARIn differences, and the OLS and NN predictions for LIDAR-SARIn differences in Figure 3. The change over 6 months (Figure 5c vs. Figure 3g) is particularly interesting for qualitative spatial agreement. One notable exception to this is the King Frederick VIII Land region, where a simple historic snowfall estimate of the snowpack state does not align well with the observations nor the OLS and NN predictions. Still, a localised summer melt (Figure 5f), leading to a dense firn layer, is a possible explanation for this. Similar qualitative spatial correlations are drawn between the snowpack density (Figure 5d) and observed LIDAR-SARIn differences (Figure 3a).

5.2. OLS and NN Models-Strengths

In light of the complex interdependencies described in Section 5.1, both OLS and NN models perform well at recreating the median observation trends per input parameter (c.f. Figure 2) and reproducing mean trends (c.f. Table 2). The ability of the NN to better reproduce non-linear trends—e.g., distance to POCA and power—is expected due to the nature of the model construction, with non-linearity enabled specifically by the multiple node layers and choice of the activation function.

Per Figure 2, the NN model outperformed the OLS model in predicting the observed spatial and temporal trends while acknowledging that neither fully predicted the full range of aggregate observations.

Strong temporal transferability is implied by the low mean and median prediction errors for unseen data, as shown in Table 2. Caution is needed with spatial transference to other regions as the complex response to inputs is expected to differ substantially between Greenland and, say, the Antarctic Ice Sheet. However, the prevalence of IS2 data would make retraining the model in other regions possible, where similar performance to GrIS would then be expected.

5.3. OLS and NN Models-Limitations

While both models reproduced broad spatial and temporal trends, the lack of ability to capture the correct aggregated magnitude of required LIDAR-SARIn adjustment implies that such a model alone is insufficient for the meaningful combination or comparison of CS2 SARIn and IS2 LIDAR datasets. Without detailed datasets for the snowpack states, including depth profiles, and a comprehensive understanding of all physical processes and sensor behaviours, analysing the reasons for the models’ shortcomings is an exercise in conjecture. There are three questions to consider: is the model lacking sufficient information, is the signal-to-noise ratio too low, and/or does the model lack complexity?

It is possible that the addition of extra datasets (including outputs of polar regional climate models) would improve a model’s ability to predict penetration. Promising candidates would include datasets describing snowpack states at different depths with high spatial and temporal resolution. One aim of this study was to find a way to meaningfully combine LIDAR and Ku-band SARIn data in such a way that was easily reproduceable and could be used to generate Level 3 (L3) CS2 products on an ongoing basis. Introducing input dependencies on datasets that do not exist for the lifespan of CS2 limits their utility in L3 product generation.

It is seen in Figure 1b that a larger than expected distribution of the LIDAR-SARIn adjustments is observed across our dataset. While it is unclear exactly how much of this distribution can be attributed to systematic, stochastic, and physical contributions, it is apparent from Section 5.1.1 that the systematic and stochastic contributions outweigh the physical ones. Given that both the OLS and NN models simply optimise their parameter space such that a loss function (least squares and Huber loss, respectively) is minimised, it can be expected that signal reproduction quality will be lost when the data has a poor signal to noise ratio as is the case in this study. For example, it would be interesting to repeat this study once a corrected CS2 L2 product is available with the changes suggested by Recchia et al. [21]. Additionally, it might be interesting to resample the dataset so that we get a fairer representation of each region and aspect instead of the north-heavy distribution, as shown in Figure 1a.

With the existing selection of inputs, we believe this study exhausted the investigation into the use of simple but deep neural networks to model the observed data, as introducing additional layers of nodes of numbers of nodes per layer did not result in an improved mean or RMSE against unseen datasets. However, given that there is a large amount of information in the entire waveform representation of power, coherence, and elevation (or elevation relative to the initial sample), the use of an autoencoder-based convolutional neural network could be employed to either perform a dimensionality reduction of each waveform parameter for the use as extra input in the standard neural network model proposed here or to directly predict a full waveform level elevation adjustment profile. Another interesting possibility is to use a waveform-level autoencoder approach with only the power to reproduce the study of slater et al. [16] covering both the low-resolution mode interior and SARIn mode periphery of GrIS.

6. Conclusions

A comparison of sub-waveform CryoSat-2 SARIn altimetry measurements to Operation IceBridge and ICESat-2 LIDAR altimetry was conducted to reproduce a required adjustment that facilitates the combination of LIDAR and SARIn datasets. Aggregation of the point level LIDAR vs. SARIn observations into a 2 km filtered product reduced the influence of stochastic and systematic effects and highlighted regions of radar signal penetration consistent with snowpack condition expectations. Comparing such products between February 2020 and August 2020 indicates changes in surface conditions consistent with accumulation, snowpack density, and melt from the regional climate MAR model forced by the ERA5 reanalysis.

A deep neural network (NN) was built to model this complex, multivariate problem. For unseen test data, the mean observed LIDAR-SARIn adjustment of −1.05 m was predicted by the NN model to be within 0.048 m, with a standard error of 0.019 m and a Pearson’s Correlation Coefficient of 0.5. The NN performance was compared by an ordinary least squares model (OLS) with the same inputs, which predicted the mean observed adjustment to within 0.175 m with a standard error of 0.020 m and a Pearson’s correlation coefficient of 0.37. Compared to the OLS, the NN more accurately recreated the complex trends of LIDAR-SARIn adjustment with respect to changes in the key input parameters.

The NN model reproduced broad, aggregated trends in spatial and temporal patterns seen in the observations. However, further investigation and enhancements would be required to capture the full magnitudes of adjustments between SARIn and LIDAR acquisitions at the point or raster level. It is postulated that one or more reasons drive the reason for poorer predictions of magnitudes in the NN: insufficient input datasets to capture the physical processes, insufficient signal-to-noise ratio, or insignificant model complexity. It is specifically recommended that these models be retrained and retested once a new CS2 product has been regenerated that incorporates the waveform level processing corrections detailed by Recchia et al. [21]. Another promising avenue for investigation is using the full information of the power, coherence, and/or elevation at the sub-waveform level in conjunction with a dimensionality reduction technique such as a neural network autoencoder.

This novel level of analysis of sub-waveform SARIn vs. LIDAR elevation measurements will be instrumental in driving and/or validating future theoretical advances into the precise nature of the differences between sensor types. Once the necessary level of accuracy is achieved in modelling and predicting the differences between measurements from different sensor types at the point level, there are many applications beyond CS2 SARIn and LIDAR. An obvious candidate is the future CRISTAL altimetry mission [54]: CRISTAL will carry a Ka/Ku radar and generate a swath and POCA elevation product. The research presented here could be used for the fusion of the swath and POCA measurement and to extract information about penetration and surface conditions from the dual frequency instrument.

Author Contributions

Conceptualization, M.E. and N.G.; methodology, A.H., N.G., M.E. and A.S.; software, A.H. and M.E.; validation, A.H., M.E. and N.G.; formal analysis, A.H., M.E. and N.G.; investigation, A.H., M.E., N.G. and A.S.; data curation, N.G. and A.H.; writing—original draft preparation, A.H.; writing—review and editing, A.H., M.E., N.G. and X.F.; visualization, A.H.; supervision, N.G.; project administration, N.G.; funding acquisition, N.G. and M.E. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the European Space Agency, grant number 4000128903/1 9/l-DT.

Data Availability Statement

The raw input data for this project is freely available from ESA (ftp://science-pds.cryosat.esa.int/TEMPO_SWATH_POINT/, accessed on 1 May 2020), NASA(https://nsidc.org/data/atl06/, accessed on 1 May 2020) and the Polar Geospatial Center (https://www.pgc.umn.edu/data/arcticdem/, accessed on 1 May 2020). Combined CryoSat-2 and ICESat-2 products are now also readily available from ESA’s CS2EO platform (https://cs2eo.org/, accessed on 1 May 2020).

Acknowledgments

A special thank you to Jerome Bouffard and Alessandro Di Bella at the European Space Agency for supporting this work and providing guidance throughout the project. An additional thank you goes to Livia Jakob, Carolyn Michael, Sarah Appleby, and George Malczyk for additional analysis and input.

Conflicts of Interest

The authors declare no conflict of interest.

References

Zwally, H.J.; Bindschadler, R.A.; Brenner, A.C.; Major, J.A.; Marsh, J.G. Growth of Greenland Ice Sheet: Measurement. Science 1989, 246, 1587–1589. [Google Scholar] [CrossRef] [PubMed]
Wingham, D.J.; Ridout, A.J.; Scharroo, R.; Arthern, R.J.; Shum, C.K. Antarctic Elevation Change from 1992 to 1996. Science 1998, 282, 456–458. [Google Scholar] [CrossRef]
Shepherd, A.; Wingham, D.; Payne, T.; Skvarca, P. Larsen Ice Shelf Has Progressively Thinned. Science 2003, 302, 856–859. [Google Scholar] [CrossRef] [Green Version]
Shepherd, A.; Wingham, D.J.; Mansley, J.A.D.; Corr, H.F.J. Inland Thinning of Pine Island Glacier, West Antarctica. Science 2001, 291, 862–864. [Google Scholar] [CrossRef]
Zwally, H.J.; Giovinetto, M.B.; Li, J.; Cornejo, H.G.; Beckley, M.A.; Brenner, A.C.; Saba, J.L.; Yi, D. Mass Changes of the Greenland and Antarctic Ice Sheets and Shelves and Contributions to Sea-Level Rise: 1992–2002. J. Glaciol. 2005, 51, 509–527. [Google Scholar] [CrossRef] [Green Version]
Wingham, D.J.; Shepherd, A.; Muir, A.; Marshall, G.J. Mass Balance of the Antarctic Ice Sheet. Philos. Trans. R. Soc. Math. Phys. Eng. Sci. 2006, 364, 1627–1635. [Google Scholar] [CrossRef]
Fricker, H.A.; Scambos, T.; Bindschadler, R.; Padman, L. An Active Subglacial Water System in West Antarctica Mapped from Space. Science 2007, 315, 1544–1548. [Google Scholar] [CrossRef] [Green Version]
Pritchard, H.D.; Arthern, R.J.; Vaughan, D.G.; Edwards, L.A. Extensive Dynamic Thinning on the Margins of the Greenland and Antarctic Ice Sheets. Nature 2009, 461, 971–975. [Google Scholar] [CrossRef]
Kääb, A.; Berthier, E.; Nuth, C.; Gardelle, J.; Arnaud, Y. Contrasting Patterns of Early Twenty-First-Century Glacier Mass Change in the Himalayas. Nature 2012, 488, 495–498. [Google Scholar] [CrossRef]
Bamber, J.L.; Griggs, J.A.; Hurkmans, R.T.W.L.; Dowdeswell, J.A.; Gogineni, S.P.; Howat, I.; Mouginot, J.; Paden, J.; Palmer, S.; Rignot, E.; et al. A New Bed Elevation Dataset for Greenland. Cryosphere 2013, 7, 499–510. [Google Scholar] [CrossRef]
McMillan, M.; Shepherd, A.; Sundal, A.; Briggs, K.; Muir, A.; Ridout, A.; Hogg, A.; Wingham, D. Increased Ice Losses from Antarctica Detected by CryoSat-2. Geophys. Res. Lett. 2014, 41, 3899–3905. [Google Scholar] [CrossRef]
Gourmelen, N.; Goldberg, D.N.; Snow, K.; Henley, S.F.; Bingham, R.G.; Kimura, S.; Hogg, A.E.; Shepherd, A.; Mouginot, J.; Lenaerts, J.T.M.; et al. Channelized Melting Drives Thinning under a Rapidly Melting Antarctic Ice Shelf. Geophys. Res. Lett. 2017, 44, 9796–9804. [Google Scholar] [CrossRef]
Gourmelen, N.; Escorihuela, M.J.; Shepherd, A.; Foresta, L.; Muir, A.; Garcia-Mondéjar, A.; Roca, M.; Baker, S.G.; Drinkwater, M.R. CryoSat-2 Swath Interferometric Altimetry for Mapping Ice Elevation and Elevation Change. Adv. Space Res. 2018, 62, 1226–1242. [Google Scholar] [CrossRef] [Green Version]
Rémy, F.; Parouty, S. Antarctic Ice Sheet and Radar Altimetry: A Review. Remote Sens. 2009, 1, 1212–1239. [Google Scholar] [CrossRef] [Green Version]
Nilsson, J.; Vallelonga, P.; Simonsen, S.B.; Sørensen, L.S.; Forsberg, R.; Dahl-Jensen, D.; Hirabayashi, M.; Goto-Azuma, K.; Hvidberg, C.S.; Kjær, H.A.; et al. Greenland 2012 Melt Event Effects on CryoSat-2 Radar Altimetry. Geophys. Res. Lett. 2015, 42, 3919–3926. [Google Scholar] [CrossRef]
Slater, T.; Shepherd, A.; Mcmillan, M.; Armitage, T.W.K.; Otosaka, I.; Arthern, R.J. Compensating Changes in the Penetration Depth of Pulse-Limited Radar Altimetry over the Greenland Ice Sheet. IEEE Trans. Geosci. Remote Sens. 2019, 57, 9633–9642. [Google Scholar] [CrossRef] [Green Version]
Gray, L. Brief Communication: Glacier Run-off Estimation Using Altimetry-Derived Basin Volume Change: Case Study at Humboldt Glacier, Northwest Greenland. Cryosphere 2021, 15, 1005–1014. [Google Scholar] [CrossRef]
Slater, T.; Shepherd, A.; McMillan, M.; Leeson, A.; Gilbert, L.; Muir, A.; Munneke, P.K.; Noël, B.; Fettweis, X.; van den Broeke, M.; et al. Increased Variability in Greenland Ice Sheet Runoff from Satellite Observations. Nat. Commun. 2021, 12, 6069. [Google Scholar] [CrossRef] [PubMed]
Arthern, R.J.; Wingham, D.J.; Ridout, A.L. Controls on ERS Altimeter Measurements over Ice Sheets: Footprint-Scale Topography, Backscatter Fluctuations, and the Dependence of Microwave Penetration Depth on Satellite Orientation. J. Geophys. Res. Atmos. 2001, 106, 33471–33484. [Google Scholar] [CrossRef] [Green Version]
Gray, L.; Burgess, D.; Copland, L.; Langley, K.; Gogineni, P.; Paden, J.; Leuschen, C.; van As, D.; Fausto, R.; Joughin, I.; et al. Measuring Height Change around the Periphery of the Greenland Ice Sheet with Radar Altimetry. Front. Earth Sci. 2019, 7, 146. [Google Scholar] [CrossRef]
Recchia, L.; Scagliola, M.; Giudici, D.; Kuschnerus, M. An Accurate Semianalytical Waveform Model for Mispointed SAR Interferometric Altimeters. IEEE Geosci. Remote Sens. Lett. 2017, 14, 1537–1541. [Google Scholar] [CrossRef]
Snauffer, A.M.; Hsieh, W.W.; Cannon, A.J.; Schnorbus, M.A. Improving Gridded Snow Water Equivalent Products in British Columbia, Canada: Multi-Source Data Fusion by Neural Network Models. Cryosphere 2018, 12, 891–905. [Google Scholar] [CrossRef] [Green Version]
Tollenaar, V.; Zekollari, H.; Lhermitte, S.; Tax, D.M.J.; Debaille, V.; Goderis, S.; Claeys, P.; Pattyn, F. Unexplored Antarctic Meteorite Collection Sites Revealed through Machine Learning. Sci. Adv. 2022, 8, eabj8138. [Google Scholar] [CrossRef] [PubMed]
Braakmann-Folgmann, A.; Donlon, C. Estimating Snow Depth on Arctic Sea Ice Using Satellite Microwave Radiometry and a Neural Network. Cryosphere 2019, 13, 2421–2438. [Google Scholar] [CrossRef] [Green Version]
Memarian Sorkhabi, O.; Asgari, J.; Amiri-Simkooei, A. Wavelet Decomposition and Deep Learning of Altimetry Waveform Retracking for Lake Urmia Water Level Survey. Mar. Georesources Geotechnol. 2022, 40, 361–369. [Google Scholar] [CrossRef]
Gray, L.; Burgess, D.; Copland, L.; Cullen, R.; Galin, N.; Hawley, R.; Helm, V. Interferometric Swath Processing of Cryosat Data for Glacial Ice Topography. Cryosphere 2013, 7, 1857–1867. [Google Scholar] [CrossRef] [Green Version]
Hawley, R.L.; Shepherd, A.; Cullen, R.; Helm, V.; Wingham, D.J. Ice-Sheet Elevations from across-Track Processing of Airborne Interferometric Radar Altimetry. Geophys. Res. Lett. 2009, 36, L22501. [Google Scholar] [CrossRef] [Green Version]
Krabill, W. IceBridge ATM L2 Icessn Elevation, Slope, and Roughness, Version 2; NASA Distribution Active Archive Center, National Snow Ice Data Center: Boulder, CO, USA, 2014.
Smith, B. ATLAS/ICESat-2 L3A Land Ice Height, Version 3; National Snow and Ice Data Center: Boulder, CO, USA, 2020. [Google Scholar]
Porter, C.; Morin, P.; Howat, I.; Noh, M.-J.; Bates, B.; Peterman, K.; Keesey, S.; Schlenk, M.; Gardiner, J.; Tomko, K.; et al. ArcticDEM, Version 3. 2019. Available online: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/OHHUKH (accessed on 20 February 2022).
Legrésy, B.; Rémy, F. Altimetric Observations of Surface Characteristics of the Antarctic Ice Sheet. J. Glaciol. 1997, 43, 265–275. [Google Scholar] [CrossRef] [Green Version]
Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods 2020, 17, 261–272. [Google Scholar] [CrossRef] [Green Version]
McNabb, R.; Nuth, C.; Kääb, A.; Girod, L. Sensitivity of Glacier Volume Change Estimation to DEM Void Interpolation. Cryosphere 2019, 13, 895–910. [Google Scholar] [CrossRef]
Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.; Chanan, G.; Killeen, T.; Lin, Z.; Gimelshein, N.; Antiga, L.; et al. PyTorch: An Imperative Style, High-Performance Deep Learning Library. In Advances in Neural Information Processing Systems 32; Curran Associates, Inc.: Nice, France, 2019; pp. 8024–8035. [Google Scholar]
Maas, A.L.; Hannun, A.Y.; Ng, A.Y. Rectifier Nonlinearities Improve Neural Network Acoustic Models. Proc. Icml. 2013, 30, 3. [Google Scholar]
Klambauer, G.; Unterthiner, T.; Mayr, A.; Hochreiter, S. Self-Normalizing Neural Networks. arXiv 2017, arXiv:1706.02515. [Google Scholar] [CrossRef]
Srivastava, N.; Hinton, G.; Krizhevsky, A.; Sutskever, I.; Salakhutdinov, R. Dropout: A Simple Way to Prevent Neural Networks from Overfitting. J. Mach. Learn. Res. 2014, 15, 1929–1958. [Google Scholar]
Ioffe, S.; Szegedy, C. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. In Proceedings of the 32nd International Conference on Machine Learning; Bach, F., Blei, D., Eds.; PMLR: Lille, France, 2015; Volume 37, pp. 448–456. [Google Scholar]
Ba, J.L.; Kiros, J.R.; Hinton, G.E. Layer Normalization. arXiv 2016, arXiv:1607.06450. [Google Scholar]
Sutskever, I.; Martens, J.; Dahl, G.; Hinton, G. On the Importance of Initialization and Momentum in Deep Learning. In Proceedings of the 30th International Conference on Machine Learning, Atlanta, GA, USA, 16–21 June 2013; Dasgupta, S., McAllester, D., Eds.; PMLR: Atlanta, GA, USA, 2013; Volume 28, pp. 1139–1147. [Google Scholar]
Kingma, D.P.; Ba, J. Adam: A Method for Stochastic Optimization. arXiv 2014, arXiv:1412.6980. [Google Scholar]
Huber, P.J. Robust Estimation of a Location Parameter. Ann. Math. Stat. 1964, 35, 73–101. [Google Scholar] [CrossRef]
Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; et al. Scikit-Learn: Machine Learning in Python. J. Mach. Learn. Res. 2011, 12, 2825–2830. [Google Scholar]
Earthwave. The University of Edinburgh. isardSAT CryoTEMPO-EOLIS—Elevation over Land Ice from Swath—Product Handbook. 2020. Available online: https://Earth.Esa.Int/Eogateway/Documents/20142/37627/CryoTEMPO-Thematic-Product-Handbook.Pdf (accessed on 20 February 2022).
Gray, L.; Burgess, D.; Copland, L.; Dunse, T.; Langley, K.; Moholdt, G. A Revised Calibration of the Interferometric Mode of the CryoSat-2 Radar Altimeter Improves Ice Height and Height Change Measurements in Western Greenland. Cryosphere 2017, 11, 1041–1058. [Google Scholar] [CrossRef] [Green Version]
Davis, C.H.; Moore, R.K. A Combined Surface-and Volume-Scattering Model for Ice-Sheet Radar Altimetry. J. Glaciol. 1993, 39, 675–686. [Google Scholar] [CrossRef] [Green Version]
Ridley, J.K.; Partington, K.C. A Model of Satellite Radar Altimeter Return from Ice Sheets. Int. J. Remote Sens. 1988, 9, 601–624. [Google Scholar] [CrossRef]
Wingham, D.J.; Francis, C.R.; Baker, S.; Bouzinac, C.; Brockley, D.; Cullen, R.; de Chateau-Thierry, P.; Laxon, S.W.; Mallow, U.; Mavrocordatos, C.; et al. CryoSat: A Mission to Determine the Fluctuations in Earth’s Land and Marine Ice Fields. Adv. Space Res. 2006, 37, 841–871. [Google Scholar] [CrossRef]
Krabill, W. Greenland Ice Sheet: Increased Coastal Thinning. Geophys. Res. Lett. 2004, 31, L24402. [Google Scholar] [CrossRef] [Green Version]
Bingham, A.W.; Drinkwater, M.R. Recent Changes in the Microwave Scattering Properties of the Antarctic Ice Sheet. IEEE Trans. Geosci. Remote Sens. 2000, 38, 1810–1820. [Google Scholar] [CrossRef]
Brunt, K.M.; Neumann, T.A.; Smith, B.E. Assessment of ICESat-2 Ice Sheet Surface Heights, Based on Comparisons over the Interior of the Antarctic Ice Sheet. Geophys. Res. Lett. 2019, 46, 13072–13078. [Google Scholar] [CrossRef] [Green Version]
Luthcke, S.B.; Thomas, T.C.; Pennington, T.A.; Rebold, T.W.; Nicholas, J.B.; Rowlands, D.D.; Gardner, A.S.; Bae, S. ICESat-2 Pointing Calibration and Geolocation Performance. Earth Space Sci. 2021, 8, e2020EA001494. [Google Scholar] [CrossRef]
Fettweis, X.; Hofer, S.; Krebs-Kanzow, U.; Amory, C.; Aoki, T.; Berends, C.J.; Born, A.; Box, J.E.; Delhasse, A.; Fujita, K.; et al. GrSMBMIP: Intercomparison of the Modelled 1980–2012 Surface Mass Balance over the Greenland Ice Sheet. Cryosphere 2020, 14, 3935–3958. [Google Scholar] [CrossRef]
Kern, M.; Cullen, R.; Berruti, B.; Bouffard, J.; Casal, T.; Drinkwater, M.R.; Gabriele, A.; Lecuyot, A.; Ludwig, M.; Midthassel, R.; et al. The Copernicus Polar Ice and Snow Topography Altimeter (CRISTAL) High-Priority Candidate Mission. Cryosphere 2020, 14, 2235–2251. [Google Scholar] [CrossRef]

Figure 1. (a) The geographical distribution of the ~29 million collocated SARIn vs. LIDAR observations in this study and (b) the statistical distribution of the same observations with a fitted normal inverse Gaussian distribution trend line. (c.f. Section 3.2 for details on the standard error calculation).

Figure 2. Summary of the observed and modelled LIDAR-SARIn differences as a function of a selection of input parameters for all 3.6 million unseen test data points across GrIS. The distribution of observation counts per parameter is shown in the grey histogram for each plot.

Figure 3. Gridded maps of observed LIDAR-SARIn differences for February 2020 (a) and August 2020 (d) and their difference (g), and equivalent model predictions from the ordinary least squares model (b,e,h), and from the neural network model (c,f,i). Note that the scale bars have been adjusted between observations and model predictions to allow for a clearer visual comparison of broad trends.

Figure 4. LIDAR-SARIn adjustment as a function of distance to POCA broken down by 6 quantile bands of coherence, power (dB), and slope across.

Figure 5. Spatial plots of MAR model outputs across GrIS for (a) cumulative snowfall in February 2020 and the prior five months; (b) cumulative snowfall in August 2020 and the prior five months; (c) the change in accumulative snowfall between August 2020 and February 2020: i.e., the difference between (b) and (a); (d) the approximate average snow density in the top 2 m of the snowpack for February 2020; (e) the average surface mass balance over 12 months prior to and including February 2020; (f) the total summer melt for June, July and August 2020.

Table 1. Summary of model input parameters.

Model Input Parameters	Source
Distance to POCA	CryoSat2-SARIn
Power	CryoSat2-SARIn
Coherence	CryoSat2-SARIn
Leading Edge Width	CryoSat2-SARIn
Relative Elevationⁱ (i = {−3, −2, −1, 1, 2, 3})	CryoSat2-SARIn
Relative Elevation Mean^j (j = {10, 20})	CryoSat2-SARIn
Across Track Slope	ArcticDEM
Along Track Slope	ArcticDEM
South/North Aspect Vector Component	ArcticDEM
West/East Aspect Vector Component	ArcticDEM

Table 2. Summary of OLS and NN accuracy on the unseen test dataset. All values relate to the difference between the model-based estimate of the LIDAR-SARIn adjustment and the observations.

Model	Mean (m)	Standard Error (m)	Median (m)	RMSE (m)	MAD (m)	r
Ordinary Least Squares (OLS)	−0.175	0.020	−0.105	3.183	1.102	0.372
Neural Network (NN)	−0.048	0.019	−0.053	2.958	0.987	0.500

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Horton, A.; Ewart, M.; Gourmelen, N.; Fettweis, X.; Storkey, A. Using Deep Learning to Model Elevation Differences between Radar and Laser Altimetry. Remote Sens. 2022, 14, 6210. https://doi.org/10.3390/rs14246210

AMA Style

Horton A, Ewart M, Gourmelen N, Fettweis X, Storkey A. Using Deep Learning to Model Elevation Differences between Radar and Laser Altimetry. Remote Sensing. 2022; 14(24):6210. https://doi.org/10.3390/rs14246210

Chicago/Turabian Style

Horton, Alex, Martin Ewart, Noel Gourmelen, Xavier Fettweis, and Amos Storkey. 2022. "Using Deep Learning to Model Elevation Differences between Radar and Laser Altimetry" Remote Sensing 14, no. 24: 6210. https://doi.org/10.3390/rs14246210

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu