<figure><img src="image_1.png" /></figure>

Article

# Modeling Small-Footprint Airborne Lidar-Derived Estimates of Gap Probability and Leaf Area Index

Tiangang Yin<sup>1,2,*</sup> (✉), Jianbo Qi<sup>3,4</sup>, Bruce D. Cook<sup>2</sup>, Douglas C. Morton<sup>2</sup>, Shanshan Wei<sup>5</sup> and Jean-Philippe Gastellu-Etchegorry<sup>4</sup>

Received: 14 November 2019; Accepted: 16 December 2019; Published: 18 December 2019

<figure><img src="image_2.png" /></figure>

check for
updates

**Abstract:** Airborne lidar point clouds of vegetation capture the 3-D distribution of its scattering elements, including leaves, branches, and ground features. Assessing the contribution from vegetation to the lidar point clouds requires an understanding of the physical interactions between the emitted laser pulses and their targets. Most of the current methods to estimate the gap probability ($P_{gap}$) or leaf area index (LAI) from small-footprint airborne laser scan (ALS) point clouds rely on either point-number-based (PNB) or intensity-based (IB) approaches, with additional empirical correlations with field measurements. However, site-specific parameterizations can limit the application of certain methods to other landscapes. The universality evaluation of these methods requires a physically based radiative transfer model that accounts for various lidar instrument specifications and environmental conditions. We conducted an extensive study to compare these approaches for various 3-D forest scenes using a point-cloud simulator developed for the latest version of the discrete anisotropic radiative transfer (DART) model. We investigated a range of variables for possible lidar point intensity, including radiometric quantities derived from Gaussian Decomposition (GD), such as the peak amplitude, standard deviation, integral of Gaussian profiles, and reflectance. The results disclosed that the PNB methods fail to capture the exact $P_{gap}$ as footprint size increases. By contrast, we verified that physical methods using lidar point intensity defined by either the distance-weighted integral of Gaussian profiles or reflectance can estimate $P_{gap}$ and LAI with higher accuracy and reliability. Additionally, the removal of certain additional empirical correlation coefficients is feasible. Routine use of small-footprint point-cloud radiometric measures to estimate $P_{gap}$ and the LAI potentially confirms a departure from previous empirical studies, but this depends on additional parameters from lidar instrument vendors.

**Keywords:** radiative transfer model; Lidar; airborne laser scan; point cloud; reflectance; leaf area index; gap probability; clumping; Gaussian decomposition; waveform

## 1. Introduction

Lidar remote sensing encompasses a broad range of technologies and applications [1]. Most remote sensing lidar devices use the time-of-flight technique to generate precise range measurements

<footnote n="1">¹ Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD 20740-3823, USA</footnote>

<footnote n="2">² NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA; bruce.cook@nasa.gov (B.D.C.); douglas.morton@nasa.gov (D.C.M.)</footnote>

<footnote n="3">³ Key Laboratory for Silviculture and Conservation of Ministry of Education, Beijing Forestry University, Beijing 100083, China; jianboqi@bjfu.edu.cn</footnote>

<footnote n="4">⁴ CESBIO - UPS, CNES, CNRS, IRD, Université de Toulouse, CEDEX 9, 31401 Toulouse, France; jean-philippe.gastellu-etchegorry@cesbio.cnes.fr</footnote>

<footnote n="5">⁵ CENSAM, Singapore-MIT Alliance for Research and Technology, Singapore S138602, Singapore; shanshan@smart.mit.edu</footnote>

<footnote>* Correspondence: tiangang.yin.85@gmail.com or tiangang.yin@nasa.gov</footnote>
based on the reflected signals of outgoing laser pulses. Lidar systems are typically categorized according to either the platform, with its corresponding laser footprint diameter (i.e., large-footprint: >25 m [2,3]; mid-footprint: 3–25 m [4,5]; small-footprint: 10 cm–3 m [6–8]; terrestrial laser scan (TLS): 0.1–3 cm [9]), or the detector system for recording return energy (e.g., waveform, discrete returns, single-photon detection, etc.). Among the current lidar technologies, waveform lidar systems store the most comprehensive information, and the recorded waveform can be converted into discrete points by deconvolution, thresholding, zero-crossings, or Gaussian Decomposition (GD) [10].

In recent years, full-waveform lidars have been explored in depth to estimate gap fraction, leaf area index (LAI) profiles, and biomass [mid-to-large footprint [11–15], small-to-mid footprint [16,17], and TLS [18,19]. Most of these physically based approaches are capable of accurate estimation of gap probability ($P_{gap}$) and effective LAI ($eLAI = \omega \cdot LAI$, where $\omega$ is the clumping index [20]) without calibration with field measurements. However, conversion from $P_{gap}$ into true LAI in most of these approaches requires a generalization of $\omega$ which was derived from globally or regionally averaged satellite products [21,22] instead of from the lidar data itself. Indeed, $\omega$ is a complex parameter that accounts for clumping induced by landscape (e.g., spatial distribution of the canopy), canopy (e.g., within-crown leaf distribution), and shoots [20]. Clumping of shoots is beyond the scope of this study and therefore not studied. For small-footprint airborne laser scan (ALS), Hu, et al. [23] computed the geometrical path length distribution within the crown to consider the clumping induced by crown shape and landscape. However, within-crown leaf area density variation was not completely addressed.

Despite the comprehensive information contained in a lidar waveform, storage of lidar data as discrete points provides an efficient alternative. Discretization of the waveform can generate uncertainties [24], and it is not necessarily straightforward to attribute precise radiometric information to points, as “intensity” values are typically assigned by proprietary onboard processing systems that usually present a black box to lidar users. Understanding the algorithms that convert waveforms into points helps to avoid misuse of ambiguous definitions of radiometric quantities that can give erroneous estimations that violate physical principles. Information and accuracy losses during waveform discretization can also undermine the use of lidar points for biophysical parameter retrieval [25,26]. Given the influence of ambiguous coefficients and residual radiometric issues, there is considerable controversy over the point-cloud inversion methods used to estimate $P_{gap}$ from a computed laser penetration index (LPI), and in a further step the effective and the true Plant/Leaf Area Index (PAI/LAI) of vegetation [27]. For example, the accuracy of TLS inversions is affected by partial hits that depend on the dimensions of the laser beam and leaves [18,28,29]. Moreover, for ALS, the point density within an area or volume is usually used to estimate the LAI [30–34] or the leaf area density [35]. However, as mentioned by Armston, et al. [16] and Chen, et al. [17], methods based on point density provide estimations of $P_{gap}$ that rely on additional empirical correlation coefficients with field measurements. Ambiguous “intensity” information from lidar points has also been used to retrieve the LAI [36]. Most point-cloud inversion approaches lack a theoretical underpinning that would enable universal application under various instrument specifications and environmental conditions. The direct computation of LPI and $\ln(LPI^{-1})$ (canopy interception described by the Beer–Lambert law by considering LPI as $P_{gap}$) from these approaches differ markedly. Furthermore, most previous empirical studies focused mainly on individual test sites (restricted environmental conditions) with a specific lidar device (fixed instrument) and a few airborne flights (limited experimental configurations). A general lack of rigorous sensitivity studies could raise questions regarding the accuracy and limitations of each approach: e.g., can any of the current methods be applied universally?

The methods developed for full-waveform data could be applied to point data if the radiometric information could be retrieved and understood for each point. The GD technique was developed to describe the waveform as a combination of discrete Gaussian profiles defined by peak amplitude, time shift, and standard deviation [37]. These quantities can be stored in point data with or without storage of the waveform [38–43]. Although these idealized quantities and the full waveform radiometric calibration [44,45] have been discussed previously, the physical interpretations under realistic conditions
have not been completely explored, such as the influences of footprint size, target size and orientation, terrain, clumped leaves, etc.

A robust radiometric linkage is required for various quantities that are involved in point-data generation and processing to understand the lidar point data and the inversion methods. Advances in lidar systems have been mirrored by development of lidar modules in radiative transfer models (RTM) [46–49], which can accurately simulate laser-target interactions under various instrument specifications (e.g., beam width, divergence and acquisition rate), acquisition and environmental conditions (e.g., flight altitude, viewing angle range, and terrain slope), and vegetation structures (e.g., leaf reflectance, leaf angle distributions, and canopy closure). The 3-D RTM framework is suitable for evaluating the sensitivity of different inversion approaches. It is capable of simulating laser pulse interactions with realistic 3-D scenes and recording the energy profile of the waveform. In recent years, many lidar waveform RTMs have been developed by extending the existing credible models, including FLIGHT [47], DIRSIG [50,51], RAYTRAN [52,53], LIBRAT [54], FLIES [55], and DART [48,56]. These RTMs have shown their capabilities in modeling lidar waveforms in the Radiation transfer Model Intercomparison (RAMI) project [57]. The DART waveform lidar module has been extended to adapt multi-pulse mode [58], which can efficiently simulate ALS and TLS waveform data in industrial lidar format, with recent bounding volume hierarchy (BVH) implementation [59] for ray tracing acceleration. Also, DART has the capability for solar signal and photon counting simulation [58]. These make DART a suitable tool for the sensitivity study and evaluation of different approaches.

This work takes advantage of the latest implementation in DART in physical modeling of point clouds from the existing multi-pulse waveform module, to conduct an in-depth study and intercomparisons of the existing gap probability and leaf area index inversion approaches of small-footprint ALS (Table 2). In Section 2, we have reviewed and built the physical basis and the modeling approaches of various radiometric quantities associated with lidar point “intensity” for investigating various laser-target interactions in radiative transfer models. Additionally, we have described the implementation in the latest DART release to convert the simulated ALS and TLS waveform data into points. Building on the theoretical basis, in Section 3, we have conducted a comprehensive review with physical interpretation about the existing $P_{\text{gap}}$ and LAI estimation approaches from either point-number-based (PNB) methods or intensity-based (IB) methods. In Section 4, we have demonstrated the utility of the DART model in evaluating these approaches from sensitivity analyses with different footprint sizes, leaf areas, leaf density variations, foliar dimensions, and homogeneous/heterogeneous scene configurations. From the results of the sensitivity studies, we have discussed the optimization of LAI estimation through small-footprint ALS, which could leverage existing and future point datasets to develop more robust LAI map products.

## 2. Theoretical Background and Implementations in DART

The radiometric theories and point-cloud modeling in DART are presented in this section. Detailed scientific and technical aspects of the waveform lidar module are described by Gastellu-Etchegorry, et al. [48] and Yin, et al. [58].

### 2.1. Lidar Pulse

DART has improved lidar modeling in realistic acquisition geometry and power distribution. For an emitted laser pulse, two Gaussian profiles are defined: the temporal convolution of the transmitted pulse and receiver response function, $S(t) = \hat{S}e^{-\frac{t^2}{2s^2}}$; and the 2-D power profile $P_1(\beta)$ within the footprint cone such that the ratio of $P_1(\beta_t)$ at the boundary (half-divergence $\beta_t$) to the central maximum $\hat{P}_{1,\beta}$ can be 0.5 (i.e., full width at half maximum (FWHM)), $1/e^2$, or $1/e$ [60]. The distribution follows: $P_1(\beta) = \hat{P}_{1,\beta}e^{-\frac{\beta^2}{2s^2}}$, where $\beta$ represents the angular offset from the pulse direction and $s_\beta$ is the standard deviation of the angular divergence. The reception is defined with telescope area $A_t = \frac{\pi D_r^2}{4}$.
The power spread within the footprint area (radius $r_{fp} = R \cdot \tan \beta_t + \frac{d_1}{2}$) is precisely defined, where R is the sensor-to-target range and $d_1$ is the beam cross-section diameter at the “exit gate” of the laser generator. $d_1$ is negligible for aircraft and spaceborne platforms, but it has a critical influence on the TLS footprint dimension.

## 2.2. Lidar Point Cloud and the Corresponding Radiometric Quantities

GD is used to extract points from simulated waveform data. According to Wagner, et al. [37], in the presence of X echoes, the temporal power profile recorded by the receiver $P_r(t)$ can be approximated as:

$$P_r(t) = \sum_{i=1}^{X} \hat{P}_i e^{-\frac{(t-t_i)^2}{2s_{p,i}^2}} \quad (1)$$

For the i-th Gaussian profile, $t_i$ is the temporal centroid; $s_{p,i} = \sqrt{s_s^2 + s_i^2}$ where $s_i$ compensates for the broadening effect from an oblique surface or a cluster of leaves that cannot be distinguished in a single return, and $s_s$ is an instrument-specific constant unless the temporal profiles of every transmitted pulse are known [61]. $s_{p,i} \approx s_s$ if the footprint size is negligible compared to the pulse duration distance. The ability of GD to distinguish different targets relies on target size, surface angle, gap size, and acquisition rate. $\hat{P}_i$ is the time-gated peak amplitude of the waveform $\hat{P}_i = \frac{\sigma_i}{C_{cal} R_i^4 s_{p,i}}$, where $C_{cal}$ is a calibration coefficient with pre-defined instrumental and experimental configurations, and detector efficiency that could be related to distance.

Wagner, et al. [37] considered the differential backscatter cross-section $\sigma = \frac{4\pi}{\hat{\Omega}} \cdot \rho \cdot A_s$ of an ideal target with an effective receiving area $A_s$, a natural reflectivity $\rho$, and a solid angle $\hat{\Omega}$ in which all reflected fluxes are assumed to be uniformly distributed. However, for an actual multi-return lidar pulse with only the target distance known and no details of the interaction, partial hits can occur and both $\rho$ and $A_s$ are unknown. Therefore, for each return, we refine this expression as:

$$\sigma_i = \frac{4\pi}{\hat{\Omega}} \cdot \rho_{a,i} \cdot A_{s,i}, \quad (2)$$

where $A_{s,i} = \frac{\pi R_i^2 \beta_i^2}{4}$ is the footprint area at a distance $R_i$ ($d_1$ is neglected). Here, we define the apparent reflectance $\rho_{a,i}$ as the ratio of reflected radiant flux from a surface, over a perpendicular infinite Lambertian surface with a reflectivity of 1, which can be considered similar to the biconical reflectance factor with the same incident and view direction [62] associated with every return to a pulse. By contrast, from the isotropic intensity given by Wagner, et al. [37], we considered isotropic scattered radiance of a Lambertian surface to define $\rho_{a,i}$ (intensity follows Lambert's cosine law). Thus, $\hat{\Omega} = \int_{2\pi} \cos \theta d\Omega = \pi$; and

$$\sigma_i = \pi \rho_{a,i} R_i^2 \beta_i^2. \quad (3)$$

For each point output, the subjected results have five types of radiometric quantities (cf. Appendix A):

1. Rapid detected peak amplitude $\hat{P}_i$ by first-derivative zero-crossing (without GD).

2. Fitted peak amplitude $\hat{P}_i$ (after GD by non-linear least-squares minimization [63]).

3. Standard deviation $s_{p,i}$ (after GD).

4. Power integral $I_i$ of a return (after GD):

$$I_i = \sqrt{2\pi} \cdot \hat{P}_i \cdot s_{p,i}. \quad (4)$$
5. Apparent reflectance $\rho_{a,i}$ (after GD):

$$ \rho_{a,i} = \frac{I_i C_{cal} R_i^2}{\pi \sqrt{2\pi} \cdot \beta_t^2} = \frac{I'_i C_{cal}}{\pi \sqrt{2\pi} \cdot \beta_t^2}. \qquad (5) $$

where $I' = I \cdot R^2$ is the distance-weighted power integral. Indeed, $\rho_{a,i}$ is physically equivalent to the lidar backscattering coefficient, which was suggested for waveform radiometric calibration due to the normalization relative to the footprint area [44,64].

By contrast, the biconical reflectance factor or apparent reflectance of an ALS return can be defined from radiation transfer theory:

$$ \rho_{a,i} = \frac{\pi I_i}{\Delta\Omega_i P_t \eta'} \qquad (6) $$

where $\Delta\Omega_i = \frac{D_i^2\pi}{4R_i^2}$ is the solid angle of return $i$ to the receiver; $P_t = \sqrt{2\pi}\hat{s}_s$ is the total pulse power, and $\eta$ is the system and atmospheric transmission factor. For various lidar devices, digital numbers that are linked to one of $\hat{P}_i$, $I_i$, or $I'_i$ could be named “intensity” in the output. In contrast to $I_i$ and $I'_i$, $\hat{P}_i$ cannot be used to directly derive $\rho_{a,i}$ without knowing $s_{p,i}$ of the Gaussian profile (unless $s_s \gg s_i$), which necessitates the recording [50] or appropriate processing [44] of the complete waveform. Note that some discrete-return lidar devices provide “reflectance” output by calibrating “intensity” values with reference targets (reflectivity of 10%, 20%, ...) at various distances [39,40]. Other information may also be available from instrument manufacturers. For example, without recording the full waveform, the online waveform processing of Riegl devices provides additional information in the extra byte associated with each point [42], but the computational details remain incomplete to the users. Correspondences between the derived radiometric quantities and the extra byte information for Riegl V-line and Q-line series are summarized in Table 1.

<table><caption>**Table 1.** Radiometric quantities provided by the Riegl V-line and Q-line extra byte associated with each point.</caption><thead><tr><th></th><th>Amplitude (P̂<sub>i</sub>)</th><th>Reflectance (ρ<sub>a,i</sub>)</th><th>Pulse Width (2.355×s<sub>p,i</sub>)</th><th>Deviation</th><th>Full Waveform Storage</th></tr></thead><tbody><tr><td>V-line</td><td>×</td><td>×</td><td></td><td>×</td><td></td></tr><tr><td>Q-line</td><td>×</td><td>×</td><td>×</td><td></td><td>×</td></tr></tbody></table>

## 2.3. Radiative Transfer under Realistic Conditions

Under realistic conditions, $\rho_{a,i}$ is not equivalent to the natural reflectivity $\rho_i$ of a target. There are usually three cases of physical interactions that could occur (Figure 1):

**Case A:** single Lambertian target perpendicular to the pulse direction, i.e., $s_i = 0$ and $s_{p,i} = s_s$ is a constant. From Equation (4), $\hat{P}_i$ is proportional to both $I_i$ and $\rho_{a,i}$, and $\rho_{a,i} = \rho_i$.

**Case B:** single oblique Lambertian target, which is the most common case for TLS and small-footprint ALS. From Lambert's cosine law, $I_i = I_{i-perp} \cdot \cos \theta$, where $\theta$ is the angle between the pulse and the target surface normal and $I_{i-perp}$ is the power integral in Case A. Therefore, $\rho_{a,i} = \rho_i \cdot \cos \theta$ if there is only a single return, which is consistent with findings by Hancock, et al. [19]. Standard deviation $s_i$ is a complex quantity that depends on the projection of the 2-D energy profile $P_1(\beta)$ onto the oblique target [41]. Unless $s_s \gg s_i$, there is no direct link between $\hat{P}_i$ and $I_i$.

**Case C:** a cluster of targets (i.e., clumped leaf volume and within-crown gap size smaller than the lidar footprint). The waveform presents a decayed profile from the peak as the laser penetrates through the cluster. The returns from several targets can merge to a single-return profile. In an ideal case:

$$ I_i = P_t \eta \cdot (1 - P_{gap}) \cdot T(\Omega_i) \cdot \frac{\rho \cdot \Delta\Omega_i}{\pi}. \qquad (7) $$

<footnote>* available only for certain devices, ** may be derived from pulse deconvolution instead of Gaussian Decomposition (GD).</footnote>
Here, $P_{gap}$ is the gap probability of the canopy. $T(\Omega_i) \cdot \frac{\rho \cdot \Delta\Omega_i}{\pi}$ is the backscatter transfer function along the direction $\Omega_i$, with $T(\Omega_i)$ defined as:

$$T(\Omega_i) = \frac{\int_{2\pi} \frac{g(\Omega_f)}{2\pi} \cos^2 \theta_{f,i} d\Omega_f}{\int_{2\pi} \frac{g(\Omega_f)}{2\pi} \cos \theta_{f,i} d\Omega_f}, \qquad (8)$$

where $\frac{g(\Omega_f)}{2\pi}$ is the leaf angle distribution of the cluster, and $\theta_{f,i}$ is the angle between the leaf normal and $\Omega_i$. Combining Equation (7) with Equation (6) gives an expression for $\rho_{a,i}$ with three unknowns:

$$\rho_{a,i} = (1 - P_{gap}) \rho \cdot T(\Omega_i). \qquad (9)$$

Hence, the use of $\rho_{a,i}$ for the inversion of any of the unknowns requires valid hypotheses on the other two unknowns.

<figure><img src="image_3.png" /><figcaption>Figure 1. The influence of three types of physical interactions between the emitted pulse and target configuration on the shape of the received pulse.</figcaption></figure>

Case A was the preliminary assumption by Wagner et al. [37]. It usually serves as an experimental setup for calibration. Actual interactions to generate small-footprint lidar points are either Case B or Case C. Both indicate that natural reflectivity $\rho$ cannot be derived from $\rho_a$ without knowing the type of interaction (number of returns, terrain slope, footprint size, leaf and gap size, etc.). Note that Case C has been extensively explored in the study for mid-to-large full-waveform lidar research [11–14,16,17]. For those cases, the waveform provides more lidar metrics than can be generated from small-footprint point clouds.

## 3. LPI Estimation and LAI/PAI Inversion from Small-Footprint ALS Point Clouds

### 3.1. Estimation of $P_{gap}$ and LAI

#### 3.1.1. Theoretical Approaches

The gap probability over a canopy that depends on zenith angle $\theta$ can be approximated as $P_{gap}(\theta) = e^{-\omega \cdot G \cdot LAI \cdot \cos \theta}$. Although the angular effect is an important factor [65,66], here we follow the majority of the past ALS work, which assumed a vertical direction. Considering only the vertical direction:

$$P_{gap} = e^{-\omega \cdot G \cdot LAI} = e^{-G \cdot eLAI}, \qquad (10)$$
where $G = \int_{2\pi} \frac{g(\Omega_f)}{2\pi} \cos \theta_f d\Omega_f$ is the unit leaf area projection along the vertical direction. For lidar, $P_{gap}$ and eLAI are typically derived from the laser penetration index (LPI) of either a single larger-footprint pulse or multiple small-footprint pulses.

For small-footprint ALS, the 3-D structure provided by lidar signal can further distinguish different components of the clumping [23]. For example, the clumping of a heterogeneous scene can be dominated by spatial distribution of tree crowns, shape of crown, and within-crown clumping. If the proportion of crown vertical projection area is defined as the vegetation canopy cover fraction $f_{vcc}$ by neglecting the within-crown gaps, then $f_{vcc}$ can be used to separate the between-crown spaces. We can derive an alternative expression for $P_{gap}$:

$$P_{gap} = f_{vcc}e^{-G\omega_{in}LAI/f_{vcc}} + (1 - f_{vcc}), \quad (11)$$

where LAI is at landscape scale, and LAI/$f_{vcc}$ is the LAI under the vegetation coverage. Note that $\omega_{in}$ accounts for both the crown shape and the within-crown clumping. In Hu, et al. [23], $f_{vcc}$ was derived from lidar first returns by neglecting the lidar incident angle. Geometrical path length distribution was used to account for the crown shape contribution in $\omega_{in}$. Equation (11) provides a parsed format of Equation (10) specifically for small-footprint ALS, relying on the accuracy of direct estimations of $f_{vcc}$ and $\omega_{in}$ from ALS data.

### 3.1.2. Practical Empirical Correlation

For ALS, the percentage of laser light penetration through a canopy [laser penetration index (LPI)] provides information about the density of foliage [67]. In practice, most of the derived values (such as LPI or eLAI) are empirically correlated with the field measurements, such as the gap probability/transmittance (e.g., [30–32], Morsdorf, et al. [33], Korhonen et al. [34]), the eLAI obtained by LAI-2000/hemispherical photos, or the LAI measured from direct methods (e.g., harvesting or destruction, etc.). The effective LAI and tree LAI of estimation can be expressed as:

$$\text{eLAI}_{\text{est}} = (G)^{-1} \ln(\text{LPI}^{-1}), \quad \text{LAI}_{\text{est}} = (\omega G)^{-1} \ln(\text{LPI}^{-1}), \quad (12)$$

where G is either known or approximated. The relationships between the reference values and the estimated values can be described as:

$$P_{gap} = \text{LPI}^{\alpha}, \quad (13)$$

$$\text{LAI} = \alpha \cdot \text{LAI}_{\text{est}} = \alpha \cdot \omega^{-1} \cdot \text{eLAI}_{\text{est}} = \alpha(\omega G)^{-1} \ln(\text{LPI}^{-1}). \quad (14)$$

Note that $\alpha$ is not only the exponential coefficient to link LPI with $P_{gap}$, but also the linear coefficient to link the estimated effective LAI (eLAI$_{est}$) with eLAI. For the ideal case $\alpha = 1$, LPI, eLAI$_{est}$, and LAI$_{est}$ are exactly equal to $P_{gap}$, eLAI, and LAI, respectively, leaving one physical parameter $\omega$ as the slope to convert from eLAI$_{est}$ into the true LAI value. The statistical description used in Case C of Figure 1 [Equations (7)–(9)] for a single large-footprint pulse is also valid for the interpretation of integrated multiple small-footprint pulses within an area.

### 3.2. Review of LPI Computation

Explicitly, the LPI can be sampled by either the area fraction or power fraction of ground returns:

$$\text{LPI} = \frac{P_{t,g}}{P_{t,g} + P_{t,v}} = \frac{A_{g,\text{proj}}}{A_{g,\text{proj}} + A_{v,\text{proj}}} = \frac{\sum I_{g}'}{\sum I_{g} + \gamma \cdot \sum I_{v}'} \quad (15)$$

where subscript indices v and g relate to vegetation and ground, respectively; $P_t$ is the incident lidar power on the targets; $A_{proj}$ is the projected area; $\gamma$ is the ratio at which to convert vegetation return reflectance to ground return reflectance as if the pulse is intercepted by the ground. Distance R, apparent reflectance $\rho_a$, and $\gamma$ are all spatially aggregated.
According to the method of sampling area or intensity as illustrated in Equation (15), the approaches for estimating LPI through ALS are categorized as either point-number-based (PNB) or intensity-based (IB) in Table 2.

<table><caption>**Table 2.** Methods for laser penetration index (LPI) estimation using airborne laser scan (ALS) points and references.</caption><thead><tr><th></th><th>Representation</th><th>Expressions</th><th>Reference Examples</th></tr></thead><tbody><tr><td rowspan="6">Point Number Based (PNB) Methods</td><td>LPI<sub>all</sub></td><td>N<sub>g</sub>/N<sub>g</sub>+N<sub>v</sub></td><td>Luo et al. [36]<br>Hu et al. [23]</td></tr><tr><td>LPI<sub>weighted</sub></td><td>(N<sub>g,1</sub>+1/2N<sub>g,2</sub>+1/4N<sub>g,3</sub>+1/4N<sub>g,4</sub>)/(N<sub>t,1</sub>+1/2N<sub>t,2</sub>+1/3N<sub>t,3</sub>+1/4N<sub>t,4</sub>)</td><td>Fleck et al. [68]<br>Schneider et al. [35]<br>Grau et al. [28]</td></tr><tr><td>LPI<sub>first</sub></td><td>(N<sub>g,single</sub>+N<sub>g,first</sub>)/(N<sub>single</sub>+N<sub>first</sub>)</td><td>Solberg et al. [30–32]<br>Riaño et al. [69]</td></tr><tr><td>LPI<sub>last</sub></td><td>(N<sub>g,single</sub>+N<sub>g last</sub>)/(N<sub>single</sub>+N<sub>last</sub>)</td><td>Lovell et al. [70]<br>Morsdorf et al. [33]</td></tr><tr><td>LPI<sub>both</sub></td><td>(N<sub>g,single</sub>+0.5(N<sub>g,first</sub>+N<sub>g last</sub>))/(N<sub>single</sub>+0.5(N<sub>first</sub>+N<sub>last</sub>))</td><td>Korhonen et al. [34]<br>Cook et al. [71]</td></tr><tr><td rowspan="3">"Intensity" Based (IB) Methods</td><td>LPI<sub>γ</sub></td><td>Σ I<sub>g</sub>' / (Σ I<sub>g</sub> + γ·Σ I<sub>v</sub>)</td><td>Empirical constant<br>γ<br>Statistically derived γ : LPI<sub>fitted</sub><br>Lefsky et al. [72]<br>Luo et al. [36]<br>Armston et al. [16]<br>Chen et al. [17]</td></tr><tr><td>LPI<sub>nearest</sub></td><td>Σ I<sub>g</sub>'<sub>pure</sub> + Σ I<sub>g,v</sub>' / (Σ I<sub>g</sub>'<sub>pure</sub> + Σ I<sub>g,v</sub>'<sub>nearest</sub> + Σ I<sub>v</sub>'<sub>nearest</sub>)</td><td>Milenković, et al. [64]</td></tr></tbody></table>

### 3.2.1. Point-Number-Based (PNB) Methods for LPI Computation

As shown in Table 2, given a small footprint size such that only a single return is retrieved from each pulse, and a high pulse density, the number of returns from leaves ($N_v$) and the ground ($N_g$) are statistical samplings of $A_{v,proj}$ and $A_{g,proj}$ in LPI<sub>all</sub> estimates. Monitoring larger regions with ALS requires a higher altitude, which makes the footprint size larger and increases the partial hit chance. LPI<sub>weighted</sub> balances the weight of each return from a pulse, where $N_{t,i} = N_{g,i} + N_{v,i}$ is the total number of returns per pulse with $i = 1, 2, 3...$ returns ($N_{t,i} = i$), and $N_{g,i}$ is the corresponding number of ground returns. In older systems, limited to 2 or 4 returns per pulse, LPI<sub>first</sub> and LPI<sub>last</sub> were used, where $N_{single}$ counts the single returns from the ground ($N_{g,single}$) or from indistinguishable vegetation elements ($N_{v,single}$). LPI<sub>first</sub> tends to provide information on between-crown spaces, whereas LPI<sub>last</sub> contains information on both the between-crown spaces and the within-crown gaps. Usually, LPI<sub>first</sub> underestimates and LPI<sub>last</sub> overestimates $P_{gap}$. LPI<sub>both</sub> balances this situation empirically, and has proven to be better sensitive to general gap sizes.

### 3.2.2. "Intensity"-Based (IB) Methods for LPI Computation

In theory, radiometric information is more suited than point density for assessing vegetation properties (e.g., LAI) [73]. Depending on the specifications of a certain device and on the user's choice, the "intensity" could be related to $\hat{P}$, $I$, $I'$, $s_{p,i}$, or $\rho_a$. However, in practice, IB methods have been used infrequently, partly due to the ambiguous radiometric quantities and the uncertain accuracies. The expression for intensity-based LPI (dependent on $\gamma$) is given in Table 2 and Equation (15). For a pulse fully intercepted by the ground, if the variation in the elevation slopes within the area is negligible and the surface is Lambertian, $\rho_{a,g} = \rho_g \cos \theta_g$, where $\theta_g$ is averaged over the reconstructed terrain model [74,75]. For a pulse fully intercepted by the vegetation, by setting $P_{gap} = 0$ in Equation (9), we get $\rho_{a,v} = \rho_v \cdot T(\Omega_i)$. An expression for $\gamma$ can be derived as:
$$ \gamma = \frac{\rho_g \cdot \cos \theta_g}{\rho_v \cdot T(\Omega_i)}. \qquad (16) $$

In previous studies, $\gamma$ was usually assumed to be a constant. For example, $\gamma = 0.5$ was an empirical value given by Lefsky, et al. [72] for lidar devices with 1064 nm wavelength, but adopting the same value at 1550 nm can be inaccurate [36]. For typical ALS assumptions on horizontal terrain, vertical pulse incidence, and spherical leaf angle distribution, $T(\Omega_i) = 2/3$, and:

$$ \gamma_{\text{ideal}} = \frac{3 \cdot \rho_g}{2 \cdot \rho_v}, \qquad (17) $$

where $\rho_g$ and $\rho_v$ can be retrieved from ground measurements as prior values. Equation (17) is an ideal expression for any wavelength with measured $\rho_g/\rho_v$. Note that in Table 2, $I'$ is used instead of $I$ for IB methods. Additionally, $\rho_a$ can replace $I'$ in Table 2 according to Equation (5), and both of them can represent the backscattering coefficient.

Without any assumptions, the four unknowns in Equation (16) can vary across different areas. Thus, $\gamma$ should not be universal. Indeed, the denominator of LPI$_\gamma$ (Table 2) represents the cumulative $I'_g$ of the ground as if there is no vegetation present. For full-waveform lidar, $\gamma$ can be estimated with linear correlation using vegetation and ground backscattering coefficients [16,17]. $\gamma$ was found to be relatively consistent with varying mid-to-large footprint size, but the correlation decreases with decreasing footprint size. In practice, the least-squares correlation is also strongly influenced by the diverse variances of $\sum I'_v$ and $\sum I'_g$.

Milenković, et al. [64] proposed another method, which utilizes the nearest pure-ground pulse (pulse that generates a single ground return). All pulses and return intensities were classified into three types: pure-ground (intensity $I'_g^{\text{pure}}$), pure-vegetation, and vegetation-ground (with ground return intensity $I'_g^{\text{v}}$). The returns of the second and third types were associated with the nearest first types (defined as $I'_v^{\text{nearest}}$ and $I'_g^{\text{nearest}}$). LPI$_{\text{nearest}}$ in Table 2 relies on the pure-ground pulses, and therefore requires large gaps to increase the probability of pure-ground pulses. Another critical assumption is that $\rho_{a,g}$ for both the pure-vegetation and vegetation-ground pulses are approximately equal to $\rho_{a,g}$ of the nearest pure-ground pulses. However, considering the $\rho_{a,g}$ calculation, $\rho_g$ could vary significantly due to understory conditions and woody debris, and $\theta_g$ could be influenced by insufficient sampling, whereas more samples (including ground returns of vegetation-ground pulses) can be used for calculating $\theta_{g,i}$ (from terrain model reconstruction) in LPI$_{\text{IB},\gamma}$ estimation.

# 4. Comparative Sensitivity Study of LPI/LAI Estimation Approaches Using DART Simulations

## 4.1. DART Simulations

The simulated point clouds of various critical parameters are studied and inter-compared with the LPI/LAI estimation approaches in Section 3. A DART scene takes 3-D cells to contain the Earth's elements for either turbid medium or facets [76]. We used homogeneous and heterogeneous forest scenes constructed from leaf facets to precisely model interactions with lidar pulses. Based on this approach, DART can directly calculate reference parameters for evaluating the inversion results.

The simulated lidar mimics the characteristics of Riegl VQ-480i (Table 3), a multi-return online waveform-processing lidar that is mounted on the G-LiHT platform [8]. The ground position of each pulse was derived from flight height, moving speed, pulse repetition rate, and scanning speed. The incident direction of each pulse was set to be vertical ($\theta_{g,i} = 0$). We focused more on the physical interactions themselves, following the majority of previous work, which neglected the angular effects, although it is an important factor for the gap fraction [65,66]. The sensor height $R$ varied from 50 m to 1000 m, increasing the footprint diameter from 0.015 m to 0.30 m. At 1550 nm, natural reflectivities $\rho_v$ and $\rho_g$ were set to 0.243 (general leaf) and 0.340 (brown moss), respectively.
<table><caption>**Table 3.** Lidar parameters used in DART simulations (Rigel VQ-480i).</caption><thead><tr><td>Parameter</td><td>Value</td><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Wavelength</td><td>1550 nm</td><td>Scanning speed</td><td>100 lines/second</td></tr><tr><td>Laser sampling interval</td><td>1 ns</td><td>Laser pulse repetition rate</td><td>200 kHz</td></tr><tr><td>Laser beam divergence</td><td>0.3 mrad</td><td></td><td></td></tr></tbody></table>

### 4.1.1. Homogeneous Scenes

The test scene consisted of a square vegetation plot and a flat ground surface (Figure 2a). The plot had a height of 12 m with a 2 m space below the canopy to facilitate the unambiguous classification of ground returns and vegetation returns. Each leaf was represented by a square (5 cm × 5 cm), and leaves within the canopy volume (22 m × 22 m × 10 m) had a spherical angle distribution. Under such ideal configuration, $\gamma$ = 2.10 was computed from Equation (17). The scene LAI ranged from 1 to 6 for the sensitivity studies.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Simulated scenes. (a) Homogeneous vegetation plot; (b) Heterogeneous forest with randomly distributed trees; (c) Heterogeneous scene with within-crown clumping; (d) and Leaf area density histogram of the clumped scene.</figcaption></figure>

### 4.1.2. Heterogeneous Scenes

A heterogeneous DART scene representing a general forest (100 m × 100 m) with randomly distributed ellipsoidal and conical crowns (327 of each) was used to test the lidar-derived LPI and LAI (Figure 2b). The crown height and diameter were set to 12 and 3.6 m, respectively. Two simulations were conducted with different leaf area densities set to 0.25 and 0.5 for each crown, resulting in LAI = 1 and 2, respectively, for the whole scene. The lidar data were simulated with 0.30 m footprint diameter. The LPI was calculated for each 10 m × 10 m area. The reference $P_{gap}$ in this case was estimated with
the irradiance intercepted by the ground surface using a large number of nadir incident rays [77,78]. DART was able to generate the map of $P_{gap}$ for each area.

Within-crown leaf density variation was added to the heterogeneous scene of LAI = 2 (Figure 2c) to study the influences of clumping and structural complexity on LAI retrieval. In DART, the non-random leaf distribution within the crown was parameterized for each species by the proportion of full cells and the vertical distribution of the leaf parameters. The simulated scene included 16 species, each with a different extent of clumping, to ensure enough variation of within-crown gaps. For each species, three layers with different vertical distributions of leaf density comprised each crown, determining the proportion over the whole canopy. The vertical distribution of leaf density and the proportion of full cells were randomly chosen for each layer of the crown. The LAI of each species were also randomly determined, with a total LAI of 2 for the whole scene. In addition, a 25% probability of empty 1 m × 1 m voxels was added into the crowns to further complicate the crown structure and clumping. The histogram of leaf area density is shown in Figure 2d.

## 4.2. Results and Analyses

### 4.2.1. Homogeneous Scene

The estimations of LPI in each 2 m × 2 m area from both the PNB and IB methods from Table 2 are illustrated in Figure 3, with various LAIs (1–6) and footprint diameters (0.015–0.30 m). The reference $P_{gap}$ is indicated with a dashed line for each subfigure.

The PNB methods to compute $P_{gap}$ across the homogenous scene have low accuracy and varied performances on LAI and footprint diameter. Persistently, $LPI_{first}$ underestimates $P_{gap}$ and $LPI_{last}$ overestimates $P_{gap}$ across all simulations. LPIs computed using other PNB approaches are in between these two values: $LPI_{last} > LPI_{both} > LPI_{weighted}$ or $LPI_{all} > LPI_{first}$, which is consistent with the previous studies referenced in Table 2. As the LAI increases, most of the PNB methods switch from underestimation to overestimation. For a fixed LAI, estimations with infinitesimal footprint size converge to $P_{gap}$ for all PNB methods, which tends to verify the capability of TLS to capture $P_{gap}$ accurately [79]. The average number of returns for each pulse approaches 1 (single return) as the footprint size approaches 0 (Figure 4a). With increasing footprint size across all simulations, $LPI_{last}$, $LPI_{first}$, and $LPI_{both}$ converge to saturated values of 1, 0, and 0.5, respectively. As the LAI increases, the saturation speed reduces due to the lower probability of energy penetration through the canopy (Figure 4b). Although $LPI_{both}$ balances the underestimation of $LPI_{first}$ and the overestimation of $LPI_{last}$, the results do not indicate that $LPI_{both}$ is more accurate. By contrast, $LPI_{weighted}$ is closer to $P_{gap}$ than $LPI_{all}$. The saturation speeds of $LPI_{weighted}$ and $LPI_{all}$ are slower than those of the other approaches. The total number of returns are illustrated in Figure 4a. The average number of returns per pulse gradually becomes saturated together with $LPI_{weighted}$ and $LPI_{all}$, as shown in Figure 3. From the sensitivity study above, $LPI_{weighted}$ and $LPI_{all}$ have broader ranges of footprint size and LAI values than the other PNB methods.

Regarding the IB methods, the constant value of nearest pure-ground intensity $I_g^{pure}$ for $LPI_{nearest}$ estimation is over-idealized for the simulations. Indeed, $I_g^{pure}$ is a diverse variable for actual data due to soil moisture, litter, and terrain slope. For $LPI_{fitted}$, the least-squares approach applied to the dataset of $\sum \rho_{a,g}$ and $\sum \rho_{a,v}$ of every single pulse is used to determine the fitted $\gamma$ (Figure 5). The scatter plots for LAI = 0.5, 2.0, and 5.0 are shown in Figure 5a–c, respectively, with increasing footprint diameter (0.06 m, 0.15 m, and 0.30 m, respectively). The fitted line determines the slope value ($-\gamma$), the y-axis intersection ($\sum \rho_{a,g}$ with a reference value of 0.340), and the x-axis intersection [$\sum \rho_{a,g}$ with a reference value of 0.162 from Equations (16) and (17)]. It should be noted that with increasing footprint diameter, the fitted $\gamma$ value and the axial intersections converge to the reference values. The data variance decreases as correlation increases in Figure 5d,e. These results are consistent with [16] and [3], in which $\gamma$ should be estimated using a larger footprint dimension with smaller variance for LAI
retrieval (e.g., the spaceborne waveform lidars [2,80]). Note that the simulation might not represent
the realistic conditions because of the constant reflectivity of the ground surface.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** LPI estimates by different methods with various LAI values [1(**a**), 3(**b**), and 6(**c**) and footprint diameters (<0.30 m) for homogeneous simulations. The dashed lines indicate the reference $P_{gap}$ for each LAI, and the shaded regions represent the variances of LPI estimations.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Average discrete returns for all pulses for the homogeneous simulation: (a) average number of returns under different footprint sizes and (b) average number (probability) of ground returns.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Fitting $\gamma$ ratio with least-squares linear regression for LPI<sub>fitted</sub> computation for homogeneous simulations. Results are illustrated for footprint diameters (a) 0.06 m, (b) 0.15 m, and (c) 0.30 m. As footprint diameter varies, (d) shows the variation of fitted $\gamma$ and correlation coefficient R, and (e) shows the variation of intersections of $\Sigma \rho_{a,g}$ and $\Sigma \rho_{a,v}$.</figcaption></figure>

The three IB methods, LPI<sub>γ=2.10</sub>, LPI<sub>nearest</sub>, and LPI<sub>fitted</sub>, generally gave accurate estimations of P<sub>gap</sub> at all footprint sizes (α ≈ 1.0), except that LPI<sub>fitted</sub> underestimated P<sub>gap</sub> at small footprint sizes as expected. Due to the over-idealized simulation conditions, LPI<sub>nearest</sub> was suitable for all footprint sizes, regardless of the structural and optical properties of the vegetation. In general, IB methods are preferable for variation of LAI values. As the LAI increases, the ground returns become weaker than the vegetation returns, but pulses with larger footprint size can still produce ground returns. Thus, the average numbers of returns for different LAI values converge to the same value in Figure 4.

The empirical coefficient (α) in Equations (13) and (14) was studied for various estimation approaches with ω = 1 (homogeneous scene). Figures 6 and 7 illustrate the correlation between the references and the estimations using different methods. The results from the homogeneous simulation
are consistent with previous assertions mentioning that LPI derived from PNB methods requires $\alpha$ to correlate with field measurements of $P_{gap}$ [16,17]. For the small footprint diameter of 0.03 m, which mimics a TLS configuration, both exponential and linear correlations are successfully established for all the approaches. For the IB methods, $\alpha \approx 1.0$ is observed for LPI and the LAI was derived for $LPI_{y=2.10}$ and $LPI_{nearest}$. $LPI_{fitted}$ approached $P_{gap}$ as the footprint size increased. For PNB methods, $LPI_{weighted}$ and $LPI_{all}$ were more accurate than the other PNB methods, which give either underestimations ($\alpha_{last} = 3.38$, $\alpha_{both} = 1.67$) or overestimations ($\alpha_{first} = 0.37$). In Figure 7, the LAI correlations are almost linear for the PNB methods with LAI > 1, and converge to the origin with varying $\alpha$ as the LAI value approaches 0, which is caused by the boundary partial hit effect due to the non-negligible ratio of footprint size to leaf size. The results also indicate that $LPI_{weighted}$ is preferable among all the PNB methods. Regarding larger footprint diameters (0.15 and 0.30 m) that mimic regular ALS configurations, $\alpha \approx 1.0$ is still maintained for the IB methods. By contrast, all the PNB methods demonstrate low sensitivities for both $P_{gap}$ and LAI correlations. $LPI_{first}$ and $LPI_{last}$ are not seen from Figure 7 because $\ln(LPI^{-1})$ is consistently equal to infinity and 0, respectively. $LPI_{both}$ shows a negligible sensitivity as the LAI varies. Unexpected non-linear correlations can be observed for $LPI_{weighted}$ and $LPI_{all}$ with LAI < 3 (0.15 m footprint) and LAI < 1.5 (0.30 m footprint). Indeed, the merging of returns with high leaf density and large LAI (Case C of Figure 1) makes the linear correlation of PNB methods distorted and inappropriate for capturing the small gaps within the crown to give accurate estimations. Indeed, the merging of returns reduces the point number, but the reflectance would be cumulated into the merged point. Therefore, IB methods are not influenced by this effect.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Plots of reference $P_{gap}$ against lidar-derived LPI for footprint diameters of (a) 0.03 m, (b) 0.15 m, and (c) 0.30 m for homogeneous simulation.</figcaption></figure>
In a forest, foliar dimensions can vary a lot with tree species and growing stages. Even for the same LAI and footprint size, variation in foliar dimensions can influence the estimated LPI. Figure 8 shows the relative bias and sensitivity of lidar-derived LPIs (LAI fixed at 3) under different footprint sizes and foliar dimensions of the homogeneous scene. LPI<sub>weighted</sub> (the most reliable PNB method, Figure 8a and c and LPI<sub>γ=2.10</sub> (the idealized IB method, Figure 8b,d are studied here. LPI<sub>γ=2.10</sub> is much less sensitive to leaf dimensions than LPI<sub>weighted</sub>. We used the ratio of footprint size to leaf size to characterize the changes in accuracy. Figure 8c,d give the relative bias against the ratios for LPI<sub>weighted</sub> and LPI<sub>γ=2.10</sub>, respectively. LPI<sub>weighted</sub> approaches the references when the ratio approaches 1.0, which can also be confirmed by the 0% difference slope in Figure 8a and the plot of LAI = 3 in Figure 3. LPI and P<sub>gap</sub> are close at 0.05 m footprint diameter but diverse for larger footprint size and for LAI = 1 (underestimation) or LAI = 5 (overestimation) as illustrated in the other plots of Figure 3. This suggests that LPI<sub>weighted</sub> is sensitive to the footprint size and foliar dimensions when it is used to estimate small within-crown gaps. The bias plot of LPI<sub>γ=2.10</sub> converges to less than 3% when the ratio of footprint diameter to leaf length becomes larger than 6, which ensures that a lidar pulse covers an adequate number of foliar elements. Indeed, the relative error significantly increases with decreasing ratio value.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Plots of reference LAI against LAI<sub>est</sub> [estimated LAI in Equation (12) with ω = 1]. The slopes are equal to the α values for footprint diameters of (a) 0.03 m, (b) 0.15 m, and (c) 0.30 m of homogeneous scenes.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The bias of lidar-derived LPI against $P_{gap}$ of homogeneous simulations for various footprint sizes and leaf dimensions: (a) $LPI_{weighted}$, (b) $LPI_{y=2.10}$, and the relative error against the ratio of footprint diameter to leaf length, (c) $LPI_{weighted}$, and (d) $LPI_{y=2.10}$.</figcaption></figure>
### 4.2.2. Heterogeneous Scene

For the heterogeneous scenes, landscape-scale areas (10 m × 10 m) were studied. Figure 9 and Table 4 show the correlation between the lidar-derived LPIs and the reference values ($P_{gap}$ and LAI), with a relatively large footprint diameter (0.3 m). For the PNB methods, the correlation depends more on the large gaps between crowns due to the sensitivity loss for the within-crown gaps at 0.3 m footprint size, as illustrated previously in Figure 6. We estimated the reference value of $f_{vcc}$ by thresholding the generated $P_{gap}$ map of 0.1 m × 0.1 m pixel size into a binary image of 1 (open space) and 0 (canopy cover), and computing the fraction within each area. Table 4 illustrates the results of correlations of various PNB and IB methods in terms of four derived parameters: $\alpha$ is the fitted exponential coefficient of LPI against $P_{gap}$ (Figure 9a-c) from Equation (13); $\alpha \cdot \omega^{-1}$ is the fitted linear correlation coefficient of reference LAI against eLAI$_{est}$ described by Equations (14) and (12) (slopes of Figure 9d-f); computed $\omega$ is generated by dividing $\alpha$ by $\alpha \cdot \omega^{-1}$; and $R^2$ is the coefficient of determination of reference LAI against eLAI$_{est}$. The last column of Table 4 provides the reference values computed by using $P_{gap}$ as LPI. Additionally, derived results of both $f_{vcc}$ and $P_{gap}$ are shown in Figure 9 for comparison.

<table><caption>**Table 4.** Derived and computed coefficients for the PNB and IB methods from Figure 9 and Equations (13) and (12). The reference derivations from $P_{gap}$ are in bold.</caption><thead><tr><th rowspan="2">Variables</th><th colspan="6">PNB Methods</th><th colspan="2">IB Methods</th><th>References</th></tr><tr><th>LPI<sub>all</sub></th><th>LPI<sub>weighted</sub></th><th>LPI<sub>first</sub></th><th>LPI<sub>last</sub></th><th>LPI<sub>both</sub></th><th>LPI<sub>y=2.10</sub></th><th>LPI<sub>nearest</sub></th><th>LPI<sub>fitted</sub></th><th>P<sub>gap</sub></th></tr></thead><tbody><tr><td rowspan="4">LAI = 2<br>(Random)</td><td>α</td><td>0.40</td><td>0.97</td><td>0.60</td><td>27.14</td><td>1.56</td><td>1.00</td><td>1.02</td><td>1.00</td></tr><tr><td>α · ω<sup>-1</sup></td><td>1.57</td><td>1.84</td><td>1.16</td><td>22.08</td><td>3.90</td><td>1.98</td><td>2.01</td><td>1.98</td></tr><tr><td>computed</td><td>0.25</td><td>0.53</td><td>0.52</td><td>1.23</td><td>0.40</td><td>0.51</td><td>0.51</td><td>0.51</td></tr><tr><td>ω</td><td>0.88</td><td>0.84</td><td>0.77</td><td>0.19</td><td>0.79</td><td>0.89</td><td>0.89</td><td>0.90</td></tr><tr><td rowspan="4">LAI = 2<br>(Clumping)</td><td>α</td><td>0.42</td><td>0.97</td><td>0.67</td><td>24.20</td><td>1.62</td><td>1.01</td><td>1.03</td><td>1.00</td></tr><tr><td>α · ω<sup>-1</sup></td><td>1.50</td><td>2.20</td><td>1.42</td><td>21.11</td><td>4.2</td><td>2.24</td><td>2.28</td><td>2.23</td></tr><tr><td>computed</td><td>0.28</td><td>0.44</td><td>0.47</td><td>1.15</td><td>0.39</td><td>0.45</td><td>0.45</td><td>0.45</td></tr><tr><td>ω</td><td>0.71</td><td>0.80</td><td>0.80</td><td>0.34</td><td>0.84</td><td>0.89</td><td>0.89</td><td>0.90</td></tr><tr><td rowspan="4">LAI = 1</td><td>α</td><td>0.27</td><td>0.61</td><td>0.42</td><td>-</td><td>1.10</td><td>1.00</td><td>1.02</td><td>1.00</td></tr><tr><td>α · ω<sup>-1</sup></td><td>0.77</td><td>0.92</td><td>0.59</td><td>-</td><td>1.95</td><td>1.49</td><td>1.49</td><td>1.48</td></tr><tr><td>computed</td><td>0.35</td><td>0.66</td><td>0.71</td><td>-</td><td>0.56</td><td>0.67</td><td>0.68</td><td>0.68</td></tr><tr><td>ω</td><td>0.90</td><td>0.84</td><td>0.78</td><td>-</td><td>0.77</td><td>0.94</td><td>0.95</td><td>0.95</td></tr></tbody></table>

For the derived $\alpha$ of Table 4 from Figure 9a-c, all of the IB methods ($LPI_{nearest}$, $LPI_{y=2.10}$, and $LPI_{fitted}$) were consistent with $P_{gap}$ with $\alpha \approx 1$ for varying LAI and within-crown clumping (bias < 3%). The varying leaf area over the entire scene benefits the building of the slope $\gamma$ for $LPI_{fitted}$, as illustrated in Figure 10. The differences between fitted $\gamma$ were less than 0.5%, but the introduction of within-crown clumping reduced the fitted $\gamma$ by 5%, which suggests that within-crown clumping increases the variance of the dataset. The correlation coefficient remained high for all conditions. For the PNB methods, all LPI estimates are strongly biased towards the reference $P_{gap}$. $LPI_{first}$ severely underestimates $P_{gap}$ because of the small gaps within the crown. $LPI_{weighted}$ tends to give a more accurate estimation than $LPI_{first}$. When the LAI varies from 2 to 1, $LPI_{weighted}$ provides an underestimation with $\alpha$ varying from 0.97 to 0.61. $LPI_{last}$ provides a strongly saturated overestimation with low sensitivity for LAI = 2, and loses its sensitivity for LAI = 1. $LPI_{both}$ averages the inaccuracy of both $LPI_{first}$ and $LPI_{last}$. $LPI_{all}$ introduces a large offset. Additionally, $LPI_{first}$ slightly underestimates $f_{vcc}$ (red diamond scatter plot in Figure 9a,c) because of the multiple returns at the boundary of the crown, but it presents a feasible approach if the lidar incident angle is neglected.

In Figure 9d-f, according to Equation (12), the eLAI$_{est}$ from various approaches were linearly correlated with the reference LAI. Unlike in the homogeneous scene, $\omega$ is a non-negligible physical unknown. Therefore, eLAI$_{est}$ is used instead of LAI$_{est}$. The slope $\alpha \cdot \omega^{-1}$ and $R^2$ values are shown in Table 4. Compared with the reference values derived from $P_{gap}$, the three IB methods give results that are close to the references (bias < 3% and $R^2 > 0.89$). The estimations using all the returns ($LPI_{all}$, $LPI_{weighted}$, $LPI_{y=2.10}$, and $LPI_{nearest}$) are closer to the references than the estimations partially using
the points (LPIfirst, LPIlast, and LPIboth). LPIfirst was insensitive to the changes in the LAI and clumping. Due to the non-linearity of Equation (11), explicit analytical references of $\alpha \cdot \omega^{-1}$ do not exist. It should be noted that although the clumping-parsed expression of Equation (11) does not show a linear relationship, most linear regression studies show high $R^2$ values in this study range, which verifies the past work in linear fitting of field measurements for LAI map generation and insensitive correlation of LPIlast in estimating the large gaps between crowns.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Plots of $P_{gap}$ against lidar-derived LPIs (a–c) and reference LAIs against eLAIest (d–f) (estimated from Equation (12)) for a simulated scene with discrete tree crowns (LAI = 1 and 2) and within-crown clumping (LAI = 2). The slope ($\alpha \cdot \omega^{-1}$) and $R^2$ values are shown in Table 4.</figcaption></figure>
With derived $\alpha$ and $\alpha \cdot \omega^{-1}$ from the $P_{gap}$ and LAI correlation studies, the scene clumping index $\omega$ is computed by dividing these two values. The reference $\omega$ computed using $P_{gap}$ and input LAI ranges from 0.45 to 0.68 in the simulations. Among all of them, the addition of within-crown clumping reduces $\omega$ to between 0.51 and 0.45. For the IB methods, the computed $\omega$ is consistent with high accuracy (bias < 2%). Indeed, $\omega$ is a physical quantity that should not vary with estimation approach. In general, the PNB methods produce various computed $\omega$ values, which suggests that most of the PNB methods have influences not only on an empirical correlation coefficient (i.e., $\alpha$), but also on other physical quantities (i.e., $\omega$). Surprisingly, among all the PNB methods, $LPI_{weighted}$ gave a close estimation of $\omega$ with bias < 4%, although the fitted $\alpha$ could be far away from 1.0. Therefore, $LPI_{weighted}$ partially interprets ALS data physically, and it is preferable over all the other PNB methods.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Fitting $\gamma$ ratio and correlation coefficient with least-squares linear regression for $LPI_{fitted}$ computation. (a) LAI = 1 and (b) LAI = 2 for discrete tree crowns; (c) within-crown clumping LAI = 2.</figcaption></figure>

In addition to the experimental and environmental variables, instrumental specifications (e.g., the reflectance threshold $\rho_a$) can also influence the LPI estimation, and this is analyzed in Appendix B. The studies in Section 4.2 assumed that every infinitesimal peak from the waveform could be detected by the sensor.
## 5. Discussion

From a homogeneous layer of randomly scattered leaves, we found the point-based LPIs computed by small-footprint ALS either overestimate or underestimate $P_{gap}$, and they become saturated relying on two factors: (1) the ratio of footprint diameter to leaf or vegetation element dimensions and (2) the LAI. The former factor determines a pulse's capability in generating multiple returns, and the latter factor determines the number of returns per pulse. To apply PNB methods over a homogeneously dense forest, the footprint size should be controlled within a range so that the sensitivity is maintained as the LAI decreases and the LPI estimations do not converge to a constant value (e.g., 1, 0.5, 0, etc.). The estimation of this range requires a deep understanding of the relationship between leaf dimensions, leaf density, and device specifications. As for increasing footprint size, the estimated $P_{gap}$ from the PNB methods becomes insensitive to LAI or leaf area density variation as well as within-crown clumping, even with the path length variation considered [23,27]. Sensitivity studies demonstrated that IB methods for physically estimating $P_{gap}$ and eLAI were accurate and less influenced by variations in footprint size, leaf area, vegetation cover, and foliar dimensions than the PNB methods. From the modeling results, the main advantage of the IB methods is that the empirical coefficient $\alpha$ can be removed from Equations (13) and (14), leaving only the physical quantities (i.e., $\omega$) for LAI estimation if the threshold of $\rho_a$ is neglected (which could introduce bias as presented in Appendix B). They also have flexibility of use in a wide variety of instrumental, experimental, and environmental configurations, which cannot be achieved by PNB methods unless the footprint size is infinitesimally small (i.e., TLS configuration).

We have shown that some of the PNB methods are only sensitive to the large gaps, and the IB methods can capture the precise $P_{gap}$ of various gap sizes. Among the radiometric quantities derived from lidar, $I'_i$ and $\rho_a$ are the only suitable quantities for inversion. In practice, lidar point intensity has more unresolved uncertainties than point number. The ground reflectance variance can be large due to moisture, litter, and low vegetation and micro-topography. Also, a LAS file stores the intensity values as 8-bit or 16-bit digital numbers, and this discretization may introduce further precision issues. Additionally, most lidar devices capture $\hat{P}$ as the intensity value, which is unsuitable for the IB methods. For example, for Riegl lidar devices, most of the past full-waveform studies used Q-line with full-waveform storage from which all the radiometric quantities can be derived. Our study supports the potential of V-line using online waveforms processed into apparent reflectance $\rho_a$ [38–40] with lower cost, smaller divergence, higher pulse repetition frequency, and more efficient data storage without waveforms. For the IB methods, the constant $\gamma$ assumption in estimating $LPI_\lambda$ heavily relies on the statistically convergent leaf angle projection, clumping index, and natural reflectivity ratio of the leaves and the ground within a region. The understory vegetation and debris can influence the ground reflectivity and the fitting of $\gamma$. The full-waveform method used by Armston, et al. [16] and Chen, et al. [17] is a practical approach without empirical correlation. However, the method needs a different least-squares approach to consider the inconsistent variances over the dataset. The $LPI_{nearest}$ could be another practical approach, but possible sources of bias exist due to insufficient single-return ground pulses and the possible large variance of ground reflectance. Furthermore, the influences of the threshold of $\rho_a$ should be noted (Appendix B).

Based on both the modeling results in this study and on previous studies (e.g., [16,17,81]), estimation of LAI from small-footprint ALS point clouds without field measurements is promising. We investigated the dependence of the linear fitting coefficient $\alpha$. From both homogeneous and heterogeneous simulations, we have shown that $P_{gap}$ can be estimated without the correlation coefficient $\alpha$. For the heterogeneous scene, the linear coefficient $\alpha \cdot \omega^{-1}$ is analyzed, with a parsed expression of $\omega$ provided in Equation (11). eLAI is divergent from the actual LAI, for which precise estimation requires the inputs of vegetation cover fraction ($f_{vcc}$), shape, and within-crown clumping index ($\omega_{in}$) and G in Equation (11). Small-footprint ALS has a great advantage in capturing the structures of tree crowns. Although our simulation shows that $f_{vcc}$ can be estimated by $LPI_{first}$ based on the vertical pulse assumption, the actual lidar incident angle would introduce non-negligible biases. Considering the influences of both crown shape and within-crown clumping, $\omega_{in}$ could be resolved with
aggregating pulses within a small area to estimate the within-crown optical depth distribution instead of only the path geometrical length distribution [25]. Another difficult configuration is the footprint size, which should be larger than the within-crown gaps but smaller than the between-crown spaces, to capture: (1) the $f_{vcc}$ through first returns to exclude the within-crown gaps and (2) the light penetration within the crown through multi-return pulses to exclude the between-crown spaces. Fortunately, most of the recent ALS data are considered satisfactory (e.g., Cook, et al. [8]). Therefore, the LAI map product could be provided more accurately with a proper method. The ray-tracing-based voxel reconstruction (e.g., AMAPVox [81]) from $\rho_a$ could potentially be a solution in resolving $f_{vcc}$ and $\omega_{in}$ together. The precise estimation of G and $\omega_{in}$ might require a multi-ALS scan from different directions.

In this study, we took Riegl lidar configurations with extra byte storage as typical examples. However, all ALS data are diverse due to the sensor specifications and environmental configurations. The acquisition methods, configurations, and point-processing algorithms from different manufacturers are also diverse. The appropriate way to understand the capability of a device requires physical modeling instead of directly applying approaches that were previously verified with only a single experiment. There are several other important effects not included in this study, such as the angle distribution of incident pulses, the woody parts of realistic trees including twigs, as well as the sensor noise, dead-time, out-of-focus effect, etc. These effects will be considered in future work in which more complex approaches (e.g., ray tracing) are involved. In this work, we specifically focused on point representation through basic and idealized experimental modeling.

## 6. Conclusions and Outlook

In this work, we conducted an in-depth study of gap probability and LAI estimates from small-footprint ALS point clouds with sensitivity studies to evaluate the methods based on either point number (PNB) or intensity (IB), using improved DART point-cloud-modeling capabilities. The simulated scene configurations include both homogeneous vegetation plots and heterogeneous forest, with clumping effects existing at both the landscape and crown scale.

Apart from the clumping index with a clear physical definition, the PNB methods require additional empirical correlation coefficients that rely on field measurements. The correlation could become saturated, with diminished sensitivity, relying on the relationship between the ratio of footprint diameter to leaf dimension and the LAI. Among all the PNB methods, the LPI<sub>weighted</sub> approach [28,35,68] is preferable because it has the broadest effective range in adapting various configurations (footprint size, LAI, and clumping), and is the most physically based, which preserves the valid clumping index without the influence of correlation (Table 4). In contrast, all the IB methods listed in Table 2 have the great advantage over the PNB methods of being able to accurately estimate the gap probability without the requirement of a correlation coefficient [$a \approx 1$ in Equations (13) and (14)], given that an appropriate radiometric quantity of lidar point is used (either the distance-weighted power integral $I'$ or the apparent reflectance $\rho_a$ of each return). The IB methods remain accurate in the presence of both landscape and within-crown clumping. However, despite the empirical correlation coefficient, the clumping index $\omega$ that is required for true LAI estimation remains a challenge to estimate merely from lidar points. Further parsed analysis of clumping to figure out the separated contributions by landscape, crown shape, and within-crown variation requires new approaches to be developed, potentially by combining voxel-based ray tracing with well-defined radiometric quantity of ALS points.

For actual data processing, the choice of PNB methods or IB methods should be analyzed case by case depending on whether field measurements are available and whether the uncertainty or noise of using point radiometric quantity is less than those of using point number. The potential of using abundant ALS point-cloud resources could provide additional validations and pre-studies of existing and future spaceborne platforms. The point-cloud module has been implemented in the latest DART release, which is also capable of TLS simulation (versions higher than 5.7.0). The recent development of the DART lidar module supports the conversion of DART-simulated waveforms into point clouds, with storage of waveforms and point-cloud data in the standard LAS format with the extra byte output.
These advances bridge the gap between simulated data and current data-processing software, to improve the usage of lidar points for assessing the accuracy of ALS and TLS inversions. Our next work will be on evaluations in 3-D point voxelization software through DART ALS and TLS simulations: e.g., AMAPVox [81], VoxelSD [28], and VoxLAD [29,82]. DART provides free licenses for scientific and educational work.

**Author Contributions:** Conceptualization, T.Y. and S.W.; methodology, T.Y., J.-P.G.-E., J.Q. and B.D.C.; software, T.Y., J.Q. and J.-P.G.-E.; validation, T.Y. and J.Q.; formal analysis, T.Y.; investigation, T.Y.; resources, B.D.C. and D.C.M.; writing—original draft preparation, T.Y., J.Q., B.D.C. and D.C.M.; writing—review and editing, T.Y., S.W., J.Q., D.C.M. B.D.C., and J.-P.G.-E.; visualization, T.Y., J.Q. and S.W.; supervision, B.D.C., D.C.M. and J.-P.G.-E.; project administration, B.D.C. and D.C.M.; funding acquisition, B.D.C. and D.C.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The research presented in this study was supported by NASA's Earth Science Division and NASA's Terrestrial Ecology Program. Shanshan Wei was supported by the project "Computerized Management of Urban Trees: Tree Inspection Using Satellite Data" funded by the Ministry of National Development Research Fund awarded to the National Parks Board and Singapore-MIT Alliance for Research and Technology, Singapore. The authors are also thankful to the French National Space Center (CNES) for supporting DART development in the frame of the TOSCA project "Fluo3D" and the computer scientists in DART group (Nicolas Lauret, Jordan Guilleux and Eric Chavanon) for the continuous supports of DART. The authors would like to thank the anonymous reviewers for their suggestions to improve the quality of this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Nomenclature

<table><tr><td>A<sub>t</sub></td><td>Receiver telescope area</td></tr><tr><td>A<sub>s</sub></td><td>Footprint area at a certain distance</td></tr><tr><td>A<sub>proj</sub></td><td>Projection area</td></tr><tr><td>C<sub>cal</sub></td><td>A calibration constant with pre-defined instrumental and experimental configurations</td></tr><tr><td>d<sub>l</sub></td><td>The beam cross-section diameter at the "exit gate" of the laser generator</td></tr><tr><td>D<sub>r</sub></td><td>Diameter of the receiver telescope</td></tr><tr><td>eLAI</td><td>Effective LAI</td></tr><tr><td>f<sub>vcc</sub></td><td>Proportion of crown vertical projection area by neglecting the within-crown gaps</td></tr><tr><td>g(Ω<sub>f</sub>)</td><td>Leaf angle distribution function</td></tr><tr><td>G</td><td>Unit leaf area projection along the pulse direction</td></tr><tr><td>I</td><td>Integral of a return Gaussian power profile, total power of a return</td></tr><tr><td>I'</td><td>The distance-weighted power integral (I · R<sup>2</sup>), proportional to ρ<sub>a</sub></td></tr><tr><td>I'<sub>g</sub><sup>pure</sup></td><td>I of pulses with only single ground return</td></tr><tr><td>I'<sub>g</sub><sup>nearest</sup></td><td>For a pure vegetation or vegetation-ground pulse, refer to I'<sub>g</sub><sup>pure</sup> For the nearest pure ground pulse</td></tr><tr><td>i</td><td>Return index of a lidar pulse</td></tr><tr><td>IB</td><td>Intensity-based</td></tr><tr><td>LAI</td><td>Leaf area index</td></tr><tr><td>LAI<sub>est</sub></td><td>Estimated LAI from lidar points without using correlation coefficient</td></tr><tr><td>LPI</td><td>Laser penetration index</td></tr><tr><td>LPI<sub>pulse</sub></td><td>Penetration index of a single pulse</td></tr><tr><td>N</td><td>Number of returns</td></tr><tr><td>N<sub>first</sub></td><td>Number of first returns</td></tr><tr><td>N<sub>last</sub></td><td>Number of last returns</td></tr><tr><td>N<sub>single</sub></td><td>Number of pulses with only a single return</td></tr><tr><td>P̂</td><td>Peak amplitude of a return Gaussian profile</td></tr><tr><td>P<sub>r</sub>(t)</td><td>Time-dependent amplitude profile of a lidar waveform recorded by the receiver</td></tr><tr><td>P<sub>t</sub></td><td>Total pulse power</td></tr><tr><td>P<sub>t,g</sub></td><td>Total incident power onto the ground</td></tr><tr><td>P<sub>t,v</sub></td><td>Total incident power onto the vegetation part</td></tr></table>
<table><tr><td>P<sub>gap</sub>(θ)</td><td>Gap probability of zenith angle θ</td></tr><tr><td>P<sub>gap</sub></td><td>The gap probability along vertical direction adapted in most airborne lidar processing</td></tr><tr><td>PNB</td><td>Point-number-based</td></tr><tr><td>R</td><td>Distance from lidar</td></tr><tr><td>r<sub>fp</sub></td><td>Footprint radius</td></tr><tr><td>S&#770;</td><td>Peak amplitude of transmitted laser pulse</td></tr><tr><td>s<sub>s</sub></td><td>The temporal standard deviation of the convolved transmitted pulse and receiver response function</td></tr><tr><td>s<sub>p,i</sub></td><td>The temporal standard deviation of a returned Gaussian profile of return index i</td></tr><tr><td>s<sub>i</sub></td><td>The temporal standard deviation broadening of a return Gaussian profile (index i) from target interaction</td></tr><tr><td>s<sub>β</sub></td><td>The standard deviation of the angular energy distribution within footprint</td></tr><tr><td>T(Ω<sub>i</sub>)</td><td>Back-scatter transfer function of Ω<sub>i</sub> of return index i</td></tr><tr><td>α</td><td>Correlation coefficient to link LPI with P<sub>gap</sub> (exponential) and eLAI<sub>est</sub> with eLAI (linear)</td></tr><tr><td>β</td><td>The angular offset from the pulse direction</td></tr><tr><td>β<sub>t</sub></td><td>Footprint divergence</td></tr><tr><td>γ</td><td>Ratio of ground apparent reflectance over vegetation apparent reflectance</td></tr><tr><td>σ</td><td>Radar back-scattering differential cross-section</td></tr><tr><td>ρ</td><td>Natural reflectivity of a target</td></tr><tr><td>ρ<sub>a</sub></td><td>Apparent reflectance of a target</td></tr><tr><td>Ω<sub>i</sub>, ΔΩ<sub>i</sub></td><td>Direction and solid angle from target to the receiver of return index i</td></tr><tr><td>Ω<sub>f</sub></td><td>Leaf normal direction</td></tr><tr><td>η</td><td>System and atmospheric transmission factor</td></tr><tr><td>ω</td><td>The overall clumping index for LAI estimation over an area</td></tr><tr><td>ω<sub>in</sub></td><td>Index of clumping that is induced by only crown shape and within-crown clumping</td></tr><tr><td>θ<sub>g</sub></td><td>The angle between ground normal and the incident direction of a pulse</td></tr></table>

## Appendix A DART Workflow of Point-Cloud Modeling

Figure A1 shows the workflow for simulating lidar point clouds. For the first step, DART reads a 3-D scene (i.e., an abstract description or field of 3-D objects) with user-defined experimental configurations for lidar acquisition (e.g., altitude, platform trajectory, pulse divergence, scan density, control points, etc.). The stages enclosed in the dashed box are DART internal processes. DART uses a quasi-Monte Carlo ray-tracing approach to model the waveform for each pulse [48]. DART-simulated waveforms can either be exported as a binary file, or processed as points with GD before being exported as a text file.

<figure><img src="image_13.png" /><figcaption>**Figure A1.** Workflow of DART outputs: discrete points in text format or discrete points with associated waveforms in LAS format.</figcaption></figure>
1. yOutput exportation to LAS:

The DART-simulated waveforms are stored as a binary file. In the next step, this binary file
can be converted into standard lidar data formats to adapt to the existing software. Previous work
documented the output conversion into sorted pulse data (SPD) format [58,83]. Here, we implement
an approach that processes the binary file into an LAS 1.3 point cloud through the LASPY library [84].
The points are stored in an LAS file, with the possibility to export the waveforms to another WDP
file that is linked to each decomposed point [Header description in [85], point type 1 or 4]. The file is
supported by existing LAS visualization and processing software (e.g., CloudCompare [86]).

2. Internal waveform processing:

Waveforms are decomposed into points in DART code that are saved as matrices in a text file. Each line represents a point. Columns store 3-D position and "intensity" information, including peak amplitude (P̂), temporal standard deviation (s<sub>p,i</sub>), integral (∫P̂ ⋅ s<sub>p,i</sub>), apparent reflectance (ρ<sub>a</sub>), return index, number of returns per pulse, etc.

The dual-output option provides great flexibility in different applications. For example, Option 1 supports the validation of algorithms that work with ALS waveform data in LAS format (e.g., Riegl's Q-line devices). Option 2 may be preferable for numerous points (e.g., Riegl's V-line devices with online waveform processing). In that case, storing waveforms in a binary file would require tremendous unnecessary computer memory and hard disk space. Additionally, users may want to investigate the sensitivity of various input instrumental variables (d<sub>1</sub>, β<sub>t</sub>, ...) or output intensity selections (P̂, σ, or P̂ ⋅ σ...) that can be stored in the extra bytes of LAS format for different applications (e.g., classification, leaf angle distribution inversion, etc.).

**Appendix B Influence of Peak Detection Threshold**

For actual devices, the parameters that determine whether a return can be detected by a lidar
device are R and ρ_a, due to the sensor capability. In this section, the points produced by GD are
filtered with threshold ρ_a > 0.009, which is a typical value retrieved from the point cloud generated by
Riegl VQ-480i.

Figure A2 illustrates the LPI calculated by the filtered point cloud under different LAIs for the homogeneous scene (LAI = 1, 3, and 6). Compared to the unfiltered LPIs illustrated in Figure 3, the filtered LPIs have a larger variance. When the LAI is small, the vegetation return usually has lower ρ_a, which may be filtered out. This causes an overestimation by LPI<sub>γ=2.10</sub>, as some energy from vegetation backscattering has been removed. LPI<sub>fitted</sub> is also influenced by the large divergence of dataset variance. The bias curve gives a parallel displacement. LPI<sub>nearest</sub> still provides accurate estimations due to the idealized ground reflectance and the adaptation with ρ_a filtering. For the PNB methods, the estimated LPIs are slightly lower than the ones from the unfiltered point cloud, because the average number of returns has been reduced (Figure A3). From Figure A3, the average number of ground returns decreases with increasing LAI, which means there are more pure-vegetation pulses. For the heterogeneous scene, the results are shown in Figure A4 and Table A1. These results correspond to a potential reason for the occlusion effect in ALS applications [87].
<table><caption>**Table A1.** Fitted $\omega$ for LPI and coefficient of determination $R^2$ for the LAI in Equation (12) for the heterogeneous scene after filtering. The reference values are bold.</caption><thead><tr><th rowspan="2">Variables</th><th colspan="6">PNB Methods</th><th colspan="2">IB Methods</th><th>References</th></tr><tr><th>LPI<sub>all</sub></th><th>LPI<sub>weighted</sub></th><th>LPI<sub>first</sub></th><th>LPI<sub>last</sub></th><th>LPI<sub>both</sub></th><th>LPI<sub>y=2.10</sub></th><th>LPI<sub>nearest</sub></th><th>LPI<sub>fitted</sub></th><th>P<sub>gap</sub></th></tr></thead><tbody><tr><td rowspan="4">LAI = 2<br>(Random)</td><td>α</td><td>0.47</td><td>0.94</td><td>0.64</td><td>6.11</td><td>1.44</td><td>1.05</td><td>1.00</td><td>1.07</td><td><strong>1.00</strong></td></tr><tr><td>α · ω<sup>-1</sup></td><td>1.70</td><td>2.01</td><td>1.26</td><td>7.26</td><td>3.35</td><td>2.10</td><td>2.10</td><td>2.10</td><td><strong>1.98</strong></td></tr><tr><td>computed</td><td>0.28</td><td>0.47</td><td>0.51</td><td>0.84</td><td>0.43</td><td>0.50</td><td>0.48</td><td>0.51</td><td><strong>0.51</strong></td></tr><tr><td>ω</td><td>0.90</td><td>0.87</td><td>0.78</td><td>0.78</td><td>0.91</td><td>0.89</td><td>0.89</td><td>0.89</td><td><strong>0.90</strong></td></tr><tr><td rowspan="4">LAI = 2<br>(Clumping)</td><td>α</td><td>0.49</td><td>1.03</td><td>0.71</td><td>4.85</td><td>1.45</td><td>1.04</td><td>1.02</td><td>1.07</td><td><strong>1.00</strong></td></tr><tr><td>α · ω<sup>-1</sup></td><td>1.66</td><td>2.36</td><td>1.51</td><td>6.99</td><td>3.51</td><td>2.29</td><td>2.27</td><td>2.34</td><td><strong>2.23</strong></td></tr><tr><td>computed</td><td>0.30</td><td>0.44</td><td>0.47</td><td>0.69</td><td>0.41</td><td>0.45</td><td>0.45</td><td>0.46</td><td><strong>0.45</strong></td></tr><tr><td>ω</td><td>0.82</td><td>0.86</td><td>0.81</td><td>0.81</td><td>0.94</td><td>0.89</td><td>0.89</td><td>0.89</td><td><strong>0.90</strong></td></tr><tr><td rowspan="4">LAI = 1</td><td>α</td><td>0.33</td><td>0.71</td><td>0.47</td><td>18.69</td><td>1.18</td><td>1.09</td><td>1.00</td><td>1.11</td><td><strong>1.00</strong></td></tr><tr><td>α · ω<sup>-1</sup></td><td>0.90</td><td>1.11</td><td>0.67</td><td>0.81</td><td>2.10</td><td>1.59</td><td>1.48</td><td>1.63</td><td><strong>1.48</strong></td></tr><tr><td>computed</td><td>0.37</td><td>0.64</td><td>0.70</td><td>1.73</td><td>0.56</td><td>0.69</td><td>0.68</td><td>0.68</td><td><strong>0.68</strong></td></tr><tr><td>ω</td><td>0.91</td><td>0.87</td><td>0.81</td><td>0.14</td><td>0.82</td><td>0.95</td><td>0.95</td><td>0.95</td><td><strong>0.95</strong></td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Figure A2.** LPI estimated by different methods with a filtered point cloud with various LAIs [1(a), 3(b), and 6(c)]. The dashed line indicates $P_{gap}$ for reference.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure A3.** Average discrete returns after filtering for the homogeneous scene: (a) average number of returns under different footprint diameters and (b) average number (probability) of ground returns.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure A4.** Reference LAI plotted against $\ln(\text{LPI}^{-1})$ for the heterogeneous scene after filtering (LAI = 1 and 2) under footprint size 0.3 m: (a,b) Correlation between reference LPI and lidar-derived LPI; (c,d) Correlation between reference LAI and eLAI$_{est}$ from Equation (12).</figcaption></figure>
## References

1. Eitel, J.U.; Höfle, B.; Vierling, L.A.; Abellán, A.; Asner, G.P.; Deems, J.S.; Glennie, C.L.; Joerg, P.C.; LeWinter, A.L.; Magney, T.S.; et al. Beyond 3-D: The new spectrum of lidar applications for earth and ecological sciences. *Remote Sens. Environ.* **2016**, *186*, 372–392. [CrossRef]

2. RDubayah, R.; Goetz, S.J.; Blair, J.B.; Fatoyinbo, T.E.; Hansen, M.; Healey, S.P.; Hofton, M.A.; Hurtt, G.C.; Kellner, J.; Luthcke, S.B.; et al. The global ecosystem dynamics investigation. In Proceedings of the American Geophysical Union, Fall Meeting (AGU, San Francisco), San Francisco, CA, USA, 15–19 December 2014.

3. Abdalati, W.; Zwally, H.J.; Bindschadler, R.; Csatho, B.; Farrell, S.L.; Fricker, H.A.; Harding, D.; Kwok, R.; Lefsky, M.; Markus, T.; et al. The ICESat-2 Laser Altimetry Mission. *Proc. IEEE* **2010**, *98*, 735–751. [CrossRef]

4. Blair, J.B.; Rabine, D.L.; Hofton, M.A. The Laser Vegetation Imaging Sensor: A medium-altitude, digitisation-only, airborne laser altimeter for mapping vegetation and topography. *ISPRS J. Photogramm. Remote Sens.* **1999**, *54*, 115–122. [CrossRef]

5. Hyde, P.; Dubayah, R.; Peterson, B.; Blair, J.B.; Hofton, M.; Hunsaker, C.; Knox, R.; Walker, W. Mapping forest structure for wildlife habitat analysis using waveform lidar: Validation of montane ecosystems. *Remote Sens. Environ.* **2005**, *96*, 427–437. [CrossRef]

6. Asner, G.P.; Boardman, J.; Field, C.B.; Knapp, D.E.; Kennedy-Bowdoin, T.; Jones, M.O. Carnegie airborne observatory: In-flight fusion of hyperspectral imaging and waveform light detection and ranging for three-dimensional studies of ecosystems. *J. Appl. Remote Sens.* **2007**, *1*, 13536. [CrossRef]

7. Asner, G.P.; Knapp, D.E.; Boardman, J.; Green, R.O.; Kennedy-Bowdoin, T.; Eastwood, M.; Martin, R.E.; Anderson, C.; Field, C.B. Carnegie Airborne Observatory-2: Increasing science data dimensionality via high-fidelity multi-sensor fusion. *Remote Sens. Environ.* **2012**, *124*, 454–465. [CrossRef]

8. Cook, B.; Nelson, R.; Middleton, E.; Morton, D.; McCorkel, J.; Masek, J.; Ranson, K.; Ly, V.; Montesano, P. NASA Goddard's LiDAR, Hyperspectral and Thermal (G-LiHT) Airborne Imager. *Remote Sens.* **2013**, *5*, 4045. [CrossRef]

9. Côté, J.-F.; Widlowski, J.-L.; Fournier, R.A.; Verstraete, M.M. The structural and radiative consistency of three-dimensional tree reconstructions from terrestrial lidar. *Remote Sens. Environ.* **2009**, *113*, 1067–1081. [CrossRef]

10. Wagner, W.; Ullrich, A.; Melzer, T.; Briese, C.; Kraus, K. *From Single-Pulse to Full-Waveform Airborne Laser Scanners: Potential and Practical Challenges*. Remote Sensing and Spatial Information Sciences 35 (Part B3). 2004, pp. 201–206. Available online: https://publik.tuwien.ac.at/files/PubDat_119591.pdf (accessed on 7 December 2019).

11. Ni-Meister, W.; Lee, S.; Strahler, A.H.; Woodcock, C.E.; Schaaf, C.; Yao, T.; Ranson, K.J.; Sun, G.; Blair, J.B. Assessing general relationships between aboveground biomass and vegetation structure parameters for improved carbon estimate from lidar remote sensing. *J. Geophys. Res. Biogeosci.* **2010**, *115*, G00E11. [CrossRef]

12. Tang, H.; Dubayah, R.; Swatantran, A.; Hofton, M.; Sheldon, S.; Clark, D.B.; Blair, B. Retrieval of vertical LAI profiles over tropical rain forests using waveform lidar at La Selva, Costa Rica. *Remote Sens. Environ.* **2012**, *124*, 242–250. [CrossRef]

13. Tang, H.; Brolly, M.; Zhao, F.; Strahler, A.H.; Schaaf, C.L.; Ganguly, S.; Zhang, G.; Dubayah, R. Deriving and validating Leaf Area Index (LAI) at multiple spatial scales through lidar remote sensing: A case study in Sierra National Forest, CA. *Remote Sens. Environ.* **2014**, *143*, 131–141. [CrossRef]

14. Tang, H.; Dubayah, R.; Brolly, M.; Ganguly, S.; Zhang, G. Large-scale retrieval of leaf area index and vertical foliage profile from the spaceborne waveform lidar (GLAS/ICESat). *Remote Sens. Environ.* **2014**, *154*, 8–18. [CrossRef]

15. Harding, D.J.; Carabajal, C.C. ICESat waveform measurements of within-footprint topographic relief and vegetation vertical structure. *Geophys. Res. Lett.* **2005**, *32*, L21S10. [CrossRef]

16. Armston, J.; Disney, M.; Lewis, P.; Scarth, P.; Phinn, S.; Lucas, R.; Bunting, P.; Goodwin, N. Direct retrieval of canopy gap probability using airborne waveform lidar. *Remote Sens. Environ.* **2013**, *134*, 24–38. [CrossRef]

17. Chen, X.T.; Disney, M.I.; Lewis, P.; Armston, J.; Han, J.T.; Li, J.C. Sensitivity of direct canopy gap fraction retrieval from airborne waveform lidar to topography and survey characteristics. *Remote Sens. Environ.* **2014**, *143*, 15–25. [CrossRef]
18. Hancock, S.; Essery, R.; Reid, T.; Carle, J.; Baxter, R.; Rutter, N.; Huntley, B. Characterising forest gap fraction with terrestrial lidar and photography: An examination of relative limitations. *Agric. For. Meteorol.* **2014**, *189–190*, 105–114. [CrossRef]

19. Hancock, S.; Gaulton, R.; Danson, F.M. Angular Reflectance of Leaves With a Dual-Wavelength Terrestrial Lidar and Its Implications for Leaf-Bark Separation and Leaf Moisture Estimation. *IEEE Trans. Geosci. Remote Sens.* **2017**, *55*, 3084–3090. [CrossRef]

20. Chen, J.M.; Rich, P.M.; Gower, S.T.; Norman, J.M.; Plummer, S. Leaf area index of boreal forests: Theory, techniques, and measurements. *J. Geophys. Res. Atmos.* **1997**, *102*, 29429–29443. [CrossRef]

21. Chen, J.M.; Menges, C.H.; Leblanc, S.G. Global mapping of foliage clumping index using multi-angular satellite data. *Remote Sens. Environ.* **2005**, *97*, 447–457. [CrossRef]

22. Wei, S.; Fang, H.; Schaaf, C.B.; He, L.; Chen, J.M. Global 500 m clumping index product derived from MODIS BRDF data (2001–2017). *Remote Sens. Environ.* **2019**, *232*, 111296. [CrossRef]

23. Hu, R.; Yan, G.; Nerry, F.; Liu, Y.; Jiang, Y.; Wang, S. Using Airborne Laser Scanner and Path Length Distribution Model to Quantify Clumping Effect and Estimate Leaf Area Index. *IEEE Trans. Geosci. Remote Sens.* **2018**, *56*, 3196–3209. [CrossRef]

24. Anderson, K.; Hancock, S.; Disney, M.; Gaston, K.J. Is waveform worth it? A comparison of LiDAR approaches for vegetation and landscape characterization. *Remote Sens. Ecol. Conserv.* **2016**, *2*, 5–15. [CrossRef]

25. Disney, M.I.; Kalogirou, V.; Lewis, P.; Prieto-Blanco, A.; Hancock, S.; Pfeifer, M. Simulating the impact of discrete-return lidar system and survey characteristics over young conifer and broadleaf forests. *Remote Sens. Environ.* **2010**, *114*, 1546–1560. [CrossRef]

26. Hancock, S.; Armston, J.; Li, Z.; Gaulton, R.; Lewis, P.; Disney, M.; Danson, F.M.; Strahler, A.; Schaaf, C.; Anderson, K.; et al. Waveform lidar over vegetation: An evaluation of inversion methods for estimating return energy. *Remote Sens. Environ.* **2015**, *164*, 208–224. [CrossRef]

27. Yan, G.; Hu, R.; Luo, J.; Weiss, M.; Jiang, H.; Mu, X.; Xie, D.; Zhang, W. Review of indirect optical measurements of leaf area index: Recent advances, challenges, and perspectives. *Agric. Forest Meteorol.* **2019**, *265*, 390–411. [CrossRef]

28. Grau, E.; Durrieu, S.; Fournier, R.; Gastellu-Etchegorry, J.-P.; Yin, T. Estimation of 3D vegetation density with Terrestrial Laser Scanning data using voxels. A sensitivity analysis of influencing parameters. *Remote Sens. Environ.* **2017**, *191*, 373–388. [CrossRef]

29. Béland, M.; Widlowski, J.-L.; Fournier, R.A.; Côté, J.-F.; Verstraete, M.M. Estimating leaf area distribution in savanna trees from terrestrial LiDAR measurements. *Agric. Forest Meteorol.* **2011**, *151*, 1252–1266. [CrossRef]

30. Solberg, S.; Næsset, E.; Hanssen, K.H.; Christiansen, E. Mapping defoliation during a severe insect attack on Scots pine using airborne laser scanning. *Remote Sens. Environ.* **2006**, *102*, 364–376. [CrossRef]

31. Solberg, S.; Brunner, A.; Hanssen, K.H.; Lange, H.; Næsset, E.; Rautiainen, M.; Stenberg, P. Mapping LAI in a Norway spruce forest using airborne laser scanning. *Remote Sens. Environ.* **2009**, *113*, 2317–2327. [CrossRef]

32. Solberg, S. Mapping gap fraction, LAI and defoliation using various ALS penetration variables. *Int. J. Remote Sens.* **2010**, *31*, 1227–1244. [CrossRef]

33. Morsdorf, F.; Kötz, B.; Meier, E.; Itten, K.I.; Allgöwer, B. Estimation of LAI and fractional cover from small footprint airborne laser scanning data based on gap fraction. *Remote Sens. Environ.* **2006**, *104*, 50–61. [CrossRef]

34. Korhonen, L.; Korpela, I.; Heiskanen, J.; Maltamo, M. Airborne discrete-return LIDAR data in the estimation of vertical canopy cover, angular canopy closure and leaf area index. *Remote Sens. Environ.* **2011**, *115*, 1065–1080. [CrossRef]

35. Schneider, F.D.; Leiterer, R.; Morsdorf, F.; Gastellu-Etchegorry, J.P.; Lauret, N.; Pfeifer, N.; Schaepman, M.E. Simulating imaging spectrometer data: 3D forest modeling based on LiDAR and in situ data. *Remote Sens. Environ.* **2014**, *152*, 235–250. [CrossRef]

36. Luo, S.-Z.; Wang, C.; Zhang, G.-B.; Xi, X.-H.; Li, G.-C. Forest Leaf Area Index (LAI) Estimation Using Airborne Discrete-Return Lidar Data. *Chin. J. Geophys.* **2013**, *56*, 233–242.

37. Wagner, W.; Ullrich, A.; Ducic, V.; Melzer, T.; Studnicka, N. Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner. *ISPRS J. Photogramm. Remote Sens.* **2006**, *60*, 100–112. [CrossRef]
38. Ullrich, A.; Pfennigbauer, M. Categorisation of Full Waveform Data Provided by Laser Scanning Devices. In *Electro-Optical Remote Sensing, Photonic Technologies, and Applications*; SPIE Security + Defence: Bellingham, WA, USA, 2011; Volume 8186.

39. Pfennigbauer, M.; Rieger, P.; Studnicka, N.; Ullrich, A. Detection of concealed objects with a mobile laser scanning system. In *Laser Radar Technology and Applications XIV*; SPIE Defense, Security, and Sensing: Orlando, FL, USA, 2009; p. 732308.

40. Pfennigbauer, M.; Ullrich, A. Improving quality of laser scanning data acquisition through calibrated amplitude and pulse deviation measurement. In *Laser Radar Technology and Applications XV*; SPIE Defense, Security, and Sensing: Orlando, FL, USA, 2010; p. 76841F.

41. Pfennigbauer, M.; Wolf, C.; Ullrich, A. Enhancing online waveform processing by adding new point attributes. In *Proc. SPIE 8731, Laser Radar Technology and Applications XVIII*; SPIE Defense, Security, and Sensing: Bellingham, WA, USA, 2013; p. 873104.

42. Riegl. LAS Extrabytes Implementation in RIEGL Software-WHITEPAPER. 2017. Available online: http://www.riegl.com/uploads/tx_pxpriegldownloads/Whitepaper_LASextrabytes_implementation_in-RIEGLSoftware_2017-12-04.pdf (accessed on 7 December 2019).

43. Schofield, L.A.; Danson, F.M.; Entwistle, N.S.; Gaulton, R.; Hancock, S. Radiometric calibration of a dual-wavelength terrestrial laser scanner using neural networks. *Remote Sens. Lett.* **2016**, *7*, 299–308.
[CrossRef]

44. Wagner, W. Radiometric calibration of small-footprint full-waveform airborne laser scanner measurements: Basic physical concepts. *ISPRS J. Photogramm. Remote Sens.* **2010**, *65*, 505–513. [CrossRef]

45. Roncat, A.; Briese, C.; Jansa, J.; Pfeifer, N. Radiometrically Calibrated Features of Full-Waveform Lidar Point Clouds Based on Statistical Moments. *IEEE Geosci. Remote Sens. Lett.* **2014**, *11*, 549–553. [CrossRef]

46. Sun, G.; Ranson, K.J. Modeling lidar returns from forest canopies. *Geosci. Remote Sens. IEEE Trans.* **2000**, *38*, 2617–2626.

47. North, P.R.J.; Rosette, J.A.B.; Suárez, J.C.; Los, S.O. A Monte Carlo radiative transfer model of satellite waveform LiDAR. *Int. J. Remote Sens.* **2010**, *31*, 1343–1358. [CrossRef]

48. Gastellu-Etchegorry, J.P.; Yin, T.; Lauret, N.; Grau, E.; Rubio, J.; Cook, B.D.; Morton, D.C.; Sun, G. Simulation of satellite, airborne and terrestrial LiDAR with DART (I): Waveform simulation with quasi-Monte Carlo ray tracing. *Remote Sens. Environ.* **2016**, *184*, 418–435. [CrossRef]

49. Ni-Meister, W.; Yang, W.; Lee, S.; Strahler, A.H.; Zhao, F. Validating modeled lidar waveforms in forest canopies with airborne laser scanning data. *Remote Sens. Environ.* **2018**, *204*, 229–243. [CrossRef]

50. Brown, S.D.; Blevins, D.D.; Schott, J.R. Time-Gated Topographic LIDAR Scene Simulation. In *SPIE Proceedings Volume 5791, Laser Radar Technology and Applications X*; SPIE Defense, Security, and Sensing: Orlando, FL, USA, 2005; pp. 342–353.

51. Wu, J.; Aardt, J.A.N.V.; Asner, G.P. A Comparison of Signal Deconvolution Algorithms Based on Small-Footprint LiDAR Waveform Simulation. *IEEE Trans. Geosci. Remote Sens.* **2011**, *49*, 2402–2414.
[CrossRef]

52. Govaerts, Y.M. *A Model of Light Scattering in Three-Dimensional Plant Canopies: A Monte Carlo Ray Tracing Approach*; JRC Catalogue No. CL-NA-16394-EN-C; Office for Official Publications of the European Communities: Luxembourg, 1996; 186p.

53. Govaerts, Y.M.; Verstraete, M.M. Raytran: A Monte Carlo ray-tracing model to compute light scattering in three-dimensional heterogeneous media. *IEEE Trans. Geosci. Remote Sens.* **1998**, *36*, 493–505. [CrossRef]

54. Disney, M.I.; Lewis, P.E.; Bouvet, M.; Prieto-Blanco, A.; Hancock, S. Quantifying Surface Reflectivity for Spaceborne Lidar via Two Independent Methods. *IEEE Trans. Geosci. Remote Sens.* **2009**, *47*, 3262–3271.
[CrossRef]

55. Kobayashi, H.; Iwabuchi, H. A coupled 1-D atmosphere and 3-D canopy radiative transfer model for canopy reflectance, light environment, and photosynthesis simulation in a heterogeneous landscape. *Remote Sens. Environ.* **2008**, *112*, 173–185. [CrossRef]

56. Gastellu-Etchegorry, J.P.; Yin, T.; Grau, E.; Lauret, N.; Rubio, J. Lidar radiative transfer modeling in the Atmosphere. In *Proceedings of the 2013 IEEE International Geoscience and Remote Sensing Symposium-IGARSS*, Melbourne, VIC, Australia, 21–26 July 2013; IEEE: Melbourne, Australia, 2013; pp. 4554–4557.
57. Widlowski, J.L.; Pinty, B.; Lopatka, M.; Atzberger, C.; Buzica, D.; Chelle, M.; Disney, M.; Gastellu-Etchegorry, J.P.; Gerboles, M.; Gobron, N.; et al. The fourth radiation transfer model intercomparison (RAMI-IV): Proficiency testing of canopy reflectance models with ISO-13528. *J. Geophys. Res. Atmos.* **2013**, *118*, 6869–6890. [CrossRef]

58. Yin, T.; Lauret, N.; Gastellu-Etchegorry, J.-P. Simulation of satellite, airborne and terrestrial LiDAR with DART (II): ALS and TLS multi-pulse acquisitions, photon counting, and solar noise. *Remote Sens. Environ.* **2016**, *184*, 454–468. [CrossRef]

59. Qi, J.; Yin, T.; Xie, D.; Gastellu-Etchegorry, J.-P. Hybrid Scene Structuring for Accelerating 3D Radiative Transfer Simulations. *Remote Sens.* **2019**, *11*, 2637. [CrossRef]

60. Roundy, C.B. *Current Technology of Laser Beam Profile Measurements*; Spiricon. Inc., 2000. Available online: http://aries.ucsd.edu/LASERLAB/TUTOR/profile-tutorial.pdf (accessed on 7 December 2019).

61. Mallet, C.; Bretar, F. Full-waveform topographic lidar: State-of-the-art. *ISPRS J. Photogramm. Remote Sens.* **2009**, *64*, 1–16. [CrossRef]

62. Schaepman-Strub, G.; Schaepman, M.E.; Painter, T.H.; Dangel, S.; Martonchik, J.V. Reflectance quantities in optical remote sensing—definitions and case studies. *Remote Sens. Environ.* **2006**, *103*, 27–42. [CrossRef]

63. Newville, M.; Stensitzki, T.; Allen, D.B.; Rawlik, M.; Ingargiola, A.; Nelson, A. *Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting For Python*; Astrophysics Source Code Library, 2016.

64. Milenković, M.; Wagner, W.; Quast, R.; Hollaus, M.; Ressl, C.; Pfeifer, N. Total canopy transmittance estimated from small-footprint, full-waveform airborne LiDAR. *ISPRS J. Photogramm. Remote Sens.* **2017**, *128*, 61–72. [CrossRef]

65. Zheng, G.; Ma, L.; Eitel, J.U.; He, W.; Magney, T.S.; Moskal, L.M.; Li, M. Retrieving Directional Gap Fraction, Extinction Coefficient, and Effective Leaf Area Index by Incorporating Scan Angle Information From Discrete Aerial Lidar Data. *IEEE Trans. Geosci. Remote Sens.* **2017**, *55*, 577–590. [CrossRef]

66. Roussel, J.-R.; Béland, M.; Caspersen, J.; Achim, A. A mathematical framework to describe the effect of beam incidence angle on metrics derived from airborne LiDAR: The case of forest canopies approaching turbid medium behaviour. *Remote Sens. Environ.* **2018**, *209*, 824–834. [CrossRef]

67. Zhao, K.; Popescu, S.; Nelson, R. Lidar remote sensing of forest biomass: A scale-invariant estimation approach using airborne lasers. *Remote Sens. Environ.* **2009**, *113*, 182–196. [CrossRef]

68. Fleck, S.; Raspe, S.; Čater, M.; Schleppi, P.; Ukonmaanaho, L.; Greve, M.; Hertel, C.; Weis, W.; Rumpf, S. Leaf area measurements. In *Manual Part XVII. United Nations Economic Commission for Europe Convention on Long-Range Transboundary Air Pollution, ICP Forests, Hamburg*; Thünen Institute of Forest Ecosystems: Eberswalde, Germany, 2012.

69. Riaño, D.; Valladares, F.; Condés, S.; Chuvieco, E. Estimation of leaf area index and covered ground from airborne laser scanner (Lidar) in two contrasting forests. *Agric. Forest Meteorol.* **2004**, *124*, 269–275. [CrossRef]

70. Lovell, J.L.; Jupp, D.L.B.; Culvenor, D.S.; Coops, N.C. Using airborne and ground-based ranging lidar to measure canopy structure in Australian forests. *Can. J. Remote Sens.* **2003**, *29*, 607–622. [CrossRef]

71. Cook, B.D.; Bolstad, P.V.; Næsset, E.; Anderson, R.S.; Garrigues, S.; Morisette, J.T.; Nickeson, J.; Davis, K.J. Using LiDAR and quickbird data to model plant production and quantify uncertainties associated with wetland detection and land cover generalizations. *Remote Sens. Environ.* **2009**, *113*, 2366–2379. [CrossRef]

72. Lefsky, M.A.; Cohen, W.B.; Acker, S.A.; Parker, G.G.; Spies, T.A.; Harding, D. Lidar Remote Sensing of the Canopy Structure and Biophysical Properties of Douglas-Fir Western Hemlock Forests. *Remote Sens. Environ.* **1999**, *70*, 339–361. [CrossRef]

73. Kashani, A.G.; Olsen, M.J.; Parrish, C.E.; Wilson, N. A Review of LIDAR Radiometric Processing: From Ad Hoc Intensity Correction to Rigorous Radiometric Calibration. *Sensors* **2015**, *15*, 28099. [CrossRef]

74. Zhang, W.; Qi, J.; Wan, P.; Wang, H.; Xie, D.; Wang, X.; Yan, G. An Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation. *Remote Sens.* **2016**, *8*, 501. [CrossRef]

75. Zhang, K.; Chen, S.-C.; Whitman, D.; Shyu, M.-L.; Yan, J.; Zhang, C. A progressive morphological filter for removing nonground measurements from airborne LIDAR data. *IEEE Trans. Geosci. Remote Sens.* **2003**, *41*, 872–882. [CrossRef]

76. Gastellu-Etchegorry, J.P.; Yin, T.; Lauret, N.; Cajgfinger, T.; Gregoire, T.; Grau, E.; Feret, J.B.; Lopes, M.; Guilleux, J.; Dedieu, G.; et al. Discrete Anisotropic Radiative Transfer (DART 5) for Modeling Airborne and Satellite Spectroradiometer and LIDAR Acquisitions of Natural and Urban Landscapes. *Remote Sens.* **2015**, *7*, 1667–1701. [CrossRef]
77. Gastellu-Etchegorry, J.P.; Martin, E.; Gascon, F. DART: A 3D model for simulating satellite images and studying surface radiation budget. *Int. J. Remote Sens.* **2004**, *25*, 73–96. [CrossRef]

78. Gastellu-Etchegorry, J.-P. 3D modeling of satellite spectral images, radiation budget and energy budget of urban landscapes. *Meteorol. Atmos. Phys.* **2008**, *102*, 187. [CrossRef]

79. Danson, F.M.; Hetherington, D.; Morsdorf, F.; Koetz, B.; Allgower, B. Forest Canopy Gap Fraction From Terrestrial Laser Scanning. *IEEE Geosci. Remote Sens. Lett.* **2007**, *4*, 157–160. [CrossRef]

80. Abshire, J.B.; Sun, X.; Riris, H.; Sirota, J.M.; McGarry, J.F.; Palm, S.; Yi, D.; Liiva, P. Geoscience Laser Altimeter System (GLAS) on the ICESat Mission: On-orbit measurement performance. *Geophys. Res. Lett.* **2005**, *32*, [CrossRef]

81. Vincent, G.; Antin, C.; Laurans, M.; Heurtebize, J.; Durrieu, S.; Lavalley, C.; Dauzat, J. Mapping plant area index of tropical evergreen forest by airborne laser scanning. A cross-validation study using LAI2200 optical sensor. *Remote Sens. Environ.* **2017**, *198*, 254–266. [CrossRef]

82. Béland, M.; Baldocchi, D.D.; Widlowski, J.-L.; Fournier, R.A.; Verstraete, M.M. On seeing the wood from the leaves and the role of voxel size in determining leaf area distribution of forests with terrestrial LiDAR. *Agric. For. Meteorol.* **2014**, *184*, 82–97. [CrossRef]

83. Bunting, P.; Armston, J.; Lucas, R.M.; Clewley, D. Sorted Pulse Data (SPD) Library. Part I: A generic file format for LiDAR data from pulsed laser systems in terrestrial environments. *Comput. Geosci.* **2013**, *56*, 197–206. [CrossRef]

84. Brown, G. Laspy: Documentation. 2018. Available online: https://github.com/grantbrown/laspy (accessed on 7 December 2019).

85. ASPRS. LAS SPECIFICATION VERSION 1.3–R10. 2009. Available online: www.asprs.org/a/society/committees/standards/asprs_las_spec_v13.pdf (accessed on 7 December 2019).

86. Girardeau-Montaut, D. *Cloudcompare-Open Source Project*; OpenSource Project; 2011; Available online: http://www.cloudcompare.org/ (accessed on 7 December 2019).

87. Kükenbrink, D.; Schneider, F.D.; Leiterer, R.; Schaepman, M.E.; Morsdorf, F. Quantification of hidden canopy volume of airborne laser scanning data using a voxel traversal algorithm. *Remote Sens. Environ.* **2017**, *194*, 424–436. [CrossRef]

<figure><img src="image_17.png" /></figure>

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).