Application and Validation of a Model for Terrain Slope Estimation Using Space-Borne LiDAR Waveform Data

Abstract

The terrain slope is one of the most important surface characteristics for quantifying the Earth surface processes. Space-borne LiDAR sensors have produced high-accuracy and large-area terrain measurement within the footprint. However, rigorous procedures are required to accurately estimate the terrain slope especially within the large footprint since the estimated slope is likely affected by footprint size, shape, orientation, and terrain aspect. Therefore, based on multiple available datasets, we explored the performance of a proposed terrain slope estimation model over several study sites and various footprint shapes. The terrain slopes were derived from the ICESAT/GLAS waveform data by the proposed method and five other methods in this study. Compared with five other methods, the proposed method considered the influence of footprint shape, orientation, and terrain aspect on the terrain slope estimation. Validation against the airborne LiDAR measurements showed that the proposed method performed better than five other methods (R² = 0.829, increased by ~0.07, RMSE = 3.596°, reduced by ~0.6°, n = 858). In addition, more statistics indicated that the proposed method significantly improved the terrain slope estimation accuracy in high-relief region (RMSE = 5.180°, reduced by ~1.8°, n = 218) or in the footprint with a great eccentricity (RMSE = 3.421°, reduced by ~1.1°, n = 313). Therefore, from these experiments, we concluded that this terrain slope estimation approach was beneficial for different terrains and various footprint shapes in practice and the improvement of estimated accuracy was distinctly related with the terrain slope and footprint eccentricity.

1. Introduction

The Earth surface survey provides fundamental and useful geo-information for terrestrial ecosystems, global climate monitoring, and landform mapping [1,2,3,4]. We extremely focus on quantifying Earth’s surface characteristics (e.g., elevation, slope, and aspect) because these characteristics provide valid indicators for the climate, precipitation, and ocean currents [5,6,7]. However, it is still a major challenge to monitor the surface topography consecutively and accurately over large areas. Development of space-borne remote sensing techniques has addressed this challenge to some extent [8,9,10]. For instance, the optical stereo photogrammetry and Synthetic Aperture Radar have generated the available global digital elevation models (DEMs) [10,11]. However, complex matching operations must be done to produce DEMs from stereos and the optical images are susceptible to the weather. Moreover, the vertical accuracies of these DEMs are not high at ~15 m [12,13,14]. In contrast, LiDAR sensors can directly acquire the surface height information by measuring the time of the transmitted and received pulses [15,16]. Its high-frequency sampling makes great promises of high vertical accuracy [17,18,19,20,21,22]. Therefore, the space-borne LiDAR has become the most promising technique for accurately measuring the terrain characteristics on a global scale [21,22].

Ice, Cloud, and land Elevation Satellite (ICESat)/Geoscience Laser Altimeter System (GLAS) is the first space-borne LiDAR mission of recording global full-waveform altimetry data within the footprint [21,23]. Effective global or region GLAS products have been applied to measure the land topography and ice sheet mass balance [24,25,26]. Terrain slope is one of the most important surface characteristics. Some previous studies have also emphasized the significance of the terrain slope estimation. For instance, the within-footprint terrain slope has an obvious impact on most waveform characteristics because its broadening effect mixes the object return and the ground return [27,28,29].

To more accurately estimate the within-footprint terrain slope and to understand the complex relation between waveform characterization and surface terrain, several proposed physical models and methods have achieved the terrain slope estimation by using the full-waveform LiDAR data [30,31,32]. Among these methods, most have been built by simplifying the footprint shape into a circle [30,31]. They ignored the fact that the true footprint shape was closer to be elliptical rather than circular [29,33,34]. In contrast, this phenomenon was taken into accounts in the study of Nie et al. [32]. They quantified the influence of footprint size, shape, orientation, and terrain aspect on terrain slope estimation by a series of theoretical physically-based analyses. A subsequent slope estimation model was proposed to reduce the impact of these various factors. This estimation approach has been proven to be theoretically feasible on estimating the terrain slope within the footprint more accurately, but it still needs to be applied and validated on a large area with a different terrain and various LiDAR footprints.

Therefore, the primary goal in this study is to practice the terrain slope estimation model of Nie et al. [32] over several study sites based on GLAS data. We validated the GLAS-estimated slopes (by the proposed method considering footprint shape, orientation, and terrain aspect and five other methods without considering them) with respect to comparisons against the airborne LiDAR data. After demonstrating the efficacy of the terrain slope estimation approach, we attempted to analyze the effects of different factors (e.g., terrain slope, footprint size, and eccentricity) on the terrain slope estimation by relating these factors to the actual absolute errors. The accuracy improvement and applicability of these terrain slope estimation methods were simply investigated based on the comparison results in different terrain and various footprint shapes. In addition, some other error sources (e.g., terrain condition and surface covering condition) were also analyzed in this work.

2. Materials

2.1. Study Area

We selected four study sites where all used datasets were available for the terrain slope estimation. Two sites are located in Greenland and the other two are located in Antarctica. These areas are barely covered by snow and vegetation. Temporal variations in topography in these areas are relatively small. The sites covered different topographies ranging from flat to high relief. Some characteristics of four study sites were summarized in

Table 1. Summary of site locations and terrain characteristics. Elevation is the study site elevation range, θ_range is the site slope range, and $\overline{\mathsf{\theta }}$ is the site mean slope (source: ASTER GDEM).

Site	Latitude	Longitude	Elevation	${\mathsf{\theta }}_{\mathbf{range}}$	$\overline{\mathsf{\theta }}$	Terrain	Location
1	69.2°N–76.7°N	46.4°W–63.3°W	611~2023 m	0~24°	4	Moderate relief	Greenland
2	65.8°N–69.3°N	32.1°W–50.9°W	28~3192 m	0~57°	11	High relief	Greenland
3	66.8°S–73.0°S	61.8°W–68.8°W	0~2058 m	0~13°	2	Low relief	Antarctica
4	67.8°S–71.0°S	65.8°W–75.6°W	0~2033 m	0~32°	6	Moderate relief	Antarctica

. Geographical locations and topographic information were shown in Figure 1.

2.2. GLAS Data

GLAS is the first space-borne laser altimeter for Earth observation, which was launched by the National Aeronautics and Space Administration (NASA) in 2003. It emitted a temporal Gaussian laser pulse for a duration of approximately 4 to 6-ns FWHM (full width half maximum) [35]. Return signals were recorded in the form of a continuous waveform reflecting multiple information (e.g., height, density, and reflectivity) of the intercepted surfaces [35,36]. The emitted laser pulses took shape of the footprints with around a 70-m diameter and a 170-m center-to-center spacing. The horizontal position accuracy of footprints is about 4.8 m and the vertical accuracy is 13.8 cm [37].

National Snow and Ice Data Center (NSIDC) provides 15 types of GLAS product (GLA01~GLA15) [38] of which GLA01, GLA05, and GLA12 were used in this study. The GLA01 product records the raw full-waveform data including the transmitted and received waveform. The GLA05 product provides the detailed description of footprint size, shape, and orientation. The GLA12 product is used to obtain laser spot geo-location and surface elevation.

Table 2. Summary of the GLAS parameters used in this study.

GLA01 Code	Description
i_rec_ndx	GLAS record index
i_tx_wf	Sampled transmit pulse waveform
i_rng_wf	Sampled received pulse waveform
i_4nsBgMean	Background mean value
i_4nsBgSDEV	Background standard deviation
GLA05 Code
i_tpazimuth	Transmit pulse azimuth
i_tpeccentricity	Transmit pulse eccentricity
i_tpmajoraxis	Transmit pulse major axis
GLA12 Code
i_lat	Coordinate data, latitude
i_lon	Coordinate data, longitude
i_elev	Ice sheet surface elevation

summarized the parameters of GLA01, GLA05, and GLA12 used in the terrain slope estimation. In this study, the GLAS datasets of Greenland sites were from Laser 3I in October 2007 and the datasets of Antarctica were from Laser 3K in October 2008. Each GLAS footprint is a variable ellipse with a slightly different orientation, size, and eccentricity. The average footprint sizes (minor and major axes) are 46 m × 57 m for 3I campaign and 41 m × 52 m for 3K campaign.

We screened the GLAS data to obtain a high-quality dataset and ensure the validity of the experiments. Invalid GLAS data without geographical coordinates or affected by clouds (the flag i_FRir_qaFlag ≠ 15 in the GLA12 product) was removed [31,33].

2.3. ASTER GDEM Data

ASTER GDEM (Advanced Space-borne Thermal Emission and Reflection Radiometer Global Digital Elevation Model) was used to calculate the prior terrain aspect, which was combined with the footprint orientation to flexibly choose the terrain slope estimation method. ASTER, which is onboard with the NASA Terra, collected the in-track stereos using nadir-looking and aft-looking near infrared cameras [39,40]. These stereos were used to produce the global digital elevation model. The GDEM data is available from 83° S to 83° N and covers most of the polar region [41]. We used Version 2 of GDEM in this study [42]. Although changes in the stereo processing have produced significant improvements in GDEM2 as compared to GDEM1, the validated horizontal accuracy of GDEM2 is still ~30 m and the vertical accuracy is ~20 m [43,44,45]. Therefore, in order to avoid introducing a new error source, the terrain aspect calculated from GDEM2 was only used as an auxiliary parameter but not an input parameter of the estimation method.

The terrain aspect was calculated as the direction of the maximum rate of change in the elevation of GDEM. The extraction of the terrain aspect was conducted by using the mapping and analytics software ArcGIS 10.2.

2.4. ATM Data

Airborne Topographic Mapper (ATM) data was used to verify the accuracy of the terrain slope estimation method in this study. It was developed by NASA to observe the Earth’s topography and, most importantly, to measure the changes of ice sheets and glaciers in North America, Greenland, and Antarctica [46,47,48]. The L1B ATM data records the spot elevation and return strength measurements over the sea ice and ice sheet surfaces where the vertical accuracy is 6.6 cm [49]. The L2 data comprises the elevation, slope, and roughness datasets that were calculated from L1B data. The Preicebridge ATM data was utilized in this paper because of the close collected times with other datasets. The acquisition dates of the ATM data, GLAS data, and ASTER GDEM are almost identical (

Table 3. Summary of ATM dataset, GLAS waveform dataset, and ASTER GDEM dataset in four sites of this study.

Site	Number of Estimates	ATM Acquisition Date	GLAS Acquisition Date	ASTER GDEM Acquisition Date
Site 1	188	14 September 2007	October 2007	October 2007
Site 2	304	23 September 2007	October 2007	October 2007
Site 3	227	14 October 2008	October 2008	October 2008
Site 4	139	30 October 2008	October 2008	October 2008

).

3. Methods

3.1. Terrain Slope Estimation Method

Below, we gave an outline of the terrain slope estimation method (referred to as the flexible method) and tested it on the GLAS waveform across four study sites.

The within-footprint terrain slope can be quantitatively inversed by the ground vertical extent and footprint diameter, which is shown in Equation (1). Five different footprint diameters were adopted for calculating the terrain slope (Equations (2)–(6)). They are the semi-major axis a, the semi-minor axis b, the arithmetic mean, the quadratic mean, and the geometrical mean of a and b. A series of theoretical analyses [32] indicated that the terrain slope error estimated by these five diameters varied regularly with the intersection angle between the footprint orientation and terrain aspect (Figure 2). Every diameter was found to be optimal only for environments with some range of intersection angles. Therefore, we chose the flexible slope estimation method based on the intersection angle once the four threshold values (${\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4$ in Figure 2) were obtained, which is shown in Equation (7).

$$\eta =\mathrm{atan}(\frac{h}{d})$$

(1)

$$\mathrm{Method}1:{\eta }_1=\mathrm{atan}(\frac{h}{2a})$$

(2)

$$\mathrm{Method}2:{\eta }_2=\mathrm{atan}(\frac{h}{2b})$$

(3)

$$\mathrm{Method}3:{\eta }_3=\mathrm{atan}(\frac{h}{a+b})$$

(4)

$$\mathrm{Method}4:{\eta }_4=\mathrm{atan}(\frac{h}{2\sqrt{ab}})$$

(5)

$$\mathrm{Method}5:{\eta }_5=\mathrm{atan}(\frac{h}{2\sqrt{\frac{a^2+b^2}{2}}})$$

(6)

$$\mathrm{Flexible}\mathrm{Method}:\eta =\{\begin{array}{l}{\eta }_1,\mathrm{if}0\le \theta <{\theta }_1\mathrm{or}\pi -{\theta }_1\le \theta <\pi \\ {\eta }_2,\mathrm{if}{\theta }_4\le \theta <\pi -{\theta }_4\\ {\eta }_3,\mathrm{if}{\theta }_2\le \theta <{\theta }_3\mathrm{or}\pi -{\theta }_3\le \theta <\pi -{\theta }_2\\ {\eta }_4,\mathrm{if}{\theta }_3\le \theta <{\theta }_4\mathrm{or}\pi -{\theta }_4\le \theta <\pi -{\theta }_3\\ {\eta }_5,\mathrm{if}{\theta }_1\le \theta <{\theta }_2\mathrm{or}\pi -{\theta }_2\le \theta <\pi -{\theta }_1\end{array}$$

(7)

where $h$ is the ground vertical extent, $d$ is the footprint diameter, $\eta$ is the estimated terrain slope, $\theta$ is the intersection angle between the terrain aspect and footprint orientation, and $a$ and $b$ are the semi-major and semi-minor axes of each footprint, respectively.

The four threshold values of the intersection angle used for choosing the estimation method can be calculated by Equations (8)–(11). The left and right sides of each equation represented the theoretical absolute errors by two methods of the left and right sides of the threshold value in Figure 2, respectively.

$${\theta }_1:\mathrm{atan}(\frac{h}{2a})-\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_1)+b^2{\sin }^2({\theta }_1)}})=\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_1)+b^2{\sin }^2({\theta }_1)}})-\mathrm{atan}(\frac{h}{2\sqrt{\frac{a^2+b^2}{2}}})$$

(8)

$${\theta }_2:\mathrm{atan}(\frac{h}{2a})-\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_2)+b^2{\sin }^2({\theta }_2)}})=\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_2)+b^2{\sin }^2({\theta }_2)}})-\mathrm{atan}(\frac{h}{2\sqrt{\frac{a^2+b^2}{2}}})$$

(9)

$${\theta }_3:\mathrm{atan}(\frac{h}{a+b})-\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_3)+b^2{\sin }^2({\theta }_3)}})=\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_3)+b^2{\sin }^2({\theta }_3)}})-\mathrm{atan}(\frac{h}{2\sqrt{ab}})$$

(10)

$${\theta }_4:\mathrm{atan}(\frac{h}{2\sqrt{ab}})-\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_4)+b^2{\sin }^2({\theta }_4)}})=\mathrm{atan}(\frac{h}{2\sqrt{a^2{\cos }^2({\theta }_4)+b^2{\sin }^2({\theta }_4)}})-\mathrm{atan}(\frac{h}{2b})$$

(11)

In general, the flexible terrain slope estimation contained two steps, which includes the calculation of ground extent and the choice of the footprint diameter. Based on these two steps, we estimated the terrain slope by the flexible method, validated the estimated accuracy, and analyzed the terrain slope errors. Our detailed scheme was shown in Figure 3.

3.1.1. Calculation of Ground Extent

Ground extent (h) represents the ground vertical height within the footprint, which can be calculated from waveform analysis [29,50,51,52]. The effect of the terrain slope on the waveform extent for the nadir-pointing laser can be expressed in Figure 4a–c. It is obvious that the ground return was stretched by introducing the terrain slope within the footprint. Therefore, the ground extent can be extracted based on the decomposed ground waveform metrics from GLAS in this study [51,52].

First, the original GLAS waveform was smoothed by using a moving mean filter to eliminate the noise signal. Then Gaussian decomposition was used to decompose the de-noised waveform into a series of Gaussian peaks [53,54]. The last decomposed waveform was considered as the ground waveform [55,56]. After obtaining the ground waveform, we extracted the ground signal beginning location and end location above 4.5 times the noise standard deviation to obtain the ground waveform extent (Figure 4d). Lastly, the ground extent was calculated from the ground waveform extent by eliminating the broadening effect of the emitted laser pulse, which is shown in Equation (12).

$$h=GroundWaveformExtent-c\cdot FWHM$$

(12)

where c is the light velocity and FWHM is the full width half maximum of the emitted pulse width.

3.1.2. Choice of Footprint Diameter

The chosen footprint diameter determines the estimation accuracy of the terrain slope. From Figure 4e, the effective footprint diameter (d) is related to the intersection angle between the terrain aspect and footprint orientation.

Five estimation methods (Equations (2)–(6)) using a certain value to represent the footprint diameter do not consider the coupling effect of the footprint shape, orientation, and the terrain aspect. In contrast, the flexible method chooses the most appropriate footprint diameter on the terrain slope estimation.

Hence, we obtained the intersection angle as the difference of the terrain aspect from ASTER GDEM2 and the footprint orientation from GLAS data. Afterward, four threshold values were calculated by substituting the footprint major axis, minor axis, and ground extent into Equations (8)–(11). The optimal footprint diameter was determined by comparing the intersection angle with four threshold angles. Lastly, the terrain slope was calculated from the ground extent divided by the chosen footprint diameter.

3.2. Validation and Analysis of the Terrain Slope Derived from GLAS

We explored the terrain slope estimation approach by GLAS data across different sites. The terrain slopes retrieved from GLAS were validated by the ATM-observed slopes. The estimation accuracy was evaluated based on the error bias, standard deviation, the root-mean-square error (RMSE), and the coefficient of determination (R²). More error analyses were done by visualizing and fitting the distribution of the slope errors. In addition, the results of the flexible method were compared with those estimated by five other slope estimation methods with the same datasets. To further evaluate the performance of the flexible method, we also compared the estimation accuracy of the flexible method and five direct methods in different terrains and various eccentricity footprints.

4. Results and Discussion

4.1. Validation Against Airborne LiDAR Slope

We compared the GLAS-derived terrain slope with the ATM-measured slope for all sites by the flexible method and five direct methods. The results were summarized in Figure 5 and

Table 4. The statistical errors of all experimental data by five methods and the flexible method.

Method	Number of Estimates	Estimated Bias (°)	Estimated Standard Deviation (°)	RMSE (°)	R2
Method 1	858	−0.556	4.142	4.147	0.760
Method 2	858	0.473	4.165	4.170	0.757
Method 3	858	−0.101	4.023	4.027	0.770
Method 4	858	−0.098	4.026	4.031	0.770
Method 5	858	−0.114	4.020	4.025	0.771
Flexible Method	858	−0.064	3.592	3.596	0.829

. Overall, the terrain slopes derived from the six methods all had moderately linear correlations with the ATM-observed terrain slopes. However, compared to the five other methods, the flexible method had a lower standard deviation (3.592°), lower RMSE (3.596°), and higher R² (0.829) in the terrain slope estimation. This indicated that the estimation of the terrain slope was improved by considering the footprint shape, orientation, and terrain aspect.

For method 1, the regression result was slightly on the left side of the 1:1 line. The estimated bias error (−0.556) also showed that choosing the major axis as the footprint diameter led to the underestimation of the terrain slope. In contrast, the regression result and bias error of method 2 represented that choosing the minor axis led to the overestimation of the terrain slope. The regression lines of methods 3, 4, 5, and the flexible method almost coincided with the 1:1 line. We also found that methods 3, 4, and 5 had slightly lower RMSE and higher R² when compared with methods 1 and 2. This indicated that choosing the value between the major and minor axes achieved slightly higher accuracy. At the same time, they corrected the underestimation of method 1 and the overestimation of method 2 to some degree.

Upon further analysis of the terrain errors of these methods, we visualized the estimated errors in the histograms and fitted them by the normal distribution curves. Relevant results were shown in Figure 6 and Figure 7. Although the errors of all methods followed the normal distribution, the means and variances of the distribution curves tell us more. The mean of error by method 1 is less than zero, by method 2 is more than zero, and, by other methods, are very close to zero. The underestimation of method 1 and the overestimation of method 2 were more clearly shown in Figure 7. However, although methods 3, 4, and 5 reduced the underestimation and overestimation of method 1 and 2, the variances hardly decreased. The flexible method not only corrected the under-or-over estimation but also made the errors more concentrated and closer to zero. These results verified the conclusions in Figure 5 and Table 4 from the perspective of the error distribution.

Although the flexible method had much higher estimation accuracy in terms of R², the RMSE of the terrain slope only improved by ~0.5° when compared with the five other methods (Table 4). This may be due to a mass of points in the low-relief region. For these points, whichever one chosen as the footprint diameter made little difference in the absolute errors of the terrain slope estimates. Thus, the estimation accuracy of the terrain slope was slightly better when using the flexible method.

4.2. Slope Estimation of Different Terrain Relief

Figure 8 showed the absolute mean difference between the ATM-observed and the GLAS-predicted terrain slopes with the gradient of the terrain slope, which also supported the conclusion that the flexible method was superior to five other methods in estimating the within-footprint terrain slope. Among the five direct methods, method 2 usually has the largest estimation error, which is followed by method 1. The estimation errors of methods 3, 4, and 5 are close with all of them being less than methods 1 and 2. This coincided with the theoretical error analyzed by Reference [32]. In addition, we also found that, as the terrain slope observed from ATM increased, the absolute mean error was reduced more by the flexible method. For instance, when the terrain slope is less than 5°, the curves of six methods in Figure 8 were almost identical. With the increment of the terrain slope, the difference between the curve of the flexible method and other curves increased. Therefore, we made quantitative statistics of 218 GLAS footprints in the high-relief region (slope observed by ATM > 5°) and 640 footprints in the low-relief region (slope observed by ATM ≤ 5°). The results of slope estimation in the high-relief and low-relief regions were summarized in

Table 5. The statistical errors of experimental data in the high-relief (slope > 5°) and low-relief (slope ≤ 5°) terrain by five methods and the flexible method.

Method	Number of Estimates (Slope > 5°)	Estimated Bias (°)	Estimated Standard Deviation (°)	RMSE (°)	R2
Method 1	218	−3.117	6.825	6.856	0.664
Method 2	218	−0.704	6.874	6.906	0.662
Method 3	218	−1.999	6.490	6.521	0.680
Method 4	218	−1.915	6.498	6.529	0.680
Method 5	218	−2.081	6.483	6.515	0.681
Flexible Method	218	−1.834	5.153	5.180	0.757
Method	Number of Estimates (Slope ≤ 5°)	Estimated Bias (°)	Estimated Standard Deviation (°)	RMSE (°)	R2
Method 1	640	0.316	1.017	1.025	0.595
Method 2	640	0.840	1.018	1.026	0.591
Method 3	640	0.545	1.017	1.024	0.594
Method 4	640	0.561	1.017	1.024	0.594
Method 5	640	0.530	1.017	1.024	0.594
Flexible Method	640	0.458	0.904	0.936	0.636

.

The results in the high-relief region (n = 218, R² = 0.757, RMSE = 5.180°) and in the low-relief region (n = 640, R² = 0.636, RMSE = 0.936°) yielded the similar conclusions as the results of all data. In addition, the flexible method significantly reduced the RMSE of the terrain slope estimation by ~1.4° when compared with methods 3, 4, and 5 and by ~1.8° when compared with methods 1 and 2 in the high-relief region. The detailed results in the high-relief region were clearly shown in Figure 9. Conversely, in the low-relief region, the RMSE of the flexible method was only reduced by ~0.1°. In theory, the coupling effect of the footprint orientation and terrain aspect will lead to the greater absolute error when the terrain slope is greater [32]. The flexible method significantly reduced this absolute error in the high-relief region both theoretically and practically. Therefore, we would use the flexible method to estimate the within-footprint terrain slope in the high-relief region.

4.3. Slope Estimation of Different Footprint Eccentricity

Figure 10 showed the absolute mean error of GLAS-derived terrain slopes with the gradient of footprint eccentricity. Overall, the terrain slope errors by all methods increased with the ellipse eccentricity. This coincided with the theoretical error analyses of footprint eccentricity [32]. Beyond that, as the footprint eccentricity increased, the flexible method decreased more error of the terrain slope. Hence, we did the statistics of 313 GLAS footprints whose ellipse eccentricity was greater than 0.6 and 545 GLAS footprints with eccentricity less than 0.6. The results of the slope estimations with different footprint eccentricity were summarized in

Table 6. The statistical errors of experimental data in the case of different eccentricity by five methods and the flexible method.

Method	Number of Estimates (e > 0.6)	Estimates Bias (°)	Estimates Standard Deviation (°)	RMSE (°)	R2
Method 1	313	−1.012	4.488	4.502	0.748
Method 2	313	0.653	4.517	4.532	0.745
Method 3	313	−0.298	4.193	4.208	0.764
Method 4	313	−0.295	4.196	4.211	0.763
Method 5	313	−0.363	4.191	4.205	0.764
Flexible Method	313	−0.230	3.409	3.421	0.838
Method	Number of Estimates (e ≤ 0.6)	Estimates Bias (°)	Estimates Standard Deviation (°)	RMSE (°)	R2
Method 1	545	−0.373	3.927	3.963	0.773
Method 2	545	0.289	3.913	3.921	0.777
Method 3	545	−0.066	3.931	3.938	0.775
Method 4	545	−0.053	3.919	3.927	0.776
Method 5	545	−0.078	3.920	3.929	0.776
Flexible Method	545	−0.041	3.589	3.596	0.813

.

From Table 6, the flexible method performed better for datasets with high-eccentricity footprints (n = 313, R² = 0.838, RMSE = 3.421°) compared with the footprints with low eccentricity (n = 545, R² = 0.813, RMSE = 3.596°). For the high-eccentricity footprints, the flexible method significantly reduced the RMSE of the terrain slope estimation by ~0.8° when compared with methods 3, 4, and 5 and by ~1.1° compared with methods 1 and 2. The RMSE of low-eccentricity footprints was only reduced by ~0.4° when compared with other methods. Figure 11 showed the detailed validation results for the dataset whose eccentricity was greater than 0.6. This verified the theoretical error analysis in Reference [32]. When the footprint eccentricity is small, there is little difference among several adopted footprint diameters. The terrain slope error is theoretically small. However, in cases of the greater eccentricity, the coupling effect of the footprint shape, orientation, and terrain aspect must be fully considered on the terrain slope estimation.

4.4. Slope Estimation of Different Footprint Size

Figure 12 showed the absolute mean error of GLAS-derived terrain slopes as a function of the footprint diameter gradients. The footprint diameters of the dataset used in this study were approximately 48 m to 63 m. In theory, the absolute terrain slope error decreases with the footprint size and the flexible method performs better for data with a small footprint size. However, for all methods, the error indicators theoretically tend to be stable when the footprint size gradually increases, especially when the major axis exceeds 50 m [32]. In practice, the footprint sizes of all data in this study were large (almost over 50 m). The difference of the absolute mean errors between the flexible method and five direct methods was almost the same even though the footprint size varied.

4.5. Limitations of the Flexible Method

All the above experimental results demonstrated that the flexible terrain slope estimation approach could be beneficial for different terrains and various footprint shapes in practice. However, there are still other error sources which this approach does not take into accounts.

First, the terrain condition containing the slope and roughness was not completely inversed. This study made an important assumption that the terrain within the footprint was only a simple slope without any roughness. However, the terrain roughness also has a broadening effect on the ground waveform [57,58]. This might induce uncertainty when estimating the terrain slope. Especially under the coupling effect of the uneven spatial distribution of the emitted laser energy and the terrain condition, the situation becomes more complicated and more uncertain. For instance, if the elevation of the footprint edge is high and the footprint center is low, the signal intercepted by the high-elevation region may be mistaken as the noise due to weak intensity. This will cause an underestimation of the terrain slope within the footprint.

A second source is the surface cover condition. In this study, we chose four sites located in Antarctica and Greenland because these areas were almost on a bare ground. The transmitted waveforms were usually expressed as a single peak. However, when there is some cover (e.g., vegetation and building) on the surface, it will vastly increase the uncertainty of extracting the ground waveform and estimating the terrain slope within the footprint. For example, the ground waveform is usually mixed with the vegetation waveform in a slope terrain covered by vegetation [35,52]. There is almost no algorithm to perfectly separate the ground waveform from the mixed waveform. In addition, the existence of surface objects leads to the reduction of the ground signal especially in the region with a dense forest. Therefore, the feasibility and accuracy of the flexible terrain estimation method should be further discussed in these regions covered by some objects.

In addition, the flexible slope estimation method was proposed based on the assumption of the nadir viewing and it relied on a rough GDEM to provide the prior terrain aspect information. In contrast, five other methods estimated the within-footprint terrain slope without any auxiliary data. Among these methods, methods 3, 4, and 5 partly corrected the underestimation of method 1 and overestimation of method 2. Therefore, in the absence of a prior terrain aspect, these three methods could be selected to estimate the terrain slope even though they hardly improved the estimation accuracy. Additionally, the flexible method should be used in the following cases: (1) high-relief region and (2) great footprint eccentricity through a series of the error analysis.

5. Conclusions

We applied and validated a slope estimation approach for large-footprint full-waveform LiDAR data to generate the accurate terrain slopes within the footprint. Building on a previous study [32] about related theoretical analyses, this experiment was conducted at several study sites with various terrain conditions, footprint sizes, and eccentricity. The results demonstrated that the flexible terrain slope approach was beneficial for improving the accuracy of GLAS slope estimation in practice. This was in agreement with the theoretical analysis results of Reference [32]. Additionally, the flexible method also improved the slope estimation results with respect to the terrain slopes derived by five other methods, which did not consider the footprint orientation and slope aspect (R² increased from ~0.757 to 0.829 and RMSE decreased from ~4.170° to 3.596°, which is supported by 858 data). In five direct estimation methods known as methods 3, 4, and 5, choosing values between the major and minor axes to estimate the terrain slope were slightly better than methods 1 and 2, which chose the major or minor axes. However, these three methods just reduced the underestimation of method 1 and overestimation of method 2 to some degree, but they did not correct the slope error caused by a footprint orientation and terrain aspect. In contrast, the flexible method took account of these factors and improved the terrain slope estimation accuracy.

Based on these comparisons, we quantified the estimation accuracy of the flexible method on GLAS terrain slope retrieval in different terrains and various footprint shapes. We obtained a more improved RMSE in the high-relief region by using the flexible method (n = 218, RMSE = 5.180° for the flexible method, reduced by ~1.8°). The flexible method also significantly improved the terrain slope estimation accuracy and reduced the RMSE for the experimental data whose footprint eccentricity was more than 0.6 (n = 313; RMSE = 3.421° for the flexible method, reduced by ~1.1°). In contrast, the slight enhancement by a flexible method was taken place in a low-relief region (n = 640, RMSE = 0.936° for the flexible method, only reduced by ~0.1°) and low-eccentricity footprints (n = 545, RMSE = 3.596° for the flexible method, only reduced by ~0.4°).

To summarize, this study practiced a flexible terrain slope estimation method by using GLAS waveform data. This method performed better when compared with those methods without considering the elliptical footprint. In particular, it corrected more error caused by the terrain aspect, footprint shape, and orientation in LiDAR footprints with high relief or great eccentricity. Future space-borne LiDAR satellites such as NASA’s Global Ecosystems Dynamics Investigation (GEDI) and GF-7 satellite in China will record the return signals in the form of the full-waveform data within the large footprint. The validated method in this study will be applied for the global terrain monitoring by future LiDAR satellites. Additionally, caution should be taken in regions with complicated terrain conditions or surface covers because these related errors increase uncertainty of the terrain slope estimation. Our next step is to discuss and analyze the feasibility of estimating the terrain slope using the space-borne waveform data over forest areas with more complex terrain. Then, this flexible method could be generalized to the global terrain slope estimation within the footprints. In addition, this paper has proven the high estimation accuracy of the terrain slope in the south and north poles. The subsequent study will pursue the change of the terrain slope and analyze the relationships between the global climate change and the terrain change of two poles.