A Pseudo-Waveform-Based Method for Grading ICESat-2 ATL08 Terrain Estimates in Forested Areas

Abstract

The ICESat-2 Land and Vegetation Height (ATL08) product is a new control point dataset for large-scale topographic mapping and geodetic surveying. However, its elevation accuracy is typically affected by multiple factors. The study aims to propose a new approach to classify ATL08 terrain estimates into different accuracy levels and extract reliable ground control points (GCPs) from ICESat-2 ATL08. Specifically, the methodology is divided into three stages. First, the ATL08 terrain estimates are matched with the raw ATL03 photon cloud data, and the ATL08 terrain estimates are used to fit a continuous terrain curve. Then, using the fitted continuous terrain curve and raw ATL03 photon cloud data, a pseudo-waveform is generated for grading the ATL08 terrain estimates. Finally, all the ATL08 terrain estimates are graded based on the peak characteristics of the generated pseudo-waveform. To validate the feasibility of the proposed method, four study areas from the National Ecological Observatory Network (NEON), characterized by various terrain features and forest types were selected. High-accuracy airborne lidar data were used to evaluate the accuracy of graded ICESat-2 terrain estimates. The results demonstrate that the method effectively classified all ATL08 terrain estimates into different accuracy levels and successfully extracted high-accuracy GCPs. The root mean square errors (RMSEs) of the first accuracy level in the four selected study areas were 0.99 m, 0.51 m, 1.88 m, and 0.65 m, representing accuracy improvement of 51.7%, 58.2%, 83.1%, and 68.8%, respectively, compared to the original ATL08 terrain estimates before classifying. Additionally, a comparison with the conventional threshold-based GCP extraction method demonstrated the superior performance of our proposed approach. This study introduces a new approach to extract high-quality elevation control points from ICESat-2 ATL08 data, particularly in forested areas.

1. Introduction

A digital elevation model (DEM) represents the bare Earth’s shape and serves as a crucial dataset for numerous applications, such as forest inventory establishment [1], climate change monitoring [2,3], and hydrology analysis [4,5,6]. Interferometric synthetic aperture radar (InSAR) is one of the most widely employed technologies for producing DEMs locally and globally, owing to its all-weather imaging capacity, extensive coverage, high resolution, and precise measurement capacity [7,8,9]. However, to calibrate systematic height errors inherent in interferometric SAR processing (e.g., unavoidable baseline errors) and enhance the accuracy of final DEM products, high-quality ground control points (GCPs) are essential. In addition, GCPs play a crucial role in calibrating translational errors and aligning elevations in DEMs derived from multiple images of the regions of interest [10,11,12], which ensures reliable and consistent topographic mapping using InSAR technologies.

To acquire high-quality GCPs, traditional approaches through direct and human measurement are expensive, time-consuming, and labor-intensive, which makes it impractical for large-scale applications and hinders timely updates. To address these challenges, a global positioning system (GPS) appears as a potential solution for GCP acquisition [13]. While GPS can provide observables at high temporal resolution over a relatively large scale with lower timing costs, its application is constrained by the high cost of equipment and challenges in deploying receivers in remote areas [10,14]. Moreover, signal reception and measurement accuracy are heavily dependent on the surrounding environment of the survey stations, further complicating receiver deployment [15,16]. Therefore, there is a need to identify alternative data sources for GCPs that can complement spaceborne InSAR remote sensing methods to generate accurate DEMs.

To date, several Earth observation lidar missions have contributed to the extraction of GCPs at the regional [17,18] and global scales [12]. The ICESat-1 mission, carrying the Geoscience Laser Altimeter System (GLAS) instrument, collected full-waveform lidar data [19]. GLAS measures the two-way travel time of laser pulses reflected off the ground from 2003 and 2009, with a ground-track repeat cycle of 183 days. It produced a series of footprints with approximately 70 m diameter, which are separated by nearly 170 m sampling intervals along the track and 30 km across the track (at the equator) and 5 km at 80º latitude. However, its single-beam laser was unable to distinguish slope effects from true elevation changes on an orbit-by-orbit basis [20]. Additionally, its lower spatial resolution in both the along-track and cross-track directions limits its ability to measure rougher terrains. To address these limitations, the National Aeronautics and Space Administration (NASA) launched the Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) in September 2018. This mission uses a single instrument, named the Advanced Topographic Laser Altimeter System (ATLAS), which uses photon-counting technology to make measurements [21]. Unlike full-waveform lidar (e.g., ICESat mission), photon-counting lidar features a smaller footprint size and a higher laser pulse repetition rate. This results in a higher point density and the ability to acquire near-continuously profiles along the flight directions [22], offering new opportunities for the fine-scale sampling of the global terrain.

The ICESat-2 ATL08 product provides users with valuable terrain and canopy height information derived from the ATL03 product. According to [23,24] studies, ATL08 has approximately 140 signal photons per 100 m segment over vegetated surfaces. Each ATL08 segment includes both canopy height and terrain height data. Additionally, the terrain parameters in the ATL08 product include both the statistical values of the ground photons within each 100 m segment and the best-fit terrain elevation at the midpoint of each segment (h_te_bestfit) [24]. However, the accuracy of these data are typically affected by numerous factors, such as terrain, vegetation coverage, satellite orbit attitude, and atmospheric conditions [25]. Since the release of the ATL08 product, numerous studies have assessed its accuracy under diverse experimental conditions. For example, Neuenschwander et al. validated the terrain and canopy height estimates from ATL08 against airborne lidar data, reporting that the vertical root mean square error (RMSE) of terrain height retrievals at the Alaska Tundra/Taiga ecotone was better than 1 m [22,26]. Similarly, Xing et al. evaluated the accuracy of ATL08 terrain estimates in forested areas, demonstrating that the estimates could acquire sub-meter accuracy and that the strong beam outperformed the weak beam [27]. Zhu et al. further conducted a quantitative assessment of ATL08 terrain accuracy and analyzed the effects of vegetation type, terrain slope, and other factors. They concluded that the accuracy of the terrain estimates was significant impacted by the presence of vegetation and terrain slope [28]. Based on the these studies and the related literature [29,30], it is evident that complex, vegetation-covered terrain inevitably introduces measurement errors (2–5 m) or outliers (tens of meters) in ATL08 terrain estimates, posing a challenge for deriving high-precision and reliable GCPs from ATL08 data.

To eliminate low-quality ICESat-2 ATL08 terrain estimates and outliers, two primary methods are typically used: threshold-based methods and evaluation labels-based methods. The core idea of the threshold-based methods is to compare ATL08 terrain estimates with an external reference dataset (e.g., a DEM) and remove data that fall outside an elevation difference threshold between the ATL08 terrain estimates and the reference data [31,32]. However, this method heavily depends on the quality of the external topographic data. In the absence of high-precision reference data, extracting high-precision control points from ICESat-2 data products becomes challenging. The second method involves extracting high-accuracy control points by utilizing attribute parameters included in the ICESat-2 dataset [30,33], such as terrain factors, terrain surface coverage data, the number of ground photons, and the cloud confidence flag. These methods aim to enhance the stability and precision of the extracted elevation control points. However, they still rely on accurate external data for evaluating the processing accuracy, and the setting of thresholds is subject to personal subjectivity.

Following this line of investigation, the objective of the present work is to propose a new approach for grading ATL08 terrain estimates, while also facilitating the extraction of high-quality GCPs. To achieve this, the proposed approach is divided into three stages. First, we combine the raw ICESat-2 photon data and ATL08 terrain estimates to fit continuous terrain. Next, the raw ICESat-2 photon data and fitted continuous terrain curve are utilized to create a photon-distribution curve, referred to as pseudo-waveform hereafter. Finally, by analyzing the positions and values of the peaks of the pseudo-waveform, the ATL08 terrain estimates are graded into various accuracy levels. To evaluate our proposed method, four typical study areas from the National Ecological Observatory Network (NEON), each characterized by various terrain features, are selected. Additionally, high-quality digital terrain models (DTMs) derived from NEON lidar measurements are employed to quantitatively validate the accuracy of the ICESat-2 ATL08 terrain estimates.

2. Study Areas and Datasets

2.1. Study Areas

To investigate the performance of our proposed approach, we selected four typical study areas within the NEON [34], including Dead Lake (DELA), Mountain Lake Biological Station (MLBS), Great Smoky Mountains National Park (GRSM), and Treehaven (TREE), as shown in Figure 1. These study areas feature diverse forest coverage types driven by local climate variations (see

Table 1. Summary of climate for each study area.

Study Area	Climate Characteristics	Mean Annual Precipitation (mm)	Average Annual Temperature (°C)
DELA	Subtropical climate.	1370	7.6
MLBS	Humid continental climate.	1227	8.6
GRSM	The variable elevations reflect the diverse climate and habitat changes characteristic of the study area.	1440 (lower valleys)~2160 (some of the highest peaks)	13.1
TREE	Bitterly cold winters and generally cool summers with brief periods of excessive heat.	797	4.8

).

Table 2. Summary of forest and topographic characteristics across the four study areas. For vegetation coverage fraction and terrain slope, the mean value is followed by the standard deviation, and the range of values is provided in parentheses.

Study Area	Forest Type	Tree Species	Average Forest Height (m)	Vegetation Coverage Fraction (%)	Terrain Slope (°)
DELA	broad-leaf forests	Bottomland hardwood forest, oak and hickory with some pine	14.84	61.14	3.45
				12.68	4.77
				(10–90)	(0–40)
MLBS	broad-leaf forests	Red maple and white oak	13.35	52.07	15.24
				9.10	8.76
				(33–67)	(0–45)
GRSM	broad-leaf forests	Yellow poplar, red maple and chestnut oak	16.27	59.68	25.41
				8.14	11.44
				(18–82)	(0–68)
TREE	coniferous forests	Black spruce and tamarack	9.09	76.27	7.82
				9.60	6.98
				(13–92)	(0–40)

The detailed tree species information is collected from https://www.neonscience.org/field-sites (accessed on 17 June 2024).

provides detailed information on the dominant forest cover types, tree species, average forest height, vegetation coverage fraction, and topographic characteristics (e.g., flat, valley, and mountain areas). Moreover, as dedicated test sites within the NEON, these four study areas benefit from frequent airborne lidar acquisitions, offering high-quality ground reference data for validating our results.

2.2. Datasets

2.2.1. ICESat-2 ATL03 and ATL08 Data

The ICESat-2 mission provides 21 data products ranging from ATL00 to ATL23, excluding ATL05 and ATL22. These products are categorized into four levels: Level-1, Level-2, Level-3A, and Level-3B. Detailed information regarding the product specifications and spatiotemporal resolution is available on the official ICESat-2 website. This study utilizes two data products: the Level-2 data product ATL03 (Global Geolocation Photon Data) and the Level-3A product ATL08 (Land and Vegetation Height Data). ATL03 provides essential geolocation and timing information for each received photon from the ATLAS instrument [35], including time, latitude, longitude, ellipsoidal height, and along-track flight distance. ATL08 is derived from ATL03 by filtering out background noises [24], delivering land and vegetation height information at a resolution of 100 m in the along-track direction and 11 m in the across-track direction. Both datasets can be freely downloaded from NSIDC (National Snow and Ice Data Center).

Table 3. The ICESat-2 data for ATL03 and ATL08 were used in this study.

Study Area	Acquisition Dates	Data Acquisition
DELA	26 November 2018	ATL03/ATL08_20181126072045_08960106_005_01.h5
	7 May 2019	ATL03/ATL08_20190507113156_05990302_005_01.h5
	4 September 2019	ATL03/ATL08_20190904054742_10410402_005_01.h5
	4 December 2019	ATL03/ATL08_20191204012732_10410502_005_01.h5
	24 May 2020	ATL03/ATL08_20200524051937_08960706_005_01.h5
	3 August 2020	ATL03/ATL08_20200803135050_05990802_005_01.h5
	23 August 2020	ATL03/ATL08_20200823005922_08960806_005_01.h5
	2 November 2020	ATL03/ATL08_20201102093039_05990902_005_01.h5
MLBS	23 March 2019	ATL03/ATL08_20190323011928_12920206_005_01.h5
	20 March 2020	ATL03/ATL08_20200320075832_12920606_005_01.h5
GRSM	6 May 2019	ATL03/ATL08_20190506232207_05910306_005_01.h5
	13 July 2019	ATL03/ATL08_20190713080125_02330402_005_01.h5
	11 August 2019	ATL03/ATL08_20190811063738_06750402_005_01.h5
	10 November 2019	ATL03/ATL08_20191110021732_06750502_005_01.h5
	10 April 2020	ATL03/ATL08_20200410190101_02330702_005_01.h5
	3 August 2020	ATL03/ATL08_20200803014101_05910806_005_01.h5
TREE	4 September 2019	ATL03/ATL08_20190904054742_10410402_005_01.h5
	5 November 2019	ATL03/ATL08_20191105025131_05990502_005_01.h5
	11 November 2019	ATL03/ATL08_20191111145042_06980506_005_01.h5
	10 February 2020	ATL03/ATL08_20200210103025_06980606_005_01.h5

summarizes the ATL03 and ATL08 data employed in this study, with the associated data strips depicted in Figure 1.

To ensure data quality and minimize signal noise, we selected only strong-beam and continuous-strip data acquisitions from ICESat-2. For the ATL08 data product, key parameters, including longitude (lon_ph), latitude (lat_ph), and the best-fit terrain elevation (h_te_best_fit), were extracted for the subsequent application. For the ATL03 data product, the parameters used in this study, along with their corresponding explanations, are listed in

Table 4. Attributes of ATL08 and ATL03 products used in this study, along with their corresponding explanations.

Dataset	Fields	Explanations
ATL08	lon_ph	Longitude and latitude of the center-most signal photon within each segment.
	lat_ph
	h_te_best_fit	The best-fit terrain elevation at the midpoint location of each 100 m segment, relative to WGS-84.
ATL03	lon_ph,	Longitude and latitude of each received photon.
	lat_ph
	h_ph	Height of each received photon, relative to WGS-84.
	ds_gt	Recordings of the ground tracks (gt1l, gt1r, gt2l, gt2r, gt3l, gt3r), which are used for recognizing strong beams.
	sc_orient	Tracking spacecraft orientations.

. These data strips were then clipped to match the boundaries of the study areas (as shown in Figure 1) based on the geolocation coordinates of ICESat-2.

2.2.2. Airborne Lidar-Derived DTM

A DTM derived from airborne lidar data was utilized as a reference for ground elevation evaluations. This DTM, freely available from the official website of the NEON foundation [36], was generated by the NEON Airborne Observation Platform (AOP) between 2017 and 2018. Additionally, the DTM is provided as a raw raster with a 1 m grid spacing, using the Universal Transverse Mercator (UTM) projection coordinate system under the ITRF2000 and NAVD88 (Geoid 12A) vertical reference. The planar error is approximately 5–15 cm, and the overall elevation accuracy ranges from 5 to 35 cm, making it suitable for accuracy assessment. To facilitate the follow-up analysis and validation, the ICESat-2 measurements, initially presented in the WGS84 horizontal coordinate system, were projected to the UTM. Additionally, the vertical reference of ICESat-2 measurements was converted from the WGS84 ellipsoidal height to the NAVD88 geoid using the VDatum 4.7 Online tool [37,38], developed by National Oceanic and Atmospheric Administration’s (NOAA) National Ocean Service.

2.2.3. Vegetation Coverage Field Data

In this study, we utilize the Global 30m Landsat Tree Canopy Cover Version 4 product (available from the NASA Land-Cover/Land-use Change Program) as a proxy of canopy density to assess the impact of vegetation coverage on the accuracy of ICESat-2 terrain estimates. This dataset, derived from the Landsat image archive and released in May 2019, provides a spatial resolution of 30 m. It shows that the tree canopy coverage per pixel is between 0% and 100%, which is also called the vegetation continuous field (VCF). Further details on the development of the VCF product can be found in [39]. To facilitate further analysis and validation, a nearest neighbor interpolation method was utilized for the VCF dataset, ensuring that values corresponding to the ICESat-2 terrain estimates were accurately extracted.

3. Methods

To extract high-accuracy GCPs from ICESat-2 ATL08 data, the proposed method consists of three primary components, as shown in Figure 2.

The first step involves continuous terrain acquisition. To reduce data quality and reduce noise, only strong beam data and continuous stripes are extracted. Subsequently, the ground points recorded in the ATL08 product are treated as prior sampled data points, and a B-spline function is used to fit the terrain along each strip. The rationale for choosing the B-spline interpolation function is twofold: (1) B-splines offer localized support for basis functions [40], and (2) B-spline functions create smooth surfaces or profiles that preserve small features while requiring fewer sampled data points [41]. This process generates continuous terrain measurements along the ground track direction.

The second part presents a generation method of pseudo-waveform. To achieve this, irregular statistical zones are established around each ATL08 terrain estimate. The height difference between ATL03 photons within the buffer zone and the fitted terrain derived from the first step is then calculated. To model the photon cloud distribution, the Gaussian kernel density estimation (KDE) is applied to fit the density distribution curve of AT03 photons. KDE is typically used to model a discrete sample of data as a continuous distribution [42] and is generally utilized for waveform decomposition [43,44]. The resulting density distribution curve represents the height difference distribution along the vertical direction and is hereafter referred to as the pseudo-waveform.

In the third part, based on the characteristics of the established pseudo-waveform, three photon distribution scenarios are identified: Scenario 1 (L1), where photons are concentrated near the ground surface; Scenario 2 (L2), where photons are either evenly distributed above the ground surface or show a distinct layered distribution; and Scenario 3 (L3), where photons located away from the ground surface are misidentified as ground points and are deemed unusable. In this study, these three scenarios are assigned to three accuracy levels: L1, L2, and L3. The subsequent sections will describe each of these steps in greater detail.

3.1. Continuous Terrain Fitting

Before fitting the continuous terrain, it is important to preprocess both the ATL08 and ATL03 products. Following the selected data strips and extracted attributes, the along-track distance for each ATL08 ground point is calculated. Then, the B-spline curve (Equation (1)) is employed to create a continuous terrain model as follows:

$$S\left(x\right)=\sum _{i=0}^n{c}_i{B}_{i,k}\left(x\right)$$

(1)

where $x$ denotes the along-track distance of each elevation point, and $c_i$ represents the fitting parameters. $B_{i,k}\left(x\right)$ stands for the B-spline basis function defined recursively.

When $k=0,$

$$B_{i,0}\left(x\right)=\left\{\begin{array}{cc}1& if t_i\le x\le t_{i+1}\\ 0& otherwise\end{array}\right.$$

(2)

When $k>0$,

$$B_{i,k}\left(x\right)={\displaystyle \frac{x-x_i}{x_{i+k-1}-x_i}}{B}_{i,k-1}\left(x\right)+{\displaystyle \frac{x_{i+k}-x}{x_{i+k}-t_{i+1}}}{B}_{i+1,k-1}\left(x\right)$$

(3)

where $x_i$ represents the knot points, $i$ is the number of knot points, and $k$ is the degree of basis functions. For example, when $k=4$, the cubic B-spline basis function is recursively calculated from lower-degree basis functions. Therefore, the interpolation function for the terrain fitting is expressed as a linear combination of several basis functions:

$$B_{ij}=B_{i,k}\left(x_j\right), with i=\mathrm{0,1},\dots ,n+k-2 and j=\mathrm{0,1},\dots ,n$$

(4)

and the coefficient vector $c$ satisfies,

$$B^{T}MBc=B^{T}y$$

(5)

where $M$ denotes a symmetric positive definite matrix, and $\mathrm{y}$ is the vertical coordinate vector of the interpolation data, i.e., elevation, expressed as $y={\left[y_0,y_1,\dots ,y_n\right]}^T$. To optimize the B-spline function fitting and preserve terrain details, the smoothness of the fitting spline is controlled using

$$\sum _i{\left(w_i\left(y_i-H\left(x_i\right)\right)\right)}^2\le s$$

(6)

where $x_i$ and $y_i$ correspond to the along-track distance and vertical height of each ATL08 ground point, respectively. $H\left(x_i\right)$ stands for the calculated elevation value from the fitted B-spline function. $w_i$ represents the weight for each ground point. Given that the ATL08 ground points are approximately equidistantly distributed and assumed to be observed with equal precision, $w_i$ is set to 1. Additionally, the smoothing parameter, denoted as $s$, which controls the smoothness of a fitted curve and plays a critical role in balancing smoothness and terrain detail preservation. A bigger $s$ results in a smoother curve, but finer terrain details may be lost [45]. In this study, the smoothing factor is set within the range of $s=0.1$ to 0.3. Alternatively, the optimal value of $s$ can be determined by splitting the sampled data into training and testing sets. Then, the training data are used to fit the B-spline functions for a series of $s$ values, and the testing data are used to assess the accuracy. The selected value for $s$ is the one that minimizes the bias between the fitted values and the ATL08 measurements.

3.2. Generation of Pseudo-Waveform

3.2.1. Establishment of Statistical Buffer Zones

After performing terrain fitting, we established a statistical buffer zone that can be adaptively adjusted based on the terrain changes at each ATL08 ground point. This buffer was defined using a window with an along-track length of 40 m and a vertical height of 40 m, and parameters were empirically determined based on the local terrain and forest height within the study area. Specifically, the vertical buffer size was designed to encompass all tree height classes within the study area, which were preliminarily estimated according to ATL08 forest height data. The along-track buffer size was chosen considering local terrain features and the number of photons available. Further, the vertical buffer size can be dynamically adjusted to reflect changes in the local terrain. A similar idea was used in [46], where it was utilized to estimate the diffuse attenuation coefficient of water using ICESat-2 data. The resulting statistical buffer zone, illustrated by the purple rectangle of Figure 3, effectively captures the majority of objects under normal terrain conditions. This preserves the critical features of photon distribution while excluding noise interference. The buffer zone is mathematically represented as:

$$B=\{\left(l,h\right)\left|L_{ATL08}-D\le l\le L_{ATL08}+D, H_{ATL08}-W\le h\le H_{ATL08}+W\right\}$$

(7)

where $L_{ATL08}$ and $H_{ATL08}$ are the along-track distance and vertical height of each ATL08 terrain estimate, respectively. $D$ and $W$ refer to the buffer size in the ground track and vertical directions, respectively, which were set to 20 m and 40 m in this study.

3.2.2. Pseudo-Waveform Generation

To generate the pseudo-waveform, the height difference between each ATL03 photon and the fitted terrain is calculated based on the previously established buffer area (Equation (7)). The height difference, denoted as $∆h$, is given by $∆h=\left[H_{ATL03}-H_{ATL08}\right]$. Next, we apply the Gaussian KDE method to fit the kernel density distribution curve of the vertical height differences, which is defined as the pseudo-waveform, as illustrated in [46]. The Gaussian KDE is mathematically expressed as

$$\widehat{f}\left(x\right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{1}{{\left(2\pi w^2\right)}^{\frac{m}{2}}}}\exp \left(-{\displaystyle \frac{|x-x_i{|}^2}{2w^2}}\right)$$

(8)

where $\widehat{f}\left(x\right)$ represents the density value at the ATL08 ground point $x$, m is the sample dimension, and $\left|x-x_i\right|$ is the distance between the ATL08 measurement point $x$ and sampling point $x_i$. w is the smoothing parameter, which controls the width of the kernel function and affects the shape of the density curve or pseudo-waveform. Specifically, a larger w leads to a smoother pseudo-waveform, reducing the features that capture ATL03 photons clustering. Conversely, a smaller $w$ results in sharp peaks, which may make it harder to distinguish the key features of the waveform. To determine the optimal parameter value, we tested various $w$ values ranging from 0.1 to 0.5. The visual inspection of the resulting curves suggested that w = 0.2 best preserved the primary features of the waveform while maintaining an appropriate number of peaks.

3.3. Grading of Terrain Estimates Based on the Pseudo-Waveform

Figure 4a shows the distribution of ATL03 photons within a statistical buffer zone, while Figure 4b the generated pseudo-waveform, which exhibits clear peaks corresponding to the concentrated area of ATL03 photons. By analyzing the characteristics of the pseudo-waveform, we can assess the similarity between ATL08 terrain estimates and ATL03 photons, as well as compare the relative heights. In this study, we focus on the three highest peaks from the pseudo-waveform, assuming the forest scenario consists of canopy, shrub, and understory layers [47]. Having said this, it is important to note that the proposed approach does not impose any constraints on the forest scenario. The corresponding peak values are labeled as LgstPeak_Value, SecPeak_Value, and ThdPeak_Value, with their respective elevation positions indicated by LgstPeak_Index, SecPeak_Index, and ThdPeak_Index, as illustrated in Figure 4b and explained in

Table 5. Parameters defined in the proposed grading strategy.

Parameters	Explanation
LgstPeak_Value	The maximum peak value in the pseudo-waveform, which is the maximum value of the function established by KDE.
LgstPeak_Index	The elevation difference between the maximum value point and ground.
SecPeak_Value	The second largest peak value in the pseudo-waveform.
SecPeak_Index	The elevation difference between the second largest value point and ground.
ThdPeak_Value	The third largest peak value in the pseudo-waveform.
ThdPeak_Index	The elevation difference between the third largest value point and ground.

.

Based on the characteristics of these peaks, we categorize the ATL08 estimates into three levels of accuracy (L1, L2, and L3). Specifically,

• L1 (highest accuracy): ATL03 photons are concentrated near the fitted terrain, indicating the highest accuracy of the ATL08 terrain estimates. These estimates provide the most precise terrain elevation information.
• L2 (moderate accuracy): Approximately 50% of ATL03 photons are concentrated near the fitted terrain, with a few photons scattered in the canopy. The accuracy of these terrain estimates is slightly lower than that of the L1 estimates.
• L3 (lowest accuracy): ATL03 photons are far away from the fitted terrain. These estimates indicate the lowest accuracy, possibly mistaking noise points as valid elevation control points. As a result, they should be excluded due to poor precision.

The specific judgment conditions for each accuracy level are summarized in

Table A1. Specific conditions and the associated descriptions for each ATL08 terrain estimation level.

Level	Sufficient Condition	Necessary Conditions	Description
Level1 (highest accuracy)	$lgstPeak\_Index\in \left[-1,1\right]$	$lgstPeak\_Value/secPeak\_Value>1.2$	The first highest peak is close to the ground surface, while the second highest peak is 10 m above the ground, and the highest peak value is 1.2 times greater than the secondary peak value.
		$secPeak\_Index-lsgtPeak\_Index>10$
		$\left|lgstPeak\_Index+secPeak\_Index\right|<2$	The first and second highest peaks are distributed on either size of the zero value, which could be caused by the fitting errors, but it still indicates a strong ground signal.
		$lgstPeak\_Value/secPeak\_Value<1.5$
	$secPeak\_Index\in \left[-1,1\right]$	$lgstPeak\_Value/secPeak\_Value<2$ $lsgtPeak\_Index-secPeak\_Index>10$	The first highest peak is close to the canopy, while the second highest peak is around the ground. However, their height difference is not significant, indicating that it still represents a strong ground signal.
Level2 (moderate accuracy)	$lgstPeak\_Index\in \left[1,10\right] \mathrm{o}\mathrm{r} lgstPeak\_Index\in \left[-3,-1\right]$	$\left|lgstPeak\_Index+secPeak\_Index\right|<2$ $lgstPeak\_Value/secPeak\_Value<1.5$	There are many branches below the canopy, and photons are concentrated near ground. In addition, two peaks are near the ground surface, which are caused by fitting errors.
		$lgstPeak\_Value/secPeak\_Value>1$& $(lsgtPeak\_Index-secPeak\_Index<5$or $lsgtPeak\_Index-secPeak\_Index>10)$	The first highest peak appears near the ground or canopy. The second highest peak is either close to or far away from the first highest peak. And their value ratio is larger than 1.5 also appear near the ground surface or canopy, which corresponds to the shrub area and forest area, respectively.
	$lgstPeak\_Index>10$	$secPeak\_Index\in \left[-4, 4\right]$ $\left(lgstPeak\_Value/secPeak\_Value\right)>1$ $\left(secPeak\_Value/thdPeak\_Value\right)>1$	The secondary peak is close to the ground surface, corresponding to relatively dense and short shrubs on the ground. In addition, there is another peak in the middle layer, corresponding to certain branches below the canopy.
		$secPea{k}_{Index}>8 \exists peak\in \left[-4, 4\right]$	The secondary peak is above 8 m, corresponding to a relatively dense sub-canopy layer. In addition, there is also a peak on the ground, corresponding to some low and dense shrub cover on the ground.
		$\left|thirdPeak\_Index+secPeak\_Index\right|<2$ $\left(secPeak\_Value/thirdPeakv\_Value\right)<1.5$	Except for the main peak in the canopy, there are two peaks near the zero value, which may be caused by fitting errors.
		$Len\left(peaks\right)=2$ $Peak\_Indexs\left[0\right]\in \left[-4, 5\right]$	There are two peaks in total, with one peak near the ground surface. This indicates that the terrain belongs to an ideal forest area with obvious stratification.
Level3 (lowest accuracy)	$All\left(peak\_Indexs\right)>5$		All peak points are greater than 5 m, and there are no peaks near the ground, indicating that the ATL08 misidentified ground noise as terrain points. In this case, ICESat-2 ATL08 ground elevation indicates a negative bias.
	$Any\left(peak\_Indexs\right)<-3$		There are peaks below 3 m, indicating that the ATL08 misidentified the photon points above the ground as terrain points. In this case, ICESat-2 ATL08 ground elevation indicates a positive bias.

.

3.4. Elevation Accuracy Assessment

To assess the elevation accuracy of ICESat-2 ATL08 terrain estimates, we employed high-quality ALS-based DTMs (described in Section 2.2.2) as a reference. Since the ICESat-2 terrain estimates are a series of discrete points, the nearest neighbor interpolation method [48] was selected to calculate the ALS DTM elevation values for each ICESat-2 data footprint. For the quantitative evaluation of elevation accuracy, several statistical matrices like root mean square error (RMSE) and the coefficient of determination $\left(R^2\right)$ were utilized [18]. Additionally, the stand deviation (STD) indicator was employed to measure the dispersion (or uncertainty) of the dataset relative to its mean, which can be computed using Equation (11).

$$RMSE={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left(X_i-Y_i\right)$$

(9)

$$R^2=1-{\displaystyle \frac{{\sum }_{I=1}^n{\left(X_i-Y_i\right)}^2}{{\sum }_{I=1}^n{\left(M-Y_i\right)}^2}}$$

(10)

$$STD=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(\left(X_i-Y_i\right)-\mu \right)}^2}$$

(11)

where $n$ is the number of ICESat-2 ALT08 ground points, $X_i$ is the elevation value of ICESat-2 measure points, and $Y_i$ is the corresponding elevation value of ALS DTM. $M$ is the mean of elevation values obtained from airborne lidar, and $\mu$ stands for the mean or expected bias value between the elevation values of ICESat-2 ATL08 terrain and ALS DTM.

4. Results and Discussion

4.1. Results

This study presents a universal method for selecting ICESat-2 ATL08 terrain estimates, which applies to both mountainous and forested areas. To evaluate the performance of the proposed method, four study areas, characterized by various forest types and terrain conditions, were selected. In each study area, the ICESat-2 ATL08 terrain estimates were divided into three groups, labeled as L1, L2, and L3, representing high, medium, and low accuracy (outliers), respectively.

Figure 5 and Figure 6 illustrate the distribution of ATL08 terrain estimates across varying terrain slopes (derived by airborne lidar DTM) and vegetation coverage conditions (derived from the VCF data). All three levels of terrain estimates are present across terrain slopes and canopy densities. Furthermore, the terrain and vegetation coverage conditions did not exhibit a systematic effect on the distribution of terrain estimates. These findings demonstrated that the proposed method is minimally influenced by the aforementioned two factors (i.e., terrain slope and vegetation coverage), highlighting its robustness in different environmental contexts.

To facilitate a visual comparison of the graded ATL08 terrain estimates within each study area, the graded terrain estimate results were superimposed onto the DEM, as shown in Figure 7. In this figure, all L1, L2, and L3 terrain estimates can be observed in regions with flatlands, valleys, and ridges. Notably, L1 terrain estimates exhibit a higher density in flatlands, whereas L3 terrain estimates are predominantly observed in valleys and ridges. This distribution aligns with the previously discussed findings, reinforcing the inference regarding the relationship between terrain features and the accuracy of ATL08 terrain estimates.

In addition to the visual assessment, we quantitatively evaluated the accuracy of the graded ATL08 terrain estimates of each level using scatter plots (Figure 8) and summarized the statistical results in

Table 6. Accuracy of the graded ATL08 terrain estimates in four study areas. RMSE_flu indicates the percentage increase (positive) or decrease (negative) in RMSE compared to the original ATL08 terrain estimates.

Study Area	Indicator	Raw	L1	L2	L3
DELA	RMSE (m)	2.05	0.99	1.75	4.85
	R2	0.92	0.98	0.92	0.26
	STD (m)	2.04	0.94	1.68	4.56
	RMSE_flu	-	−51.7%	−14.63%	+136.6%
MLBS	RMSE (m)	1.22	0.51	1.46	4.96
	R2	1.0	1.0	1.0	1.0
	STD (m)	1.17	0.46	1.36	4.48
	RMSE_flu	-	−58.2%	+19.7%	+306.5%
GRSM	RMSE (m)	11.13	1.88	6.69	14.4
	R2	1.0	1.0	1.0	1.0
	STD (m)	11.03	1.88	6.66	13.98
	RMSE_flu	-	−83.1%	−39.9%	+29.4%
TREE	RMSE (m)	2.08	0.65	1.54	6.57
	R2	0.99	1.0	1.0	0.93
	STD (m)	2.03	0.62	1.45	6.07
	RMSE_flu	-	−68.8%	−26.0%	+215.9%

, which includes the RMSE, R², STD, and growth rate of RMSE. The scatter plot (see Figure 8) shows that the L1 terrain estimates across all four areas demonstrate a high level of consistency. The RMSE of the L1 estimates is at least 50% smaller than that of the raw ATL08 terrain estimates, with the highest accuracy reaching 0.51 m. Specifically, in DELA, MLBS, and TREE, the accuracy of the L1 estimates achieves sub-meter levels. In GRSM, the RMSE of L1 is reduced by 83% compared to the raw data before grading (1.88 m vs. 11.13 m). These findings indicate that the L1 terrain estimates extracted by the proposed method exhibit significantly higher accuracy. Furthermore, when focusing on the STD indicator, it appears that all L1 terrain estimates across the four study areas show a lower STD value compared to the original ATL08 terrain estimates. This result suggests that the L1 terrain estimates derived by the proposed method exhibit reduced uncertainty, offering more reliable terrain data across varying study areas.

Additionally, by comparing the accuracy of L1 terrain estimates across various forest coverage types, we observed that the TREE study area, characterized by coniferous forests, exhibits overall higher accuracy (RMSE = 0.65 m) than the two broadleaf-forest study areas, DELA and GRSM, which had RMSE values of 0.99 m and 1.88 m, respectively. However, in the MLBS study area, which is also characterized by broadleaf forests, the L1 terrain estimates showed higher accuracy than those in the coniferous forests study area (i.e., TREE). This result could be attributed to the nighttime acquisition of ICESat-2 data in the MLBS area, which can help mitigate the impact of broadleaf forests on terrain estimation accuracy. This finding is consistent with previous studies conducted in different forest ecosystems [22,28,49]. Based on these results, we recommend using strong beam data acquired at night for extracting high-accuracy elevation control points.

In the case of L3 estimates, the RMSE in three study areas (DELA, MLBS, and TREE), is more than double that of the pre-graded terrain estimates. In particular, the RMSE of MLBS exhibits the largest increase, more than three times that of the original terrain estimates. This result indicates that the proposed method is effective in identifying outliers. Regarding the L2 terrain estimates, their accuracy lies between the L1 and L3 estimates, supporting the validity of the grading strategy proposed in Section 3.3. These findings can be summarized as follows: (1) in regions with relatively flat or uniform slope variations, such as MLBS, the proposed method can effectively remove outliers (i.e., L3 terrain estimates); (2) in regions with significant terrain slope variations, like GRSM, the grading strategy successfully identifies high-accuracy elevation control points (i.e., L1 terrain estimates).

4.2. Discussion

4.2.1. Comparison Between the Proposed Method and Existing Methods

By testing the proposed method in various study areas, its effectiveness has been demonstrated. To further evaluate its suitability and advantages for selecting high-precision GCPs, we conducted a comparative analysis between the proposed control points selection strategy and the conventional threshold-based method [30,50,51]. It is noteworthy that the threshold-based method divides the terrain estimates into two levels, namely Level 1 and Level 2, whereas our proposed method grades the estimates into three levels.

To quantitatively assess the performances of both methods, we calculate the RMSE using Equation (9), as well as the retention rate for L1 estimates alone and for the combination of L1 and L2 estimates. These results are then compared with those of the L1 estimates obtained through the threshold-based method, considering its L2 estimates as outliers (Figure 9). From Figure 9a, we observed that, in terms of RMSE, the accuracy of L1 estimates achieved by our proposed method is significantly superior to that of the threshold-based method. Next, we analyzed the retention rate, which indicates the proportion of specific points of interest, with the results shown in Figure 9b. The result demonstrates that our proposed method outperforms the threshold-based method in the MLBS and TREE study areas and performs similarly in the DELA and GRSM study areas.

Furthermore, the accuracy of the combination of L1 and L2 estimates achieved by our proposed method is higher than that of the L1 estimates obtained through the threshold-based method in the DELA, MLBS, and TREE study areas, as evidenced by the lower RMSE values. Additionally, the data retention rate of the proposed method far exceeds that of the threshold-based method in these three study areas. However, in the GRSM study area, the accuracy of the L1 and L2 combination is lower than that of the L1 estimate achieved by the threshold-based method, although it has a higher data retention rate. This discrepancy could be attributed to the steep and undulating terrain conditions in this study area (see Table 2). On the other hand, when fitting continuous terrain using the method described in Section 3.1, the best-fit terrain elevation (h_te_best_fit) from the ATL08 data product is utilized. In this case, the 100 m interval along the track may be insufficient to accurately represent the complex terrain characteristics [29,52], potentially leading to deviations in the generated pseudo-waveform from the actual terrain shapes. As a result, this can influence the effectiveness of the grading strategy (outlined in Section 3.3).

In general, the comparisons indicate that our proposed method outperforms the threshold-based method in selecting high-accuracy terrain estimates from ICESat-2 ATL08 data. Whether considering the L1 terrain estimates or the combination of L1 and L2, our method demonstrates superior performance compared to the L1 terrain estimates obtained by the threshold-based method, especially in regions with relatively uniform terrain changes. Furthermore, the proposed method significantly improves the data retention rate of ICESat-2 ATL08 terrain estimates, ensuring the better preservation of valuable data points for related applications.

4.2.2. Limitations and Improvement for Future Work

The proposed method enables the extraction of graded ATL08 terrain estimates, but there are still some limitations that need to be addressed. Firstly, in the area with steep and undulating terrain, the established pseudo-waveform often exhibits significant variance, which introduces uncertainty in the grading of ATL08 terrain estimates. Figure 10 displays the statistical distribution of height errors (between the ATL08 terrain estimates and lidar-derived DTM) across different terrain slope classes. Overall, the dispersion in the samples’ height errors was associated with terrain slopes in all four study areas, showing a positive correlation. In particular, a higher dispersion was observed for the L3 terrain estimates, indicating higher uncertainty in the corresponding terrain estimates. Additionally, Figure 11 shows the statistical distribution of height errors across different VCF classes. Similarly to the results observed with the terrain slopes, a general positive correlation was found between the dispersion of height errors and VCF values. This result can partly be attributed to the fact that, as the canopy density increases, the overlap between canopy signals and ground signals becomes more significant. This overlap poses challenging in accurately identifying ground features from the established pseudo-waveform, thereby leading to higher uncertainty in the terrain estimates.

Secondly, the choice of the statistical buffer zone window directly impacts the ability of the generated pseudo-waveforms to accurately reflect the signal distributions of the canopy and ground. Specifically, a smaller window size may obscure the statistical characteristics of the pseudo-waveform due to a reduced number of photons, while a larger one may result in the inaccurate extraction of pseudo-waveform features in complex terrain conditions, thereby influencing the outcomes. In this study, the corresponding parameters were empirically determined based on the forest height bins and local terrain conditions.

To further enhance the performance and application capability of our proposed method, future research will focus on accurately extracting pseudo-waveforms. A key improvement will involve computing the slope of ATL08 photon points along the ground track, which will allow the method to adapt more effectively to complex and changing terrain. Building upon this foundation, an adaptive selection of the window for establishing the statistic buffer zone will enhance the method’s ability to capture terrain variations, thereby enabling the extraction of more accurate pseudo-waveform characteristics.

By addressing and managing the aforementioned limitations, the graded ATL08 terrain estimates extracted by our method can further enhance its applicability in various scenarios. For example, high-precision ATL08 terrain data could help improve the accuracy of InSAR-derived DEMs or global topographic mapping [53,54]. Additionally, it can provide essential ground reference data for identifying and correcting errors in existing globally available DEMs [55,56,57], such as SRTM DEM and TanDEM-X 90m DEM.

5. Conclusions

In this study, we proposed using ICESat-2 ATL03 photon point cloud data to establish pseudo-waveform functions for ICESat-2 ATL08 terrain points. This approach enables the grading of ATL08 terrain estimates based on pseudo-waveform characteristics and their spatial positions. The evaluation tests demonstrated that the proposed method effectively classified the ATL08 terrain estimates into different accuracy levels, which is crucial for extracting reliable ground control points (GCPs). Additionally, compared to existing methods for extracting GCPs from ATL08 products, our proposed method does not rely on external reference data, eliminating uncertainties related to the accuracy and accessibility of such data. Furthermore, it is well-suited for batch processing of ATL08 terrain estimates. These findings demonstrate that pseudo-waveform features derived from ATL03 photons can effectively aid in grading ATL08 terrain estimates.

Due to the simplicity and efficiency of the approach described in this paper, the entire process has been automated. This method also holds great potential for the large-scale extraction of elevation control points. As the ICESat-2 satellite continues to observe the Earth and generates abundant data, the proposed method will play a crucial role in ensuring the consistency and accuracy of terrain elevation control points in geodetic applications such as global mapping.