<figure><img src="image_1.png" /></figure>

HAL
open science

# Improving Satellite Waveform Altimetry Measurements With a Probabilistic Relaxation Algorithm

Song Shu, Hongxing Liu, Frédéric Frappart, Emily Lei Kang, Bo Yang, Min Xu, Yan Huang, Bin Wu, Bailang Yu, Shujie Wang, et al.

<figure><img src="image_2.png" /></figure>

## To cite this version:

Song Shu, Hongxing Liu, Frédéric Frappart, Emily Lei Kang, Bo Yang, et al. Improving Satellite Waveform Altimetry Measurements With a Probabilistic Relaxation Algorithm. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59 (6), pp.4733-4748. 10.1109/TGRS.2020.3010184. hal-04720543

HAL Id: hal-04720543

https://hal.science/hal-04720543v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Improving Satellite Waveform Altimetry
Measurements With a Probabilistic
Relaxation Algorithm

Song Shu<sup>D</sup>, Hongxing Liu<sup>D</sup>, Member, IEEE, Frédéric Frappart, Emily Lei Kang, Bo Yang<sup>D</sup>, Min Xu<sup>D</sup>,
Yan Huang, Bin Wu, Bailang Yu<sup>D</sup>, Senior Member, IEEE, Shujie Wang,
Richard Beck, and Kenneth Hinkel

**Abstract**—The Geoscience Laser Altimeter System onboard the NASA Ice, Cloud, and land Elevation Satellite (ICESat/GLAS) provided elevation measurements of Earth's surface between 2003 and 2009. The centroid and maximum-amplitude-peak (MAP) retracking methods have been designed and applied to process the returned laser waveforms for elevation measurements. Although these two methods work well in general, they may generate erroneous measurements when the returned waveform was complicated by adverse atmospheric conditions (clouds, ice fogs, blowing snow, and dust storms). The centroid retracking method is often more severely affected when compared with the MAP retracking method. In this study, we present a new retracking method that exploits the spatial contextual information from neighboring footprints along the satellite ground track, in addition to the single return waveform shape information. Our method uses a probabilistic relaxation (PR) algorithm to integrate the spatial contextual information and the waveform shape information to identify the waveform

Manuscript received April 19, 2019; revised September 25, 2019, February 26, 2020, and May 9, 2020; accepted July 9, 2020. This work was supported in part by the National Science Foundation (NSF) under Grant ARC-1107607 and Grant 0713813 and in part by the National Natural Science Foundation of China under Grant 41701502, Grant 41771461, and Grant 41471449. (Corresponding authors: Hongxing Liu; Song Shu.)

Song Shu is with the Department of Geography and Planning, Appalachian State University, Boone, NC 28608 USA (e-mail: shus@appstate.edu).

Hongxing Liu is with the Department of Geography, The University of Alabama, Tuscaloosa, AL 35487 USA, and also with the School of Geographical Science, East China Normal University, Shanghai 200241, China (e-mail: hongxing.liu@ua.edu).

Frédéric Frappart is with the Laboratoire d'Etudes en Géophysique et Océanographie Spatiales (LEGOS), UMR 5566, CNRS/IRD/UPS, OMP, 31400 Toulouse, France.

Emily Lei Kang, Richard Beck, and Kenneth Hinkel are with the Department of Mathematical Sciences, University of Cincinnati, Cincinnati, OH 45221 USA, and also with the Department of Geography and Geographic Information Science, University of Cincinnati, Cincinnati, OH 45221 USA.

Bo Yang is with the Department of Sociology, University of Central Florida, Orlando, FL 32816 USA.

Min Xu is with the College of Marine Science, University of South Florida at St. Petersburg, St. Petersburg, FL 33701 USA.

Yan Huang, Bin Wu, and Bailang Yu are with the Key Laboratory of Geographical Information Science, Ministry of Education, School of Geographical Science, East China Normal University, Shanghai 200241, China.

Shujie Wang is with the Department of Geography, Earth and Environmental Systems Institute (EESI), The Pennsylvania State University, University Park, PA 16802 USA, and also with the Institute for Computational and Data Sciences, The Pennsylvania State University, University Park, PA 16802 USA.

Color versions of one or more of the figures in this article are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2020.3010184

peak that most likely represents the true surface elevation, rather than simply detecting the peak with the maximum magnitude. For different types of land surfaces, such as inland lakes, polar tundra, ice sheet, and sand deserts, we demonstrate that our new PR retracking method is able to produce more reliable, consistent, and accurate elevation measurements than the standard NASA ICESat/GLAS data products. The root mean squares error (RMSE) is reduced from 0.85 to 0.17 m for inland lake, from 0.81 to 0.23 m for polar tundra, from 1.25 to 0.33 m for ice sheet, and from 2.48 to 2.34 m for sand desert.

**Index Terms**—ICESat/GLAS, probabilistic relaxation (PR), retracking algorithm, spatial contextual information, waveforms.

## I. INTRODUCTION

SATELLITE altimetry has been widely used for global elevation measurements of various Earth surface types on a regular basis since the mid-1980s [1]–[4]. The spaceborne altimeter sensors emit a series of electromagnetic pulses in the nadir direction toward Earth’s surface and accurately measure the two-way travel time of the signal. Depending on the frequency of the emitted pulses, altimeters can be classified into two types: radar altimeters and laser altimeters. Radar altimeters employ microwave electromagnetic waves at Ku-band [e.g., 13.5 GHz for the Radar Altimeter 2 onboard the Environmental Satellite (ENVISAT/RA-2)] or Ka-band [e.g., 35.5 GHz for the Ka-band altimeter onboard the Satellite with ARgos and ALTiKa (SARAL/AltiKa)]. For most of the missions, a secondary frequency channel at C-band (e.g., 5.3 GHz for Jason-1/Poseidon-2) or S-band (e.g., 3.2 GHz for ENVISAT/RA-2) is also used for the correction of the range delay caused by the ionosphere over the ocean [5]. Laser altimeters operate at visible and near-infrared wavelengths. ICESat/GLAS was launched by NASA in January 2003 with an orbit of 600-km altitude and an inclination angle of 94°. The Geoscience Laser Altimeter System (GLAS) was the sole payload for the mission and worked at two wavelengths, 532 and 1064 nm [6]. It was the first satellite laser altimetry mission that provided elevation measurements of Earth surface between 86°S and 86°N latitudes from 2003 to 2009 [6]. GLAS transmitted 40 laser pulses per second toward Earth’s surface. The illuminated footprints of laser pulses on the ground surface have an elliptical shape with a diameter of
about 70 m, and these laser footprints were spaced at 170-m intervals along the satellite tracks. The primary objective of ICESat/GLAS mission was to measure ice sheet elevation changes in Antarctica and Greenland [6], [7]. The data have also been widely used to estimate sea-ice freeboard and sea ice thickness [8], forest canopy height [9], snow depth [10]–[12], as well as lake water-level changes [13], [14]. The accuracy of ICESat/GLAS elevation measurements could reach 2 cm under clear-sky conditions over the Antarctic Ice Sheet [15]. In normal situations, the expected vertical accuracy is about 14 cm [16], [17].

The elevation of each laser footprint is jointly determined by satellite orbital position, laser beam direction, and the range between satellite and the reflecting surface. The laser energy intensity of transmitted and returned pulses are densely sampled and recorded over time as a line plot, which is referred to as laser waveform. The range is calculated by tracking the elapsed time between the emission and the reception of the laser pulse. The elapsed time between the emission and the reception is precisely determined by identifying the reference points on the transmitted and the returned altimeter waveforms. The transmitted waveform is approximately a simple Gaussian function [6]. The waveform returned from a flat and homogeneous surface should have a single dominant peak with a Gaussian shape. For a rough and heterogeneous surface, the returned waveform may exhibit several peaks with a complex shape [18], [19].

In the creation of the ICESat/GLAS data products, Earth's surface was classified into four major land cover types (land, ocean, sea ice, and ice sheet) to determine the reference points on the returned waveforms with different retracking methods. The surfaces of ocean, sea ice and ice sheet, are relatively homogeneous and flat, while the land surfaces are usually heterogeneous and complex. Two standard retracking methods were designed specifically to cope with these two categories of surfaces according to the characteristics of their waveform shapes. The maximum-amplitude-peak (MAP) retracking method is suitable for relatively smooth and gradually changing surfaces [19], while the centroid retracking method is designed for all different types of terrestrial land surfaces. The terrestrial land surface is usually a complex mosaic of different geomorphic features. The land surface elevation may vary significantly within the laser footprint, and the elevation measurement corresponding to the geometric centroid reference point represents the average surface height within the illuminated footprint. In contrast, the surfaces of ice sheets, sea ice, and oceans are relatively smooth and homogeneous, and the returned waveforms from such surfaces are mostly characterized by a single dominant peak with a Gaussian shape. The MAP retracking method identifies individual peaks on the returned waveform by fitting a series of Gaussian functions and then selects the Gaussian peak with a maximum amplitude to represent the surface height within the footprint.

Although the centroid method is the choice for land surfaces with drastic elevation variations, the MAP method could be a better choice for some terrestrial geographical features with a relatively small elevation variation under good weather

conditions, such as, inland lakes, sand deserts, coastal zones, and so on. Under bad weather conditions, the MAP retracking method may also fail to produce reliable elevation measurements for these relatively flat surfaces due to the atmospheric influence [15], [20]–[23]. In fact, atmospheric events may lead to the low energy intensity and the non-Gaussian shape of the returned waveform [6]. These affected waveforms may significantly deviate from the theoretical waveform model and hence fail the two standard retracking methods, leading to centimeter to meter level of error in elevation retrievals [15], [20], [23]. In some extreme situations, the error could be up to 8 m in Antarctica [23].

The ICESat/GLAS measurements contaminated by atmospheric scattering usually exhibit high detector gain value and low surface reflectivity [24]–[29]. The detector gain value ranges from 1 to 255. It is set to a higher level when the returned laser pulse energy becomes weaker. The reflectivity is the ratio of the returned laser pulse energy to the emitted energy, ranging from 0 to 1. The higher gain value and the lower reflectivity indicate the stronger level of atmospheric scattering. The Arctic region is covered by cloud over half the time of the year [30]. Blowing snow is also a common phenomenon in the winter season in this region [31]. ICESat/GLAS measurements in the polar regions are more likely subject to the scattering than in other regions. By using the two variables, surface reflectivity ($i$_refl UC in GLA12–15 products) and detector gain value ($i$_gval_rcv in GLA12–15 products), we found out that in the Arctic coastal plain of northern Alaska, about 19% and 13% of the total observations collected by ICESat/GLAS during 2003 and 2009 have the detector gain value higher than 250 or a reflectivity value less than 0.1 and hence may have been influenced by the atmospheric scattering due to bad weathers. The relatively large percentages indicate the widespread influences of atmospheric scattering in the Arctic regions.

Some heuristic approaches have been used to identify and then exclude the ICESat/GLAS observations that were heavily affected by atmospheric scattering. Yi et al. [24] used the gain threshold value of 30, Kwok et al. [25] used 50, and Felikson et al. [26] used 150 to filter potentially contaminated footprints over Arctic sea ice and Greenland Ice Sheet, respectively. Later, Yi et al. [27] recommended to use different thresholds of gain values for the measurements collected in different campaigns. Siegfried et al. [28] suggested that the measurements with surface reflectivity less than 0.24 most likely do not meet the data quality requirement of ICESat mission. Yi et al. [27] and Zwally et al. [29] used 0.05 as the surface reflectivity to identify the affected measurements. However, no robust numerical algorithms have been reported to treat the waveforms contaminated by the atmospheric scattering for the purpose of reliable elevation retrieval.

In this study, we present a new robust retracking method that exploits the spatial contextual information from neighboring footprints along the satellite ground track, in addition to the waveform shape information of an individual footprint as in the standard centroid and MAP retracking methods. Our method utilizes a probabilistic relaxation (PR) algorithm to integrate the spatial contextual information with the individual
<figure><img src="image_4.png" /><figcaption>Fig. 1. Transmitted and returned waveforms and the retracking methods. This figure is adapted from Fig. 3 in the ICESat Algorithm Theoretical Basis Document [19].</figcaption></figure>

waveform shape information to identify the reference peak that most likely corresponds to the true surface elevation.

In the rest of the paper, we will first introduce the standard NASA ICESat/GLAS data products and two associated standard waveform retracking methods. Then, the conditions that would fail these two retracking methods will be discussed. Next, we will present the mathematical formulation of the PR retracking algorithm. This is followed by a number of application examples to demonstrate its effectiveness. Finally, we summarize key research findings and draw conclusions.

## II. DATA SETS

ICESat/GLAS altimetry system conducted a total of 18 measurement campaigns during 2003 and 2009. Each campaign lasted 12–55 days, as a subcycle of the more densely spaced 91-day exact repeat orbit [29]. There are 15 different types of ICESat/GLAS data products available at National Snow and Ice Data Center (NSIDC) (http://nsidc.org/data/icesat). GLA01 provides the original transmitted and returned waveform data. GLA12, GLA13, GLA14, and GLA15 provide surface elevations for ice sheets, sea ice, land, and oceans, respectively. In this study, GLA01, GLA12, and GLA14 data products are used.

We used three independent high-resolution topographical data sets to validate and evaluate the surface elevation estimates produced by our PR retracking method. Those data sets include IceBridge LVIS elevation data over the Greenland Ice Sheet [32], the Arctic DEM over the tundra in the Arctic Coastal Plain, and the Shuttle Radar Topography Mission (SRTM) 1 arc-second elevation data over the Taklimakan desert in China. IceBridge LVIS elevation data were acquired by an airborne laser altimeter with a footprint diameter of 20 m and a spacing between footprints of 10 m. The vertical accuracy of LVIS data is better than 0.12 m [33]. The LVIS data set used in this study was acquired on 2 May 2012, following ICESat/GLAS track 0174. The Arctic DEM data set provided by the Polar Geospatial Center (http://pgc.umn.edu/arcticdem) at the University of Minnesota has a high spatial resolution of 5 m, which was created from high-resolution satellite

stereo images acquired by DigitalGlobe Inc. The SRTM DEM, provided by U.S. Geological Survey (USGS) Earth Explorer (https://earthexplorer.usgs.gov/), was created using the interferometry radar data acquired by the SRTM in 2000. It has a spatial resolution of 30 m and variable vertical accuracy depending on the topography and the presence of vegetation. The accuracy over vegetated area was estimated to be 4.07 ± 0.47 m [34]. The accuracy of SRTM measurements is expected to be even better for sand deserts where no vegetation cover exists at all. To ensure fair and valid comparison, all elevation data sets have been vertically referenced to the same geoid—EGM2008.

## III. ICESAT/GLAS STANDARD RETRACKING METHODS

The transmitted and returned laser pulse waveforms represent the pulse energy over the elapsed time in nanosecond (ns), as shown in Fig. 1. The elevation of the reflecting surface is given as follows [6]:

$$h = H - (R + \Sigma \Delta R_c) - h_{\text{geoid}} \quad (1)$$

where $h$ is the elevation of the laser footprint, $H$ is the height of the mass center of the satellite above the reference ellipsoid, $R$ is the nadir range between the mass center of the satellite and the laser footprint, $\Sigma \Delta R_c$ represents the sum of the instrumental and geophysical corrections applied to the range $R$ to account for instrumental biases, propagation delays of the electromagnetic waves in the atmosphere and geophysical effects, and $h_{\text{geoid}}$ is the difference between the reference ellipsoid and a specific geoid model (EGM2008 in this study).

The range $R$ is given as follows [19]:

$$R = (T_{\text{return}} - T_{\text{trans}}) \times c/2 \quad (2)$$

where $c$ is the speed of light; $T_{\text{trans}}$ and $T_{\text{return}}$ represent the times of the transmission and the reception of the laser pulse, respectively, which are precisely determined by the two reference points on the transmitted and the returned waveforms, as shown in Fig. 1.
The standard centroid retracking method defines the geo-
metric centroid of the waveform above the background noise
level as the reference point to calculate the range and surface
elevation [35]. The standard MAP retracking method identifies
individual peaks (up to two peaks) on the returned waveform
by fitting multiple Gaussian functions to the peaks. The
midpoint of the highest Gaussian peak is then identified as
the time reference point to calculate the range and surface
elevation [19]. The transmitted waveform is fit by a single
Gaussian peak. For the standard MAP retracking method
(GLA12, 13, and 15 products), the returned waveform was
fit by up to two Gaussian peaks. For GLA14 product, even
though the range was determined by the geometric centroid
of the waveform, the returned waveform was also modeled by
up to six Gaussian peaks. The information of these peaks, such
as the centroid, amplitude, area, and sigma of each Gaussian
peak, were also included in the data product for each laser
shot. The number of Gaussian peaks on the returned waveform
was controlled by merging the adjacent peaks during the
fitting process. For GLA12, 13, and 15 products, the adjacent
peaks with an interval of less than 30 ns are merged, while
for GLA14 product, the interval between adjacent peaks was
15 ns [19]. This merging procedure often results in a wide
Gaussian peak for the returned waveform, as shown by the
green line (W<sub>FR</sub>) in Fig. 1. Both the standard centroid and the
MAP method try to capture the main shape of returned wave-
form to represent the surface elevation within the illuminated
area. As a result, they often omit the small signal peak on the
returned waveform.

In general, the centroid retracking method is a good choice
for terrestrial land surfaces due to its relatively high relief
and heterogeneous vegetation cover. However, many terres-
trial geographical features, such as inland lakes, large rivers,
sand deserts, polar tundra, and coastal zones, are relatively
homogeneous with a small elevation variation. For these
types of terrestrial surfaces, the MAP retracking method may
be more suitable and reliable than the centroid retracking
method [19], [21]. Also, it should be pointed out that some
terrestrial landscapes exhibit strong seasonal variations. In the
wet season, wetlands are likely to be filled with water during
the flood period and effectively become flat surfaces. During
the long winter season, the polar tundra regions are more
homogeneous and smooth with a full, cover of snow and ice,
compared with the brief summer season when mosses, herba-
ceous sedges and grasses prevail [36]. For those seasonally
changing terrestrial surfaces, the MAP retracking method may
be a better choice to retrieve their surface elevations for some
specific seasons. However, ICESat/GLAS GLA14 products for
all types of terrestrial land surfaces were produced solely
by the centroid retracking method, without distinguishing
between permanent and seasonal homogeneous terrestrial land
surfaces for the MAP retracking method.

No matter which retracking method is adopted, producing
a reliable measurement is a challenge when the returned
waveforms are affected by bad weather conditions. The
suspended particles of clouds, ice fog, blowing snow, and
dust could cause multiple particle scattering on laser pulse
photons and increase the photon travel path, referred to as

"atmospheric path delay" [21]. The lower the cloud is, the greater the scattering effect would be [22]. The near-surface blowing snow and ice fog can have a much stronger scattering effects than clouds. Even a thin layer of blowing snow can result in a large path delay [20]. The scattering changes the returned waveform to non-Gaussian shape and could introduce many spurious peak features as shown in Fig. 1. Another typical manifestation of the weather contamination is the long tail on the right side of the returned waveform. Both the MAP and the centroid retracking methods try to capture the main shape of the returned waveform to represent the surface elevation, but are often misled by spurious peak features and the long tail. Under the influences of the long tail and the spurious peaks, the time reference points defined by them are often shifted to the right on the time axis, resulting in a longer range and hence biased, lower elevation measurements. We recognize that the shape information of a single waveform alone is not adequate to address the problem associated with the spurious peaks and the long tail. In this research, we attempt to introduce the spatial contextual information to tackle this problem.

IV. NEW PR RETRACKING METHOD

A. Spatial Contextual Information Along the Satellite Track

On a flat or gradually changing terrain surface, the adjacent
footprints along the satellite track would generate waveforms
with similar shape characteristics due to the geographical
proximity and spatial autocorrelation. The spatial contextual
information among the adjacent footprints may offer additional
clues to help retrieve the true surface elevation.

Fig. 2(a) shows a five-footprint neighborhood of the
satellite track on a horizontal continuous smooth surface.
These footprints should produce gradually varying elevation
measurements lacking a dramatic change. Without signal
contamination, the returned waveforms of these footprints
should be dominated by a single Gaussian peak, and the
dominant peaks of these waveforms could align on the time
axis with close time proximity. Assuming that the highlighted
red footprint in Fig. 2(a) is a footprint contaminated by thin
cloud cover and that no dominant peak can be clearly iden-
tified for its waveform, this contaminated waveform probably
fails by the MAP retracking method. However, when this
waveform is examined together with its neighboring wave-
forms in Fig. 2(b), a small peak feature on this contaminated
waveform [as shown in Fig. 2(c)] stands out, which is close
to the dominant peaks of the neighboring waveforms. With
reference to these neighboring waveforms, we can infer that
this small peak on the contaminated waveform of the central
red footprint most likely corresponds to the true surface
elevation. This small signal peak is often missed by the two
standard retracking methods that rely solely on the major shape
information of the single waveform returned from the red
footprint as shown in Fig. 1. For our application examples
in Section V, we adopted a three-footprint or a five-footprint
neighborhood, depending on the surface roughness at each
example site.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Spatial contextual information in the neighboring footprints along a satellite track. (a) Definition of the neighborhood. (b) Central contaminated waveform and its neighboring waveforms. (c) Gaussian peaks and the decomposition of the central contaminated waveform. The centroid and the MAP retracking points on time axis are denoted as a blue point and a black point, respectively, in (c).</figcaption></figure>

By integrating the waveform shape information of the central contaminated footprint with the spatial contextual information from the neighboring footprints, our new retracking method is capable of identifying the signal peak for elevation measurement. It consists of several computational steps as described below.

## B. Decomposition of the Returned Waveform

The returned waveform can be modeled by a combination of multiple Gaussian functions [6], as defined in the following equations:

$$W(t) = \varepsilon + \sum_{j=1}^{M} W_j(t) \quad (3)$$

$$W_j(t) = A_j e^{\frac{(t-t_j)^2}{2\sigma_j^2}} \quad (4)$$

where $W(t)$ is the energy intensity of returned waveform at time $t$, $W_j(t)$ is the intensity of Gaussian peak $j$ ($j = 1, \dots, M$) at time $t$, $M$ is the total number of Gaussian peaks on the waveform, $A_j$ is the amplitude of Gaussian peak $j$, $t_j$ and $\sigma_j$ are the center location on the time axis and the standard deviation (bandwidth) of Gaussian peak $j$, respectively, and $\varepsilon$ is the background noise level of the waveform. Following the method by Brenner et al. [19], the initial peaks are identified using the second derivative of the waveform and a nonlinear least square fit is then implemented to refine the parameters for each Gaussian peak. The fitting minimizes the difference between the modeled waveform and the observed waveform. To retain the small peaks on the returned waveform, the adjacent peaks with an interval of less than 3 ns (instead of 15 or 30 ns) are merged during the fitting process. This interval value is adopted since 3 ns is the least number of bins on a waveform to form a Gaussian peak, one bin as the

starting point, one bin as the crest point, and one bin as the end point. The number of Gaussian peaks is determined by the complexity of the returned waveform shape. As shown in Fig. 2(c), the central contaminated waveform is modeled by four Gaussian peaks, instead of a single wide Gaussian peak given by MAP in Fig. 1. Each of the four neighboring waveforms in Fig. 2(b) is also decomposed into individual Gaussian peaks using the same procedure.

Our new retracking method tries to identify the peak with the largest probability of being the signal peak for the surface elevation retrieval. Each Gaussian peak of the central contaminated waveform is assigned an initial probability of being the signal peak according to its amplitude, which is calculated using (5) in Section IV-C. Then, the initial probability of each peak is iteratively modified through its compatibility with the peaks of the neighboring waveforms. Once the process has converged or reached the predefined maximum number of iterations, the peak with the maximum posterior probability is selected as the signal peak for the elevation retrieval. As shown in Fig. 2(b) and (c), the first peak (the red dotted line) is not the one with the largest amplitude, but its probability of being the signal peak should be very high when the neighboring waveforms are taken into consideration. We adopt the PR method to incorporate the spatial contextual information from the neighboring waveforms.

## C. Identification of the Signal Peak Using the PR Method

The PR method was first introduced by Rosenfeld [37]. It is an iterative technique that allows the incorporation of contextual information in the labeling (classification) problem. The basic idea is to make a tentative, rather than firm, initial classification of an object based on the attributes of the object, and then repeatedly refine the classification decision based on the local contextual information at each iteration [38].
The initial classification decision will be positively or negatively reinforced by the contextual information, depending on the compatibility between the object being processed and its neighboring objects. The mathematical formulation of the PR method is in line with Bayes' theorem, which updates the prior probability of an event to a posterior probability by incorporating the conditions related to the event. In this study, we implemented the PR algorithm to compute the posterior probability of a Gaussian peak to be the signal peak by introducing the spatial contextual information. The initial prior probability of a peak is determined according to its amplitude, then the prior probability will be positively or negatively updated at each iteration by evaluating its compatibility with the neighboring waveforms.

After the returned waveforms, including the central waveform and the neighboring waveforms, are decomposed into several Gaussian peaks [see Fig. 2(c)], an initial prior probability of being the signal peak is assigned for each peak according to its peak amplitude (the shape information). The higher the peak, the larger the initial prior probability will be. Numerically, the initial prior probability of a peak is defined by the following equation:

$$P_{ij}^{0} = \frac{A_{ij}}{\sum_{j=1}^{M_i} A_{ij}} \qquad (5)$$

$$\sum_{j=1}^{M_i} P_{ij}^{0} = 1 \qquad (6)$$

where $P_{ij}^{0}$ is the initial prior probability of peak j on waveform i, $A_{ij}$ is the amplitude of the peak j on waveform i, and $M_i$ is the number of peaks on the waveform i, which varies among different waveforms. $P_{ij}^{0}$ is in the range from 0 to 1, and the sum of the probabilities of all peaks on the waveform i is equal to unity as shown in (6).

The spatial contextual information (conditional information) is quantified by the compatibility between each candidate peak of the central waveform and the peaks of the neighboring waveforms. The total number of footprints along the satellite track is S. Each footprint along the satellite track will be treated as a central footprint and is denoted by s (s = 1, 2, ..., S). The candidate peak of the central waveform (s) under examination is denoted by index k (k = 1, 2, ..., L). The compatibility ($Q_k^r$) of candidate peak k with all the peaks on the neighboring waveforms at iteration r (r = 1, 2, ..., R) is given by the following equations:

$$Q_k^r = \sum_{i=1}^{N} \sum_{j=1}^{M_i} w_i C_{k,ij} P_{ij}^{r-1} \qquad (7)$$

$$C_{k,ij} = \frac{\frac{1}{|t_k - t_{ij}|}}{\sum_{j=1}^{M_i} \frac{1}{|t_k - t_{ij}|}} \qquad (8)$$

where N is the number of neighboring footprints, $w_i$ is the weight of the neighboring waveform i, $M_i$ is the number of peaks on the neighboring waveform i, $t_k$ is the midpoint location of the candidate peak k on the central waveform under examination, $t_{ij}$ is the midpoint location of the peak j of the neighboring waveform i on the time axis, and $C_{k,ij}$ is the

normalized temporal proximity between the candidate peak k on the central waveform under examination and the peak j on the neighboring waveform i. The closer the two peaks, the higher $C_{k,ij}$ is. $P_{ij}^{r-1}$ is the probability of the peak j to be signal peak for the neighboring waveform i at (r - 1)th iteration.

At each iteration r, the probability $P_k^r$ of peak k on the waveform under examination will be updated using the prior probability $P_k^{r-1}$ in the preceding iteration (r - 1) and the compatibility $Q_k^r$ at current iteration r as in the following equation:

$$P_k^r = \frac{P_k^{r-1} Q_k^r}{\sum_{k=1}^{L} P_k^{r-1} Q_k^r} \qquad (9)$$

where the compatibility $Q_k^r$ at iteration r is calculated based on the prior probability of the peaks on all neighboring footprints in the preceding iteration (r - 1) according to (7). The posterior probability of each peak on the central waveform under examination (k = 1, 2, ..., L) is updated using the computational steps given in (7)–(9).

For each iteration r, every footprint along the satellite track will be processed as a central footprint using the procedure described above. We will then obtain the posterior probability $P_{ij}^r$ for all waveform peaks of all footprints along the satellite track, which will be used as the prior probabilities for the computation at the next iteration (r + 1).

The calculated probabilities of the peaks on every waveform along the satellite track at current iteration (r) are compared with their probabilities at the preceding iteration (r - 1). We consider that the convergence state of the waveform under examination is achieved if the magnitude of the probability changes is smaller than a predefined threshold $\alpha$ as in the following equation:

$$\frac{1}{M_i} \sum_{j=1}^{M_i} |P_j^r - P_j^{r-1}| < \alpha \qquad (10)$$

where $P_{ij}^r$ and $P_{ij}^{r-1}$ are the probabilities of peak j on the waveform i at the current and preceding iterations. The convergence criterion of the iterative process is set by the threshold value of $\alpha$, which is set to 0.005 in our case studies. The whole iteration process will cease, when the convergence criteria are satisfied for all footprints along the satellite track or when the predefined maximum number of iterations (100 for our case studies) is reached. The threshold value and the maximum iteration are determined through trial and error. When $\alpha$ reaches 0.005, the probability of each peak usually becomes stable without significant change. In most cases, the iteration process converges after a few iterations. The maximum iteration (100) allows for the convergence of almost all footprints. For the four application examples in Section V, there were only five footprints in the sand desert profile that reached the maximum iteration, due to their very complex waveform shape.

In current iteration, if the convergence state of the footprint under examination is achieved, then this footprint will be skipped in all the following iterations and the probabilities of the peaks on this waveform will retain the same.
Therefore, the whole iterative computation process can be summarized as Algorithm 1 below.

**Algorithm 1**

For iteration $r = 1, ..., R$
For the central footprint $s = 1, ..., S$
Check the convergence state of footprint $s$
If not satisfied
Update the probabilities of each peak on waveform $s$
Else
Move to the next footprint $s + 1$
End
End

After the iteration process ceases, the peak with the highest probability on each waveform will be used as the signal peak to determine the reference point ($t_{new}$) and hence to calculate the corresponding surface elevation using the following equation:

$$h_{new} = h_{default} + (t_{default} - t_{new}) \times \frac{c}{2} + \Delta h_{GC} \quad (11)$$

where $h_{new}$ is the elevation value obtained from the reference point determined by our PR method, $h_{default}$ is the surface elevation given in ICESat/GLAS products, $t_{default}$ is the time reference point given by either centroid or MAP retracking method, $t_{new}$ is the reference point (midpoint of the signal peak) derived from our new PR method, $c$ is the speed of light, and $\Delta h_{GC}$ is the Gaussian-centroid (G-C) offset ($i_GmC$ in GLA12/14 product) on the transmitted waveform. The inclusion of this term aims to maintain the consistency of the types of reference points on both transmitted and returned waveforms and can correct a bias up to 6 cm in the elevation measurement [39].

In our PR retracking method, a 1-D neighborhood needs to be defined for the utilization of the spatial contextual information. For our case studies, we adopted a three-footprint or five-footprint neighborhood. For a five-footprint neighborhood, when a footprint is processed, the two adjacent footprints before it and the two adjacent footprints after it along the satellite track are used to define the size of the neighborhood. The immediate neighbor footprint received a weight twice as large as that of the distal neighbor footprint in the iterative posterior probability computation. For the same type of terrain surface, the smoothness and continuity condition may be easily satisfied for a geographically small neighborhood but may not be fulfilled for a large neighborhood. Therefore, an appropriate choice of the neighborhood size based on the surface smoothness is important for applying our retracking method. A larger neighborhood usually facilitates faster convergence of the iteration, due to the inclusion of more conditional information from neighboring footprints. But at the same time, it poses higher requirement of surface smoothness. The three-footprint and the five-footprint neighborhood used here were determined by considering both the convergence of iteration and the requirement of surface roughness. The footprint size of ICESat/GLAS laser pulse is about 70 m in diameter, and the spacing between adjacent footprints is

about 170 m. In our case studies, we adopted a five-footprint neighborhood (~680 m) for relatively smooth lake surface and the Greenland Ice Sheet and a three-footprint neighborhood (~340 m) for relatively bumpy tundra and desert surfaces.

The incorporation of the spatial contextual information from footprints within the neighborhood is only applicable to the types of terrain surfaces that are relatively smooth and have gradually changing elevation. If the majority of the neighboring footprints were extremely contaminated by atmospheric factors or there is abrupt and dramatic elevation change (e.g., cliffs, trees) within the neighborhood, the spatial contextual information given by the neighboring footprints would not be useful and hence our PR method would not be applicable.

To ensure the applicability of our PR method, we first identify the extremely contaminated footprints along the satellite track and then check the surface continuity and smoothness condition before applying our PR method. We utilize the reflectivity ($i$-reflectUC) and the gain value ($i$-gval_rcv) in GLA12/14 data product (https://nsidc.org/data/glas/data-dictionary-glah014) and the signal-to-noise ratio (SNR) as indicator variables to detect extremely contaminated footprints. The SNR value for each footprint is calculated based on the maximum intensity of the waveform and the mean background noise following the method in [12]. The footprints with ($i$-reflectUC < 0.05 && $i$-gval_rcv >= 250 && SNR < 12) are considered extremely contaminated. They are not used as neighboring footprints for the central footprint under processing. The surface smoothness condition is examined through checking the topographic slope on the two sides of the central footprint. The slope on each side of the central footprint is calculated based on the elevation value of each footprint determined by the highest peak in (3). A line is fit on each side for the neighboring footprints and the central footprint to compute the slope. The neighboring footprints will be incorporated in our PR method, only if the slope on that side is less than 0.25°. If the slope on both sides is larger than 0.25°, the terrain surface is considered unsatisfying the smoothness and continuity condition. In the case that a central footprint under processing has no neighboring footprint, the elevation measurement determined by the highest peak in (3) and (4) will be adopted.

## V. APPLICATION EXAMPLES

We applied our PR retracking method to four different types of terrestrial surfaces including frozen lakes, polar tundra, ice sheet, and sand deserts. The locations of ICESat/GLAS ground tracks used in our application examples are shown in Fig. 3. The ground track number, geographic coordinates, and the acquisition dates of these data are shown in Table I. The surface elevation profiles collected by ICESat/GLAS on these four dates were all affected by bad weather and were processed by our new method. For each application site, we compare the topographic profile derived from our PR retracking method with the profiles from original NASA ICESat/GLAS data product and from the high-resolution reference elevation data to demonstrate the effectiveness of our method.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Locations of the four application sites.</figcaption></figure>

<table><caption>TABLE I
GEOGRAPHIC COORDINATES AND ACQUISITION DATES
OF THE ICESAT/GLAS TRACKS IN APPLICATIONS</caption><thead><tr><th>Site</th><th>Surface</th><th>Lon<sup>a</sup></th><th>Lat<sup>a</sup></th><th>Track</th><th>Acquisition Date<sup>b</sup></th></tr></thead><tbody><tr><td>1</td><td>Arctic lake</td><td>-153.84°</td><td>70.72°</td><td>0134</td><td>11/04/2005</td></tr><tr><td>2</td><td>Arctic tundra</td><td>-153.35°</td><td>70.89°</td><td>0334</td><td>03/16/2004</td></tr><tr><td>3</td><td>Greenland Ice Sheet</td><td>-45.01°</td><td>64.79°</td><td>0174</td><td>11/07/2005</td></tr><tr><td>4</td><td>Xinjiang desert</td><td>86.14°</td><td>40.13°</td><td>1128</td><td>10/06/2003</td></tr></tbody></table>

## A. Snow Surface of Lake Teshekpuk in Arctic Coastal Plain

First, we evaluated the applicability of our PR method to the frozen lake surface. A lake normally contains liquid water, when the temperature is higher than the freezing point. It is covered by solid ice, when the temperature is below the freezing point. It may also be covered by snow when snow falls on the frozen ice surface. The water level of a lake may vary at seasonal and inter-annual time scales in response to hydrological processes (e.g., rainfall, evaporation, runoff, and exchanges with groundwater). No matter what the phase of a lake surface is (water, ice, or snow), it is in general very flat, and there is no drastic (e.g., several meters) surface elevation changes within the lake extent. Therefore, the prospect of applying our method to lake surfaces is good for all seasons.

Our case study site is Lake Teshekpuk, located on the Arctic Coastal Plain of northern Alaska (see Fig. 3). It is the largest lake in this region with an area over 800 km². We examined ICESat/GLAS elevation measurements along the satellite track 0134 (segment A-B in Fig. 3) at different dates from the ICESat/GLAS repeat campaigns during 2003 and 2009. Because this site was classified as part of terrestrial land, in the creation of ICESat/GLAS data products, the elevations measurements on the lake surface were retrieved using the standard centroid retracking method. We plotted several elevation profiles acquired at different dates, including March 2 and October 18 in 2004, November 4, 2005, and March 2 and November 18 in 2008. On all these five dates, the lake surface was frozen up and covered by snow. We observed that the elevation profiles acquired on March 2 and October 18 in 2004

<figure><img src="image_7.png" /><figcaption>Fig. 4. Elevation profiles of the frozen snow surface over Lake Teshekpuk on the Arctic coastal plain along the transect from footprint A to footprint B. (a) Problematic profile collected at November 4, 2005 and the normal profiles collected at other dates. (b) Profiles are retrieved using the standard centroid method, the PR methods, and the highest peak given in (3), respectively.</figcaption></figure>

and March 2 and November 18 in 2008 are very similar and consistent, and meet our expectation for a flat lake surface. However, the profile (the blue line) acquired on November 4, 2005 is very different from other repeat profiles [see Fig. 4(a)]. The drastic elevation variation up to 4 m along the track was counter to our expectation of a relatively smooth surface of the lake. This raised our suspicion that the ICESat/GLAS observations on that date were contaminated and that the contaminated waveforms failed the default centroid retracking method, leading to the spurious elevation measurements.

We plotted and examined the waveforms of the footprints along the track AB (see Fig. 3) acquired on November 4, 2005. We found out that these waveforms are characterized with multiple spurious peaks and a long tail. The shapes of the waveforms for five footprints in a neighborhood are plotted in Fig. 5(a). These five footprints were highlighted in different colors in Fig. 4(b), and the elevation measurements given by the standard ICESat/GLAS data product are apparently spurious and highly likely to be erroneous by visually inspecting the topographic profile. As the first step to tackle the spurious measurements produced by the centroid retracking method, we retrieved the elevation measurements using the highest peak given by (3) and (4). After this retrieval, the topographic profile becomes much more reasonable than that of the original ICESat/GLAS data product, as shown by the purple triangles in Fig. 4(b).

Furthermore, we applied our PR retracking method to enhance the computational results using a five-footprint neighborhood. First, ten footprints in this profile were identified as extremely contaminated footprints and were not included as eligible neighborhood footprints for the PR computation.

<footnote><sup>a</sup>Geographic coordinates for the northernmost footprint of each track.</footnote>

<footnote><sup>b</sup>The date format is MM/DD/YYYY.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 5. Illustration of the iteration process. (a) Waveforms within the neighborhood highlighted in Fig. 4. (b) Contaminated waveform and its decomposed Gaussian peaks. (c) Probability of each peak of being the signal peak at each iteration. (d) Average probability change at each iteration. The black vertical dashed line in (a) and (b) indicates the position on time axis at which the dominant peaks of the four neighboring waveforms line up. The retracking point given by our PR and by the centroid method have been denoted in (b).</figcaption></figure>

Our new retracking method and its iteration process is shown in Fig. 5. The waveform of the central footprint [red diamond in Fig. 4(b)] was under examination, with reference to the four footprints [gray diamond in Fig. 4(b)] in its neighborhood. Fig. 5(a) shows the waveforms of the footprints within this neighborhood. The waveform of the central footprint was decomposed into eight separate Gaussian peaks. As shown in Fig. 5(b), Peak 2 is not consistent with the dominant waveform peaks of the four neighboring footprints.

To start our iterative process, we first assigned the initial probabilities for eight separate peaks according to (5). Since Peak 1 has a lower amplitude than Peak 2, its initial prior probability of being the signal peak is lower than Peak 2. According to (7), we computed the compatibilities between the eight peaks of the central footprint and all the peaks of its four neighboring footprints. Since Peak 1 is closer to the dominant signal peaks of four neighboring waveforms than Peak 2 and other six peaks, the compatibility of Peak 1 with four neighboring footprints is the highest among the eight peaks of the central waveform. Therefore, the posterior probability of Peak 1 calculated according to (9) will increase. As shown in Fig. 5(c), after four iterations, the posterior probability of Peak 1 exceeds that of Peak 2 as well as other six peaks. As shown in Fig. 5(d), the process converges quickly. The average change in the posterior probabilities of these peaks calculated at tenth iteration compared with those at the ninth iteration, defined in (10), is less than the predefined threshold value (0.005). Peak 1 has the largest posterior probability when the convergence criteria for the iteration process are satisfied; therefore, it was used as the reference point to calculate the elevation for this central footprint in our PR retracking method. Apparently, our new retracking method has effectively identified the signal peak, Peak 1,

by incorporating the contextual information from neighboring footprints and increased the elevation measurement for central footprint by 2.35 m.

After applying our new PR retracking method to all footprints along the satellite track, some of the erroneous measurements that are not consistent with neighboring footprints are corrected using the spatial contextual information. The new elevation profile [the green line in Fig. 4(b)] derived from our PR method is more consistent with the profiles acquired on other dates in normal condition. As shown in Fig. 4, the elevation measurements from original ICESat/GLAS data products on November 18, 2008 with the default centroid retracking method are consistently lower than the new profile with our probabilistic retracking method.

Since there is no high-resolution reference data available, we use the original ICESat/GLAS data acquired on October 18, 2004 as the reference data for our assessment of our results in comparison with original ICESat/GLAS data product on November 4, 2005. This is because the profile acquired on October 18, 2004 is consistent with that acquired on other three dates (March 2 in 2004, March 2 and November 18 in 2008) and plausibly correct for the lake surface. Using the ICESat/GLAS data acquired on October 18, 2004 as the reference, the root mean squares error (RMSE) of the elevation measurements derived from our PR retracking method is 0.17 m, much smaller than that (0.85 m) of original ICESat/GLAS elevation measurements acquired on November 4, 2005 (see Table II). The standard deviation of the elevation measurements along the satellite track AB has been reduced from 0.81 m for original ICESat/GLAS data products to 0.18 m for new elevation retrievals from our PR method, indicating a significant improvement in the consistency and accuracy of the elevation measurements over the frozen lake snow surface.

## B. Tundra Surface in the Arctic Coastal Plain

The application site for tundra surface is also located on the Arctic Coastal Plain of northern Alaska (see Fig. 3), near Lake Teshekpuk. The Arctic tundra is underlain by thick permafrost, and there are two main seasons, a long winter and a brief summer. During the summer, the top layer of seasonally frozen soil melts, and the tundra is covered by marshes, lakes, bogs streams, and low growing tundra plants, such as sedges, grasses, dwarf willows, moss, crowberry, black bearberry, and lichen [40]. In most cases, these plants can only grow up to several centimeters above ground [41], e.g., dwarf willows, moss, black bearberry, and so on. During the winter, it is very cold and dark, and covered by snow. The ICESat/GLAS track 0334 (segment C-D in Fig. 3) traverses several creeks and vegetated tundra. The ICESat/GLAS data under examination were acquired on March 16, 2004, and at that time the site was covered by snow. Although the tundra surface here is expected to be more bumpy and rougher than the snow-covered lake surface, it still represents a relatively low-relief terrain with homogeneous surface, particularly when it is covered by snow in winter.
<table><caption>TABLE II
EVALUATION OF THE PR DERIVED MEASUREMENTS VERSUS ORIGINAL ICESAT/GLAS DATA PRODUCTS WITH REFERENCE TO HIGH-RESOLUTION TOPOGRAPHIC DATA</caption><thead><tr><th rowspan="2">Site</th><th rowspan="2">Surface</th><th rowspan="2">Number of Footprints</th><th colspan="4">Original ICESat/GLAS Products (m)</th><th colspan="4">Probabilistic Relaxation Results (m)</th></tr><tr><th>Max Diff</th><th>Mean Diff</th><th>RMSE</th><th>r<sup>a</sup></th><th>Max Diff</th><th>Mean Diff</th><th>RMSE</th><th>r<sup>a</sup></th></tr></thead><tbody><tr><td>1</td><td>Arctic lake</td><td>89</td><td>3.79</td><td>1.49</td><td>0.85</td><td></td><td>0.63</td><td>0.18</td><td>0.17</td><td></td></tr><tr><td>2</td><td>Polar tundra</td><td>69</td><td>4.02</td><td>2.03</td><td>0.81</td><td>0.94</td><td>0.95</td><td>0.25</td><td>0.23</td><td>0.99</td></tr><tr><td>3</td><td>Ice sheet</td><td>73</td><td>5.34</td><td>1.04</td><td>1.25</td><td>0.98</td><td>0.73</td><td>0.19</td><td>0.34</td><td>0.99</td></tr><tr><td>4</td><td>Sand desert</td><td>94</td><td>8.14</td><td>2.22</td><td>2.48</td><td>0.99</td><td>8.05</td><td>1.78</td><td>2.34</td><td>0.99</td></tr></tbody></table>

<sup>a</sup>r is the Pearson's correlation coefficient between the ICESat/GLAS data and the reference DEM data.

Max Diff and Mean Diff are calculated by subtracting original ICESat/GLAS measurement or our probabilistic relaxation results from the reference topographic data. For Arctic Lake case, the ICESat/GLAS data acquired on October 18, 2004 (similar season) was used as the reference for the comparison.

<figure><img src="image_9.png" /><figcaption>Fig. 6. Comparison of the PR derived profiles with the original ICESat/GLAS data product and the high-resolution reference data. (a) Over tundra surface in the Arctic coastal plain. (b) Over ice-sheet surface of Greenland. (c) Over sand desert surface of Xinjiang. The footprints barely affected by atmospheric events are denoted for the ice sheet and the desert profiles.</figcaption></figure>

For the ICESat/GLAS track segment C-D, the elevation measurements, which were retrieved by the standard centroid retracking method as this site was classified as land in the original ICESat/GLAS data product, were plotted as a topographic profile (in blue) in Fig. 6(a). This original ICESat/GLAS profile was compared with the reference profile from the 5-m spatial resolution Arctic DEM (http://pgc.umn.edu/arcticdem), which was created using high-resolution satellite stereo images

acquired during 2013 and 2015. A time difference exists between the acquisition of the ICESat/GLAS data and the Arctic DEM data. Since this remote tundra site is situated far away from human activities, the general topographic undulation is very similar and without significant change during the time period. The gray line in Fig. 6(a) is the elevation profile derived from the interferometric synthetic aperture radar (IFSAR) Alaska digital surface model (DSM) which was collected on July 24, 2002 and can be downloaded from the Division of Geological and Geophysical Surveys in the State of Alaska (https://elevation.alaska.gov/). The spatial resolution of IFSAR DSM is 5 m. A comparison between the IFSAR and the Arctic DEM profiles indicates that the general topographic undulation has no significant change from 2002 to 2013 at this site. Seasonal surface elevation change could occur due to the vegetation growth and the snow accumulation. The vegetation height in summer and the snow accumulation in winter in this region are averaged 17.3 ± 8.1 cm [41] and 32.5 ± 8.8 cm [42], respectively. The rate of surface uplift due to glacial isostatic adjustment is about 0.01 m per year [43].

Since the IFSAR DSM profile is much noisier than the Arctic DEM profile, we used the latter one for the validation. The difference between the original ICESat/GLAS data product and the Arctic DEM along the satellite track segment C-D ranges from 0.5 to 4.02 m, with an RMSE of 0.81 m (see Table II). It appears that the original ICESat/GLAS measurements from the centroid retracking method significantly underestimated the surface elevation, in comparison with the high-resolution Arctic DEM as the reference. We applied our PR retracking method to this portion of ICESat/GLAS track, namely, decomposing each waveform into separate peaks and then incorporating the spatial contextual information through an iterative process to find the peak with the maximum posterior probability for elevation retrieval. The elevation profile derived from our PR retracking method is in much better agreement with the reference Arctic DEM than the original ICESat/GLAS data product. Pearson's correlation coefficient *r* between our elevation retrievals and the measurements of the Arctic DEM is as high as 0.98, much larger than that (0.88) of the original ICESat/GLAS data product [see Fig. 7(a)]. The RMSE of our elevation retrievals with the reference to the Arctic DEM is reduced to 0.23 m, a significant improvement over the original NASA ICESat/GLAS data product with an RMSE of 0.81 m (see Table II).
<figure><img src="image_10.png" /><figcaption>Fig. 7. Scatter plots of the original ICESat/GLAS data products and the measurements derived from the PR retracking method with the reference elevation data. PR in the graph is the abbreviation of PR. (a) Over tundra surface in the Arctic coastal plain. (b) Over ice-sheet surface of Greenland. (c) Over sand desert surface of Xinjiang.</figcaption></figure>

### C. Ice Sheet Surface of Greenland

The application site for the ice sheet surface is in southern Greenland (see Fig. 3). Greenland is the second largest mass of glacial ice in the world after the Antarctic Ice Sheet. The thickness of Greenland Ice Sheet is generally more than 2 km and over 3 km at its thickest point [44]. Our application site is located near the central region of southern Greenland. The distance from this site to the eastern, the western, and the southern coast is about 250, 300, and 500 km. The ice sheet surface at our application site is generally smooth with a relatively small surface slope.

The elevation measurements from the original ICESat/GLAS data product were plotted as a topographic profile along the satellite track 0174 (segment E-F in Fig. 3), which was acquired on November 7, 2005. The ICESat/GLAS data product for this site was created using the standard MAP retracking method. Although the profile is smooth in general, there are some spurious undulations [see Fig. 6(b)]. Both MAP- and PR-derived profiles are compared with the high-resolution reference topographic data (green line), which was acquired by an airborne wide-swath laser altimeter system LVIS through the IceBridge program on May 2, 2012. The LVIS measurements within 70 m of each ICESat/GLAS footprint were averaged. Then, the mean value was compared with the ICESat/GLAS measurement of this footprint given by the two methods. There is a 6.5-year span between the acquisition dates of the ICESat/GLAS data and IceBridge LVIS data. However, the surface elevation change in this region is very small. As reported in a previous study, the ice sheet in the southern central Greenland thickens by 0.05 ± 0.02 m each year [45]. The surface uplift due to glacial isostatic adjustment is 0.012 m every year [46]. The insignificant change in this region can also be observed through the comparison between the IceBridge profile [the green line in Fig. 6(b)] and the Arctic DEM profile [the gray line in Fig. 6(b)]. The Arctic DEM profile was generated using the WorldView-3 stereo images acquired on August 29, 2016, which was four years later than the acquisition of the IceBridge profile. These temporal changes in the ice sheet elevation are negligible as compared with the magnitude of errors observed in this profile below.

As shown in Table II, the RMSE of the original ICESat/GLAS data product is 1.25 m, with a maximum error up to 5.34 m. After applying our new retracking method, the derived elevation measurements (red) are almost the same as the reference IceBridge LVIS measurements. The RMSE of the elevation measurements from our new retracking method is reduced to 0.34 m from 1.25 m, and the Pearson correlation coefficient *r* is as high as 0.99 as compared with the 0.98 of the original ICESat/GLAS product. This application example indicates that our new PR method is capable of producing more accurate and reliable elevation measurements than the standard MAP retracking method.

### D. Sand Dune Surface in an Arid Desert

The application site is located in the northeast Taklimakan Desert (see Fig. 3), the largest desert in China. It is made up of continuous sand dunes, which are usually over 100 m high and some even higher than 300 m. In this arid region, a thin layer of blowing particles near the surface due to winds often occurs, which may lead to the forward scattering effect on the ICESat/GLAS waveforms. The ICESat/GLAS satellite track 1128 (segment G-H in Fig. 3) extends across four major sand dunes from north to south, as shown in Fig. 6(c). The original ICESat/GLAS data along this track segment under examination were acquired on October 6, 2003. We use the 30-m SRTM DEM data collected on February 11, 2000 as the reference. The northeast Taklimakan desert is an area of low wind energy [47].

Even though there might be changes for the small sand dunes, comparison of the Landsat images acquired on February 6, 2000 and October 4, 2003 indicates that there was no clear horizontal shift in the major dunes in our application site.

The original elevation measurements from ICESat/GLAS data products for this site were retrieved by the standard centroid retracking method. Although the original ICESat/GLAS elevation measurements correlate well with the SRTM data in general, the elevation measurements at some points are clearly underestimated (e.g., footprint 57) as shown by the enlarged inset in Fig. 6(c). Fig. 8(a) shows the waveforms of
<figure><img src="image_11.png" /><figcaption>Fig. 8. Waveforms of the footprints highlighted in the inset of Fig. 6(c). (a) Waveforms of the three footprints within the neighborhood. (b) Reference points determined by the centroid and the PR methods for the waveform of footprint 57. The vertical dashed line in (a) indicates the position on time axis at which the dominant peaks of the three waveforms line up. The upper and the lower horizontal axes in (a) and (b) represent the location of waveform bins in terms of distance and time, respectively.</figcaption></figure>

the three adjacent footprints 55, 56, and 57 in Fig. 6(c) inset. These waveforms line up with each other as indicated by the vertical dashed line in Fig. 8(a). In Fig. 8(b), the waveform of footprint 57 is decomposed into six peaks. The first red peak is identified by our PR method as the signal-peak for elevation retrieval. This elevation value is more consistent than the original elevation measurement with the SRTM data. It improves the elevation measurement of footprint 57 by 3.63 m as shown by the inset in Fig. 6(c). Applying our new PR retracking method, we derived the elevation measurements for all footprints along the satellite track segment G-H in Fig. 3. As shown in Fig. 6(c), the elevation measurements from our method are in a good agreement with the reference SRTM data with Pearson's correlation coefficient of 0.99. The RMSE of our new elevation measurements is reduced to 2.34 m, in comparison with that (2.48 m) of the original ICESat/GLAS data product.

## VI. DISCUSSION

In the creation of ICESat/GLAS data products, two standard retracking methods were used to process the laser waveforms returned from Earth's surface for the elevation retrievals, the centroid retracking method for the complex terrestrial surfaces, and the MAP retracking method for smooth and homogeneous surfaces (sea ice, ocean water, and ice sheet). Our analysis suggests that a more detailed terrestrial land cover classification should be adopted to distinguish relatively smooth and gradually changing land surfaces from complex land surfaces for the selection of an appropriate waveform retracking method for elevation retrievals. In particular, when the returned waveform is contaminated by bad weather conditions through multiple particle scattering, both methods tend to generate erroneous measurements that are significantly lower than true surface elevation.

For regions with more frequent clouds, fog, blowing snow, or sand dust, the laser waveforms are more prone to atmospheric contamination. Clouds prevail in the Arctic over half of the year [24], [48]. Snow is the primary feature of the Arctic, covering land surfaces for about eight–ten months of a year. The frequent and strong winds in this region can disturb the loose and recent snow and carry it long distances

across the low-relief tundra or ice sheet surfaces [25], forming a blowing snow layer. As reported, over half the wind events in the Arctic region have speed greater than 5 m/s, and the highest speed could reach to 25 m/s [49]. Blowing snow has been demonstrated to have a more significant influence than cloud cover on the estimation of surface elevation [20], [22]. In our application examples, the measurement error of the original ICESat/GLAS data products can be as large as about 4 m for the Arctic frozen lake, about 5 m for the Greenland Ice Sheet, and about 4 m in the Arctic tundra. In the desert areas, blowing sand and dusts often occur due to lack of vegetation cover. The dust particles may cause the forward scattering effect on the laser waveforms, leading to the failure of the standard centroid retracking method. In our desert application example, the error of the original ICESat/GLAS data product can be as large as 8.56 m. Therefore, the regions prone to clouds, ice fog, and blowing snow and dust may require a more robust retracking method that is resistant to atmospheric contamination.

Our PR retracking method combines spatial contextual information from adjacent footprints along the satellite track with the waveform shape information in the determination of the reference point for the elevation retrieval. This is in contrast to the standard centroid and MAP retracking methods that solely rely on the waveform shape information. Therefore, our method is more resistant to the atmospheric forward scattering effect and can more reliably identify the signal peak on the waveform for elevation determination.

Our PR retracking method is an iterative computational process. Its first step is to assign the initial prior probabilities to the decomposed peaks in terms of their peak amplitudes. The highest peak has the largest prior probability to be the signal peak. Clearly, the first step of our PR method is virtually the same as the standard MAP retracking method. The difference of our retracking method from MAP retracking method is that we select the peak with the maximum posterior probability as the signal peak after the spatial contextual information is iteratively incorporated through evaluating the neighborhood compatibility, instead of simply selecting the peak with the maximum amplitude in the MAP retracking method. From this perspective, PR retracking method extends the MAP retracking method and is less susceptible to the non-Gaussian waveform caused by the atmospheric forward scattering. It should be noted that our new method does not decrease the along-track resolution since laser footprints are processed one-by-one sequentially in the probabilistic sense along the satellite track and no weighted averaging of neighboring footprints is involved.

As demonstrated in the four application examples, our new retracking method is applicable to different types of terrain surfaces, including the ice sheets, inland lake surfaces, tundra, and sand deserts. With independent/separate reference data sets, we have done a detailed quantitative assessment of the performance of our new method in terms of absolute accuracy (RMSE) and relative accuracy (Pearson's *r*). As reported in Table II and Fig. 7, our new method reduced the RMSE by 80%, 73%, 76%, and 5% for frozen lake, polar tundra, ice sheet, and sand desert, respectively. The significantly increased Pearson's *r* values also indicate the enhanced relative accuracy
and consistency of the elevation measurements derived from our new method as compared with the original NASA ICESat/GLAS data products produced either by the centroid or MAP retracking method. The relatively damped performance of our new method over sand desert could partly attribute to the lower accuracy of the reference SRTM data as compared with the reference data used in other two sites. These four types of terrain surface are abundant and widely distributed on Earth's surface, and the geophysical processes underlying these terrain types are unique. For instance, about 10% of Earth's land surface is covered by ice sheet [50], 3.7% by inland lakes [51], 5.5% by polar tundra [52], and 30% by sand desert [53]. Our new method could be very valuable to derive improved elevation measurements for these surfaces if the footprints were affected by adverse atmospheric events and be beneficial to science and application communities.

The biases between our new method and the two standard retracking methods were evaluated using the footprints along the Greenland Ice Sheet and the sand desert profiles that were barely affected by atmospheric multiple scattering. These footprints were identified through checking the gain and the reflectivity values associated with each footprint and have been denoted on the two profiles in Fig. 6. The bias between our new method and the centroid method is $-7.0 \pm 2.3$ cm. The bias between our new method and the MAP method is $15.6 \pm 9.8$ cm. For these high-quality ICESat footprints, the elevation measurements produced by our new method are slightly lower than the centroid method and slightly higher than the MAP method. For the footprints affected by atmospheric events, the biases (estimated based on high-quality footprints) could be added to or subtracted from the measurements produced by our new method in order to align with the data sets produced by the two standard methods. In such a way, we would be able to obtain improved and consistent ICESat data sets, instead of simply discarding the low-quality footprints. A more comprehensive and rigorous evaluation of the biases between our new method and the two standard methods will be conducted in our future studies.

Although our PR retracking method is more robust than the standard MAP retracking method for some landscapes, it is still not suitable for complex terrain surface with strong topographic relief within the laser footprint or for satellite tracks where all the neighboring footprints are extremely contaminated. Particularly, it is limited for the footprints over complex surface that were simultaneously affected by adverse weather conditions. The spatial contextual information utilized in our retracking method stems from the smooth and gradually varying neighborhood, in which the footprints have similar surface roughness and elevation. The relative smoothness and topographic continuity of the terrain surface are the basic condition for applying our PR retracking method. As demonstrated by application examples, the inland lake surfaces, tundra, ice sheets, and deserts with large sand dunes satisfy the smoothness and continuity condition, making our PR retracking method applicable to such Earth surfaces. Also, it should be noted that if the majority of the neighboring waveforms have been extremely contaminated by the atmospheric forward scattering effect, the spatial contextual information

becomes invalid and unhelpful for identifying the signal peak and hence our retracking method would not work. In our case studies, we used the gain and reflectivity thresholds to identify footprints that were extremely contaminated in the evaluation of the neighborhood compatibility for the iterative calculation of posterior probability.

Our PR retracking method is specially designed to process ICESat/GLAS laser altimetry data. We believe that the basic idea and algorithm can be applied to process the airborne waveform laser and radar altimetry data. Since the footprint of airborne laser is much smaller than ICESat/GLAS laser altimetry system, and the density of footprints is much higher, a much geographically smaller neighborhood can be used to incorporate the spatial contextual information for the improved elevation retrieval. Within a smaller neighborhood, the requirement of surface smoothness and continuity conditions can be easily satisfied. We anticipate that our PR retracking method can be applied to wider, even more complicated terrestrial land surfaces with airborne waveform altimeters. The PR retracking algorithm is likely to provide realistic estimates for the multipeaked echoes sometimes observed over inland water bodies using radar altimetry when classical retrackers such as Offset Center of Gravity (OCOG) [54], commonly used over rivers, lakes, floodplains, and wetlands [55], fail to provide accurate height estimates.

The source code of this algorithm and the waveform data of the four sites in this article are publicly available on GitHub via https://github.com/ssshusong/Probabilistic-Relaxation-for-ICESat.git.

## VII. CONCLUSION

The ICESat/GLAS provided elevation measurements of Earth's surface during 2003 and 2009, which have been widely used in various scientific and practical applications. In the creation of ICESat/GLAS data products, two standard retracking methods have been employed to process the laser waveforms for elevation retrieval. The MAP retracking method was used for ice sheets, oceans, and sea ice, while the centroid retracking method was used for all different types of terrestrial lands on Earth's surface. Both the MAP and centroid retracking methods rely on the shape information of the returned laser waveform. Although these two standard methods work well in general, they may generate erroneous measurements when the returned waveform had multiple spurious peaks and a long tail generated by bad weather (cloud cover, ice fog, blowing snow, and dust storms). Some types of terrestrial land cover, such as inland waters, polar tundra, sand deserts, wetlands, and coastal zones, have relatively smooth and gradually changing surfaces all year-around or seasonally. Our study suggests that these types of land covers should have been identified for the selection of an appropriate waveform retracking method for elevation retrievals. ICESat/GLAS data acquired over these areas for all seasons or for a specific season should have been processed with the MAP retracking method, rather than using the centroid retracking method for all terrestrial lands. The magnitude and frequency of the atmospheric forward scattering effect on the ICESat/GLAS
waveforms may vary from place to place. Our case studies show that the original ICESat/GLAS data products over the frozen Arctic lake, the Arctic tundra, arid sand desert, and the ice sheets may contain erroneous elevation measurements due to the use of standard retracking methods and the forward scattering effect by frequent clouds, ice fogs, blowing snows, and dusts.

We have incorporated spatial contextual information in waveform analysis for the first time as part of a novel PR retracking method. This new method initially assigns a prior probability value for each candidate peak on a returned waveform to be the signal peak in terms of its amplitude, and then iteratively computes the posterior probability by incorporating the spatial contextual information through evaluating each peak's compatibility with the footprints in a predefined neighborhood. After the iterative process is converged, the peak with the maximum posterior probability is selected as the signal peak for the elevation retrieval.

Our application examples demonstrate that our PR retracking method is more resistant to atmospheric contaminations of the laser waveforms and is able to create more reliable and accurate elevation measurements than the original NASA ICESat/GLAS data products for different types of land surfaces, including inland lakes, polar tundra, sand deserts, and ice sheets. Our method is developed specifically for ICESat/GLAS laser waveform altimetry data, and the case studies have been conducted with ICESat/GLAS data. Nevertheless, we believe that the basic idea of incorporating spatial contextual information and the corresponding numerical algorithm can find even wider application in processing the waveform data from airborne laser and radar altimetry systems.

## REFERENCES

[1] F. Frappart et al., “Satellite altimetry: Principles and applications in Earth sciences,” in *Wiley Encyclopedia of Electrical and Electronics Engineering*. Hoboken, NJ, USA: Wiley, 2017.

[2] C. H. Davis, “Satellite radar altimetry,” *IEEE Trans. Microw. Theory Techn.*, vol. 40, no. 6, pp. 1070–1076, Jun. 1992.

[3] L. L. Fu, and A. Cazenave, *Satellite Altimetry and Earth Sciences: A Handbook of Techniques and Applications*. Amsterdam, The Netherlands: Elsevier, 2000.

[4] Y. I. Troitskaya et al., “Satellite altimetry of inland water bodies,” *Water Resour.*, vol. 39, no. 2, pp. 184–199, Mar. 2012.

[5] J. Benveniste, “Radar altimetry: Past, present and future,” in *Coastal Altimetry*. Berlin, Germany: Springer, 2011, pp. 1–17.

[6] H. J. Zwally et al., “ICESat’s laser measurements of polar ice, atmosphere, ocean, and land,” *J. Geodyn.*, vol. 34, nos. 3–4, pp. 405–445, Oct. 2002.

[7] U. C. Herzfeld and B. Wallin, “Spatio-temporal analysis of surface elevation changes in pine island glacier, antarctica, from ICESat GLAS data and ERS-1 radar altimeter data,” *Ann. Glaciology*, vol. 55, no. 66, pp. 248–258, 2014.

[8] R. Kwok, “ICESat observations of arctic sea ice: A first look,” *Geophys. Res. Lett.*, vol. 31, no. 16, 2004, Art. no. L16401.

[9] M. Keller, “Revised method for forest canopy height estimation from geoscience laser altimeter system waveforms,” *J. Appl. Remote Sens.*, vol. 1, no. 1, Sep. 2007, Art. no. 013537.

[10] R. Bindschadler, H. Choi, C. Shuman, and T. Markus, “Detecting and measuring new snow accumulation on ice sheets by satellite remote sensing,” *Remote Sens. Environ.*, vol. 98, no. 4, pp. 388–402, Oct. 2005.

[11] D. Treichler and A. Kääb, “Snow depth from ICESat laser altimetry— A test study in southern norway,” *Remote Sens. Environ.*, vol. 191, pp. 389–401, Mar. 2017.

[12] S. Shu et al., “Estimation of snow accumulation over frozen arctic lakes using repeat ICESat laser altimetry observations—A case study in northern alaska,” *Remote Sens. Environ.*, vol. 216, pp. 529–543, Oct. 2018.

[13] X. Wang et al., “Water-level changes in China’s large lakes determined from ICESat/GLAS data,” *Remote Sens. Environ.*, vol. 132, pp. 131–144, May 2013.

[14] C. Song, B. Huang, and L. Ke, “Inter-annual changes of Alpine inland lake water storage on the Tibetan plateau: Detection and analysis by integrating satellite altimetry and optical imagery,” *Hydrolog. Processes*, vol. 28, no. 4, pp. 2411–2418, Feb. 2014.

[15] H. A. Fricker, “Assessment of ICESat performance at the Salar de Uyuni, Bolivia,” *Geophys. Res. Lett.*, vol. 32, no. 21, pp. 1–5, 2005.

[16] N. T. Kurtz, T. Markus, D. J. Cavalieri, W. Krabill, J. G. Sonntag, and J. Miller, “Comparison of ICESat data with airborne laser altimeter measurements over Arctic sea ice,” *IEEE Trans. Geosci. Remote Sens.*, vol. 46, no. 7, pp. 1913–1924, Jul. 2008.

[17] C. A. Shuman et al., “ICESat antarctic elevation data: Preliminary precision and accuracy assessment,” *Geophys. Res. Lett.*, vol. 33, no. 7, pp. 1–4, 2006.

[18] F. Pirotti, “IceSAT/GLAS waveform signal processing for ground cover classification: State of the art,” *Italian J. Remote Sens.*, vol. 42, no. 2, pp. 13–26, Jun. 2010.

[19] A. C. Brenner et al., “Derivation of range and range distributions from laser pulse waveform analysis for surface elevations, roughness, slope, and vegetation heights,” *Geosci. Laser Altimeter Syst. (GLAS), Algorithm Theor. Basis Document, Version*, vol. 5, pp. 1–125, Sep. 2011.

[20] Y. Yang, A. Marshak, T. Varnai, W. Wiscombe, and P. Yang, “Uncertainties in ice-sheet altimetry from a spaceborne 1064-nm single-channel lidar due to undetected thin clouds,” *IEEE Trans. Geosci. Remote Sens.*, vol. 48, no. 1, pp. 250–259, Jan. 2010.

[21] D. P. Duda, J. D. Spinhirne, and E. W. Eloranta, “Atmospheric multiple scattering effects on GLAS altimetry. I. Calculations of single pulse bias,” *IEEE Trans. Geosci. Remote Sens.*, vol. 39, no. 1, pp. 92–101, Jan. 2001.

[22] A. Mahesh, J. D. Spinhirne, D. P. Duda, and E. W. Eloranta, “Atmospheric multiple scattering effects on GLAS altimetry. II. Analysis of expected errors in antarctic altitude measurements,” *IEEE Trans. Geosci. Remote Sens.*, vol. 40, no. 11, pp. 2353–2362, Nov. 2002.

[23] S. P. Palm, Y. Yang, J. D. Spinhirne, and A. Marshak, “Satellite remote sensing of blowing snow properties over antarctica,” *J. Geophys. Res.*, vol. 116, no. D16, pp. 1–16, 2011.

[24] D. Yi, H. Jay Zwally, and X. Sun, “ICESat measurement of greenland ice sheet surface slope and roughness,” *Ann. Glaciology*, vol. 42, no. 1, pp. 83–89, 2005.

[25] R. Kwok, G. F. Cunningham, H. J. Zwally, and D. Yi, “ICESat over arctic sea ice: Interpretation of altimetric and reflectivity profiles,” *J. Geophys. Res.*, vol. 111, no. C6, pp. 1–20, 2006.

[26] D. Felikson et al., “Comparison of elevation change detection methods from ICESat altimetry over the Greenland ice sheet,” *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 10, pp. 5494–5505, Oct. 2017.

[27] D. Yi, H. J. Zwally, and J. W. Robbins, “ICESat observations of seasonal and interannual variations of sea-ice freeboard and estimated thickness in the weddell sea, antarctica (2003–2009),” *Ann. Glaciology*, vol. 52, no. 57, pp. 43–51, 2011.

[28] M. R. Siegfried, R. L. Hawley, and J. F. Burkhart, “High-resolution ground-based GPS measurements show intercampaign bias in ICESat elevation data near Summit, Greenland,” *IEEE Trans. Geosci. Remote Sens.*, vol. 49, no. 6, pp. 3393–3400, Sep. 2011.

[29] H. J. Zwally, D. Yi, R. Kwok, and Y. Zhao, “ICESat measurements of sea ice freeboard and estimates of sea ice thickness in the Weddell sea,” *J. Geophys. Res.*, vol. 113, no. C2, pp. 1–17, 2008.

[30] M. D. Shupe et al., “Clouds at arctic atmospheric observatories.—Part I: Occurrence and macrophysical properties,” *J. Appl. Meteorol. Climatol.*, vol. 50, no. 3, pp. 626–644, Mar. 2011.

[31] G. E. Liston and M. Sturm, “Winter precipitation patterns in Arctic Alaska determined from a blowing-snow model and snow-depth observations,” *J. Hydrometeorology*, vol. 3, no. 6, pp. 646–659, Dec. 2002.

[32] J. B. Blair and M. Hofton, “IceBridge LVIS L2 geolocated surface elevation product, version 1,” NASA Nat. Snow Ice Data Center Distrib. Act. Arch. Center, Boulder, CO, USA, 2012, doi: 10.5067/OIKFGJNBM6OO.
[33] K. M. Brunt et al., "Assessment of NASA airborne laser altimetry data using ground-based GPS data near Summit station, Greenland," *Cryosphere*, vol. 11, no. 2, pp. 681-692, Mar. 2017.

[34] Y. Gorokhovich and A. Voustianiouk, "Accuracy assessment of the processed SRTM-based elevation data by CGIAR using field data from USA and thailand and its relation to the terrain characteristics," *Remote Sens. Environ.*, vol. 104, no. 4, pp. 409-415, Oct. 2006.

[35] C. S. Gardner, "Ranging performance of satellite laser altimeters," *IEEE Trans. Geosci. Remote Sens.*, vol. 30, no. 5, pp. 1061-1072, Sep. 1992.

[36] K. M. Hinkel, Z. Lin, Y. Sheng, and E. A. Lyons, "Regional lake ice meltout patterns near Barrow, Alaska," *Polar Geography*, vol. 35, no. 1, pp. 1-18, Mar. 2012.

[37] A. Rosenfeld, R. A. Hummel, and S. W. Zucker, "Scene labeling by relaxation operations," *IEEE Trans. Syst., Man, Cybern.*, vol. SMC-6, no. 6, pp. 420-433, Jun. 1976.

[38] A. Rosenfeld, and A. C. Kak, *Digital Picture Processing*. New York, NY, USA: Academic, 1982.

[39] A. A. Borsa, G. Moholdt, H. A. Fricker, and K. M. Brunt, "A range correction for ICESat and its potential impact on ice-sheet mass balance studies," *Cryosphere*, vol. 8, no. 2, pp. 345-357, Mar. 2014.

[40] C. J. Markon and D. V. Derksen, "Identification of tundra land cover near Teshekpuk lake, Alaska using SPOT satellite data," *Arctic*, vol. 47, no. 3, pp. 222-231, Jan. 1994.

[41] J. C. von Fischer, R. C. Rhew, G. M. Ames, B. K. Fosdick, and P. E. von Fischer, "Vegetation height and other controls of spatial variability in methane emissions from the arctic coastal tundra at Barrow, Alaska," *J. Geophys. Res.*, vol. 115, no. 3, pp. 1-11, Sep. 2010.

[42] M. Sturm and G. E. Liston, "The snow cover on lakes of the Arctic coastal plain of Alaska, U.S.A.," *J. Glaciology*, vol. 49, no. 166, pp. 370-380, 2003.

[43] G. F. Sella et al., "Observation of glacial isostatic adjustment in 'stable' North America with GPS," *Geophys. Res. Lett.*, vol. 34, no. 2, pp. 1-6, 2007.

[44] J. L. Bamber, R. L. Layberry, and S. P. Gogineni, "A new ice thickness and bed data set for the greenland ice sheet: 1. Measurement, data reduction, and errors," *J. Geophys. Res., Atmos.*, vol. 106, no. D24, pp. 33773-33780, Dec. 2001.

[45] K. C. Jezek, "Surface elevation and velocity changes on the South-central Greenland ice sheet: 1980-2011," *J. Glaciology*, vol. 58, no. 212, pp. 1201-1211, 2012.

[46] S. A. Khan et al., "Geodetic measurements reveal similarities between post-last glacial maximum and present-day mass loss from the Greenland ice sheet," *Sci. Adv.*, vol. 2, no. 9, Sep. 2016, Art. no. e1600931.

[47] X. Wang, Z. Dong, J. Zhang, and G. Chen, "Geomorphology of sand dunes in the Northeast Taklimakan desert," *Geomorphology*, vol. 42, nos. 3-4, pp. 183-195, Jan. 2002.

[48] X. Wang and J. R. Key, "Arctic surface, cloud, and radiation properties based on the AVHRR polar pathfinder dataset—Part I: Spatial and temporal characteristics," *J. Climate*, vol. 18, no. 14, pp. 2558-2574, Jul. 2005.

[49] D. Small, E. Atallah, and J. Gyakum, "Wind regimes along the Beaufort sea coast favorable for strong wind events at Tuktoyaktuk," *J. Appl. Meteorol. Climatol.*, vol. 50, no. 6, pp. 1291-1306, Jun. 2011.

[50] R. Massom, and D. Lubin, *Polar Remote Sensing: Ice Sheets*, vol. 2. Chichester, U.K.: Springer, 2006, pp. 1-511.

[51] C. Verpoorter, T. Kutser, D. A. Seekell, and L. J. Tranvik, "A global inventory of lakes based on high-resolution satellite imagery," *Geophys. Res. Lett.*, vol. 41, no. 18, pp. 6396-6402, Sep. 2014.

[52] P. D. Moore, *Tundra*. New York, NY, USA: Infobase Publishing, 2009.

[53] M. Allaby, and R. Garratt, *Deserts*. New York, NY, USA: Infobase Publishing, 2009.

[54] D. J. Wingham, C. G. Rapley, and H. Griffiths, "New techniques in satellite altimeter tracking systems," in *Proc. IGARSS*, Sep. 1986, pp. 1339-1344.

[55] F. Frappart, S. Calmant, M. Cauhope, F. Seyler, and A. Cazenave, "Preliminary results of ENVISAT RA-2-derived water levels validation over the Amazon basin," *Remote Sens. Environ.*, vol. 100, no. 2, pp. 252-264, Jan. 2006.

<figure><img src="image_12.png" /></figure>

**Song Shu** received the B.S. degree in physical geography and the M.S. degree in cartography and geographic information science from East China Normal University, Shanghai, China, in 2010 and 2013, respectively, and the Ph.D. degree in geography and geographic information science from the University of Cincinnati, Cincinnati, OH, USA, in 2019.

He is an Assistant Professor with the Department of Geography and Planning, Appalachian State University, Boone, NC, USA. His research interests include the geographic information science (GIS) and remote sensing applications on lake hydrology, water resources, Arctic snow, cryosphere, global climate change, and urbanization.

<figure><img src="image_13.png" /></figure>

**Hongxing Liu** (Member, IEEE) received the Ph.D. degree in geography from The Ohio State University, Columbus, OH, USA, in 1999.

He is a Professor with the Department of Geography, The University of Alabama, Tuscaloosa, AL, USA. His research interests include remote sensing, geographical information science, hydrological modeling, terrain analysis, coastal mapping, and change studies. His recent research projects have been funded by the National Aeronautics and Space Administration (NASA), the National Science Foundation (NSF), and the National Oceanic and Atmospheric Administration (NOAA) Sea Grant Program.

Dr. Liu has received a number of prestigious awards, including the NASA Group Achievement Award, the Best Paper Award from Computers and Geoscience, and the Distinguished Achievement Award in Teaching from the College of Geoscience, Texas A&M University.

<figure><img src="image_14.png" /></figure>

**Frédéric Frappart** received the master's degree in detection and information system from ENSTA Bretagne (formerly ENSIETA) Engineer's School, Brest, France, in 2001, and the Ph.D. and Habilitation degrees in geophysics from Paul Sabatier University, Toulouse, France, in 2006 and 2015, respectively.

He is an Associate Professor at Observatoire Midi-Pyrénées, Federal University of Toulouse Midi-Pyrénées, Toulouse. His research interest includes the monitoring of the terrestrial water cycle using remote sensing.

<figure><img src="image_15.png" /></figure>

**Emily Lei Kang** received the B.S. degree in applied mathematics from Tianjin University, Tianjin, China, in 2004, the B.A. degree in finance from Nankai University, Tianjin, China, in 2004, and the M.S. and Ph.D. degrees in statistics from The Ohio State University, Columbus, OH, USA, in 2006 and 2009, respectively.

From 2009 to 2011, she was a Post-Doctoral Fellow with the Statistical and Applied Mathematical Sciences Institute, Research Triangle Park, NC, USA, and North Carolina State University, Raleigh, NC, USA. She is an Associate Professor with the Division of Statistics and Data Science and the Department of Mathematical Sciences, University of Cincinnati, Cincinnati, OH, USA. Her research interests include statistical methodology and algorithms for large data sets with complex dependence structures, uncertainty quantification, and their applications in remote sensing, climate, and environmental sciences.
<figure><img src="image_16.png" /></figure>

**Bo Yang** received the B.S. degree in applied mathematics from Shaanxi Normal University, Shaanxi, China, in 2008, and the M.S. degree in computer science from Capital Normal University, Beijing, China, in 2011.

He is a Post-Doctoral Research Scholar with the University of Central Florida (UCF), Orlando, OH, USA. His research interests include spatio-temporal modeling in geographic information science (GIS), unmanned aerial vehicle (UAV) and satellite remote sensing, spatial statistics, environmental, and sociological modeling.

<figure><img src="image_17.png" /></figure>

**Bailang Yu** (Senior Member, IEEE) received the B.S. and Ph.D. degrees in cartography and geographic information systems from East China Normal University, Shanghai, China, in 2002 and 2009, respectively.

He is a Professor with the Key Laboratory of Geographic Information Science, Ministry of Education, East China Normal University, where he is also with the School of Geographic Sciences. His research interests include urban remote sensing, nighttime light remote sensing, LiDAR, and object-based methods.

<figure><img src="image_18.png" /></figure>

**Min Xu** received the B.S. degree in computer science and technology from Central South University, Changsha, Hunan, China, in 2013, the M.A. degree in geography from the University of Cincinnati, Cincinnati, OH, USA, in 2015, and the Ph.D. degree in geography from The University of Alabama, Tuscaloosa, AL, USA, in 2020.

She is a Post-Doctoral Scholar with the College of Marine Science, University of South Florida at St. Petersburg, St. Petersburg, FL, USA. Her research interests include remote sensing, water quality, harmful algal blooms, and spatio-temporal analysis of remote sensing images.

<figure><img src="image_19.png" /></figure>

**Shujie Wang** received the B.S. degree in geographic information system and the M.A degree in cartography and geographic information science from Sun Yat-sen University, Guangzhou, China, in 2010 and 2012, respectively, and the Ph.D. degree in geography and geographic information science from the University of Cincinnati, Cincinnati, OH, USA, in 2018.

She is an Assistant Professor with the Department of Geography, The Pennsylvania State University, University Park, PA, USA. Her research interests include remote sensing, machine learning, numerical modeling, ice sheet dynamics, snow/ice albedo, and glacial microbes.

<figure><img src="image_20.png" /></figure>

**Yan Huang** received the B.S. degree in geographical information system from the Chengdu University of Technology, Chengdu, China, in 2010, and the Ph.D. degree in cartography and geographic information systems from East China Normal University, Shanghai, China, in 2016.

She is an Associate Professor with the Key Laboratory of Geographic Information Science, Ministry of Education, East China Normal University, where she is also with the School of Geographic Sciences. Her research interests include remote sensing of the cryosphere and spatio-temporal modeling of remote sensing images.

<figure><img src="image_21.png" /></figure>

**Richard Beck** received the Ph.D. degree in earth sciences from the University of Southern California, Los Angeles, CA, USA, in 1995.

He is an Associate Professor of geography and geographic information science with the University of Cincinnati, Cincinnati, OH, USA. His research interests include geographic information networks, geographical information science (GIS), remote sensing, climate change, South Asia, and water resources.

<figure><img src="image_22.png" /></figure>

**Bin Wu** received the Ph.D. degree in cartography and geographic information systems from East China Normal University, Shanghai, China, in 2018.

He is a Post-Doctoral Researcher with the Key Laboratory of Geographic Information Science, Ministry of Education, and the School of Geographic Sciences, East China Normal University. His research interests include urban remote sensing, LiDAR, and spatio-temporal analysis.

<figure><img src="image_23.png" /></figure>

**Kenneth Hinkel** received the B.S. degree in physical geography from Michigan State University, East Lansing, MI, USA, in 1976, and the M.S. degree in physical geography and the Ph.D. degree in geology from the University of Michigan, Ann Arbor, MI, USA, in 1982 and 1986, respectively.

He is currently a Professor Emeritus and a McMicken Scholar with the University of Cincinnati, Cincinnati, OH, USA, and a Research Professor in geological sciences and mining engineering and sciences with Michigan Technological University, Houghton, MI, USA. His research interests include climatic and human impacts on permafrost stability and Arctic lake hydrology.