Estimating Terrain Slope from ICESat-2 Data in Forest Environments

Abstract

The global digital elevation measurement (DEM) products such as SRTM DEM and GDEM have been widely used for terrain slope retrieval in forests. However, the slope estimation accuracy is generally limited due to the DEMs’ low vertical accuracy over complex forest environments. The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) mission shows excellent potential for slope estimation because of the high elevation accuracy and unique design of beam pairs. This study aimed to explore the possibility of ICESat-2 data for terrain slope retrieval in the United States forests. First, raw ICESat-2 data were processed to obtain accurate ground surfaces. Second, two different methods based on beam pairs were proposed to derive terrain slopes from the ground surfaces. Third, the estimated slopes were validated by airborne LiDAR-derived slopes and compared with SRTM-derived slopes and GDEM-derived slopes. Finally, we further explored the influence of surface topography and ground elevation error on slope estimation from ICESat-2 data. The results show that the ground surface can be accurately extracted from all scenarios of ICESat-2 data, even weak beams in the daytime, which provides the basis for terrain slope retrieval from ICESat-2 beam pairs. The estimated slope has a strong correlation with airborne LiDAR-derived slopes regardless of slope estimation methods, which demonstrates that the ICESat-2 data are appropriate for terrain slope estimation in complex forest environments. Compared with the method based on along- and across-track analysis (method 1), the method based on plane fitting of beam pairs (method 2) has a high estimation accuracy of terrain slopes, which indicates that method 2 is more suitable for slope estimation because it takes full advantage of more ground surface information. Additionally, the results also indicate that ICESat-2 performs much better than SRTM DEMs and GDEMs in estimating terrain slopes. Both ground elevation error and surface topography have a significant impact on terrain slope retrieval from ICESat-2 data, and ground surface extraction should be improved to ensure the accuracy of terrain slope retrieval over extremely complex environments. This study demonstrates for the first time that ICESat-2 has a strong capability in terrain slope retrieval. Additionally, this paper also provides effective solutions to accurately estimate terrain slopes from ICESat-2 data. The ICESat-2 slopes have many potential applications, including the generation of global slope products, the improvement of terrain slopes derived from the existing global DEM products, and the correction of vegetation biophysical parameters retrieved from space-borne LiDAR waveform data.

1. Introduction

The accurate retrieval of Earth’s surface topography is essential for modeling coastal flooding, quantifying glacier elevation change, monitoring global climate change, mapping forest structure, and estimating global vegetation biomass [1,2,3,4,5,6,7,8]. Although several previous studies have realized the importance of terrain slope retrieval and successfully derived accurate terrain slopes in small areas [9,10], it is still challenging to consecutively and accurately estimate the terrain slopes over large areas.

The development of space-borne remote sensing techniques provides an opportunity for the global measurement of terrain slopes. Both optical stereo-photogrammetry and synthetic aperture radar (SAR) have been adopted to generate global digital elevation models (DEMs) (e.g., the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation model (GDEM), TerraSAR-X add-on for digital elevation measurement (TanDEM-X) DEM, Shuttle Radar Topography Mission (SRTM) DEM) [2,5,11,12,13,14], which can be used to calculate terrain slopes at a global scale. However, the vertical accuracies of these DEM products are generally less than 15 m [15,16,17], which significantly limits their application in the accurate retrieval of terrain slopes. In contrast, space-borne LiDAR can provide effective samplings with nearly a global extent for accurately measuring the elevation of the Earth’s surface, making it a promising technology for terrain slope retrieval at a global scale [18,19,20].

As the first earth observation LiDAR system, the Ice, Cloud, and Land Elevation Satellite (ICESat)/Geoscience Laser Altimeter System (GLAS) has been successfully applied for terrain slope estimation [21,22,23,24]. According to previous studies, there are two methods capable of deriving terrain slopes from large-footprint LiDAR waveform data [20,21,25,26,27]. The first type of method can simultaneously calculate along-track slopes and across-track slopes based on repeated tracks, but they fail to obtain accurate small-scale terrain slopes due to the limitation of sparse footprint distribution. The other kind of method can estimate within-footprint terrain slopes by analyzing the received waveforms [20,21,22,28]. Specifically, the within-footprint slopes were obtained by examining the relationship between laser characteristics, surface topography, and waveform feature parameters. This kind of method can receive small-scale terrain slopes to some extent; however, the estimation accuracy is generally limited due to the coupling effect of laser-emitted pulses and surface roughness on slope estimation. Additionally, the canopy returns and ground returns are mixed and not easy to distinguish in complex forest environments, which also induces great uncertainty in slope estimation.

The follow-on mission ICESat-2 was launched in September 2018 to measure the Earth’s surface characteristics [18,29,30,31,32,33]. Different from ICESat, ICESat-2 offers an opportunity to accurately estimate small-scale terrain slopes due to the unique sensor design of dense along-track sampling and beam pairs [18,29,34,35]. Specifically, the along-track distance of GLAS adjacent footprints is approximately 170 m, while that of ICESat-2 is only 70 cm. The dense sampling of ICESat-2 can provide more detailed information about along-track terrain slopes. Additionally, in contrast to ICESat, which only emitted a single beam, ICESat-2 carries a photon-counting LiDAR named the Advanced Topographic Laser Altimeter System (ATLAS), which adopted three pairs of beams. Each pair includes a strong beam and a weak beam, and the beam distance within a pair is approximately 90 m [18,36,37]. Due to the design of beam pairs, ICESat-2 can simultaneously provide abundant information for along-track and across-track slope estimation. In fact, a previous study successfully estimated the local slope of the ice sheet using beam pairs of simulated ICESat-2 data [38]. However, their method has not yet been applied to real ICESat-2 data for testing the algorithm availability. Additionally, no research has been conducted to explore the possibility of ICESat-2 data in terrain slope retrieval over complex forest environments.

In this study, we aimed to propose effective methods to estimate terrain slope from real ICESat-2 data in forest environments. To fulfill this objective, the ICESat-2 data were first processed to obtain accurate ground surfaces. Then, two different methods were proposed to derive the terrain slopes from ICESat-2 beam pairs. Subsequently, the comparisons among ICESat-2 slopes, airborne LiDAR-derived slopes, SRTM-derived slopes, and GDEM-derived slopes were conducted to evaluate our methods’ performance. Finally, we investigated the effects of surface topography and ground elevation error on terrain slope retrieval from ICESat-2 data.

2. Materials

2.1. Study Area

Our study area includes four experimental forest sites: Harvard Forest (HARV), Smithsonian Conservation Biology Institute (SCBI), San Joaquin Experimental Range (SJER), and Soaproot Saddle (SOAP). Two of the experimental sites are located in the western United States, and the other two are located in the eastern United States. The relevant information about the study area, including the site location and ICESat-2 data distribution, is shown in Figure 1. These experimental sites were selected according to the following three criteria: (1) all the data used for slope estimation and validation must be available; (2) these sites covered different topographies ranging from flat to high relief; and (3) the experimental sites encompassed a variety of forest types and coverage.

2.2. ICESat-2 Data

As a follow-on mission of ICESat, ICESat-2 carries a photon-counting LiDAR named ATLAS to measure the elevation of a changing Earth surface. ATLAS employs three-beam pairs to obtain surface elevation profiles, and the pairs are separated by approximately 3.3 km in the across-track direction. Each pair consists of a strong beam and a weak beam, and their energy ratio is approximately 4:1. The pulse repetition rate of each beam is 10 kHz. The ground footprint diameter is 17 m, and the adjacent footprint distance in the along-track direction is ~0.7 m. There are 20 data products (ATL01 to ATL21) for ICESat-2/ATLAS, among which ATL03 is the global geolocated photon product containing the position and elevation information about each photon [18,39]. In this study, we collected ATL03 data from October 2018 to September 2019 on the National Snow and Ice Data Center (NSIDC) website (https://nsidc.org/data/ATL03/versions/3).

2.3. Airborne LiDAR Data

Airborne LiDAR data were used to generate reference slopes for accuracy validation and error analysis. The airborne LiDAR data were collected through the National Ecological Observatory Network (NEON) Airborne Observation Platform [40,41]. According to the NEON LiDAR Algorithm Theoretical Basis Document (ATBD), the maximum horizontal and vertical accuracies of NEON LiDAR data are less than 0.4 m and 0.36 m, respectively, and the elevation datum is referenced to the mean sea level. Given that ICESat-2 data were collected from October 2018 to September 2019, only airborne LiDAR data acquired in 2019 were used in this study to ensure time consistency with ICESat-2 data.

After collecting airborne LiDAR data, a number of data processing steps, including LiDAR filtering, DTM generation, and reference slope estimation, were taken to generate reference slopes. The details are presented as follows:

(1) LiDAR filtering. The ground returns were separated from non-ground returns using a revised progressive triangular irregular network (TIN) densification (PTD) method [42]. Different from the classic PTD method, the revised method builds an improved TIN and changes the original iterative judgment criterions, thus it has an advantage in preserving key ground points and removing non-ground points that belong to low vegetation.

(2) DTM generation. Once ground returns were obtained, a 1 m-resolution digital terrain model (DTM) was generated by interpolating the ground returns into grids.

(3) Reference slope estimation. Reference slopes were calculated based on the DTMs by the plane fitting for accuracy validation.

2.4. Other Data

In this study, both ASTER GDEM and SRTM DEM data were adopted to further evaluate our slope estimation methods’ performance, and the terrain slopes of ASTER GDEM and SRTM DEM data were calculated by the ArcGIS software [43]. As the terrain slope retrieval in forests over large areas is generally based on GDEM and SRTM DEM, we compared the estimated slopes with GDEM and SRTM DEM derived slopes to demonstrate that the ICESat-2 data is superior to the existing global DEMs for retrieving terrain slopes over complex forest environments.

ASTER GDEM is a product from METI and NASA, available from 83°N to 83°S, with a spatial resolution of 30 m [44]. In this study, the second version of GDEM (GDEM2) was used to calculate terrain slopes [2,44]. Compared with the first version of GDEM, GDEM2 has been significantly improved by altering the stereo processing algorithm. However, the horizontal and vertical accuracies of GDEM2 are still very low at ~30 m and ~20 m, respectively.

The 1-arcsecond (30 m) SRTM DEM is the product generated by a collaborative effort of NASA, NGA, DLR, and ASI. In contrast to ASTER GDEM, SRTM DEM only covers 80% of the Earth’s land surface between 60°N and 56°S latitude [2,45,46]. In this study, we downloaded the third version of the SRTM DEM from the United States Geological Survey’s EarthExplorer site (http://earthexplorer.usgs.gov).

3. Methodologies

3.1. ICESat-2 Data Processing and Ground Surface Extraction

A list of ICESat-2 data processing steps, including data filtering, noise removal, ground surface extraction, and ground elevation validation, was conducted in this study to obtain accurate and reliable ground surfaces, which are prerequisites for terrain slope estimation. The details are presented as follows.

(1)

ICESat-2 data filtering. The aim of ICESat-2 data filtering is to obtain photons with high quality by eliminating the cloud effect. Given that the clouds are far above the land surface, ICESat-2 data whose elevation values are more than 60 m above the ASTER DEM were regarded as photons affected by clouds and excluded.

(2)

Noise photon removal. In this study, the signal photons were separated from noise photons using the noise removal algorithm proposed by our previous studies [37,47]. This algorithm includes two steps. First, elevation frequency histograms were built to remove visible noise photons. Second, density distribution histograms were built to eliminate the remaining noise photons. In prior to establishing density distribution histograms, the density of each photon should be calculated based on photon number in the elliptical neighborhood. A more detailed description of the noise removal algorithm can be found in Zhu et al. [37,47].

(3)

Ground surface extraction. There are four key steps to implement ground surface retrieval according to our previous studies [37,47], including initial ground photon extraction, iterative detection of ground photons, ground photon densification, and ground surface fitting. First, a local elevation frequency histogram was built, and the photon at the lowest elevation peak with the maximum density was regarded as the initial ground photon. Second, the empirical mode decomposition (EMD) algorithm was utilized to eliminate wrong ground photons such as canopy photons and noise photons. Third, the number of accurate ground photons was increased based on progressive densification. Finally, all ground photons and corresponding ground surfaces were retrieved utilizing cubic spline interpolation. A more detailed description of the ground surface extraction algorithm can be found in Zhu et al. [47].

(4)

Ground elevation validation. Once the ground surface was obtained, the airborne LiDAR-derived DTM was adopted to assess the accuracy of the ground elevation. Compared with strong beams, fewer signal photons in weak beams may reduce the accuracy of the estimated ground elevation. Additionally, photons acquired in the daytime may also introduce more solar background noises than that in the nighttime. To examine the capability of ICESat-2 data in the ground elevation estimation, the ICESat-2 data in four different scenarios were used for ground elevation estimation and validation: 1) ICESat-2 strong beams collected in the daytime, 2) ICESat-2 weak beams collected in the daytime, 3) ICESat-2 strong beams collected in the nighttime, and 4) ICESat-2 weak beams collected in the nighttime. The validation results will guide the slope estimation of ICESat-2 data.

3.2. Slope Estimation

In this study, two methods based on beam pairs were proposed to derive terrain slopes from ground surfaces. One method derives the terrain slopes based on along- and across-track analysis and the other method calculates the terrain slopes by the plane fitting.

3.2.1. Slope Estimation Method Based on Along- and Across-Track Analysis

In the method based on along- and across-track analysis, the terrain slope can be divided into the along-track slope ${\phi }_a$ and across-track slope ${\phi }_c$. Only the along-track ground photons are used to calculate the along-track slope ${\phi }_a$ by linearly fitting the photon elevations, while the retrieval of the across-track slope ${\phi }_c$ relies on both the strong and weak beams. Once ${\phi }_a$ and ${\phi }_c$ are obtained, the terrain slope $\phi$ can be subsequently calculated using the following Equation.

$$\phi =\arctan (\sqrt{{\tan }^2({\phi }_a)+{\tan }^2({\phi }_c)})$$

(1)

As shown in Figure 2, there are a pair of beams (red for beam A and green for beam B), and photon p is located in beam A. To determine the terrain slope of photon p, we should first calculate both the along-track slope ${\phi }_a(p)$ and across-track slope ${\phi }_c(p)$ of photon p. There are three steps to derive the along-track slope of photon p. First, we obtain all the photons that are located in beam A and in the 45-m neighborhood of photon p. Second, the elevations of all the obtained photons are linearly fitted by line 1, as shown in Equation (2). Finally, ${\phi }_a(p)$ is calculated according to Equation (3).

To calculate the across-track slope of photon p, we should first determine the position of the interpolated photon p′ in beam B on the basis of the criterion that line p-p′ is perpendicular to beam A and beam B (see Figure 2c). Then, the elevation of photon p′ is obtained by inverse distance weighted (IDW) interpolation. Once the position and elevation of photon p′ is determined, the across-track slope can finally be calculated based on the elevations of photons p and p′, as shown in Equation (4).

In this study, the slope method based on along- and across-track analysis (method 1) can be further divided into three different methods: method1_All, method1_Strong, and method1_Weak. The method1_Strong means that only the terrain slopes of strong beams were calculated by method 1, while the method1_Weak estimates the terrain slopes of weak beams. In contrast to method1_Strong and method1_Weak, method1_All extracted terrain slopes from both the strong and weak beams. For each slope estimation method, we calculated the corresponding slope value at an interval of ~90 m in the along-track direction given that the distance between beam pairs is approximately 90 m.

$$z=f(x)=ax+b$$

(2)

$${\phi }_a(p)=\arctan (a)$$

(3)

$${\phi }_c(p)=\arctan \left[\frac{H(p^{\prime })-H(p)}{L}\right]$$

(4)

where z is the photon elevation, x is the along-track distance, a and b are the slope coefficient and the intercept, respectively, H(p′) refers to the elevation of interpolated photon p′, and H(p) refers to the elevation of photon p, and L is the distance between beam A and beam B, which is approximately 90 m.

3.2.2. Slope Estimation Method Based on Plane Fitting of Ground Track Pairs

This method (named method 2) calculates the terrain slopes from ground track pairs by the plane fitting. It includes four steps. First, the centerline of the ground track pairs is extracted based on the criterion that the distances from the centerline to each ground track are equal. Second, a list of center points are sampled from the centerline by setting the point interval of 90 m, and then we obtain all photons within the grid of 90 m for each center point. Third, the photon elevations are expressed as a function of the along-track distance and across-track distance and fitted to a plane. Finally, the terrain slope of each center point is calculated based on the fitted plane. As shown in Figure 3, there are a pair of ground beams (red for beam A and green for beam B). Line 2, which is regarded as the centerline, is extracted based on the beam pairs. To estimate the terrain slope of center point O, we first obtain all the photons within the dotted rectangle; then a plane (see Equation (5)) is used to fit the photon elevations; the terrain slope of point O is finally calculated based on Equation (6).

$$z=g(x,y)=a_{x}x+a_{y}y+b$$

(5)

$$\phi =\arctan (\sqrt{a_x{}^2+a_y{}^2})$$

(6)

where z is the photon elevation, x is the along-track distance, y is the across-track distance, b is the intercept, and a_x and a_y are the slope coefficients.

3.3. Accuracy Validation

In this study, airborne LiDAR-derived slopes were adopted to check the validity of our slope estimation methods. Four indicators, including the coefficient of determination (R²), error bias, standard deviation (SD), and root-mean-square error (RMSE), were used to quantitatively validate the accuracies of terrain slope retrievals from ICESat-2 data. Additionally, the comparison among ICESat-2 slope estimates, GDEM-derived slopes, and SRTM-derived slopes was also conducted further to evaluate the performance of our slope estimation methods.

3.4. Error Analysis

Both surface topography and ground elevation error have a significant impact on terrain slope retrieval from ICESat-2 data. Therefore, it is essential to explore the effects of the factors mentioned above on slope estimation. In this study, error analysis was conducted by investigating the changes in terrain slope errors concerning terrain slope and ground elevation error.

4. Results

4.1. Validation of the Ground Elevation Accuracy

In this study, both qualitative and quantitative assessments were conducted to evaluate the ground surface extraction from ICESat-2 data. The qualitative analysis was made by visual inspection. Figure 4 shows the qualitative results of ground surface retrievals from four different scenarios of ICESat-2 data. From Figure 4, we can see that the ground surfaces can be extracted from all ICESat-2 data, even the weak beams in the daytime.

To further evaluate the performance of ICESat-2 data in extracting the ground surface, we validated the accuracy of the ground elevation with airborne LiDAR-derived DTM. The quantitative statistics of ground elevation errors are summarized in

Table 1. The statistical errors of ground elevation using four different scenarios of ICESat-2 data.

Datasets	R2	Bias (m)	SD (m)	RMSE (m)
Strong beams in the daytime	1.000	−2.563	3.112	4.031
Weak beams in the daytime	1.000	−2.477	4.464	5.105
Strong beams in the nighttime	1.000	−1.388	2.037	2.465
Weak beams in the nighttime	1.000	−1.638	3.211	3.605
All datasets	1.000	−1.816	3.056	3.555

, and Figure 5 shows the frequency histogram of ground elevation errors. According to Table 1 and Figure 5, strong consistency was observed between the estimated ground elevation with the airborne LiDAR-derived DTM, as indicated by high R² and low RMSE for all scenarios of ICESat-2 data. Additionally, Table 1 and Figure 5 also indicate that the strong beams have higher estimation accuracy of the ground elevation than that of the weak beams, and ICESat-2 data at night perform much better than those in the daytime for retrieving ground elevation due to higher R² and lower SD and RMSE.

4.2. The Assessment of ICESat-2-Derived SLOPES

To evaluate our methods’ performance, we compared the estimated terrain slopes with the airborne LiDAR-derived terrain slopes, and the relevant results are shown in

Table 2. The statistical slope errors using different methods across various experimental sites.

Sites	Methods	Number	Bias (°)	SD (°)	RMSE (°)
HARV	method1_All	391	0.891	2.174	2.347
	method1_Strong	201	0.630	2.123	2.210
	method1_Weak	190	1.168	2.199	2.485
	method 2	201	0.805	1.332	1.553
SCBI	method1_All	1204	−0.086	2.870	2.870
	method1_Strong	602	−0.116	2.879	2.879
	method1_Weak	602	−0.056	2.863	2.862
	method 2	602	−0.270	2.108	2.124
SJER	method1_All	1752	−0.399	3.304	3.327
	method1_Strong	876	−0.558	3.146	3.193
	method1_Weak	876	−0.240	3.450	3.456
	method 2	876	−0.780	2.128	2.265
SOAP	method1_All	1564	−0.773	4.327	4.394
	method1_Strong	782	−0.617	4.201	4.243
	method1_Weak	782	−0.929	4.446	4.539
	method 2	782	−0.769	2.828	2.929
All	method1_All	4911	−0.339	3.528	3.543
	method1_Strong	2461	−0.371	3.412	3.432
	method1_Weak	2450	−0.306	3.640	3.652
	method 2	2461	−0.522	2.363	2.419

and Figure 6. According to Table 2 and Figure 6, a strong positive correlation was observed between the estimated slopes and LiDAR-derived slopes due to the high R² and low RMSE for all slope estimation methods, and the slope estimation accuracies differ across various experimental sites. Additionally, there is a slight underestimation of ICESat-2-derived slopes compared with LiDAR-derived slopes regardless of estimation methods. Among all methods, method 2 performed best, followed by method1_Strong, and method1_Weak performed worst in terrain slope retrieval as indicated by Table 2 and Figure 6.

Additionally, we also compared the slope estimation accuracy of ICESat-2 data in the daytime with that of ICESat-2 data in the nighttime, as shown in

Table 3. The statistical slope errors of our methods using different ICESat-2 data.

Data Time	Methods	Number	Bias (°)	SD (°)	RMSE (°)
Daytime	method1_All	1462	−0.752	4.311	4.374
	method1_Strong	731	−0.582	4.198	4.236
	method1_Weak	731	−0.922	4.417	4.509
	method 2	731	−0.769	2.831	2.932
Nighttime	method1_All	3449	−0.163	3.121	3.125
	method1_Strong	1730	−0.282	3.016	3.028
	method1_Weak	1719	−0.043	3.220	3.220
	method 2	1730	−0.418	2.126	2.166

. From Table 3, we know that the accuracies of terrain slopes extracted by ICESat-2 data in the daytime are lower than those in the nighttime regardless of slope estimation methods.

To further evaluate the performance of ICESat-2 in estimating terrain slopes, we also compared the estimated slopes with SRTM-derived slopes and GDEM-derived slopes. The relevant results are shown in

Table 4. The comparison results among ICESat-2 estimated slopes, SRTM-derived slopes, and GDEM-derived slopes.

Methods	Slope Errors	ICESat-2 Slopes	SRTM-Derived Slopes	GDEM-Derived Slopes
method1_All	R2	0.736	0.676	0.592
	Bias	−0.339	−0.928	1.125
	SD	3.528	3.883	4.822
	RMSE	3.543	3.992	4.951
method1_Strong	R2	0.754	0.692	0.607
	Bias	−0.371	−0.935	1.062
	SD	3.412	3.770	4.705
	RMSE	3.432	3.884	4.822
method1_Weak	R2	0.717	0.660	0.577
	Bias	−0.306	−0.92	1.194
	SD	3.640	3.994	4.942
	RMSE	3.652	4.097	5.083
method 2	R2	0.872	0.681	0.573
	Bias	−0.522	−0.868	1.064
	SD	2.363	3.766	4.837
	RMSE	2.419	3.864	4.952

and Figure 7. According to Table 4 and Figure 7, our new proposed methods perform much better than ASTER GDEM and SRTM DEM in estimating terrain slopes due to the higher R², lower SD, and RMSE.

4.3. Error Analysis of ICESat-2 Derived Slopes

4.3.1. The Effect of Surface Topography on Terrain Slope Retrieval

To further explore the influence of surface topography on terrain slope estimation from ICESat-2 data, we investigated the relationship between slope errors and terrain slope. In this paper, ICESat-2 data were divided into four groups according to the slope values with intervals of 5°. The statistical slope errors of the different methods for each group are summarized in

Table 5. Statistical errors of different methods, classified by slope.

Methods	Slope Errors (°)	Slope (°)
		0–5	5–10	10–15	>15
method1_All	Number	1217	1533	1097	1064
	Bias	1.208	−0.076	−0.768	−2.043
	SD	2.340	2.842	3.451	4.621
	RMSE	2.633	2.842	3.534	5.051
method1_Strong	Number	573	817	553	518
	Bias	1.002	−0.172	−0.672	−1.885
	SD	2.330	2.957	3.436	4.294
	RMSE	2.535	2.960	3.498	4.685
method1_Weak	Number	644	716	544	546
	Bias	1.392	0.033	−0.866	−2.194
	SD	2.335	2.702	3.467	4.911
	RMSE	2.717	2.700	3.570	5.374
method 2	Number	617	782	559	503
	Bias	0.465	−0.405	−0.884	−1.513
	SD	1.462	1.937	2.417	3.171
	RMSE	1.533	1.977	2.572	3.511

. Additionally, the relationship between the RMSE of slope errors and terrain slope is shown in Figure 8. According to Table 5 and Figure 8, the slope errors show an increasing trend with terrain slope regardless of slope estimation methods. This means that the surface topography has a significant impact on terrain slope retrievals from ICESat-2 data.

4.3.2. The Effect of Ground Elevation Error on Terrain Slope Estimation

Table 6. Statistical errors of different methods, classified by ground elevation error.

Methods	Slope Errors (°)	RMSE of Ground Elevation Errors (m)
		0–1	1–2	2–3	>3
method1_All	Number	490	1556	1133	1732
	Bias	0.313	−0.320	−0.543	−0.406
	SD	2.120	3.032	3.597	4.141
	RMSE	2.141	3.048	3.637	4.160
method1_Strong	Number	260	877	628	696
	Bias	0.010	−0.366	−0.383	−0.511
	SD	2.225	3.114	3.536	3.974
	RMSE	2.220	3.133	3.554	4.004
method1_Weak	Number	230	679	505	1036
	Bias	0.656	−0.261	−0.743	−0.335
	SD	1.943	2.925	3.665	4.251
	RMSE	2.047	2.935	3.736	4.262
method 2	Number	286	664	462	801
	Bias	−0.096	−0.392	−0.663	−0.614
	SD	1.276	1.885	2.370	2.899
	RMSE	1.278	1.924	2.459	2.962

shows the statistical errors of different methods for each group classified by ground elevation error. Figure 9 shows the relationship between the RMSE of slope errors and ground elevation error. From Table 6 and Figure 9, we can see that the slope errors increase with an increase in the ground elevation error for all slope estimation methods.

5. Discussion

5.1. Ground Surface Retrieval

Our study demonstrated that both strong and weak beams have a strong capability of accurately extracting ground surfaces, regardless of ICESat-2 data in the daytime or nighttime. This is in agreement with the findings of [48]. In their study, all the ICESat-2 datasets perform well in extracting ground surfaces, which lays the foundation for terrain slope estimation based on beam pairs of ICESat-2 data. Additionally, results also indicated that the strong beams have higher estimation accuracy of the ground elevation than that of the weak beams, and ICESat-2 data at night perform much better than those in the daytime for ground surface retrieval. It is well known that the strong beams with higher energy can obtain denser signal photons than weak beams, resulting in a more detailed and accurate description of the ground surface. Compared with ICESat-2 data in the night, more solar background noises would be induced to the ICESat-2 data in the daytime, thereby significantly affecting the accuracy of ground elevation estimation.

5.2. The Slope Estimation from ICESat-2 Data

The ICESat-2 data show strong capability in deriving terrain slopes regardless of estimation methods, which demonstrated that terrain slopes could be accurately estimated from ICESat-2 data even in complex forest environments. Besides, the slope estimation accuracies differ across various experimental sites, which can be explained by the fact that both forest coverage and surface topography may be quite different across various experimental sites.

Compared with LiDAR-derived slopes, a slight underestimation of ICESat-2-derived slopes was observed regardless of estimation methods. This is maybe due to the large footprint size of ICESat-2 data. In this study, the resolution of LiDAR DTM is 1 m, while the footprint size of ICESat-2 data is large of 17 m. The large size of ICESat-2 data smooths or flattens surface topography, hence resulting in slope underestimation, whereas the smaller DTM cells are more sensitive to surface micro-topography.

Results indicated that terrain slopes could be better extracted from ICESat-2 data in the nighttime compared to ICESat-2 data in the daytime. The reason may be that ICESat-2 data in the daytime introduced abundant solar background noises, which significantly affected the extraction of the ground surface, finally inducing more substantial errors in the terrain slope retrievals. In practical applications, we can only select ICESat-2 nighttime datasets for slope estimation to ensure the terrain slopes’ overall accuracies.

Among all the slope estimation methods, method 2 performed best, followed by method1_Strong, and method1_Weak performed worst in terrain slope retrieval. This may be because the slope estimation accuracy can be improved by using more detailed ground surface information. Specifically, method 1 mainly utilizes one ground track to obtain terrain slopes, while method 2 takes full advantage of both strong and weak ground tracks to derive terrain slopes; compared with method1_Weak, method1_Strong can obtain more accurate along-track slopes because strong beams can obtain more effective ground photons to derive along-track slopes.

5.3. Error Analysis of ICESat-2 Derived Slopes

In this study, we explored the effect of surface topography on slope estimation using ICESat-2 data and drew a conclusion that the slope errors show a strong correlation with the terrain parameter. This means that the surface topography has a significant impact on terrain slope retrievals from ICESat-2 data. A possible reason is that ground surface extraction is significantly affected by surface topography, as demonstrated by our previous study [37], particularly in forest environments with complex terrain. Given that the slopes were estimated based on ICESat-2-derived ground surfaces, larger uncertainty will be induced with increasing terrain slopes. In future studies, we should focus on ground surface extraction over extreme terrains to improve the accuracy of terrain slope retrieval from ICESat-2 data. Additionally, the results also indicated that the slope estimation from ICESat-2 data is highly dependent on ground elevation accuracy. To ensure the reliability of ICESat-2 estimated slopes, we should propose effective methods to improve the estimation accuracy of ground elevation in the case of complex forests or low signal-to-noise ratio (SNR).

5.4. The Potential Applications of ICESat-2 Sopes

As demonstrated above, ICESat-2 slopes are accurate and reliable, thus showing great potential in many applications. Below, we presented the potential applications of ICESat-2 slopes.

(1) The generation of the global terrain slope product.

Although ICESat-2 data can provide accurate terrain slope samplings with global coverage, ICESat-2 slopes are not continuously distributed, limiting their applications in many fields. To generate global terrain slope products, ICESat-2 slopes can be directly interpolated into grids or combined with other continuously distributed remote sensing data using regression analysis.

(2) The improvement of terrain slopes derived from existing global DEM products.

The existing global DEM products, such as GDEM, SRTM DEM, and TanDEM-X DEM, can be used to calculate terrain slopes. However, their estimation accuracies of terrain slopes are generally limited due to the low vertical accuracies of global DEM products. Additionally, these global DEM products generally represent canopy-top surfaces rather than true ground surfaces in forests. Hence, these existing DEM products cannot provide reliable terrain slopes. Given that ICESat-2 can provide effective sampling data with nearly global spatial coverage for accurate retrieval of terrain slopes, ICESat-2 slopes can be adopted as reference slopes to improve the accuracy of terrain slope retrieval from existing global DEM products.

(3) The correction of space-borne LiDAR-derived vegetation biophysical parameters.

The retrieval of vegetation biophysical parameters, including forest height, leaf area index (LAI), and biomass, may be significantly affected by surface topography. To reduce the influence of surface topography on the retrieval of vegetation biophysical parameters from space-borne LiDAR waveform data, terrain slopes should be provided as an input parameter of terrain correction models or auxiliary information. In previous studies, terrain slopes are generally derived from existing global DEM products, even though they are inaccurate. In future studies, the ICESat-2 slopes can be used to correct the vegetation biophysical parameters that are derived from space-borne LiDAR waveform data, such as GLAS and GEDI.

5.5. The Innovations and Limitations

Although a previous study has made an initial attempt to derive the local slope of ice sheet from simulated ICESat-2 data [38], no research has been conducted to explore the potential of real ICESat-2 data in terrain slope retrieval, especially over complex forest environments. Compared with the previous study, our paper has three main innovations:

(1) This paper is the first attempt to derive terrain slopes in forests using ICESat-2 data, providing a new technique for terrain slope retrieval over forest environments.

(2) Two effective methods, which make full use of the unique design of beam pairs, were proposed to accurately estimate terrain slopes from ICESat-2 data. Given ICESat-2 can provide effective samplings at a global extent, the ICESat-2 slopes will lay the foundation for the generation of global terrain slope product and the correction of vegetation biophysical parameters retrieved from space-borne LiDAR waveform data.

(3) This paper not only proposed effective solutions, but also carried out experiments and analysis, illuminating that the proposed methods are capable of obtaining accurate terrain slopes, and perform much better than the existing global DEM products (ASTER GDEM, SRTM DEM) in estimating terrain slopes over forest environments. Additionally, we also quantitatively investigated the effects of surface topography and ground elevation error on terrain slope retrieval from ICESat-2 data.

Although several innovations have been made in this study, there is also a limitation that we fail to accurately extract ground surface form low SNR ICESat-2 data or over extremely complex forest environments, significantly affecting the accuracy of terrain slope estimation. In future studies, we should propose effective methods to improve the estimation accuracy of ground elevation in the case of extremely complex forests or low SNR.

6. Conclusions

In this paper, we explored the possibility of ICESat-2 data for estimating terrain slopes over complex forest environments. Two methods based on beam pairs were first proposed to derive terrain slopes from ICESat-2 data. Then, the slope estimates were compared with airborne LiDAR-derived slopes, GDEM-derived slopes, and SRTM-derived slopes to evaluate the performance of our methods. Finally, we investigated the effects of surface topography and ground elevation error on terrain slope retrieval from ICESat-2 data. Based on the experimental results, we drew the following conclusions: (1) both strong and weak beams are capable of accurately extracting ground surfaces regardless of using ICESat-2 daytime or nighttime data; (2) the slope estimation accuracies of our methods are high, demonstrating that the ICESat-2 data are suitable for slope estimation in forests; (3) compared with method 1, method 2 can achieve terrain slopes with higher accuracy because it takes full advantage of ground surface information; (4) ICESat-2 data perform better than existing global DEM products in estimating terrain slopes over complex forest environments; and (5) both ground elevation error and surface topography have a significant impact on terrain slope retrieval from ICESat-2 data.

Overall, this paper was the first attempt to derive terrain slopes from ICESat-2 data over complex forest environments, which demonstrated the capability of ICESat-2 data in estimating terrain slope. Additionally, this study provided an effective solution for terrain slope retrieval from ICESat-2 data, which will lay the foundation for the generation of global terrain slope products. As terrain slope estimation is affected by ground elevation error and surface topography, how to propose effective methods for ground surface extraction over extremely complex forest environments or from low-SNR ICESat-2 data has become an emphasis on future research. In addition, our next plan will also focus on the potential applications of ICESat-2 slopes.