Estimation of Forest Canopy Height from Spaceborne Full-Waveform LiDAR Data Using a Bisection Approximation Decomposition Method

Abstract

Forest canopy height (FCH) is a vital indicator for assessing forest health and ecosystem service capacity. Over the past two decades, full-waveform (FW) LiDAR has been widely employed for estimating forest biophysical variables due to its high precision in measuring vertical forest structures. However, the impact of terrain undulations on forest parameter estimation remains challenging. To address this issue, this study proposes a bisection approximation decomposition (BAD) method for processing GEDI L1B data and FCH estimation. The BAD method analyzes the energy composition of simplified echo signals and determines the fitting parameters by integrating overall signal energy, the differences in unresolved signals, and the similarity of inter-forest signal characteristics. FCH is subsequently estimated based on waveform peak positions. By dynamically adjusting segmentation points and Gaussian fitting parameters, the BAD method achieved precise separation of mixed canopy and ground signals, substantially enhancing the physical realism and applicability of decomposition results. The effectiveness and robustness of the BAD method for FCH estimation were evaluated using 2049 footprints across varying slope conditions in the Harvard Forest region of Petersham, Massachusetts. The results demonstrated that digital terrain models (DTMs) extracted using the GEDI data and the BAD method exhibited high consistency with the DTMs derived using airborne laser scanning (ALS) data (coefficient of determination R² > 0.99). Compared with traditional Gaussian decomposition (GD), wavelet decomposition (WD), and deconvolution decomposition (DD) methods, the BAD method showed significant advantages in FCH estimation, achieved the smallest relative root mean square error (rRMSE) of 17.19% and greatest mean estimation accuracy of 84.57%, and reduced the rRMSE by 10.74%, 21.49%, and 28.93% compared to GD, WD, and DD methods, respectively. Moreover, the BAD method exhibited a significantly stronger correlation with ALS-derived canopy height mode data than the relative height metrics from GEDI L2A products (r = 0.84, p < 0.01). The robustness and adaptability of the BAD method to complex terrain conditions provide great potential for forest parameters using GEDI data.

1. Introduction

Forests play a crucial role in terrestrial ecosystems as the largest carbon sink on land. Forest canopy height (FCH), a key indicator of forest growth, is closely linked to forest health and directly affects ecosystem service capacities [1,2]. Rapid and accurate acquisition of FCH is essential for assessing forest carbon sequestration and carbon emission, providing a solid scientific basis for formulating effective climate change mitigation and adaptation measures [3,4].

Compared with traditional ground-based measurements, remote sensing technologies greatly improve the efficiency and spatial details of information in forest monitoring [5]. Based on different data sources, remote sensing methods are categorized into optical, microwave, and Light Detection and Ranging (LiDAR) remote sensing. Optical remote sensing, however, is limited by spectral saturation, making it challenging to accurately capture the three-dimensional (3D) structures of forests [6,7]. Microwave remote sensing offers penetration capability, but its accuracy in estimating forest parameters remains limited [8,9]. LiDAR, on the other hand, emits laser pulses and measures the time difference of returned signals, enabling high-precision acquisition of 3D structural information of terrestrial features. Compared with optical and microwave remote sensing, LiDAR can penetrate canopy gaps and is unaffected by weather or cloud cover, providing detailed 3D structural information about forests [10,11]. Among LiDAR platforms, spaceborne LiDAR, with its high orbital altitude and extensive observation coverage, surpasses ground-based and airborne systems in monitoring forest parameters on a regional scale [12,13].

The Global Ecosystem Dynamics Investigation (GEDI) mission, launched in April 2019, is mounted on the International Space Station (ISS), covering latitudes from 51.6° N to 51.6° S. GEDI’s primary goal is to collect global 3D forest ecosystem information for studying forest growth, species diversity, tree structure, and vegetation dynamics [14,15]. GEDI data feature high resolution and dense sampling, capturing continuous vertical forest structures, making them uniquely advantageous for forest parameter estimation [16,17]. As the latest generation of full-waveform (FW) LiDAR platforms, GEDI provides valuable resources for studying vertical structural parameters such as FCH.

In FW LiDAR data, the first peak often represents the canopy top, while the last peak corresponds to the ground. The relative distance between these two peaks is used to estimate FCH [18,19]. The waveform decomposition method enables more accurate identification of peak positions, facilitating the precise extraction of forest canopy height [19,20]. Unlike optical and SAR data sources, which often rely on regression models and ground-truth data [12], waveform decomposition directly estimates canopy height using only LiDAR data, offering significant advantages for this application. Existing waveform decomposition methods primarily include Gaussian decomposition (GD), wavelet decomposition (WD), and deconvolution decomposition (DD). GD fits the waveform with multiple Gaussian components, offering simplicity and computational efficiency [20]. WD leverages wavelet functions for local feature extraction and multi-scale analysis, suitable for complex signal decomposition [21]. DD removes system noise, yielding waveform features closer to the actual vertical distribution of forest structures [22].

Although these methods perform well in flat terrains, challenges arise in non-flat areas where complex forest structures and topographic variations broaden and mix waveform signals. In rugged terrains, echoes from the ground, shrubs, and trees of the lower layer are mixed, increasing uncertainty in FCH estimation [23,24]. Notably, slopes significantly affect FCH estimation accuracy, with root mean square error (RMSE) increasing by 14% when the slope exceeds 20° compared to slopes between 10° and 20° [25].

Slope data can be used to correct FCH estimation. Lefsky et al. [24] corrected waveform decomposition results using slope data from a digital elevation model (DEM), achieving explainable tree heights with an RMSE of 5 m. Similarly, Nie et al. [26] improved terrain correction methods with DEM data, enhancing FCH estimation accuracy (R² = 0.77).

Machine learning models can be used to estimate FCH by combining GEDI-extracted relative height metrics with field data. Zhu et al. [27] used GEDI L2A relative height metrics and Landsat imagery to model FCH in the Harvard Forest, achieving high accuracy with an RMSE of 3.11 m. Similarly, Chen et al. [28] extracted GEDI waveform and terrain features to establish regression models for estimating forest canopy height in the mountainous regions of northern Henan, obtaining satisfactory accuracy with an RMSE of 2.23 m. However, this approach relies on external data sources, increasing model complexity and uncertainty. Moreover, errors in intermediate FCH products are difficult to quantify during subsequent modeling processes, affecting overall accuracy.

To address these challenges, we developed a bisection approximation decomposition (BAD) method for processing spaceborne LiDAR full-waveform (FW) data to estimate FCH. This method analyzed the internal energy composition of waveform signals, determined fitting coefficients by integrating overall energy, differences in unresolved signals, and canopy signal characteristics, and ultimately estimated FCH. GEDI L1B waveform data were decomposed using this method, and the results were validated against airborne-laser-scanning-derived canopy height model (ALS-CHM) data from the National Ecological Observatory Network (NEON) Airborne Observation Platform (AOP) [29]. The proposed method was compared with traditional methods including GD, WD, and DD based on the FCH values derived from ALS-CHM. This study aimed to provide an accurate FCH estimation method using spaceborne FW LiDAR while mitigating the influence of terrain slopes on FCH estimates.

2. Materials and Methods

2.1. Study Area

The study area is Harvard Forest (Figure 1), located in Petersham, MA, USA. Harvard Forest has an area of approximately 1500 hectares and features a hilly terrain with elevations ranging from 90 to 410 m above sea level. The region is characterized by a temperate continental climate, with cold winters, warm and humid summers, and an average annual precipitation of approximately 1000 mm, evenly distributed throughout the year. The forest vegetation represents a typical New England mixed forest ecosystem. Dominant coniferous species include white pine (Pinus strobus) and eastern hemlock (Tsuga canadensis). Broadleaf species are represented by red maple (Acer rubrum), sugar maple (Acer saccharum), paper birch (Betula papyrifera), and American beech (Fagus grandifolia).

2.2. GEDI Data and Processing

The GEDI tracks of the study area are shown in Figure 1b. We obtained the GEDI L1B footprint data from the Land Processes Distributed Active Archive Center (LP DAAC, https://cmr.earthdata.nasa.gov/search (accessed on 7 November 2024)). The GEDI data products are categorized into four levels: Level 1 provides geographic location waveform data, Level 2 offers footprint-level canopy height and profile metrics data, Level 3 includes gridded canopy height metrics and variability data, and Level 4 comprises footprint and gridded aboveground carbon estimation data [30]. Our dataset consists of version 2.0 GEDI L1B footprint data collected between February and December 2022, covering the Harvard Forest area. These data record geolocated ground return waveforms, essential for studying forest canopy structures and dynamics.

The GEDI L1B footprint data were filtered based on the geographical boundaries of the study area to exclude footprints falling outside the target region. To ensure data quality, footprints were selected based on the criteria of “degrade_flag = 0” and “solar_elevation < 0”. The “degrade_flag” set to 0 indicates that the data were not degraded, preserving their original quality and precision. Additionally, a “solar_elevation” value of less than 0 corresponds to nighttime acquisitions, which reduce interference between solar radiation and laser pulses [30,31].

Waveform signals during propagation generally approximate a Gaussian distribution. However, atmospheric factors inevitably introduce background noise, while interactions with foliage and atmospheric water vapor may result in sharp spikes in the waveform [32,33]. Background noise was estimated using Equation (1), where the constant v was set to 4 [32]. Figure 2a presents a comparison of the echo signals before and after denoising. The waveform amplitude reflects the intensity of the LiDAR signal, with the green curve representing the average background noise. The main oscillatory component (blue curve) is clearly emphasized. To mitigate sudden disturbances and high-frequency components, a Gaussian smoothing filter was applied (Figure 2b), which effectively extracted the underlying trends and fundamental features of the signals [34]. According to the Algorithm Theoretical Basis Document (ATBD), a smoothing width of 6.5 ns was chosen to ensure the integrity of waveform characteristics [35,36].

$$RXwave=RXwave-\left(nois{e}_{mean}+v\times nois{e}_{\theta }\right)$$

(1)

where $RXwave$ is the receiver echo waveform, $nois{e}_{mean}$ is the mean background noise, $nois{e}_{\theta }$ is the standard deviation of background noise, and $v$ is a constant.

Although the majority of the Harvard Forest area is covered by forest, certain GEDI footprints inevitably fell on non-forest areas. After denoising and smoothing the waveform signals, the footprints clearly represent the echo characteristics of ground objects at their respective locations. Non-forest footprints were excluded based on the following criteria: (a) chaotic signals lacking distinct main components or recognizable features, indicating interference at the footprint location and rendering the result unreliable [24,37]; (b) a single peak observed, whereas waveforms in forested areas typically exhibit multiple peaks corresponding to different reflective surfaces, such as the ground, vegetation, and canopy [24,37]; and (c) the heights of the first and last peaks being less than 4 m, as these peaks reflect canopy height, which is generally above this threshold in our study area.

Ultimately, 2049 high-quality GEDI footprints were retained (Figure 2c) for subsequent waveform decomposition and analysis.

2.3. ALS Data and Processing

The ALS data used in this study were acquired from the NEON field observatory (https://data.neonscience.org/data-products/ (accessed on 7 November 2024)) in August 2022. Funded by the U.S. National Science Foundation and operated by Battelle, NEON is a continental-scale observation facility designed to collect long-term, open access ecological data to better understand ecosystem changes across the United States. NEON employs a multi-scale sampling strategy that integrates a ground-based sensor network and high-resolution AOP to link terrestrial and spaceborne remote sensing scales [28].

The LiDAR instrument used in the NEON AOP is the Optech Gemini system, which has a laser footprint diameter of 0.25 m (0.8 m in wide beam divergence mode) and an elevation accuracy of <5–35 cm (1σ) [29,38]. For this study, we collected CHM, digital terrain model (DTM), and slope raster data provided by NEON AOP, all with a spatial resolution of 1.0 m. The ALS-CHM served as the reference dataset for evaluating the FCH extracted from waveform decomposition.

To account for the GEDI footprint diameter of 25 m, the NEON AOP raster data (ALS-CHM, ALS-derived DTM (ALS-DTM), and ALS-derived slope (ALS-slope)) were resampled to a resolution of 25 m using an averaging method. Forest canopy height (ALS-CHM) values were extracted for the 2049 retained footprints. The overall ALS-CHM distribution across all footprint locations ranged from 4.45 m to 31.15 m, with an average value of 19.50 m. The statistical distribution of heights within specific intervals is shown in

Table 1. Statistical distribution of forest canopy heights for footprint locations.

Height Range (m)	Number of Footprints	Minimum Height (m)	Maximum Height (m)	Average Height (m)	Coefficient of Variation (%)
<10	92	4.45	9.97	7.76	20.85
10–15	233	10.01	14.99	12.92	11.00
15–20	682	15.06	19.99	17.98	7.50
20–25	855	20.00	24.98	22.24	6.32
>25	187	25.02	31.15	26.50	4.61

.

2.4. Bisection Approximation Decomposition

During the transmission of forest echo pulses, the pulse signal initially interacts with the top of the forest canopy, followed by the ground, as illustrated in Figure 3a. The signal composition of the echo pulses during transmission is characterized by

$$Signa{l}_{total}=Signa{l}_{canopy}+Signa{l}_{ground}$$

(2)

where $Signa{l}_{total}$ is the total echo signal, $Signa{l}_{canpoy}$ is the canopy echo signal, and $Signa{l}_{ground}$ is the ground echo signal.

However, in areas with steep terrain slopes and complex forest stands, mixed waveforms can occur at the same horizontal elevation position. Moreover, the terminal peak in the waveform within the footprint represents the lowest point on the ground, leading to an overestimation of FCH estimates, as depicted in Figure 3b. In such cases, the signal composition of the echo pulses during transmission is influenced by

$$Signa{l}_{total}=Signa{l}_{canopy}+{\displaystyle \sum }Signa{l}_{c\&t}+{\displaystyle \sum }Signa{l}_{trunk}+{\displaystyle \sum }Signa{l}_{t\&g}+Signa{l}_{ground}$$

(3)

where $Signa{l}_{c\&t}$ is canopy–trunk mixed echoes and $Signa{l}_{t\&g}$ is trunk–ground mixed echoes.

However, when FCH is estimated, only the positions of the first and last wave peaks are required. Although these positions may shift due to slope effects and deviate from the actual locations, the broadening of the waveform ensures that the true first and last wave peaks are located within the range of the original peaks. To simplify the waveform decomposition process, we divided the signal into three components: ground signal, canopy signal, and intermediate mixed signal. We assume that the total energy of the entire echo signal is the sum of these three components, even though slope-induced effects may cause mixing between sub-canopy signals and the ground or canopy signals.

$$Signa{l}_{total}=Signa{l}_{canpy}+\sum Signa{l}_{midst}+Signa{l}_{ground}$$

(4)

where $Signa{l}_{midst}$ is intermediate mixed echoes.

To accurately determine the locations of these three signal components, a bisection search-based decomposition method was adopted. This approach iteratively approximates the positions of the ground, canopy, and mixed signals within the waveform. The specific steps of the proposed BAD method are described as follows:

(1) The number of peaks in the pre-processed echo signal is first determined, and the signal is divided into two segments based on the median peak count. This segmentation ensures that one segment predominantly contains canopy-mixed signals and the other contains ground-mixed signals, providing a foundation for subsequent decomposition.

For the canopy-mixed signal segment, the peak count (k) is determined:

• If k = 1, the entire segment is considered a canopy signal and is directly fitted using a Gaussian model.
• Otherwise, the peak position (M₁) nearest to the split point (O) is identified, and a tangent point (T) is defined at M₁. A tangent (L_t) is drawn perpendicular to the x-axis, centered at T. The width is set to L_t O, and the average height of the peaks is used as the amplitude to fit the sub-canopy signal (Figure 4a).

(2) For the remaining portion of the signal, half the horizontal distance between the tangent point (T) and the first peak (F) is used as the width, with the midpoint as the center and the average height of the peaks as the amplitude to fit the canopy signal (Figure 4b,c).

Throughout this process, key features are recorded, including the energy ratio of the canopy-mixed signals and the amplitude-to-width ratio of sub-canopy signals. If the distance (M₁F − TF) exceeds a predefined threshold (Threshold), the tangent point (T) is adjusted by moving it toward F by a distance proportional to k/(M₁F − TF). This adjustment is repeated until the distance falls below the threshold.

(3) The ground-mixed signal is processed using the same method as the canopy-mixed signal to extract features for both sub-canopy and ground signals.

(4) After completing all iterations, the decomposition results are evaluated. The optimal result is selected based on a combination of the overall energy difference between the decomposed and raw signals and the similarity of sub-canopy signal features. As shown in Equation (5), the evaluation $E_{score}$ is used to assess the performance of waveform decomposition. When $E_{score}$ approaches 1, it indicates accurate decomposition results.

$$E_{score}={\displaystyle \frac{\left|E_C-E_c\right|}{E_C}}+{\displaystyle \frac{\left|E_G-E_g\right|}{E_G}}+{\displaystyle \frac{d_2{a}_1{w}_1}{d_1{a}_2{w}_2}}$$

(5)

where $E_{score}$ is the decomposition score; $E_C$ and $E_c$ represent the total energy in the canopy region before and after decomposition, respectively; $E_G$ and $E_g$ represent the total energy in the ground region before and after decomposition, respectively; $d_1$ and $d_2$ denote the total distance in the canopy and ground regions, respectively; $a_1$ and $w_1$ represent the amplitude and width of the inter-layer signal in the canopy region, respectively; and $a_2$ and $w_2$ represent the amplitude and width of the inter-layer signal in the ground region, respectively.

2.5. Forest Canopy Height Estimation and Accuracy Evaluation

We utilized the proposed BAD method to analyze waveform signals obtained within forest footprints and estimated FCH based on the relative positions of decomposed sub-wave signals. Additionally, we compared this method with traditional methods including GD, WD, and DD. The FCH values derived from ALS were used as references for accuracy validation using metrics including coefficient of determination (R²), root mean square error (RMSE), mean absolute error (MAE), relative RMSE (rRMSE), and estimation accuracy (EA) [13,39].

$${\mathrm{R}}^2=1-{\displaystyle \frac{{{\displaystyle \sum }}_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{{{\displaystyle \sum }}_{i=1}^n{(y_i-\overline{y})}^2}}$$

(6)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{{\displaystyle \sum }}_{i=1}^n{({\widehat{y}}_i-y_i)}^2}{n}}}$$

(7)

$$\mathrm{MAE}={\displaystyle \frac{1}{n}}{{\displaystyle \sum }}_{i=1}^n\left|{\widehat{y}}_i-y_i\right|$$

(8)

$$\mathrm{rRMSE}={\displaystyle \frac{\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{({\widehat{y}}_i-y_i)}^2}{n}}}{\overline{y}}}\times 100\%$$

(9)

$$\mathrm{EA}=1-{\displaystyle \frac{1}{n}}{{\displaystyle \sum }}_{i=1}^n{\displaystyle \frac{\left|{\widehat{y}}_i-y_i\right|}{y_i}}\times 100\%$$

(10)

where $y_i$ is the ALS-derived FCH values, ${\widehat{y}}_i$ is the estimated FCH using GEDI L1B footprint data, and $n$ is the number of samples. In addition, we conducted the tests of significant differences of mean residuals from zero by comparing the FCH estimates with the airborne-LiDAR-derived FCH reference values and the tests of significant differences of mean absolute residuals from the paired methods based on Student-t distribution and the significance level of smaller than 0.05.

3. Results

3.1. Waveform Decomposition

Figure 5 illustrates the decomposition results of four waveform decomposition methods applied to forest footprint waveforms under varying slope conditions. The figure highlights notable differences in the ability of each method to handle slope-induced waveform distortions. The GD method maintains the original positions of waveform peaks and relies on Gaussian fitting based on the identified peak locations. While effective in scenarios with minimal slope influence, this method’s performance heavily depends on the accuracy of background noise estimation and signal smoothing. In highly variable signal conditions, GD often introduces spurious peaks. Additionally, the method does not account for slope-induced waveform broadening, making it particularly susceptible to significant errors in peak identification under steep slope conditions. The DD and WD methods exhibit robust smoothing capabilities, akin to secondary smoothing of the raw signal. This process effectively reduces spurious peaks and partially mitigates slope-induced waveform broadening. However, the secondary processing of the signal may lead to shifts in peak positions, particularly under steep slope conditions. The limited controllability of these shifts can compromise the stability and accuracy of the decomposition results. The BAD method redefines peak positions by decomposing the waveform signal into three components: ground, canopy, and branches (upper and lower). By integrating overall energy distribution, residual differences from undecomposed signals, and intra-forest signal similarities, the BAD method achieves a more precise and stable decomposition. Notably, the BAD method excels in addressing slope-induced waveform broadening while minimally altering the initial and final peak positions in flat terrain. This dual capability enhances the physical realism of the decomposition results and extends their applicability to diverse terrain conditions.

3.2. Ground Elevation Estimation with Waveform Decomposition

The estimation of ground elevation beneath each GEDI footprint is a critical component of waveform decomposition, as its accuracy directly reflects the effectiveness of the decomposition algorithm [26,40]. In this study, ground elevation data were extracted from 2049 pre-processed footprints within the study area using four different waveform decomposition methods. The extracted digital terrain models (DTMs) were then compared with the ALS-DTM to assess the accuracy of the decomposition methods. The results are shown in Figure 6.

DTMs obtained from all four decomposition methods exhibited high consistency with the ALS-DTM, as indicated by strong coefficients of determination (R² > 0.99). This consistency demonstrates the reliability of ground elevation results regardless of the chosen decomposition method, confirming the robustness of waveform decomposition for this application. Despite the overall high accuracy, differences in error levels were observed among the methods. Although the GD method shows the best performance in MAE (2.93 m) and EA (98.83%), it is evident from Figure 6b that the majority of data points fall below the Y=X line, indicating that the GD method tends to underestimate the DTM. This is primarily due to the waveform broadening effect, which causes the GD method to misclassify lower elevation points as ground level. In contrast, the DTM estimated using the BAD method is more evenly distributed around the Y=X line, with lower RMSE and rRMSE values of 4.33 m and 1.61%, respectively, demonstrating better overall consistency. These results are closely aligned with those of Nie et al. [26] and Wang et al. [41], but our method stands out by not requiring the introduction of additional data sources. These findings demonstrate the superiority of the BAD method in accurately estimating ground elevation.

3.3. Forest Canopy Height Estimation with Waveform Decomposition

In

Table 2. Accuracy comparison of FCH estimates among four decomposition methods.

Methods	R2	RMSE/m	MAE/m	rRMSE/%	EA/%
BAD	0.67	3.33	2.83	17.19	84.57
GD	0.37	5.66	4.20	29.21	75.82
WD	0.25	7.49	5.60	38.68	67.47
DD	0.23	8.94	7.42	46.12	57.19

, the results of FCH estimation using four distinct methods are compared. The proposed BAD demonstrated the highest estimation accuracy by achieving a greatest R² of 0.67 and EA value of 84.57%, with the smallest RMSE of 3.33 m, MAE of 2.83 m, and rRMSE of 17.19%. Based on error and accuracy metrics, the best performance was achieved by BAD, followed by GD and DD, and WD yielded the poorest performance. Compared to GD, DD, and WD, BAD reduced the rRMSE by 10.74%, 21.49%, and 28.93%, respectively. Moreover, the fitting plots illustrated that BAD closely matches the measured values, demonstrating superior estimation accuracy with a close alignment along the Y=X line (Figure 7). Although all the methods led to approximate normal distributions of the predicted values similar to those of the reference values, BAD resulted in estimates close to the line of Y=X, overall, while the three traditional methods led to the overestimation being obvious.

In Figure 8, the distributions of residual errors from four methods are compared. The residuals of the predicted values of FCH from BAD were predominantly within ±6 m and evenly spread around the zero line, indicating robust performance. Conversely, GD, DD, and WD exhibited residuals mainly falling in the range of 10 m to −10 m, suggesting less accurate and notable overestimations.

Additionally, as shown in

Table 3. Comparison of the four decomposition methods based on their correlation with ALS-CHM and the significance difference test of the absolute residual differences between methods (significance level ≤ 0.01).

Methods	Correlation Difference (R-Value) Between ALS-CHM and Estimated Results	Significant Difference (t-Value) Among Methods of Mean Absolute Residuals
	FCH	BAD	GD	WD	DD
BAD	0.82(0.00)	-	39.62	47.11	77.95
GD	0.61(0.00)	39.62	-	26.40	64.26
WD	0.50(0.00)	47.11	26.40	-	55.82
DD	0.48(0.00)	77.95	64.26	55.82	-

, a comparison of the FCH values estimated by the four methods with the ALS-CHM reference values revealed significant correlations (p < 0.01) for all methods. Among them, the BAD-FCH method exhibited the highest correlation (0.82) with ALS-CHM, while the correlations for GD-FCH, WD-FCH, and DD-FCH with ALS-CHM were 0.61, 0.50, and 0.48, respectively. Furthermore, a comparison of the absolute residuals of FCH extraction results among the four methods and the obtained t-values showed greater than the critical t-value at the significance level of less than 0.01, indicating that the mean absolute residuals of the methods were statistically significantly different from each other.

Compared to the study by Chen et al. [28], which utilized GEDI waveform feature data and machine learning models to estimate forest canopy height in the mountainous regions of northern Henan, our approach achieved comparable estimation accuracy (rRMSE of 17.19% and 17%, respectively). Similarly, when compared to the work by Zhu et al. [27], which estimated forest canopy height in the Harvard Forest, the accuracy of the BAD method also proves to be comparable, with RMSE values of 3.33 m and 3.11 m, respectively. However, the key difference is that we did not employ a regression model; instead, we directly analyzed the waveform data, which enhanced the interpretability, practicality, and physical realism of the results.

4. Discussion

4.1. The Impact of Terrain Slope

After the laser pulse from spaceborne LiDAR reaches the ground, it forms a large footprint, with diameters ranging from several tens of meters to hundreds of meters. Terrain slope variations significantly affect the shape of these footprints, thereby influencing the estimation of forest parameters [39]. Currently, there is no universally effective method to mitigate the impact of slope on full-FW LiDAR data. To further investigate the effect of slope on forest canopy height (FCH) extraction, we extracted ALS-slope data from NEON AOP and obtained the corresponding slope values for each footprint based on its location. We classified the slope values into six intervals: I [0–5°], II [5–10°], III [10–15°], IV [15–20°], V [20–25°], and VI [>25°], and analyzed the robustness of the four methods within each slope interval. The slope gradations were determined based on natural breaks observed in the topographic slope distribution of the study area, ensuring that the classification reflects the terrain’s inherent variability. This approach is consistent with the classification thresholds adopted by Fayad et al. [25]. Moreover, the use of 5° intervals provides a balanced resolution to analyze the effect of slope on FCH extraction results, offering clear insights into the trends within each range. Figure 9 shows the trends of RMSE and MAE for FCH estimation using the four decomposition methods across different slope intervals.

The results showed that for the GD, WD, and DD methods, RMSE gradually increased with the slope interval, a trend consistent with the findings of Lee et al. [42] and Liu et al. [29], indicating that the accuracy of FCH extraction decreases as the slope increases. In contrast, the RMSE for the BAD method remained relatively stable across the slope intervals, with no significant correlation to slope changes. In all slope intervals, the RMSE for the BAD method was consistently lower than that of the other three methods. Specifically, the order of RMSE from lowest to highest was BAD < GD < WD < DD. For each slope interval, the RMSE for the BAD method was always lower than the RMSE for any of the other three methods across all slope intervals. The trend for MAE (Figure 9e) followed a similar pattern.

The BAD method we proposed effectively redefines the sub-wave positions within the valid signal range by integrating the overall energy of the signal, the difference from the undecomposed signal, and the similarity of the in-between signals. This process adjusts the peak positions and reduces the broadening effect caused by the slope. Overall, the BAD method significantly outperforms the GD, WD, and DD methods in terms of FCH extraction accuracy across different slope intervals, showing stronger slope adaptability. In contrast, the GD, WD, and DD methods are more sensitive to slope changes. Compared to the studies of Lefsky et al., Liu et al., and others, the BAD method demonstrates stronger slope adaptability and achieves the smallest estimation error (rRMSE = 17.19%) without relying on additional data sources.

4.2. The Impact of ALS Canopy Height Resampling Method

This study evaluated the effects of different sampling methods on the accuracy of FCH estimation by comparing waveform-decomposition-estimated FCH with CHM data generated from airborne LiDAR. During the ALS-CHM generation process, point cloud data were gridded to match the GEDI footprint diameter (25 m), and window mean resampling was applied to obtain the average stand height within each footprint. However, the window mean resampling does not fully align waveform-decomposition-extracted canopy height with the average FCH across the stand.

To explore the suitability of various sampling methods, this study further resampled the ALS-CHM data at 1 m resolution using three distinct approaches (Figure 10):

(1)

Maximum Value Sampling (ALS-CHM(max)): extracts the maximum canopy height within the footprint, representing the vertical distance from the highest treetop to the lowest ground point. This metric corresponds to the distance between the start of the first echo and the peak of the last echo, providing an approximation of the total vertical range of the canopy.

(2)

95% Relative Height Sampling (ALS-CHM(95%)): extracts the 95% relative canopy height within the footprint. This approach accounts for potential deviations of the waveform’s first peak due to canopy structural complexity or environmental factors such as wind or water vapor, correcting for these influences.

(3)

Interval Average Sampling (ALS-CHM(sec)): calculates the average height of pixels within the 45%–95% relative height interval. This approach better reflects the upper canopy height while mitigating the effects of gaps and low ground-level values.

When using ALS-CHM(max) and ALS-CHM(95%) as reference variables, the DD method achieved the highest accuracy, with estimation accuracies of 81.33% (ALS-CHM(max)) and 80.75% (ALS-CHM(95%)). Both WD and GD methods yielded comparable results, with accuracies of 76.88% and 77.06% (ALS-CHM(max)) and 76.20% and 77.43% (ALS-CHM(95%)), respectively. However, both methods showed underestimation in some regions. The BAD method demonstrated the lowest accuracy, with extraction accuracies of 66.49% (ALS-CHM(max)) and 68.67% (ALS-CHM(95%)).

When using ALS-CHM(sec) as the reference variable, the BAD method outperformed other methods, achieving an accuracy of 84.31%. The GD method ranked second, with an accuracy of 76.04%. The WD and DD methods showed lower accuracies of 67.83% and 57.84%, respectively. Furthermore, the GD, WD, and DD methods generally overestimated canopy height compared to reference values.

Across different CHM sampling strategies, the WD and DD methods are more suited for extracting maximum canopy height but are prone to underestimating high values in ALS-CHM(max) and ALS-CHM(95%), while overestimating low values in ALS-CHM(95%) and ALS-CHM(sec). The BAD method consistently demonstrated superior performance in ALS-CHM(sec) and ALS-CHM, achieving significantly higher accuracy, particularly for estimating average canopy height.

The results highlighted the distinct advantage of the BAD method in estimating average FCH, which is a critical metric for assessing forest quality and estimating carbon sequestration potential. Average FCH serves as a key indicator in forestry applications, providing essential insights for forest ecosystem management and decision-making.

4.3. Comparison with GEDI L2A Products

The GEDI Level 2A (L2A) dataset, a Level 2 product released by NASA, includes height metric products that extend beyond the waveform data provided in GEDI Level 1B (L1B). GEDI L2A consists of 100 relative height (RH) metrics, which capture vertical height variations at footprint locations. These RH metrics have been widely used for estimating forest parameters, including FCH [43,44].

In this study, we obtained corresponding GEDI L2A data together with the downloaded GEDI L1B waveform data. By matching footprint identifiers, we extracted 11 default relative height metrics (RH90–RH100) for the same footprints. GEDI L2A (version 2.0) includes RH metrics extracted using six different algorithms, with the officially recommended results set as default [17,45]. To evaluate the effectiveness of the BAD method, we analyzed the correlation between FCH extracted from these 11 relative height metrics and the BAD-extracted canopy heights, using four sampling strategies for CHM data: ALS-CHM, ALS-CHM(sec), ALS-CHM(95%), and ALS-CHM(max). The results are presented in Figure 11.

The overall variation in RH90–RH100 metrics was small, with RH95, RH96, RH97, and RH98 showing the highest correlations, which is consistent with the findings of Zhu et al. [27] and Potapov et al. [46]. However, the FCH extracted from the BAD method demonstrated significantly stronger correlations with the four CHM sampling strategies compared to GEDI L2A metrics. The correlation coefficients between BAD-extracted FCH and ALS-CHM were 0.816 (ALS-CHM(mean)), 0.816 (ALS-CHM(sec)), 0.704 (ALS-CHM(95%)), and 0.692 (ALS-CHM(max)). In contrast, GEDI L2A RH90–RH100 metrics did not achieve comparable levels of correlation in the same comparisons.

The results highlighted the superior performance of the BAD method in FCH estimation. Its extracted canopy heights demonstrated significantly higher correlations with CHM-extracted canopy heights compared to the default RH90–RH100 metrics in GEDI L2A. This indicated that the BAD method offered greater adaptability and potential for accurately estimating FCH.

5. Conclusions

This study utilized the GEDI full-waveform spaceborne LiDAR data to estimate forest canopy height (FCH) in the Harvard Forest region of Petersham, Massachusetts, and validated the results using airborne-laser-scanning-derived canopy height model (ALS-CHM) data from the National Ecological Observatory Network. A novel bisection approximation decomposition (BAD) method was developed, which outperformed Gaussian decomposition (GD), wavelet decomposition (WD), and deconvolution decomposition (DD) approaches based on estimation accuracy (rRMSE: 17.19%, 29.21%, 38.68%, and 46.12%; EA: 84.57%, 75.82%, 67.47%, and 57.19%, respectively) and adaptability across varying slope conditions. At the same time, the digital terrain model (DTM) obtained using the BAD method showed good consistency with the airborne-laser-scanning-derived digital terrain model (ALS-DTM) data (R² > 0.99), further confirming the robustness of its waveform decomposition. The BAD method demonstrated a significantly stronger correlation with ALS-CHM data compared to the relative height metrics from GEDI L2A products (r = 0.84, p < 0.01). The BAD method’s ability to dynamically adjust segmentation points and Gaussian fitting parameters enables precise separation of mixed signals from the canopy and ground, enhancing the physical realism and applicability of decomposition results. These results highlighted the BAD method’s potential for fine-scale LiDAR data analysis and its value in providing robust technical and theoretical support for forest parameter estimation, particularly under complex terrain conditions.