Slope Estimation from ICESat/GLAS

Abstract

We present a novel technique to infer ground slope angle from waveform LiDAR, known as the independent slope method (ISM). The technique is applied to large footprint waveforms (∼60 m mean diameter) from the Ice, Cloud and Land Elevation Satellite (ICESat) Geoscience Laser Altimeter System (GLAS) to produce a slope dataset of near-global coverage at 0.5° × 0.5° resolution. ISM slope estimates are compared against high resolution airborne LiDAR slope measurements for nine sites across three continents. ISM slope estimates compare better with the aircraft data (R² = 0.87 and RMSE = 5.16°) than the Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) inferred slopes (R² = 0.71 and RMSE = 8.69°). ISM slope estimates are concurrent with GLAS waveforms and can be used to correct biophysical parameters, such as tree height and biomass. They can also be fused with other DEMs, such as SRTM, to improve slope estimates.

1. Introduction

Imaging the Earth’s surface using satellite Light Detection and Ranging (LiDAR) has shown great potential for mapping surface topography [1–3] and the distribution of vegetation parameters, e.g., forest height, from regional to global scales [4–10]. In particular, waveform LiDAR has proven its potential on a global scale via the Ice, Cloud and Land Elevation Satellite (ICESat) Geoscience Laser Altimeter System (GLAS). ICESat/GLAS was operational from 2003–2009 for about one month at a time, 1–3 times per year [11]. To date, this is the only mission to have recorded near global waveform LiDAR measurements of the Earth’s surface, providing data for studies of the oceans, cryosphere, atmosphere and biosphere.

A number of vegetation parameters, such as vegetation height, cover fraction, timber volume and biomass, can be estimated from GLAS waveforms [6,9,10,12–14]. However, accurate estimation of these parameters is compromised by the presence of sloped terrain. The LiDAR waveform reflected from a slope mixes with that of the overlying vegetation [1,15–17]. This dictates sloped terrain be a major factor affecting waveform signal-to-noise ratios [9] and, hence, induces difficulties in identifying ground returns. The failure to distinguish between signals from terrain and from vegetation can lead to waveform misinterpretations and inaccurate vegetation parameter retrieval [17]. Moreover, the waveform simulations of Rosette et al. [17] suggest that vegetation height retrievals for slopes ≥ 15° and high density canopy cover conditions (≥98%) are severely underestimated. In addition, footprint dimensionality can induce a directional bias and associated inconsistencies in waveform retrieved vegetation and terrain parameters; the dimensions of the elliptical GLAS footprints vary as a function of data acquisition period [18]. For instance, the first acquisition of Laser 1 is highly elliptical, and measurements are subject to an element of directionality, whereas the majority of Laser 3 acquisitions are more circular and will exhibit less directionality [18]. Accurate quantification of LiDAR inferred vegetation parameters is paramount in understanding the carbon cycle and forest ecosystem dynamics and in adequately assessing forest productivity and carbon sequestration rates [15,19]. Here, we develop a novel technique to retrieve slope estimates directly from LiDAR waveform returns, employing ICESat/GLAS waveforms for demonstration and validation purposes. The technique is expected to be applicable to all other waveform technologies, such as, for example, the Scanning LiDAR Imager of Canopies by Echo Recovery (SLICER) [20] and the Laser Vegetation Imaging Sensor (LVIS) [21].

In recent research, externally sourced slope information has been employed to filter suspected spurious waveform data from analysis [9,10,22]. For example, Simard et al. [9] and Los et al. [10] produced global forest height and global vegetation height distributions from GLAS data, employing the SRTM [23,24] data to derive close to spatially concurrent slope information to filter spurious data. However, as SRTM has a limited latitudinal extent between ±60°, reliable (direct inference) slope information is restricted between these boundaries. Extension beyond this latitudinal limit is possible by use of the Global 30 Arc-Second Elevation (GTOPO30) DEM product, which offers full global coverage between ±90°. However, the (global) vertical accuracy of this latter product (±18 m at 90% confidence) [25] is far less than that associated with SRTM equivalent measurements (±8m at 95% confidence) [23].

The uncertainties associated with the use of secondary DEM products depend on spatial resolution, location error and the amount of overlap. For example, the global SRTM data have a spatial resolution of 90 m (at the equator), and slope is estimated from two pixels on either side of a central pixel; this effectively reduces the resolution to 270 m, whereas the mean GLAS footprint size varies between 51 m and 102 m. This, combined with directional effects in the shuttle radar and location errors will introduce spatial discrepancies between GLAS footprints and corresponding SRTM DEM grid cells. By contrast, the method for retrieving ground slope in the present study utilises GLAS data, itself eradicating spatial inconsistencies, therefore minimising spatial uncertainties. Predicted within-footprint slopes should yield superior accuracy with respect to slope obtained from other sources (noted above), mainly due to identical spatial agreement between data.

2. Data

2.1. Waveform Data

In this study, we use the ICESat/GLAS land data (GLA14) product, release 33 [26,27]. Whilst in operation (2003–2009), GLAS emitted a Gaussian shaped laser pulse in the 532 and 1064 nm wavebands, approximately 0.75–0.90 m in width (corresponding to a 5–6-ns duration) at the full width half maximum (FWHM). The returned waveform energy profile (in volts, V) is utilised as a proxy to inform of intercepted surfaces, such as vegetation and/or the ground. The returned waveform was measured for a duration equivalent range (i.e., two-way travel) window of 81.6 m (544 ns) for earlier campaigns or 150 m (1000 ns) for later campaigns at 0.15-m (0.5 ns) intervals; the capture duration is expressed relative to a reference range time from the spacecraft [18]. GLAS was the first spaceborne waveform LiDAR instrument designed to offer almost continuous, near global coverage between ±86° latitude [12]. However, technical difficulties dictated that data were collected intermittently for its operational lifetime. This consisted of typically 2–3 periods of approximately one month-long data acquisitions per year [28,29], known as laser campaigns. GLAS footprints are variable size ellipses, with average dimensions (minor and major axes) of 52 m × 95 m for campaigns L1–L2 C and 47 m × 61 m for other campaigns, equivalent to an average circular diameter of approximately 64 m [30]. The returned waveforms contain various peaks according to intercepted surfaces and are fitted with up to six Gaussians (Figure 1) as described by Duong et al. [31], the sum of which define the “model alternate fit” return pulse. Note, not all waveforms exhibit such distinct ground returns, as in Figure 1 [16,17,32], which can infer difficulties in identifying the ground from waveform return profiles.

Table 1. Summary of the GLAS ancillary data parameters used in this study.

GLA14 Code	Description
i_rec_ndx	Unique waveform batch number
i_UTCTime	Time code in reference to 01/01/2000 00:00:00 UTC
i_lat	Latitude of waveform centroid
i_lon	Longitude of waveform centroid
i_SigBegOff	Signal begin range increment
i_ldRngOff	Land range offset
i_SigEndOff	Signal end range offset
i_gpCntRngOff	Centroid range increment for up to six peaks
i_FRir_qaFlag	Cloud detection flag
GLA05 Code
i_parmTr	Gaussian parameters of emitted pulse

provides GLA14 parameters used in this study in methodological processing (Section 3). Data without geographical coordinates or that exhibit cloud contamination (i_FRir_qaFlag ≠ 15) have been removed.

2.2. Validation Data

We tested ISM slope estimates on data from nine study sites where high resolution (≤1 mpixels) airborne LiDAR (ALS) (raster) data were available. These sites are located in Canada, Sweden, Germany, the Netherlands and Australia. Five sites are located within the northern boreal forests, three in temperate forests and one in a sub-alpine mountainous region.

The sites exhibit topographic variability, ranging from flat terrain to areas of ≥ 35° slope and cover major forest biomes. The ISM slope estimates were tested predominantly within forested regions (chosen to test the capability of retrieving slope estimates in the presence of vegetation). Test regions encompassed a diverse range of conditions, both in terms of terrain and vegetation cover, including conditions known to be problematic for the large footprint data of GLAS, e.g., sparse and dense canopy.

Three of the four Canadian sites are of the former southern BOREAS (Boreal Ecosystem–Atmosphere Study) study sites in Saskatchewan (CA); these exhibit a flat topography with relatively homogeneous (land-cover types) forested areas of Trembling aspen (SOA), black spruce (SOBS) and jack pine (SOJP) stands [33,34]. Airborne LiDAR data were collected and processed by the authors in August, 2008. Data at Mosquito Creek (MC), near Banff (CA), a mountainous region with a mean elevation of 2500 m a.s.l, were collected and processed by the authors in August 2011. Site topography is relatively complex with localised karst surface features, where forest coverage (20%) is generally restricted to valleys below 2100 m [35]. The European boreal site of Norunda (NR), near Uppsala, Sweden (S), exhibits a relatively homogeneous forest with flat topography and a high boulder presence, inducing greater terrain roughness [36,37]. Airborne LiDAR data were collected and processed by the authors in June 2011. The second Swedish boreal site near Edsbyn (ED) exhibits a relatively homogeneous forest with moderate topography and scattered wetlands throughout; here, airborne LiDAR data were processed and provided by Lantm¨ateriet (the Swedish land survey) as part of the Swedish national LiDAR dataset. At Loobos (LB), near Apeldoorn, the Netherlands (NL), the topography is generally flat with some local undulation and open areas throughout the forest [38]. At Tharandt (T), near Dresden, Germany (D), a relatively dense, multi-layered heterogeneous forest stand is situated on undulating terrain [39]. Data from these two sites were collected and processed by the authors in June 2010. The Australian (AUS) site is located near the Tumbarumba (TM) research station within the Bago State Forest, New South Wales, exhibiting a quite open canopy with relatively complex terrain [40]. Airborne LiDAR data were collected and processed by the authors during November 2009. Airborne LiDAR survey information is summarised in

Table 2. Airborne LiDAR meta data. Date is the acquisition date of the airborne LiDAR data, and P_D is the mean point density of the airborne LiDAR data per square meter (pm⁻²). See the text for the site code reference.

Site	Date	Instrument	PD
SOA	August 2008	Optech ALTM3100	4 pm−2
SOBS	August 2008	Optech ALTM 3100	4 pm−2
SOJP	August 2008	Optech ALTM 3100	4 pm−2
MC	August 2011	Optech ALTM 3100C	4 pm−2
EB	2009–2013	Multiple	0.5 − 1 pm−2
NR	June 2011	Leica ALS50-II	6 pm−2
LB	June 2010	Leica ALS50-II	5 pm−2
T	June 2010	Leica ALS50-II	5 pm−2
TM	November 2009	Riegl LMS-Q560	5 pm−2

and study site characteristics in

Table 3. Summary of site locations and terrain characteristics. ψ, λ are latitude and longitude, respectively. N_G is the number of usable (post filtering) GLAS footprints present at each site. E_Range is the study site elevation range (minimum to maximum). θ̄ is site mean slope, and H_m is the site maximum vegetation height. These measures are sourced from airborne LiDAR data. See the text for site code reference. MC, Mosquito Creek; T, Tharandt; LB, Loobos; NR, Norunda; TM, Tumbarumba.

Site	ψ, λ (°)	NG	ERange (m)	θ̄ (°)	Hm (m)	Terrain	References
SOA	53.63, −106.20	0	524–572	2	21	Very low relief	[41–43]
SOBS	53.99, −105.12	43	551–593	2	16	Very low relief	[41–43]
SOJP	53.92, −104.69	21	468–495	1	13	Very low relief	[41–43]
MC	51.66, −116.31	120	2059–2958	20	20	High relief, complex	[35]
EB	61.72, 15.16,	913	129–520	6	28	Moderate relief, wetlands	-
NR	60.50, 17.29	36	34–83	2	28	Very low relief, roughness	[36,37,44]
LB	52.10, 5.44	66	56–101	3	17	Low relief	[38]
T	50.57, 13.34	119	141–653	7	26	Moderate relief	[39]
TM	−35.39, 148.09	531	812–1724	15	40	High relief, complex	[40,45]

. Study site geographic locations and terrain information are illustrated in Figure 2.

The interpolated SRTM DEM version 4.1, distributed by the Consultative Group for International Agriculture Research-Consortium for Spatial Information (CGIAR-CSI) [23,24], is employed to infer slope estimates at seven of the nine airborne LiDAR sites for comparison purposes. SRTM data are not available for > 60°N latitude, and therefore, the Edsbyn and Norunda sites (Sweden) are excluded from this comparison. Data from these two sites are useful as a demonstration of the proof of concept.

3. Method

3.1. Independent Slope Model (ISM)

Below, we derive the independent slope model (referred to as ISM throughout) and test it on GLAS waveforms. We use the GLAS model alternate fit return waveform, derived from the GLA14 Gaussian parameters (i_gpCntRngOff; see Figure 1), to infer slope. Estimating slope requires the accurate identification of the waveform ground return. In the case of GLAS, this is defined as the least elevated peak (greatest range) in the waveform profile between the i_SigBegOff (signal start) and i_SigEndOff (signal end) GLA14 elevation parameters (defined as the waveform extent; Figure 3a). We only analyse waveforms whose ground return maximum amplitude exceeds an imposed threshold of ≥ 0.2V. This threshold provided the best distinction (empirically) between returns representing ground only and mixed signal returns from ground and vegetation. A single Gaussian (G_f; Figure 3a) is fitted to the isolated waveform ground return, ensuring a single Gaussian representation, implemented to provide consistency in ground return shape. Ground return shape can be skewed, particularly when it represents mixed signal returns (as described above) and/or consists of multiple Gaussian fits.

It is assumed that at an intensity threshold of 0.001 V, the width of G_f (Figure 3b) globally (for GLAS) represents the vertical displacement between the ground’s least and most elevated positions within the footprint. This intensity threshold was selected to exceed Gaussian distribution “tail” values, noted on average in the order of 10⁻⁹ here. The selected threshold is the maximum functional value before significant changes are noted in the recorded width of G_f. Lower thresholds yield no significant change in the width of G_f, but are more susceptible to being masked in Gaussian tail values, which can result in spurious measurements. The distance measure of G_f is used in combination with the mean footprint diameter to obtain slope (θ) by Equation (1).

$$\theta =\text{atan}\left(\frac{{\text{W}}_{{\text{G}}_{\text{f}}}}{\overline{\text{D}}}\right)$$

(1)

where W_Gf is the width of G_f at 0.001 V and D̄ is the mean footprint diameter for individual waveforms (available from NSIDC[18] on a unique laser campaign basis; values are in the order of 10’s m).

3.1.1. Minimum Measurable Slope

Returned waveforms will always exhibit a finite value of ▿_E for (distinguishable) ground returns, because of the duration of the emitted signal and associated atmospheric attenuation. The implication is that even for completely flat surfaces, the lowest and highest within-footprint ground elevations differ, which translates to a finite slope measurement. An approach developed by Los et al. [10] for correcting vegetation height estimates can be applied here for similar purpose (with modification). The original method related vegetation height to the 5% values of the area under Gaussian 1. This yields a linear model that provides minimum measurable vegetation height estimates as a function of the area under Gaussian 1.

Based on this method, similar refinements can be applied in this study, where 1% FWHM values of G_f (in nanoseconds, ns) are linearly related to the returned waveform maximum amplitude (A_max), illustrated in Figure 4. This allows W_Gf (and subsequently, θ estimates from Equation (1) to be adjusted according to:

$${\text{W}}_{\text{m}}=\text{a}+\text{bA}$$

(2)

where W_m is the minimum observable slope, A is the returned waveform maximum amplitude (V) and a and b are fitting coefficients. Using ∼ 51, 000 waveforms from all study sites in addition to four predominantly flat desert regions, the fitting coefficients can be estimated as a = 4.689 and b = 0.759; flat regions allow slope overestimation to be reviewed where little slope is expected. The value of W_m is subtracted from all ISM estimates of W_Gf (Equation (1)), such that:

$$\theta =\text{atan}\left(\frac{{\text{W}}_{{\text{G}}_{\text{f}}}-{\text{W}}_{\text{m}}}{\overline{\text{D}}}\right)$$

(3)

3.1.2. Gaussian Fit Filter

Quality control was performed on G_f, ensuring ground returns that are well represented by G_f are analysed for slope. To test this, the intensity profiles of the waveform ground return and G_f were isolated above an intensity threshold ≥ 0.001 V and compared. The implication being that these intensity profiles should be identical (i.e., R² = 1.00) for perfect overlays of G_f and the waveform ground return, whereas poor fits are expected to yield reduced values of R². Where resulting values of R² ≤ 0.90, these waveforms are removed from further analysis, as results will unlikely represent true measurements.

3.2. Testing ISM Slope with ALS and SRTM

Both ISM and SRTM slope estimates were compared with slopes from high resolution ALS data for the 9 study sites. The comparison with SRTM data serves as a reference against which the performance of ISM is tested. ALS data were matched to unique GLAS footprints by overlaying the GLAS footprint perimeter (respecting size, azimuth rotation and location) over the ALS data and extracting data cells that were > 50% contained within the perimeter; slopes were calculated from the difference between minimum and maximum within-footprint elevations, divided by the mean footprint diameter (Equation (1)). SRTM slopes were calculated as the maximum difference between the centre SRTM cell located within the GLAS footprint and its eight neighbour cells [10].

4. Results and Discussion

4.1. Validation against Airborne LiDAR Slope

The comparisons of ISM retrieved slopes for all nine sites with airborne LiDAR (ALS) slopes are shown in Figure 5 and in

Table 4. Model performance statistics for ALS-, ISM- and SRTM-derived slope information. N = sample size, R² = square of Pearson correlation coefficient, p_R² = fractional probability of no correlation (significance), D = Kolmogorov–Smirnov statistic [46], F₂ = fraction of predictions within a factor of two of observed values, F_B = fractional bias, and RMSE = root mean square error.

	N	R2	pR2	D	F2	FB	RMSE (°)
ALS-ISM (All sites)	1849	0.81	< 10−10	0.06	0.61	0.06	5.08
ALS-ISM (SRTM available sites)	900	0.87	< 10−10	0.06	0.74	−0.02	5.16
ALS-SRTM (SRTM available sites)	900	0.71	< 10−10	0.38	0.40	−0.66	8.69

. The equivalent comparisons for SRTM (and ISM) slopes with ALS for seven sites are shown in Figure 6 and Table 4. The results show a closer similarity between ISM slopes and ALS slopes (R² = 0.87 and RMSE = 5.16°) compared to the SRTM slopes (R² = 0.71 and RMSE = 8.69°). Additionally noted statistics in Table 4 are the: Kolmogorov-Smirnov statistic (D) [46], the fraction of (ISM/SRTM) predictions within a factor of two of the (ALS) observations (F₂) and fractional bias (F_B). Each corroborate that ISM slopes relate to ALS equivalents more closely than SRTM slopes.

Figure 5a illustrates a comparison of ALS and ISM slopes for raw and filtered (according to Section 3.1.2) data; Figure 5b illustrates filtered data only, indicating the density of points per coloured pixel. Results suggest that the highest density of points occur near the 1:1 line (marked in both figures), up to an approximate maximum slope of 40°. Here, data filtering removed approximately 15% of all waveforms, and it is possible that some spurious data still remain in the current selection. Adjusting the acceptance threshold beyond R² = 0.9 (see Section 3.1.2) will remove a larger proportion of spurious data, but this will be at the cost of removing too many reliable data.

Variability in results can arise from factors, such as ground roughness, terrain complexity and variations in the density of multi-layered canopies [1,16,47], all of which are unfavourable for the clear identification of waveform ground returns, particularly in large footprint systems.

4.2. Comparison with SRTM Slope

Figure 6a,b and Figure 6c,d illustrate the same as Figure 5a,b, but for ISM and SRTM and a corresponding subset of ISM and ALS data, respectively. It is highlighted that SRTM-derived slopes tend to underestimate ALS equivalents, an effect that becomes more prominent with increasing slope.

Figure 7a illustrates the absolute mean difference between observed (ALS) and predicted (ISM and SRTM) slope values as a function of 1° slope intervals, supporting that SRTM more severely underestimates slope than ISM equivalents with respect to ALS observations.

The statistical significance of the difference between ISM and SRTM slope estimates (at 1° intervals) was quantified with a non-parametric Mann–Whitney–Wilcoxon test [48]. Statistically significant differences are highlighted in Figure 7b as black points within the grey region, illustrated with respect to sample population. In particular, SRTM slopes are not statistically different from ALS equivalents up to 5° (Figure 7a), as well as for larger slopes (>28°); in the latter case, likely due to the small sample population (see Figure 7b). From 5°–28°, the majority of results exhibit statistically significant differences, as supported by Figure 7b. For slopes up to ∼28°, ISM slopes compare well with ALS equivalents; for slopes larger than 5°, they compare better than SRTM estimates. However, beyond 28°, further data would be required to yield a valid test agreement, as the sample is small and the population poorly represented.

The mean differences between ALS and ISM slopes exhibit less variability than the mean differences between ALS and SRTM slopes; this is an artefact of the cell size used to compare data, i.e., the number of samples in the ALS tiles is larger, and the variability as measured by the standard deviation is smaller. The large cell size of SRTM data smooths or flattens surface features found on the natural landscape, hence causing slope underestimates, whereas the smaller GLAS footprints/cells are more sensitive to surface features.

4.3. Global ISM Coverage

In accordance with GLAS coverage, a near-global ISM product was realised on a 0.5° × 0.5° spatial grid (Figure 8). Note that Greenland and Antarctica are not part of the GLA14 land surface altimetry data product [18].

ISM slope from unique footprints are binned between 0°–70° at 0.5° intervals within each 0.5° × 0.5° grid cell, a method adopted from Los et al. [10]; values > 70° are deemed unrealistic and excluded. The value of each grid cell is the mean of all slopes between 0°–70° weighted by the frequency distribution of individual slopes (from the created histogram).

Visual inspection of the ISM product indicates most severe slopes occur in known mountainous regions, such as: the North American Rocky Mountains, the South American Andes and the Central Asian Himalayas. However, some indication is given to suggest that areas of the tropics exhibit significant slope, particularly within the Amazon basin and central Africa. Some features in tropical regions are real; however, ISM appears to often overestimate slope in these regions. This is suspected to be a legacy of dense, multi-layered vegetation canopies (cf. Section 4.4). Known areas with little terrain relief appear to be well represented by ISM inferred slope, particularly desert regions (e.g., the Sahara) and prairies (e.g., central Canada).

The ISM product exhibits a mean and standard deviation of 4.07° and 5.15°, respectively. It is estimated that approximately 73% of the observed land surface exhibits a slope ≤ 5°, where slopes from 5° to 25° cover approximately 27% of the land surface. Beyond 25°, slopes are observed to a maximum of 54°; however, these slopes constitute < 1% of the observed land surface. Even greater slopes are expected with finer spatial resolutions over mountainous regions. The coarse resolution of the current product is subject to high within-cell variability, the effect of which is diminished in taking the weighted mean of all slopes within a given grid cell. Whilst this is true for all cells, this effect is further exacerbated by the complex terrain that co-relates to mountainous regions, as a result a greater standard deviation, is noted for cells in such regions.

4.4. Benefits and Limitations of ISM

The ISM method for slope retrieval can be applied to any waveform LiDAR data source (not only GLAS, as employed here). The principles of the ISM essentially allow the retrieval of within-footprint vertical displacement. Provided the footprint’s mean diameter is known, slope can be retrieved from any footprint (dependent on waveform interpretation) from any waveform LiDAR data source with some minor changes required to accommodate instrument specific technicalities (e.g., emitted pulse length). If footprint dimensions are unavailable, they can be estimated using laser optical and physical (orientation) parameters [49]. Additionally, when used in conjunction with global GLAS data, the ISM is capable of providing slope estimates beyond the limits of other data sources. For example SRTM, used in most cases, offers limited coverage of ±60° latitude. The use of ISM with GLAS data allows slope retrievals to be made up to ±86° (limited by GLAS’ orbit). In addition, for dense canopies, SRTM penetrates to an inconsistent depth [50]. Hence, it neither represents a digital elevation model nor a digital surface model and cannot be assumed to be accurate for derivations of slope information.

Based on ALS and ISM comparisons, the ISM product is capable of reducing slope uncertainties within GLAS footprints based on a simplistic, robust methodology. The results presented in this study are subject to large within-footprint terrain variability, which scales as a function of increasing footprint mean diameter. Hence, for waveform LiDAR systems yielding smaller footprints than those used in this study, it is not unreasonable to expect improved results with reduced uncertainty. This has been noted in SRTM data also [23,24]. ISM results are limited according to the quality of the returned waveform profile, which can be disturbed by high density and/or multi-layered vegetation canopies [22]. Such an effect has been noted to some extent in this study, particularly at the German site of Tharandt. Increased variability noted at this site (RMSE = 5.92°; Figure 9a) is a legacy of vegetation structure, not terrain complexity. This conclusion can be drawn as sites with low canopy density and high terrain complexity, for example Tumbarumba (Figure 9b), produced more coherent results (RMSE = 4.32°) within the same slope range (0° ≤ θ ≤ 25°). Here, in order to ensure that a like-for-like comparison was performed between sites without any sample bias, the same number of (filtered) data present at Tharandt were randomly sampled within the same ALS slope range from Tumbarumba.

5. Conclusions

This study has demonstrated the use of waveform LiDAR as a means of estimating ground slope with high accuracy when compared to high resolution airborne LiDAR equivalent results. Slope estimates were realised from Geoscience Laser Altimeter System (GLAS) waveforms via the implementation of the specifically developed Independent Slope Model (ISM), a novel technique for retrieving slope from waveform profiles. The current study has achieved the first quantitative retrieval of ground slope estimates from any waveform LiDAR instrument, offering improved accuracies on equivalent estimates from another spaceborne sensor and has potential to be applied to multiple instruments at a near global scale.

ISM slopes were compared with Shuttle Radar Topology Mission (SRTM)-derived slope information, as this has been employed in the published literature [9,10], such as to filter global ICESat/GLAS continuous waveform LiDAR data of spurious returns according to an imposed slope threshold. The current study has demonstrated that based on high resolution airborne LiDAR equivalent results, slope retrievals from SRTM data tend to underestimate surface slope across such footprints. The findings note that the ISM improves slope prediction results, with respect to equivalent SRTM results, when compared to airborne LiDAR-derived equivalents (R² increases from 0.71 to 0.87, and RMSE decreases from 8.69° to 5.16°). Moreover, with the use of ISM in conjunction with global GLAS data, slope information can be retrieved from LiDAR waveforms outside the SRTM coverage region (±60°), up to ±86°, the latitudinal extent of GLAS coverage. The implications of such observations are that: (1) spatially concurrent slope information can be derived directly from all LiDAR waveform sources without the need of calibration and/or supplementary data; (2) accurate slope information can be derived beyond 60°, an SRTM coverage limitation; and (3) slope filtering can be applied to data with improved accuracy, as slope predictions from ISM appear more representative of the natural landscape than SRTM equivalent results.