A Satellite Full-Waveform Laser Decomposition Method for Forested Areas Based on Hidden Peak Detection and Adaptive Genetic Optimization

Abstract

Laser waveform data that contain rich three-dimensional structural object information hold significant value in forest resource monitoring. However, traditional waveform decomposition algorithms are often constrained by complex waveform structures and depend on the initial parameter selections, which affect the accuracy and robustness of the results. To address the issues of the strong dependence on initial parameters, susceptibility to local optima, and difficulty in detecting hidden peaks during waveform overlap in the traditional satellite laser waveform decomposition algorithms, this study proposes a waveform decomposition method that combines hidden peak detection and an adaptive genetic algorithm (HAGA). This method uses hidden peak detection algorithms to improve the accurate extraction of the Gaussian components from the original waveform and provides the initial parameters. The high-precision extraction of waveform parameters is achieved through the adaptive genetic algorithm (AGA) combined with Levenberg–Marquardt (LM) optimization. In the experimental validation, the proposed method outperformed the traditional methods in both waveform decomposition fitting accuracy and tree height extraction. The average waveform decomposition accuracy $R_{mean}^2$ for more than 2000 laser spots reaches 0.955, whereas the $RMSE$ of the tree height extractions is better than 2 m, demonstrating strong robustness and applicability.

1. Introduction

Forests are an integral part of terrestrial ecosystems and play crucial roles in carbon storage and ecological functions. The detection of forest structures aids in understanding the global carbon cycle [1,2,3]. LiDAR (light detection and ranging) is a system capable of high-precision distance measurements, providing highly accurate and high-quality data [4,5]. Among its types, full-waveform LiDAR captures the complete information of a laser pulse from emission to reception [6], enabling users to extract richer target information from the scanning data [7,8,9,10,11], and accurately represents the vertical structures of the forests [12,13].

Full-waveform LiDAR data can be considered the superposition of echoes generated by multiple targets at different heights within a laser footprint [6,14]. By decomposing the full-waveform data, it is possible to obtain single pulses generated by the scattering of various targets. Analyzing these pulses provides information about the distance to the target and other characteristic features [15]. Therefore, decomposing full-waveform data is a critical step. Traditional waveform decomposition methods include deconvolution [16] and Gaussian decomposition [17]. The deconvolution method considers that the received echo is actually the convolution of the transmitted pulse and the ground response, and the ground surface response function data can be obtained by deconvolution; however, the method needs to be carried out on the basis of the known transmitted waveform and the received echo waveform. However, because the deconvolution method itself has a certain morbidity problem, it has strict requirements on the quality of the preceding noise reduction calculation. Moreover, the implementation of the deconvolution method is more complicated, and different practical problems often require the use of different deconvolution methods; thus, its generality is weak [16,18,19]. The Gaussian decomposition method regards the echo waveform data as functions of the reception time, and the distribution adheres to the Gaussian distribution. The processing needs to involve Gaussian decomposition of the original echo data and fitting the decomposed waveforms, but the background noise and the waveform aliasing caused by the superposition of individual waveforms affect the accuracy of the extraction of the waveform feature’s parameters [17,20,21]. Traditional waveform decomposition methods can extract waveform features to a certain extent, but their accuracy and adaptability are limited, particularly when dealing with complex waveforms. To address these issues, researchers in recent years have proposed various improved methods. These methods incorporate more advanced models and algorithms, aiming to overcome the shortcomings of the traditional approaches and provide more precise waveform decomposition capabilities. For example, the decomposition method for full-waveform LiDAR data, which is based on a variable-component mixture Gaussian model and reversible jump Markov chain Monte Carlo (RJMCMC), as well as an optimized decomposition method (ODM) for full-waveform LiDAR data [22,23], can extract more effective information from the original waveform. These methods achieve better fitting of the original waveform, and effectively decompose closely overlapping echo components. However, they cannot accurately fit asymmetric or trailing waveforms. Xie proposed a full-waveform data processing method based on the genetic algorithm (GA) to address the sensitivity to the initial parameters during waveform decomposition [24]. However, this method cannot accurately extract the number of waveform components in complex waveforms with hidden peaks. Fayad et al. used principal component analysis and random forest regression to process GLAS data and estimate canopy height in French Guiana [25]. Azadbakht proposed a sparse-constrained deconvolution regularization method, which demonstrated excellent performance in eliminating the influence of system waveforms [26]. A decomposition method based on long short-term memory (LSTM) and the sparrow search algorithm (SSA) directly determines the number of Gaussian components in an echo. However, a large amount of accurate data are required to train the model [27].

In summary, traditional decomposition methods rely on the accuracy of the initial parameters [28,29], and if the accuracy of the initial values is poor, they fall into the local optimum or even fit dispersion, whereas the GA can perform preliminary optimization of the initial parameters to obtain globally better values; however, the number of Gaussian components is identified incorrectly when aliasing of the waveform occurs, resulting in less-than-optimal waveform fitting accuracy. Therefore, to overcome these problems, this paper proposes a decomposition method based on hidden peak identification with an adaptive genetic algorithm (HAGA). One of these algorithms, the adaptive genetic algorithm (AGA), is an improved genetic algorithm in which higher adaptability and efficiency are achieved by adjusting the crossover probability and mutation probability in the genetic algorithm. While the traditional GA runs with these parameters fixed, the adaptive genetic algorithm adjusts the parameters according to the fitness of the individual to better balance the ability of the global and local searches and to perform the initial optimization of the initial parameters [30]. Waveform parameters serve as an indicator of whether the identified waveform components were correctly extracted. Since waveform components are echoes generated by the scattering of ground targets, they reflect the true characteristics of the targets. Therefore, accurate extraction of waveform parameters is essential to ensure that the echo components reliably represent the physical properties of the ground targets, improving measurement accuracy and providing a solid basis for subsequent applications and analysis. The HAGA decomposition method consists of the following steps: (1) Use a hidden peak recognition method to extract all the initial waveform parameters. (2) Apply the AGA to optimize the waveform parameters, which provides more accurate parameters for fitting optimization. The proposed method helps improve the accuracy of waveform fitting, thereby providing more precise forest-sensitive waveform parameters.

2. Materials and Methods

2.1. Study Area

To validate the effectiveness of the waveform processing method proposed in this paper, this research area included the Chaocha Forest Farm and Shangyanggeqi Forest Farm in Genhe city, Inner Mongolia, China (50°45′42″N–51°7′29″N, 121°20′39″E–121°41′12″E), as shown in Figure 1. This study area has a cold temperate continental monsoon climate and is located in a low mountainous region. The mountain range runs in a northeast-southwest direction, with steeper sunny slopes and gentler shady slopes. The average slope is 13°, and the forest cover is approximately 75%. The dominant forest type is Xing’an larch (Larix gmelinii), and it has the largest and most widespread distribution, while the other major vegetation types include white birch (Betula platyphylla) and camphor pine (Pinus sylvestris var. mongolica).

2.2. Experimental Satellite Information

This study used experimental data from the Land Ecosystem Carbon Monitoring Satellite (“Goumang” satellite) with multibeam LiDAR as the data source. The “Goumang” satellite primarily targets terrestrial forest ecosystems and is equipped with a variety of payloads, including multibeam LiDAR, aerosol detection LiDAR, multiangle multispectral cameras, hyperspectral imagers, and multiangle polarization imagers. These payloads are used for the vertical structure detection of forests, vegetation fluorescence retrieval, atmospheric parameter inversion, and other tasks. The satellite’s observation range covers latitudes from 80°S to 80°N, with its multibeam LiDAR system comprising five beams. According to the specifications, the “Goumang” satellite can provide both L1- and L2-level data. L1-level data comprise laser data derived from waveform preprocessing of the raw strip data and preliminary localization of the laser spots (including the main waveform, delayed main waveform, echo waveform, and coordinate data for each laser channel in each payload observation). L2-level data are obtained by applying atmospheric tidal corrections, precise attitude and orbit determination, and three-dimensional coordinate calculation of the footprint points to the L1-level data. The satellite parameters are shown in

Table 1. Main parameters of the “Goumang” satellite’s vegetation LiDAR.

Index	Design Value
Laser beams	5
Laser repetition rate	40 Hz
Laser operation wavelength	1064 nm
Laser waveform frequency	1.2 GHz
Number of bits of laser emission and echo quantization	12 bit
Camera resolution of optical axis monitoring	≤8 m
Number of bits quantized by camera in optical axis monitors	12 bit

. This study utilized L2-level data, with the selected L2-level data information from the “Goumang” satellite presented in

Table 2. Information on the “Goumang” satellite/L2 data used.

“Goumang” Satellite/L2	Data Acquisition Time
CM1_CASAL_A1_20220814_0000000155_L20000234185.h5	20220814
CM1_CASAL_A1_20220902_0000000443_L20000234183.h5	20220902
CM1_CASAL_A1_20220902_0000000443_L20000234184.h5	20220902
CM1_CASAL_A1_20220907_0000000520_L20000234190.h5	20220907
CM1_CASAL_A1_20221002_0000000906_L20000234188.h5	20221002

. As the corresponding airborne LiDAR data were acquired between 20 August and 13 September 2022, it was necessary to select data covering this study area during the period from August to October 2022. The distribution of the laser spots within this study area is illustrated in Figure 2.

2.3. Real Validation Data

To further verify the practicality of the proposed method in extracting waveform parameters, the average tree height of all the trees within the sample plots, obtained from field-collected airborne LiDAR data, was used for validation. The airborne LiDAR data were acquired using a CAF-LiCHy system, with a Riegl VQ580-∏ laser scanner (manufactured by RIEGL Laser Measurement Systems GmbH, located in Horn, Austria) mounted on a UAV. The relevant parameters are shown in

Table 3. Riegl VQ580-∏ scanner parameters.

Parameter	Value
Wavelength	1550 nm
Pulse Divergence Angle	0.25 mrad
Transmit Pulse Width	3 ns
Scanning Angle Range	±37.5°
Maximum Re-Frequency	2000 kHz
Minimum Re-Frequency	150 kHz
Maximum Scan Rate	300 lines/s
Point Cloud Density	1–3 pts/m2
Echo Count	Up to 15 per beam
Ranging Accuracy (m)	0.02 m

. The data had an average point density of 3 points/m², with a reported horizontal and vertical accuracy of 30 cm.

Figure 3 shows the distribution of the laser spots corresponding to the location of the “Goumang” satellite data within the study area. This distribution is based on the acquisition of airborne point cloud data. The location information and average tree heights within the sample plots are presented in

Table 4. Sample point locations and airborne LiDAR data used to extract the mean tree heights.

Index	Lat/(°)	Lon/(°)	Hmean/(m)	Index	Lat/(°)	Lon/(°)	Hmean/(m)
0	50.976	121.506	13.724	24	50.867	121.572	13.018
1	50.843	121.434	13.241	25	50.947	121.64	11.576
2	50.842	121.63	7.522	26	50.875	121.533	14.070
3	50.766	121.423	12.542	27	50.967	121.569	7.959
4	50.823	121.371	12.689	28	51.116	121.672	6.403
5	50.791	121.488	13.709	29	50.938	121.417	10.601
6	50.784	121.52	3.854	30	50.804	121.505	14.301
7	50.967	121.537	17.105	31	51.016	121.632	9.176
8	50.808	121.507	16.445	32	50.836	121.485	19.227
9	50.842	121.672	10.380	33	50.941	121.583	8.219
10	50.897	121.51	12.905	34	50.84	121.684	15.023
11	50.927	121.663	8.797	35	50.889	121.397	13.299
12	50.78	121.537	5.747	36	50.84	121.433	12.229
13	50.937	121.667	13.961	37	51.016	121.577	11.145
14	50.881	121.469	15.552	38	50.942	121.602	13.748
15	50.798	121.667	14.784	39	50.996	121.549	15.553
16	50.785	121.54	10.863	40	50.967	121.593	8.317
17	51.062	121.631	4.073	41	50.938	121.546	10.759
18	50.792	121.555	17.389	42	50.954	121.607	13.482
19	51.07	121.599	8.605	43	50.981	121.599	11.562
20	50.769	121.546	6.032	44	50.842	121.52	14.556
21	50.82	121.609	7.873	45	50.925	121.595	8.120
22	50.774	121.548	3.510	46	50.816	121.619	10.819
23	50.828	121.569	14.178	47	50.78	121.429	4.607

.

2.4. Research Methods

To address the limitations of the traditional spaceborne laser waveform decomposition algorithms, such as the strong dependency on the initial parameters, susceptibility to local optima, and difficulty in resolving hidden peaks during waveform overlap—resulting in low decomposition accuracy—this paper proposes the HAGA method. This method provides more accurate and reliable waveform parameters for the application of full-waveform data. As shown in Figure 4, the HAGA method consists of four main steps. First, the waveform is preprocessed to remove background noise and Gaussian white noise by thresholding and Gaussian filtering. Second, a hidden peak recognition method is applied to identify all the potential waveform components and extract their corresponding initial parameters, ensuring that all the information within the waveform is captured. Third, the AGA is employed to optimize the initial parameters, preventing the waveform fitting divergence caused by large discrepancies between the initial values and the actual parameters. Finally, the Levenberg-Marquardt (LM) algorithm is used for final optimization and fitting, yielding more precise waveform indices.

2.4.1. Full-Waveform Data Preprocessing

Full-waveform LiDAR echoes capture multiple reflections from various targets, including the ground surface, vegetation, and buildings. Each reflection can be approximately regarded as a Gaussian signal, meaning that the echo waveform can be considered the superposition of multiple Gaussian signals. Therefore, it can be modeled and fitted using a Gaussian model [5], which can be expressed in the form of Equation (1).

$$F(t)={\displaystyle \sum _{i=1}^N{A}_i}\times e^{{\displaystyle \frac{-{(t-t_i)}^2}{2{\sigma }_i^2}}}+{\mu }_{noise\hspace{0.33em}}$$

(1)

where $F(t)$ is the full-waveform echo; $N$ is the number of Gaussian components in the echo; $A_i$, $t_i$, and ${\sigma }_i$ are the amplitude, position, and half-width of the i-th Gaussian component, respectively; and ${\mu }_{noise\hspace{0.33em}}$ is the noise in the original full-waveform echo. Each Gaussian component is determined by the scattering characteristics of the corresponding target. Waveform fitting involves iteratively searching for the optimal solution of $N$ × 3 Gaussian parameters and one noise parameter through n sampling points, such that the fitting waveform minimizes the discrepancy with the original waveform.

During the acquisition of full-waveform echoes, noise inevitably arises from factors such as the background light and the sensor itself. Therefore, the denoising and filtering of the raw waveform are necessary before fitting and decomposition.

To ensure that no data are missed, the recorded full-waveform echo signal often exceeds the effective range. As a result, the waveform can generally be divided into a background-noise-only section and an effective section mixed with noise. By estimating the noise component, its influence can be eliminated. The beginning and ending portions of the full-waveform echo signal have a low probability of containing effective signals. Therefore, these portions are selected to calculate the noise mean ${\mu }_{noise\hspace{0.33em}}$ and standard deviation ${\sigma }_{noise}$ [31]. The noise threshold, $t{h}_{noise}$, is calculated using the following formula:

$$t{h}_{noise}={\mu }_{noise}+n{\sigma }_{noise}$$

(2)

where n ranges from 3 to 4.5.

The noise in the effective portion of the full-waveform echo signal can be considered Gaussian white noise. Gaussian filtering is applied to smooth the raw waveform, retain the effective waveform, and improve the signal-to-noise ratio. The Gaussian filtering process involves the convolution of the signal $s(t)$ with the filtering function $g(t)$ and is expressed as

$$s(t)\ast g(t)={\displaystyle \sum _{i=1}^n\frac{A_i{\sigma }_i}{\sqrt{{\sigma }^2+{\sigma }_i^2}}}{e}^{{\displaystyle \frac{-{(t-t_i)}^2}{2({\sigma }^2+{\sigma }_i^2)}}}$$

(3)

where $\sigma$ and $t$ are the half-wave width and position of the filter function, respectively.

2.4.2. Hidden Peak Recognition

For subsequent waveform fittings, the initial parameters of the waveform need to be estimated, and the peak points and inflection points in the echo data need to be extracted through the principle of numerical differentiation. Although most of the random noise present in the waveform is eliminated from the original waveform data after data preprocessing, some of the residual noise may lead to errors in the peak and inflection point extractions, and the error points need to be excluded according to the noise threshold $t{h}_{noise}$.

For simple full-waveform echoes, where each Gaussian component is independent, such as G1 and G3 in Figure 5, the number and positions of the waveform components can be determined by identifying the peak points and inflection points. However, for overlapping components—where the distance between two Gaussian components is too close, causing the peak of one Gaussian component to be hidden by another, as shown with G2 and G3 in Figure 5—an alternative method is adopted. First, the peak point of a simple waveform is identified in the signal, such as $A1(t_{A1},A_{A1})$. From the peak point, both directions are searched for the first encountered inflection points, $I1(t_{I1},A_{I1})$ and $I2(t_{I2},A_{I2})$. The amplitude $A_{A1}$, position $t_{A1}$, and half-width ${\displaystyle \frac{(t_{I1}+t_{I2})}{2}}$ can then be used as the initial values for the Gaussian component G1. Next, the presence of hidden peaks, such as G2, is checked. Within the range around A2, which exceeds the threshold, the inflection points are searched on both sides. If the number of inflection points on the same side is greater than or equal to three, such as $I3(t_{I3},A_{I3})$, $I4(t_{I4},A_{I4})$, and $I5(t_{I5},A_{I5})$, it is preliminarily determined that a hidden peak exists.

If the distance between two adjacent inflection points exceeds the half-width of the transmitted waveform, it indicates the presence of a hidden peak within this range. As this is a rough estimation of the initial parameters, high precision is not needed. Take G2 as an example

$$A_{G2}=\left\{\begin{array}{cc}{A}_{I4}& A_{I4}>A_{I3}\\ A_{I3}& A_{I3}>A_{I4}\end{array}\right.$$

(4)

$$t_{G2}={\displaystyle \frac{(t_{I3}+t_{I4})}{2}}$$

(5)

$${\sigma }_{G2}=\left|\begin{array}{c}{\displaystyle \frac{t_{I4}-t_{I3}}{2}}\end{array}\right|$$

(6)

where $A_{G2}$, $t_{G2}$, and ${\sigma }_{G2}$ are the initial parameters of the Gaussian component G2. If $A_{G2}$ > $t{h}_{noise}$, G2 is retained.

2.4.3. Initial Parameter Preliminary Optimization

Owing to the significant deviation between the initial parameters estimated in the previous step, especially for the hidden peaks and expected values, the fitting process may converge to a local optimum or even diverge. Therefore, the adaptive genetic algorithm (AGA) is used for preliminary optimization to bring the initial parameters as close as possible to the expected values, ensuring the quality of the final fit. The AGA optimization process is illustrated in Figure 6.

The fitness $f$ of all the individuals is calculated by calculating the sum of the fitted residuals, and those with high fitness values are selected to enter the next generation; the crossover probability $P$ and the variance probability $g$ of the individuals are determined via a sigmoid adaptive strategy [30].

$$P=\left\{\begin{array}{cc}{P}_{max}& f<f_{avg}\\ P_{min}+{\displaystyle \frac{P_{max}-P_{min}}{1+e^{a\times \frac{x\times \left(f-f_{avg}\right)}{f_{max}-f_{avg}}}}}& f\ge f_{avg}\end{array}\right.$$

(7)

$$g=\left\{\begin{array}{cc}{g}_{max}& f<f_{avg}\\ g_{min}+{\displaystyle \frac{g_{max}-g_{min}}{1+e^{a\times \frac{x\times \left(f-f_{avg}\right)}{f_{max}-f_{avg}}}}}& f\ge f_{avg}\end{array}\right.$$

(8)

where $f$ is the larger fitness of the two individuals; $f_{avg}$ and $f_{\max }$ are the average and maximum fitness of the entire individual, respectively; $P_{max}$, $P_{\min }$, $g_{\max }$, and $g_{\min }$ are the maximum and minimum cross-variance probabilities, respectively; and $a$ = 9.903438.

The calculated probability is used to determine whether to perform the cross-mutation operation, and the next generation is obtained in accordance with Equations (9) and (10) in the case of cross-mutation.

$$\begin{array}{c}{x}_A^{t+1}=\alpha x_A^t+\left(1-\alpha \right)x_B^t\\ x_B^{t+1}=\alpha x_B^t+\left(1-\alpha \right)x_A^t\end{array}$$

(9)

$$x_C^{t+1}=\left\{\begin{array}{c}{x}_C^t+\beta \times \left(x_{max}-x_c^t\right)\\ x_C^t-\beta \times \left(x_c^t-x_{min}\right)\end{array}\right.$$

(10)

where $x_A^{t+1}$ and $x_B^{t+1}$ are individuals after the crossover operation; $x_C^{t+1}$ are individuals after the mutation operation; $x_A^t$, $x_B^t$, and $x_C^t$ are randomly selected parent individuals; and $\alpha$ and $\beta$ are crossover variation constants, with values ranging from [0, 1].

When the number of iterations exceeds n (where n is the maximum number of iterations set, with n = 50 in this study) or the required accuracy is achieved (waveform fitting residuals < three times the standard deviation of the waveform), the AGA terminates and outputs the individual with the highest fitness in the population, completing the preliminary optimization of the initial parameters.

2.4.4. Optimal Fit

After the initial parameters are extracted and preliminarily optimized, the Gaussian parameters that are closer to the expected values are obtained. These parameters are then used for nonlinear least squares fitting to derive more precise waveform parameters.

The Levenberg-Marquardt (LM) algorithm [32,33,34] is a widely used method for optimizing waveform data fitting. It combines the advantages of the Gauss-Newton (GN) method and the gradient descent method. When close to the convergence point, it exhibits fast convergence characteristics; when far from the convergence point, it demonstrates strong adaptability. With the fast convergence of the GN method and the robustness of the gradient descent method, they enable more efficient and rapid fitting optimization.

Using the initial parameters optimized by the adaptive genetic algorithm (AGA), LM optimization fitting is performed. The iteration stops when the fitting residual between the fitted waveform and the original waveform meets the threshold or the maximum number of iterations is reached. This process yields Gaussian parameters with higher accuracy and derives the expression in the form of Equation (1), thereby completing the waveform decomposition.

2.4.5. Evaluation Indicators

The accuracy of waveform decomposition fitting was evaluated using the coefficient of determination $R^2$, the root mean square error ($RMSE$), and its standard deviation $\sigma$ to assess the fitting performance. Both $R^2$ and $RMSE$ adopted data self-validation, where the fitted echo waveform was compared with the original echo data for calculation [23]. The mean tree height within the footprint area was calculated using Equation (13) [35], and the accuracy of the extracted mean tree height was assessed using the mean absolute error ($MAE$), standard deviation ${\sigma }_h$, $RMSE$, and 95% confidence intervals (CIs).

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

(11)

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{\left(A_i-a_i\right)}^2}}$$

(12)

$$H=(Ground\_Peak-Canopy\_Peak)\times BinSize$$

(13)

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

(14)

$${\sigma }_h=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^n{(\Delta h_i-\overline{\Delta h})}^2}}{n}}}$$

(15)

$$CI=\overline{\Delta h}\pm 1.96\times {\displaystyle \frac{{\sigma }_h}{\sqrt{n}}}$$

(16)

where $f_i$ and $y_i$ denote the intensity values of the fitted waveform and the original waveform at the i-th point, respectively; and $\overline{y}$ represents the mean intensity value of the original waveform. The $RMSE$ used to evaluate the accuracy of waveform decomposition fitting is denoted as $RMS{E}_y$, which was calculated by setting $A_i$ = $f_i$ and $a_i$ = $y_i$. $H$ represents the average tree height within the laser footprint, $Ground\_Peak$ denotes the peak position of the ground return, and $Canopy\_Peak$ represents the peak position of the canopy return. $BinSize$ corresponds to the distance represented by each data bin. The sampling interval of each frame of the “Goumang” satellite is 0.833 ns, so the bin size was calculated as $BinSize={\displaystyle \frac{1}{2}}\mathrm{c}\times \mathrm{t}=0.125 \mathrm{m}$, where c represents the speed of light. When evaluating the accuracy of the average tree height, $RMSE$ was denoted as $RMS{E}_h$, which was calculated by setting $A_i$ = ${\widehat{h}}_i$ and $a_i$ = $h_i$, where ${\widehat{h}}_i$ is the average tree height extracted from the waveform, and $h_i$ is the average tree height derived from airborne LiDAR point cloud data. The average tree height error is defined as $\Delta h_i={\widehat{h}}_i-h_i$, with its mean value represented as $\overline{\Delta h}$.

3. Results

In this study, based on the selected full-waveform data of the “Goumang” LiDAR system from the Genhe study area (a total of 2533 laser points from August to October 2022, as shown in Figure 2) and 49 airborne LiDAR plot s(as shown in Figure 3), the HAGA method was employed to decompose waveforms and extract parameters, verifying the reliability and practicality of this waveform decomposition approach.

3.1. Waveform Decomposition Results and Accuracy Evaluation

We used the HAGA method to decompose the “Goumang” full-waveform data from the Genhe study area. The same dataset was also decomposed using the ODM and GA methods for comparison. Figure 7 shows the results of decomposing the full-waveform data containing different numbers of waveform components using the various methods. Figure 8 illustrates the results when the hidden peaks were not identified.

For simple waveforms such as those in Figure 7a, the Gaussian components and initial parameters contained in the original waveform were accurately extracted. At this stage, all methods demonstrated good decomposition and fitting performance, with the similarity $R^2$ between the fitted waveform and the original waveform exceeding 0.99. However, as the number of Gaussian components in the original waveform increased, the waveform became more complex. For waveforms with consistent Gaussian component identification, as shown in Figure 7b–f, discrepancies began to appear between the fitted waveforms and the original waveforms across the three methods, and the differences in the $R^2$ values became more pronounced. Nonetheless, the HAGA consistently achieved the highest $R^2$ values. However, as shown in Figure 7d, the $R^2$ values of all three methods were less than 0.99. This was because some waveform components with lower peak values were too similar to noise and were thus discarded. When all the waveform components were identified but the parameters extracted from the waveform deviated significantly from the true values, as shown in Figure 7e,f, the ODM exhibited the poorest fitting performance, with $R^2$ dropping below 0.9. The GA performed better in waveform fitting. The HAGA, on the other hand, achieved the best fitting performance, with its fitted waveform almost completely overlapping the original waveform. As shown in Figure 8, when the hidden peaks were not identified, the $R^2$ of the waveform decreased. The HAGA successfully identified all the hidden peaks, achieving the highest decomposition and fitting accuracy. However, since the ODM and GA each missed one waveform component, resulting in significant differences between the extracted waveform parameters and the true values, both the ODM and GA showed poorer decomposition and fitting performance. The decomposition and fitting of all the waveform data are shown below. Figure 7 and Figure 8 show a random selection of the waveforms from the entire dataset. From the figures, it can be seen that the final fitted waveform closely matched the original waveform, which was the goal of waveform fitting—to make the two as similar as possible, ensuring that no overfitting occurred. The results for all the waveforms are presented in

Table 5. Statistics of the initial position fitting results.

Method	${\mathit{R}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	${\mathit{\sigma }}_{{\mathit{R}}^2}$	${\mathit{\sigma }}_{\mathit{R}\mathit{M}\mathit{S}\mathit{E}}$
HAGA	0.955	8.071	0.065	7.94
ODM	0.916	10.049	0.151	12.1
GA	0.862	17.293	0.092	11.84

. $R_{mean}^2$ and $RMS{E}_{mean}^2$ represent the mean $R^2$ and $RMS{E}_y$ across 2533 laser spots, respectively, whereas ${\sigma }_{R^2}$ and ${\sigma }_{RMSE}$ represent their standard deviations.

Table 5 shows that the decomposition and fitting results obtained via the HAGA yielded an $R_{mean}^2$ = 0.955 and an $RMS{E}_{mean}^2$ = 8.071. Compared with ODM, these values were improvements of 0.093 in $R_{mean}^2$ and 9.222 in $RMS{E}_{mean}^2$. Compared with those of the GA, the improvements were 0.039 for $R_{mean}^2$ and 1.978 for $RMS{E}_{mean}^2$. Additionally, the standard deviations of both $R_{mean}^2$ and $RMS{E}_{mean}^2$ under the HAGA method were the smallest, indicating that the HAGA decomposition method provided a more stable and reliable fitting performance.

To further verify the stability of the HAGA method under inaccurate initial parameters, the peak position $t_i$ in the extracted initial parameters was changed, and it was shifted left and right by 5, 10, and 15 frames to observe the fitting of the waveforms and evaluate the decomposition fitting results.

Table 6. Statistics of the fitting results of $t_i$ offset by 5 frames.

Method	${\mathit{R}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	${\mathit{\sigma }}_{{\mathit{R}}^2}$	${\mathit{\sigma }}_{\mathit{R}\mathit{M}\mathit{S}\mathit{E}}$
HAGA	0.936	9.74	0.1	10.135
ODM	0.852	14.315	0.232	17.851
GA	0.794	22.087	0.107	14.614

,

Table 7. Statistics of the fitting results of $t_i$ offset by 10 frames.

Method	${\mathit{R}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	${\mathit{\sigma }}_{{\mathit{R}}^2}$	${\mathit{\sigma }}_{\mathit{R}\mathit{M}\mathit{S}\mathit{E}}$
HAGA	0.928	9.87	0.123	11.334
ODM	0.454	34.16	0.413	31.571
GA	0.706	25.987	0.143	17.295

, and

Table 8. Statistics of the fitting results of $t_i$ offset by 15 frames.

Method	${\mathit{R}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}^2$	${\mathit{\sigma }}_{{\mathit{R}}^2}$	${\mathit{\sigma }}_{\mathit{R}\mathit{M}\mathit{S}\mathit{E}}$
HAGA	0.913	10.663	0.149	12.662
ODM	0.064	46.806	0.403	33.818
GA	0.573	31.12	0.182	20.67

show the results when $t_i$ was shifted by 5, 10 and 15 frames, respectively, and Figure 9 shows the comparison of the results after shifting the initial position.

As shown in Figure 9, when the initial positions of the peak values were shifted, the fitting accuracy of all three methods decreased. Among them, the HAGA resulted in the smallest decline, demonstrating superior robustness. In particular, the ODM method experienced a rapid decrease in fitting accuracy when the initial position shift exceeded five frames. When the shift reached 15 frames, $R_{mean}^2$ approached 0, indicating that the fitting results were completely unreliable. Table 5, Table 6, Table 7 and Table 8 show that the HAGA was insensitive to the initial values. When the initial peak position was shifted by five frames, the $R_{mean}^2$ of the HAGA method decreased by 0.019, and the $R_{mean}^2$ increased by 1.669. In contrast, the $R_{mean}^2$ values of the ODM and GA methods decreased by 0.064 and 0.068, respectively, whereas their $RMS{E}_{mean}^2$ values increased by 4.266 and 4.794, respectively. When the peak position shift reached 10 or 15 frames, the changes in $R_{mean}^2$ and $RMS{E}_{mean}^2$ for the HAGA remained relatively small compared with the initial position, whereas significant changes were observed in the $R_{mean}^2$ and $RMS{E}_{mean}^2$ of the ODM and GA. Furthermore, the HAGA method resulted in smaller standard deviations regardless of whether the parameters were original or shifted, indicating that its results were more stable for data processing.

3.2. Evaluation of Tree Height Extraction Accuracy

To further validate the reliability of the HAGA method in practical applications, “Goumang” full-waveform data from the study area were used to calculate the average tree height. Moreover, the ODM and GA methods were employed to extract the average tree height, which was then compared with the tree height extracted from airborne LiDAR data. The reason for using mean tree height for validation is that due to the complexity of the environment, the ground echo in the raw waveform tends to merge with echoes from low-lying shrubs, making it difficult to identify the ground echo and consequently affecting the accuracy of the mean tree height estimate. The HAGA method can identify hidden peaks and preserve more of the information in the raw waveform, allowing the ground echo to be identified, resulting in a more accurate estimate of average tree height [36,37]. The results and errors of the average tree height calculated by different decomposition methods compared with the average tree height extracted from airborne LiDAR data are shown in Figure 10.

Figure 10b shows that over 60% of the errors across the three decomposition methods were negative, indicating an underestimation of the average tree height. Figure 10a shows that, among the three methods (ODM, GA, and HAGA), the HAGA exhibited the best overall performance, with the smallest difference between the estimated Hmean and the LiDAR-measured tree heights. This improvement was particularly significant in sample plots 26 and 31, where the estimation errors were greatly reduced. The reason lies in the HAGA’s ability to fit the original waveform accurately, thereby providing more precise forest-sensitive waveform parameters.

To evaluate the overall practical applicability of the different methods, the $MAE$, MAX, MIN, ${\sigma }_h$, $RMS{E}_h$, and CI for each method were calculated, as shown in

Table 9. Comparison of the average tree height extraction accuracies of the three methods.

Method	$\mathit{M}\mathit{A}\mathit{E}$/(m)	MAX/(m)	MIN/(m)	${\mathit{\sigma }}_{\mathit{h}}$/(m)	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{h}}$/(m)	CI
ODM	2.03	9.008	−4.277	2.77	2.74	[−0.73, 0.837]
GA	1.795	9.574	−5.292	2.621	2.595	[−0.683, 0.801]
HAGA	1.499	3.472	−4.171	1.816	1.899	[−1.126, −0.099]

.

As shown in Table 9, the MAE of the average tree height extracted by the HAGA method was 1.499 m, which was 0.531 m and 0.296 m lower than that of the ODM and GA methods, respectively. The ${\sigma }_h$ and $RMS{E}_h$ values of the three methods showed exactly the same trend, indicating that the HAGA method outperformed the other methods in terms of estimation accuracy and stability. The HAGA method had a maximum value of 3.472 m, a minimum value of −4.171 m, and a CI of [−1.126, −0.099]. Compared to ODM and GA, it had a smaller maximum value, the largest minimum value, and the smallest confidence interval. These results indicate that the HAGA method has a lower dispersion and uncertainty in the estimation of mean tree height. For these evaluation metrics, the HAGA consistently had the smallest values, a narrower range between the maximum and minimum values, and a narrower confidence interval, demonstrating greater consistency and stability. This indicates that HAGA achieves the best accuracy and robustness in tree height estimation, making it suitable for most scenarios.

4. Discussion

Through the analysis in Section 3, it was demonstrated that the proposed HAGA method outperforms the ODM and GA in terms of accuracy. However, in practical applications, the success rate of waveform decomposition is also a prerequisite for evaluating the applicability of a method. This section further discusses the success rates of the different methods for waveform decomposition, combined with accuracy, to assess the practicality of these methods. Waveform decomposition is considered successful when the $R^2$ of the fitted waveform exceeds 0.8 [38]. The number of successfully decomposed waveforms is shown in Figure 11.

The experimental results show that the number of converged waveform decompositions fits is close to 100% on the same test set, but the number of converged decompositions reflects the mathematical properties of the algorithm rather than the ability to effectively extract the real physical features of the waveform and, more importantly, the success rate of waveform decomposition. The success rate of waveform decomposition directly determines the applicability and reliability of the method in practical applications. Figure 11 shows that the HAGA successfully decomposed 2383, and its decomposition success rate was 94.1%, which was significantly better than those of ODM (80.8%) and the GA (79.7%). The HAGA was optimal in terms of both decomposition accuracy and decomposition success rate, whereas the ODM had a relatively high number of successes and accuracy, despite a low number of decomposition convergences; GA had the highest number of decomposition convergences, but the number of decomposition successes and its overall accuracy were the lowest. Since the accuracy of waveform decomposition is affected by the accuracy of the initial parameters and whether all the waveform components are identified, the HAGA method has a strong ability to recognize waveform components and is less sensitive to the initial parameters.

Considering both the decomposition accuracy and the number of successful decompositions, the HAGA method achieves a good balance between accuracy and the success rate. This is attributed to the algorithm’s strong ability to identify waveform components, and its insensitivity to the initial parameters, making the HAGA highly practical for processing waveform data. This contrast between the number of successful decompositions and the accuracy actually reflects the differences in processing between the different methods. The HAGA method focuses on the combined optimization of accuracy and decomposition success, whereas the GA method focuses more on the wide coverage of decomposition, which often comes at the cost of a loss of accuracy. In contrast, the ODM method sacrifices the decomposition range for higher accuracy, but its applicability is relatively limited, such as in flat terrain and forests with low forest cover.

To further consider the relationship between the decomposition accuracy, success rate, and waveform noise, a signal-to-noise ratio grading evaluation system was employed (as shown in

Table 10. Relationship between waveform decomposition success rate and waveform SNR.

Method	10 < SNR < 11	11 < SNR < 12	12 < SNR < 13	13 < SNR < 14	14 < SNR < 15
HAGA	91%	91%	95%	94%	96%
ODM	64%	70%	82%	87%	91%
GA	74%	61%	76%	91%	86%

/Figure 12).

In the 10–15 dB range, the decomposition success rates of all three methods decreased with decreasing SNR. When the SNR was in the 14–15 dB range, the decomposition success rate of the HAGA reached 96%, significantly higher than those of ODM and GA, which were 91% and 86%, respectively. When the SNR dropped to the 10–11 dB range, the HAGA still maintained a success rate of 91%, while ODM and GA dropped to 64% and 74%, respectively. In terms of trends, the decomposition success rate of the ODM method decreased monotonically with increasing noise, while GA showed fluctuations. Although the HAGA also decreased, it showed greater stability. Overall, the variation in R² further confirms this result. As the noise increased, the R² values decreased for all methods. However, the HAGA maintained values of 0.93, 0.93, 0.94, 0.95, and 0.96 over different intervals, while ODM and GA remained below 0.9. This advantage is due to the HAGA’s ability to identify hidden peaks within the waveform and detect potential initial parameters. In addition, it uses probabilistic modeling rather than deterministic rules, reducing the dependence on initial parameters and improving overall robustness.

Overall, the HAGA method outperforms the other methods in terms of its general performance, ensuring high decomposition accuracy while maintaining a reasonable decomposition success rate. This approach strikes a balance between accuracy and applicability. However, under complex environmental conditions, where the waveform signal-to-noise ratio decreases, the decomposition accuracy of the HAGA method deteriorates. Future improvements should focus on noise handling to enable the method to maintain high accuracy while adapting to more complex scenarios.

5. Conclusions

To address the challenges of hidden peak identification and sensitivity to initial parameters in the processing of full-waveform LiDAR echo data, this paper proposes the HAGA decomposition method. This method first uses the inflection point method to identify hidden peaks in the waveform, thereby extracting the initial parameters of all potential waveform components. AGA is then used to optimize these initial parameters, reducing the loss in fitting accuracy caused by inaccurate initial parameters of traditional methods. Finally, the LM algorithm is combined with further optimized fitting to achieve the high-precision decomposition of complex waveforms. According to the comparison of the tree height extraction results obtained from “Goumang” multibeam LiDAR data and airborne LiDAR data, the HAGA method outperforms ODM and GA in terms of decomposition success rate and R², and its extracted mean tree heights have higher correlations with the measured data and a lower mean square error. The HAGA method overcomes the strong dependence of the traditional optimization algorithm on the initial parameters and, at the same time, has higher decomposition accuracy for complex waveforms (e.g., superimposed waveforms or asymmetric waveforms), which fully proves its effectiveness in practical applications. Future research will further improve the noise suppression strategy and further optimize the method so that it can adapt to more complex waveforms and environmental conditions while maintaining high accuracy.