A Method for Quantifying Understory Leaf Area Index in a Temperate Forest through Combining Small Footprint Full-Waveform and Point Cloud LiDAR Data

Abstract

Understory vegetation plays an important role in the structure and function of forest ecosystems. Light detection and ranging (LiDAR) can provide understory information in the form of either point cloud or full-waveform data. Point cloud data have a remarkable ability to represent the three-dimensional structures of vegetation, while full-waveform data contain more detailed information on the interactions between laser pulses and vegetation; both types have been widely used to estimate various forest canopy structural parameters, including leaf area index (LAI). Here, we present a new method for quantifying understory LAI in a temperate forest by combining the advantages of both types of LiDAR data. To achieve this, we first estimated the vertical distribution of the gap probability using point cloud data to automatically determine the height boundary between overstory and understory vegetation at the plot level. We then deconvolved the full-waveform data to remove the blurring effect caused by the system pulse to restore the vertical resolution of the LiDAR system. Subsequently, we decomposed the deconvolved data and integrated the plot-level boundary height to differentiate the waveform components returned from the overstory, understory, and soil layers. Finally, we modified the basic LiDAR equations introducing understory leaf spectral information to quantify the understory LAI. Our results, which were validated against ground-based measurements, show that the new method produced a good estimation of the understory LAI with an R² of 0.54 and a root-mean-square error (RMSE) of 0.21. Our study demonstrates that the understory LAI can be successfully quantified through the combined use of point cloud and full-waveform LiDAR data.

1. Introduction

Vegetation in forests often displays distinct layering [1], which can be divided into overstory and understory. The overstory refers to the uppermost layer made up of woody plants, while the understory is located beneath the overstory, typically consisting of trees stunted by lack of light, other small trees with low light requirements, saplings, vines, and undergrowth [2]. Between the overstory and understory, a gap height interval often exists (Figure 1).

In view of the great financial value and critical ecological functions, the retrieval of overstory vegetation attributes has been advanced considerably with the help of remote sensing techniques [3,4,5]. In contrast, the characteristics of understory vegetation remain poorly explored to date, despite its important impacts on a wide range of processes in forest ecosystems. For example, it aids in soil nutrient cycling and soil structure maintenance, protecting soil and water from erosion [1]. It provides wildlife with habitats and forage, contributing to species diversity [6]. It affects overstory regeneration and thus changes the rate and direction of forest succession [7]. It also contributes to carbon stocks, affecting Aboveground biomass (AGB) and net primary production (NPP) estimation [8]. In addition, it plays an important role in driving the transformation from a surface fire to a crown fire, which can result in enormous economic losses [9]. Thus, quantifying understory characteristics is of great significance.

The leaf area index (LAI) has been an extensively studied indicator for the quantification of vegetation [10,11,12,13,14,15]. Watson first defined it as the total one side leaf area per unit ground surface area [10]. Later, Chen and Black proposed a more common definition that LAI is half of the total leaf area per unit ground surface area [11]. LAI is a key vegetation structure parameter controlling the biological and physical processes of forests. It regulates primary leaf physiological processes such as photosynthesis, respiration, and transpiration [12]. LAI is highly correlated with carbon content, NPP, and the fraction of absorbed photosynthetically active radiation (FPAR) [16]. Furthermore, LAI is critical to the processes of carbon, water, and energy exchange between land and atmosphere [12,17]. It has been utilized as an important input or output parameter in many ecological, hydrological, and climate models [18,19,20]. However, failure to consider the understory LAI contribution can lead to large uncertainties in assessing forest properties and modeling surface processes [13,21]. Moreover, quantifying the understory LAI is vital for biodiversity conservation since, motivated by the significant role of the LAI in ecosystem functioning, the remote sensing and ecology communities have deemed it one of the crucial biodiversity variables for tracking from space the progress toward meeting the Aichi Biodiversity Targets [22]. For these reasons, it is critical to be able to accurately quantify the understory LAI of forest ecosystems.

Despite its importance, quantifying the understory LAI has conventionally been difficult. Traditional field surveys are very costly and limited to small areas. Passive remote sensing techniques are also insufficient due to their inherently limited ability to penetrate forest canopies and provide vertical structure information [23]. In contrast, more recently developed active remote sensing techniques, such as light detection and ranging (LiDAR), can penetrate forest canopies to varying degrees and capture a wide range of information on the three-dimensional (3D) distribution of forest structural components, such as canopy height, timber volume, and aboveground biomass [24,25].

Generally, LiDAR data can be classified into point clouds and waveforms, depending on how the reflected signals are recorded and processed. Discrete LiDAR systems (or LiDAR system in discrete mode) record the exact location and intensity information of the detected peaks, and the resulting point cloud LiDAR data have been widely used to estimate forest structural properties, including mapping understory vegetation height [26], discriminating vegetation strata [27], estimating understory cover [28], and characterizing understory vegetation density [29]. Although these studies did not provide a direct method for quantifying understory LAI, they strongly indicated that point cloud data can be utilized to distinguish overstory and understory vegetation.

Compared with discrete LiDAR, full-waveform LiDAR usually records intensity at every nanosecond, instead of the discrete peaks, of the returned pulse, which provides more detailed information on the interaction between LiDAR pulse and canopy components. Therefore, the resulting full-waveform data provide new possibilities for a comprehensive understanding of the mechanism of pulse–vegetation interaction and for further inversion of the vegetation parameters, provided that sophisticated models are established, such as DART [30], HELIOS [31], and a waveform simulation model developed by Sun et al. [32]. Based on these models, more detailed forest structural properties, such as the vertical distribution of the foliage area volume density (FAVD), can be estimated from full-waveform data [33]. The vertical distribution of gap probability of heterogeneous and discrete canopies is the key concept modeling LiDAR waveforms, as it links the canopy structure with the LiDAR waveforms through the lidar equation, which provides us with the means to estimate the vertical gap probabilities from LiDAR full-waveform data, with some knowledge of the foliage and ground reflectance ratio [34]. The LAI and the cumulative vertical LAI profile over a tropical rainforest were retrieved by the vertical gap probabilities using Laser Vegetation Imaging Sensor (LVIS) data [35]. Then the vertical foliage profile (VFP) over the United States was further characterized on the basis of Geoscience Laser Altimeter System (GLAS) data [36]. However, none of these studies separated the understory LAI from the total LAI.

Notably, Li et al. [37] applied the basic LiDAR equations to retrieve both the understory LAI and the overstory LAI of an orchard by combining full-waveform and point cloud LiDAR data. The two layers were separated by acquiring the height boundary in a field survey, which is time-consuming and laborious work. In addition, the differences in leaf optical properties between the two layers were not considered in that study.

The aim of this study was to develop a new method for quantifying the understory LAI in temperate forests by combining the advantages of both point cloud and full-waveform LiDAR data. Specifically, we set out to (1) determine the height boundary between the overstory and understory vegetation using point cloud data at the plot level; (2) integrate the boundary information derived from point cloud data with full-waveform data to distinguish the waveform components returned from the overstory, understory, and soil layers; and (3) modify the basic LiDAR equations proposed by Ni-Meister for understory LAI estimation using the discriminated waveform components.

2. Materials

2.1. Study Area and Field Data

The study area is located at the Saihanba National Forest Park (42.35° N, 117.32° E), Chengde, Hebei Province, China (Figure 2). The park covers a total area of 94,000 hectares, and its elevation ranges from 1100 m to 1900 m. The main soil types are aeolian sandy soil, meadow soil, and swamp soil [38]. The park has a mean annual precipitation of 450 mm and an average temperature of −1.4 °C, with a short growing season (between May and September) and a long winter (from November to March of the next year) [39]. The dominant overstory species are Mongolian pine (Pinus sylvestris var. mongolica Litv.), Northern China’s larch (Larix principis-rupprechtii Mayr), and birch (Betula platyphylla). The understory layer is dominated by Carex rigescens, Thalictrum aquilegifolium, Galium verum, Geum aleppicum, Artemisia tanacetifolia, and Agrimonia pilosa [39].

Field surveys, including LAI and spectrum measurements, were carried out from late August to early September during the peak vegetation growing seasons in 2018 and 2019. Here, we neglect the slight differences in forest properties between the two consecutive years. LAI was measured in 16 plots (in Figure 1), each with the same area of 25 m × 25 m. The sampling scheme within each plot was a cross formed by 5 measurement points: two perpendicular 3-point transects, with the central one overlapping. The smartphone positioning software Ovital Map (version 7.5.8) (https://www.ovital.com/download/ (accessed on 4 May 2018)) was used to record these points (with an accuracy of between 1 and 5 m). Considering the low accuracy of the positioning software, all the plots we chose were relatively homogeneous within the range of error.

A Nikon fisheye camera was used to collect downward digital hemispherical photographs (DHPs) for the acquisition of the understory LAI. The camera was operated by hand at 1.5 m above ground during overcast sky conditions or, alternatively, in the very early morning or at dusk to avoid direct illumination from the sun. The hemispherical photographs were processed using CAN_EYE software (version 6.4.91) (http://www6.paca.inra.fr/can_eye (accessed on 10 June 2018)). When using this software, users can discriminate image pixels as belonging to either vegetation elements or background. The background corresponds to soil pixels in downward images. The undesired parts of the images, such as the feet of the operator, were interactively masked out before discrimination. Based on the classification results, gap fractions in different zenith and azimuth directions were calculated to generate understory LAI values. The plot-level understory LAI was calculated as the mean value derived from DHPs within each plot, holding a distribution range of between 0.32 and 1.47 (

Table 1. Statistics of in situ understory LAI measurements at the plot level (plot number = 16).

	Understory LAI
Min	0.32
Max	1.47
Mean	0.86
SD	0.31
Range	1.15

). The mean value of the plot-level understory LAI was 0.86, with a stand deviation value of 0.31.

Understory vegetation, overstory vegetation, and soil spectra were measured by using a handheld SVC HR-1024 Spectroradiometer (Spectra Vista Corporation, Poughkeepsie, NY, USA) during the survey. For the understory vegetation, seven samples of the dominant species of Carex rigescens mixed with a few other species were measured, and the reflectance was calculated to have an average value of 0.21 (s.d. = 0.04) at a wavelength of 1550 nm (equal to the wavelength of the LiDAR system introduced in Section 2.2.). For the soil, only three measurements were available (0.38, 0.40, and 0.33 at 1550 nm), and we used their average value of 0.37 as the ground reflectance. The overstory reflectance spectrum was measured using an ASD FieldSpec spectrometer (Analytical Spectral Devices, Boulder, CO, USA) in combination with an ASD integrating sphere. Leaves from three dominant overstory species were picked by an experimenter standing on a crane and were immediately stored in zipper-close plastic bags before measurement in the laboratory; for each species, nine leaves were collected separately from the upper, middle, and lower crowns of three trees. The reflectance of each overstory species was calculated as the mean value of the nine samples at 1550 nm; the value for Mongolian pine was 0.16 (s.d. = 0.02), the value for Northern China’s larch was 0.25 (s.d. = 0.02), and the value for birch was 0.31 (s.d. = 0.02).

2.2. LiDAR Data

LiDAR data were obtained in the study area in early September 2018 using the LiCHy system [40] by the Chinese Academy of Forestry (CAF). The full-waveform sensor (Riegl LMS-Q680i) included in this system has a wavelength of 1550 nm and a laser pulse width of 3 ns. The waveform data were recorded at intervals of 1 ns, resulting in approximately 15 cm per waveform bin. The flying altitude was approximately 1000 m above ground, and the diameter of the footprint was approximately 0.25 m. Point cloud data were acquired at a density of 4.8 pts/m² without considering multiple returns. The corresponding processing steps are described in the next section.

3. Methods

The main methodological workflow for our study is shown in Figure 3. It consisted of the determination of the height boundary between the overstory and understory vegetation at the plot level, the preprocessing of the full-waveform data, and the modification of basic LiDAR equations for quantifying the understory LAI. Further details of each step are described in the next subsections.

3.1. Determination of the Overstory–Understory Height Boundary at the Plot Level

Since vegetation in forests often displays distinct layering, and a gap height interval can often be found between the overstory and understory layers, the overstory and understory vegetation can be distinguished through the height gap interval. Assuming that the vertical gap probability remains the same within the height interval but changes in intervals with leaves, the profile of the vertical gap probability can be used to separate the overstory and understory vegetation. However, when there is no distinct boundary between the overstory and understory vegetation, a cutoff threshold can be set to separate them. In our study, we set this threshold to 2 m according to [41].

The vertical gap probabilities at the plot level were estimated using the point cloud data normalized with respect to the height above ground. First, the point cloud data were classified as ground returns and nonground returns using an improved progressive triangulated irregular network (TIN) densification filtering algorithm [42]. Next, the elevation of each return was normalized by subtracting the elevation of the ground returns. Then, a method of using the first returns to calculate the vertical gap probabilities [4] was adopted:

$$P_{gap}\left(z\right)=1-\frac{{\sum }_z^{z_{max}\sum }n}{N}$$

(1)

where $P_{gap}\left(z\right)$ denotes the gap probability at height $z$, $N$ is the total number of emitted pulses, and $n$ is the total number of first returns that fall in the height strata between the top of the canopy and the height $z$.

To determine the optimal boundary height at the plot level, we took the derivative with respect to the gap probability along the vertical direction, $\mathrm{dp}/\mathrm{dz}$. Since we assumed that the vertical gap probability remained the same only at the height boundary between overstory and understory vegetation, the values of $\mathrm{dp}/\mathrm{dz}$ were supposed to be zero only at the height gap interval.

Let $H_b$ denote the height boundary between the overstory and understory vegetation, and let $H_b$ have a distribution range of between a and b (shown in Equation (2)). Near the ground height, ground returns and low nonground returns are likely to be misclassified, leading to the relatively low accuracy of the estimated the vertical gap probabilities. To mitigate this effect, we set the value of a as 1 m. Except for the boundary height, gap probabilities along the vertical direction are likely to remain the same near the height of the canopy top (equal to 1), so we set the value of b as 4 m (according to [43]).

$$\begin{array}{cc}a& m\end{array}\le H_b\le \begin{array}{cc}b& m\end{array}$$

(2)

3.2. Preprocessing of Full-Waveform Data

3.2.1. Waveform Deconvolution

The observed waveform is the sum of the convolutions of terms related to the system and environment contributions [44]:

$$P_r\left(t\right)={\displaystyle \sum }_{i=1}^N\frac{D^2}{4\pi {\lambda }^2{R}_i{}^4}\underset{\begin{array}{cc}system& contribution\end{array}}{\underbrace{P_t\left(t\right)\ast {\eta }_{sys}\left(t\right)}}\ast \underset{\begin{array}{cc}environment& contribution\end{array}}{\underbrace{{\eta }_{atm}\left(t\right)\ast {\sigma }_i\left(t\right)}}$$

(3)

where “$\ast$” denotes the convolution operation, $t$ is the travel time of the transmitted laser pulse, $P_r\left(t\right)$ is the observed signal as a function of time, $N$ is the number of targets that contribute to the returned waveform, $D$ is the aperture diameter of the receiver optics, $\lambda$ is the wavelength, and $R$ is the distance from the LiDAR system to the target. $P_t\left(t\right)$ is the outgoing signal, which has a certain distribution function, while ${\eta }_{sys}$ is the receiver impulse response, which is typically assumed to be Gaussian (for the Riegl LMS-Q680i, the outgoing pulse and a model of the system impulse response as calibrated in [45] are shown in Figure 4). ${\eta }_{atm}$ is the atmospheric factor, and $\sigma \left(t\right)$ is the cross section of the illuminated area. The system contribution is not a perfect delta function but rather an approximate Gaussian function with a certain width. Consequently, the convolution of the system and environmental contributions can result in a reduction of the temporal resolution of the observed signal. The signal resolution can be restored after such a decrease through waveform deconvolution algorithms.

The Richardson–Lucy (RL) algorithm [46] was employed to deconvolve the waveform data. The RL algorithm has the following advantages: a small root-mean- square error (RMSE) for recovering the true cross section, a low false discovery rate for detecting the unobservable local peaks in the stretched raw waveforms, and high classification accuracy for differentiating herbaceous biomass levels [47]. The RL algorithm was initially designed for restoring an image by iteratively searching for solutions to the deconvolution problem. Based on Bayes’ theorem, it aims to calculate the most likely value $f_t\left(x\right)$ given the observed $d\left(x\right)$ and the known point-spread function $g\left(x\right)$. A LiDAR waveform can be seen as an image with dimensions of 1 × N, and the solution from the $t\mathrm{th}$ iteration can be written as follows in convolutional form:

$$f_{t+1}\left(x\right)=f_t\left(x\right)·\left(\frac{d\left(x\right)}{f_t\left(x\right)\ast g\left(x\right)}\ast g\left(-x\right)\right)$$

(4)

where $d\left(x\right)$ is the observed value at location $x$, $f_t\left(x\right)$ is the most likely value at location $x$, $g\left(x\right)$ is the known point-spread function, and $t$ is the number of iterations. One can solve for $f\left(x\right)$ by iterating Equation (4) until convergence is reached.

3.2.2. Waveform Decomposition

Before waveform decomposition, the noise contained in the deconvolved waveform should be removed. In our study, first, a noise threshold was defined as the mean value of the last 5% of the raw waveform intensities. We subtracted this threshold from the raw waveform. Then, Savitzky–Golay (S-G) filtering [48] was applied to smooth the waveform and remove the remaining noise.

To discriminate the waveform components originating from the overstory, understory, and ground, it is essential to decompose the deconvolved waveform. Gaussian decomposition is an effective way to decompose the waveform into a sum of sub-Gaussian waveforms [49]. The corresponding analytical expression can be written as:

$$f\left(x\right)={\displaystyle \sum }_{i=1}^n{A}_i·exp(-\frac{\left(x-{\mu }_i{}^2\right)}{2{\sigma }_i{}^2})$$

(5)

where n is the number of sub-Gaussian waveforms, $A_i$ is the peak amplitude of the $i\mathrm{th}$ sub-Gaussian waveform, ${\sigma }_i$ is the standard deviation of the $i\mathrm{th}$ sub-Gaussian waveform, and ${\mu }_i$ is the location of the peak in the $i\mathrm{th}$ sub-Gaussian waveform.

The deconvolved waveform was fitted as a mixture of sub-Gaussian waveforms using the nonlinear least squares (NLS) method and the sequential quadratic programming (SQP) as the optimization algorithm [50]. To ensure the quality of the decomposition, suitable initial values of n, $A_i$, ${\sigma }_i$, and ${\mu }_i$ must be obtained. For this purpose, we calculated these parameters by first finding the local maxima of the deconvolved waveform that exceeded a certain threshold and then used the corresponding peak values as the initial amplitudes, the half widths of adjacent peaks as the initial standard deviations, and the peak locations as the corresponding initial locations.

After waveform decomposition, we took the last sub-Gaussian waveform as the ground return. Except for the ground return, the waveform components for which the peak locations were below the boundary height were taken as understory returns, and the rest were taken as overstory returns.

3.3. Retrieval of the Understory LAI

The understory LAI and overstory LAI were retrieved based on LiDAR equations [34] and gap theory [51,52]. The LiDAR equations were used to build a link between the LiDAR waveforms and the vertical gap probability distribution, while gap theory was used to quantify the relationship between the gap probability and the LAI.

The LiDAR equations define the energy of a waveform received from the forest environment as a function of the emitted pulse energy, the vertical distribution of the gap probability, and the leaf and ground backscattering coefficients, as shown in the following equations:

$$R_v\left(z\right)=J_0·{\rho }_v·\left(1-P\left(z\right)\right)$$

(6)

$$R_v\left(0\right)=J_0·{\rho }_v·\left(1-P\left(0\right)\right)$$

(7)

$$R_g=J_0·{\rho }_g·P\left(0\right)$$

(8)

where $R_v\left(z\right)$ is the integrated laser energy returned from the top of the canopy to a height $z$, with $R_v\left(0\right)$ denoting the integrated laser energy returned from the top of the canopy to the bottom of the vegetation layer; $P\left(z\right)$ represents the gap probability at height z, with $P\left(0\right)$ representing the gap probability at the ground; $J_0$ is the emitted pulse energy after atmospheric correction; ${\rho }_v$ is the backscattering coefficient of the vegetation, which is related to the leaf spectral properties, the leaf angular orientation, and the phase function of leaf scattering; $R_g$ is the laser energy returned from the ground; and ${\rho }_g$ is the backscattering coefficient of the ground.

From Equation (8), we obtain:

$$P\left(0\right)=\frac{R_g}{J_0\cdot {\rho }_g}$$

(9)

Substituting Equation (9) into Equation (7), we obtain:

$$R_v\left(0\right)=J_0\cdot {\rho }_v\cdot \left(1-\frac{R_g}{J_0\cdot {\rho }_g}\right)$$

(10)

$$J_0=\frac{R_v\left(0\right)}{{\rho }_v}+\frac{R_g}{{\rho }_g}$$

(11)

Substituting Equation (11) into Equation (6), we obtain:

$$P\left(z\right)=1-\frac{R_v\left(z\right)}{R_v\left(0\right)}·\frac{1}{1+\frac{{\rho }_v}{{\rho }_g}\cdot \frac{R_g}{R_v\left(0\right)}}$$

(12)

From the gap probability distribution in the vertical direction, the LAI can be calculated using the method of Macarthur and Horn [15]. The LAI can be determined as the integrated FAVD value between the desired height intervals if it is not too difficult to obtain the continuous FAVD between these intervals. However, this process can be replaced with the logarithmic transformation of the gap probabilities, which is easy to calculate as follows:

$$LA{I}_{cum}\left(z\right)={\int }_{z0}^z{F}_a\left(z\right)dz=-\frac{1}{G}·{\int }_{z0}^z\frac{dlnP\left(z\right)}{dz}dz$$

(13)

where $LA{I}_{cum}\left(z\right)$ is the cumulative LAI as a function of the height $z$; $z_0$ is the bottom height; $F_a\left(z\right)$ is the FAVD, in units of m²/m³; and $G$ is the projection coefficient, which was set to 0.5 here under the assumption of a random foliage distribution.

Equation (12) does not consider the backscattering coefficient differences between the overstory and understory vegetation, and from Equation (13), we cannot separate the understory LAI from the total LAI.

To calculate the understory LAI, we introduced the understory leaf backscattering coefficient and emphasized $P\left(b\right)$, which is the gap probability at the boundary height between the overstory and understory layers, as shown in the following three equations:

$$R_o=J_0·{\rho }_o·\left(1-P\left(b\right)\right)$$

(14)

$$R_u=J_0·{\rho }_u·P\left(b\right)·(1-P\left(0{)}^{\prime }\right)$$

(15)

$$R_g=J_0·{\rho }_g·P\left(b\right)·P{(0)}^{\prime }$$

(16)

Combining Equations (14)–(16), we obtain:

$$P{(0)}^{\prime }=\frac{1}{\frac{R_u}{R_g}\cdot \frac{{\rho }_g}{{\rho }_u}+1}$$

(17)

$$P\left(b\right)=\frac{1}{\frac{R_o}{R_g}\cdot \frac{{\rho }_g}{{\rho }_o}\cdot P{(0)}^{\prime }+1}$$

(18)

$$P\left(0\right)=P{(0)}^{\prime }·P\left(b\right)$$

(19)

$R_o$, $R_u,$ and $R_g$ are the laser energy returns from the overstory, understory, and ground, respectively. Their ratios, $R_u/R_g$ and $R_o/R_g$, can be easily calculated from the decomposed waveforms. $J_0$ is the atmospherically corrected laser pulse energy, as previously defined. $P{\left(0\right)}^{\prime }$ denotes the gap probability that a laser pulse can penetrate from the boundary height to the ground without hitting vegetation. $P\left(0\right)$ denotes the gap probability that a pulse can penetrate from the canopy top to the ground without hitting vegetation. ${\rho }_o$, ${\rho }_u,$ and ${\rho }_g$ are the backscattering coefficients of the overstory, understory, and ground, respectively. Here, we treated these coefficients in a simplified way, as has been done in previous studies [14,35,53], by replacing them with the reflectance values of the overstory leaves, understory leaves, and ground.

The understory LAI is the integrated value between the boundary height and the ground height:

$$ULAI=-\frac{1}{G}·{\int }_0^b\frac{dlnP\left(z\right)}{dz}dz=-\frac{1}{G}\left(ln(P\left(b\right)\right)-ln(P\left(0\right)))$$

(20)

3.4. Sensitivity Analysis of the Spectral Parameters

Leaf and soil spectral properties of the same species often differ in different periods, and they also vary even within a small area [54,55,56]. In addition, leaf spectral properties are affected by the vertical canopy position [57]. In our study, two spectral parameters, ${\rho }_g/{\rho }_o$ and ${\rho }_g/{\rho }_u$, were simply treated as the mean values of the measured reflectance, which may cause large uncertainties in the quantification of the understory LAI. Thus, the influence of these uncertainties on the modified LiDAR equation performance should be tested. On the basis of our field measurement, the mean values of ${\rho }_g/{\rho }_o$ and ${\rho }_g/{\rho }_u$ were set to 1.48 and 1.76, respectively. Here we changed their corresponding values to range from 1 to 2 and 1.2 to 2.2, respectively, to evaluate their effects on the understory LAI estimation.

4. Results

4.1. Overstory–Understory Height Boundary

Two patterns of the gap probability distribution in vertical direction were observed, representing with obvious gap stratum and without obvious gap stratum between the overstory and understory vegetation, respectively (Figure 5). In both cases, the gap probability increased from the ground to the top of the canopy. For the pattern with an obvious gap stratum, the gap probability remained the same over a height stratum between 2.1 m and 5 m, except in the height interval near the top. For the pattern without an obvious height gap, the gap probability continue to increase until it reached the top height.

Figure 6 displays the distribution of $\mathrm{dp}/\mathrm{dz}$ along the vertical direction. For a plot with obvious gap stratum, we determined the height boundary as the height corresponding to the value of $\mathrm{dp}/\mathrm{dz}$ that was equal to zero in the height interval between 1 and 4 m. For a plot without an obvious gap stratum, we set the height boundary to the aforementioned 2 m.

Table 2. Height boundaries at the plot level (unit: m).

Plot ID	Height Boundary	Plot ID	Height Boundary
1	1.35	9	1.50
2	3.90	10	2.40
3	3.30	11	2.40
4	2.10	12	3.60
5	3.90	13	2.00
6	3.90	14	2.10
7	2.85	15	1.05
8	1.50	16	2.00

lists the height boundary determined for each plot. Six gap height values of the 16 plots fell between 1 and 2 m, while the remaining 10 plots had height gap values greater than 2 m. The highest boundary value was 3.90 m, belonging to plots 2, 5, and 6. In contrast, the lowest boundary value was 1.05 m, corresponding to plot 15. There was no obvious gap stratum between the overstory and understory vegetation at plots 13 and 16.

4.2. Waveforms after Deconvolution and Decomposition

A received waveform and its corresponding deconvolved and decomposed waveforms are displayed in Figure 7. The results show that after deconvolution, the intensities become higher, and the widths become narrower. Compared with those in the raw received waveform, the peaks in the deconvolved waveform are much easier to find. This verifies our assumption that deconvolution can improve the vertical resolution of the waveform and help to more effectively distinguish waveform components from different targets. After decomposition, the waveform components from the overstory, understory, and ground can be explicitly discriminated.

4.3. Sensitivity of the Method to Spectral Parameters

Figure 8 shows the estimated understory LAI values for a randomly selected plot. We can see that the predicted understory LAI is not sensitive to ${\rho }_g/{\rho }_o$ but that it increases with increasing ${\rho }_g/{\rho }_u$; however, the variation is slight.

4.4. Retrieved Understory LAI

The plot-level understory LAI and overstory LAI values were estimated by averaging the corresponding retrieved values at the footprint level. Figure 9 shows the relation between the estimated understory and overstory LAIs and the corresponding field-measured LAIs. We can see that the estimated understory LAI values were in good agreement with the measured understory LAI values ($R^2$= 0.54, RMSE = 0.21, bias = 0.02).

5. Discussion

Prior to this study, some understory structural parameters, e.g., height, cover, and density, were estimated using discrete return LiDAR data [26,41]. Our study demonstrates that the understory LAI can also be quantified through the combined use of discrete return and full-waveform LiDAR data.

The agreement between the estimated understory LAI and the measured LAI is very good, but still not ideal. A primary reason may be that the range of our measured understory LAI values is relatively small. With this limitation, a bias of 0.5 in a plot could lower the consistency between the predicted and measured values. Thus, experiments should be developed to acquire larger-range understory values to further validate the effectiveness of our method.

The use of inaccurate input parameters in the modified LiDAR equations could also impact the estimated results. In the modified LiDAR equations, there are three kinds of parameters. The first is the overstory–understory height boundary parameter. This parameter is crucial not only because of its key role in the LiDAR equations but also because of its critical influence on waveform component discrimination. In our study, we set the plot size used to determine the height boundary equal to the size of our investigated plots (25 m). However, this value may not be appropriate since larger plots are more likely to be heterogeneous. Smaller, more homogeneous plots are more suitable for use in determining the boundary height. Smaller plot sizes, such as 10 m and 5 m, could have been tested to determine a more accurate height boundary if the accuracy of our positioning software had been sufficient. The second parameter type corresponds to the waveform-related parameters, which were calculated in accordance with the full-waveform data. The third parameter type corresponds to the spectral parameters, which are subject to uncertainty. A sensitivity analysis of the spectral parameters showed that only the uncertainties of ${\rho }_g/{\rho }_u$ slightly affected the predicted understory LAI, which indicated the robustness of our method in quantifying understory LAI.

The good agreement between the estimated and measured values is encouraging, but there are still some limitations in this paper. First, the understory LAI estimation was based on the extraction of the waveform components associated with the understory vegetation, which may not be effective for certain understory vegetation types. For example, it seems impossible to extract waveform returns from moss and low herbs. There are two main factors contributing to this. One is that understory vegetation is rooted in the soil. Unlike in the overstory–understory spatial distribution pattern, there is no gap stratum between the understory vegetation and the soil, which leads to unavoidable mixing of the understory and soil waveform components. The other factor is that the heights of these understory species are often lower than the vertical resolution of the LiDAR system. This increases the degree of mixing of the understory and soil waveform components, making it more difficult to separate them. This limitation comes from the height-based distinction strategy. In fact, except for height, the intensity of discrete LiDAR data—especially for the single, first, and last return LiDAR data—is also a good indicator for the distinction of forest components [27]. In the future, LiDAR intensity should be considered to advance the research. Second, our method was developed for application to a temperate artificial forest, and it may not be applicable to more complex forests. Both the overstory and understory vegetation types are relatively simple in this man-made forest, causing the spectral parameters in the modified LiDAR equations to have relatively low uncertainty. Therefore, theoretically, the uncertainty of the estimated understory LAI values is relatively low. For more complicated forests, such as tropical forests, the uncertainty of these parameters will be larger, which will lead to higher uncertainty in the estimated results. As analyzed in the previous paragraph, ${\rho }_g/{\rho }_o$ and ${\rho }_g/{\rho }_u$ can influence the estimated LAI values. Thus, our method should be tested on more complex forests.

6. Conclusions

In this study, we proposed a new method for quantifying understory LAI in a temperate forest, which combined the advantages of both point cloud and full-waveform LiDAR data. The results of this study demonstrate that the gap probability distribution derived from point cloud data is a good indicator of overstory and understory height boundaries, and therefore, it is useful in overstory and understory waveform component discrimination. Waveform deconvolution can restore the vertical resolution of the full-waveform data and thus make it easier to find waveform components belonging to understory vegetation. Our modified LiDAR equations provide a sensible solution for estimating the understory LAI on the basis of the discriminated waveform components. This study indicates the feasibility of our method for quantifying the understory LAI in temperate forests. Further studies should be developed to test its effectiveness in complex forests, where the understory is composed of many types of vegetation, such as understory tree, shrub, herb, and moss. In addition, the next study should check whether this method is suitable for spaceborne LiDAR data.