Direct Estimation of Forest Leaf Area Index based on Spectrally Corrected Airborne LiDAR Pulse Penetration Ratio

: The leaf area index (LAI) is a crucial structural parameter of forest canopies. Light Detection and Ranging (LiDAR) provides an alternative to passive optical sensors in the estimation of LAI from remotely sensed data. However, LiDAR-based LAI estimation typically relies on empirical models, and such methods can only be applied when the ﬁeld-based LAI data are available. Compared with an empirical model, a physically-based model—e.g., the Beer–Lambert law based light extinction model—is more attractive due to its independent dataset with training. However, two challenges are encountered when applying the physically-based model to estimate LAI from discrete LiDAR data: i.e., deriving the gap fraction and the extinction coe ﬃ cient from the LiDAR data. We solved the ﬁrst problem by integrating LiDAR and hyperspectral data to transfer the LiDAR penetration ratio to the forest gap fraction. For the second problem, the extinction coe ﬃ cient was estimated from tiled (1 km × 1 km) LiDAR data by nonlinearly optimizing the cost function of the angular LiDAR gap fraction and simulated gap fraction from the Beer–Lambert law model. A validation against LAI-2000 measurements showed that the estimates were signiﬁcantly correlated to the reference LAI with an R 2 of 0.66, a root mean square error (RMSE) of 0.60 and a relative RMSE of 0.15. We conclude that forest LAI can be directly estimated by the nonlinear optimization method utilizing the Beer–Lambert model and a spectrally corrected LiDAR penetration ratio. The signiﬁcance of the proposed method is that it can produce reliable remotely sensed forest LAI from discrete LiDAR and spectral data when ﬁeld-measured LAI are unavailable.


Introduction
The vegetation leaf area index (LAI), defined as one half of the total leaf area per unit ground surface area [1], is one of the important variables related to many ecological applications [2]. Remote sensing technology offers a cost-effective method for surveying LAI in a large coverage area of land. Passive remote sensing data from sensors such as MODIS are routinely used to estimate vegetation LAI for large areas [3]. Active Light Detection and Ranging (LiDAR) sensors are alternatives in estimating forest LAI [4]. The most commonly used method to estimate LAI from LiDAR data is based on fitting empirical models linking field-measured LAI to LiDAR-derived metrics. These metrics include cover fractions [5,6], height [7] or height percentiles [8], and varied penetration metrics [9,10]. However, in the above empirical methods, the field-measured LAI data are required as training data, and such data sets are not always available.
Physically-based methods are attractive because they do not rely on empirical relations but attempt to model the underlying physical processes, which can be extended across space and time. Most physically-based methods that link LAI with LiDAR-measured variables rely on the light extinction models, such as the Beer-Lambert Law [11,12], in which the vegetation gap fraction in a given direction is proportional to the exponential function of canopy LAI and the light extinction coefficient.
The canopy gap fraction is defined as the amount of open area within a canopy in a certain viewing angle [13]. A commonly used approach is to estimate the forest canopy gap fraction from the LiDAR pulse penetration ratio, which typically describes the fraction of LiDAR pulses with returns from ground level [14]. Although the percentage of ground returns is related to the canopy gap fraction, it is not equivalent to the actual gap fraction, so numerous studies have linked different penetration ratios with the gap fraction or gap fraction-like variables (e.g., cover fraction) using linear [15] or nonlinear [16] regression methods. The difference between the penetration ratio and gap fraction also hampers the use of physically-based models in the direct estimation of vegetation LAI at the regional scale [17].
To fill the gap between the LiDAR-derived penetration ratio and canopy gap fraction, Lefsky et al. [5] suggested that the penetration ratio could be adjusted as a function of spectral properties of soil and vegetation. However, they did not provide a physically-based approach to obtain a spectral calibration coefficient. Instead, they used an empirical constant. Ni-Meister et al. [18] proposed a physically-based model for the simulation of waveform LiDAR data, and Tang et al. [4] inverted this model to estimate the LAI in a tropical forest using waveform LiDAR data. In this model, the LiDAR penetration ratio derived from accumulated power of returns from vegetation and ground was corrected with a spectral coefficient. Armston et al. [17] and Chen et al. [19] further investigated the sensitivity of the model while they retrieved the gap fraction from the waveform LiDAR data. However, their methods were applied with waveform LiDAR data, which is not available as frequently as discrete return data. There are no reports considering the application of their methods in the direct estimation of LAI using discrete LiDAR data.
Another factor that prevents the LiDAR data from being directly used in the physically-based estimation of LAI is the difficulty of estimating the extinction coefficient, which is defined as the mean projection of the unit leaf area on the plane perpendicular to the direction of light beam [20,21]. For simplicity, the extinction coefficient is often assigned a value of 0.5 following an assumption of a spherical leaf distribution and a near-vertical viewing angle [22]. In practice, spherical leaf angle distribution is not always valid for all the tree species [23]. For example, Korhonen and Morsdorf [24] suggested that the average value of 0.38 for a vertical viewing angle is suitable for the boreal forest.
The estimation of the extinction coefficient requires information about the leaf angle distribution, which can be estimated by measuring the angular gap fractions. The difficulty of obtaining the angular gap fraction has resulted in fewer studies considering the estimation extinction coefficient from airborne LiDAR data. Zheng et al. [25] used overlapping airborne LiDAR flight lines to calculate the angular gap fraction, and then derived the plot-level extinction coefficient and corresponding LAI. However, most airborne LiDAR systems are designed for measuring surface elevation rather than vegetation structure, so the range of the scan zenith angles observed in plot level is usually limited due to the small area of the plot.
The main overall objective of this work is to develop a physically-based method that is capable of the direct estimation of LAI from airborne LiDAR data. The specific objective is to link discrete a LiDAR penetration ratio with the canopy gap fraction with the help of the target spectral information, and to derive the canopy extinction coefficient from angular gap fraction data. This is obtained by using a Beer-Lambert law-based light extinction model. To calculate the forest gap fraction, we follow the theory of Lefsky et al. [5] and Ni-Meister et al. [18], but extend their work from waveform to discrete LiDAR data, and apply hyperspectral data to calculate the spectral reflectivity of vegetation and background. To calculate the extinction coefficient, we use tiled (1 km × 1 km) LiDAR data, where a larger range of scan zenith angles is available in the estimation of the angular distribution of gap fractions. A nonlinear optimization method is then employed to construct a reliable estimate of leaf angle distribution for the research area. Finally, we average the extinction coefficient from tiled LiDAR data, and use it as input to a Beer-Lambert law-based model in the estimation of plot level LAI.

Study Area
The study area (34 • 36 40"N, 81 • 42 50"W) is located in the state of South Carolina, USA (Figure 1a,b) and the field experiment was conducted in the Calhoun Experimental Forest over an area of 13 km × 15 km [26] (Figure 1c). This investigation was conducted as part of the Calhoun Critical Zone Observatory (CZO). The dominant tree species is Loblolly pine (Pinus taeda L.), and there is a varying degree of hardwood species (e.g., Quercus). The mean tree height is 21.4 ± 6.4 m (mean ± standard variance), and the mean basal area is 0.08 ± 0.07 m 2 [27].

Datasets and Preprocessing
All the data in this study were downloaded from the Calhoun CZO program website (http: //criticalzone.org). The dataset includes the field LAI of LAI-2000 measurement, airborne discrete return LiDAR point cloud data and hyperspectral imagery.

Field Measured Data
Field LAI data were collected in July-August 2014 using the LAI-2000 (LI-COR Inc., Lincoln, NE) instrument. Six repeats with different directions were made for each plot. There were 35 LAI plots deployed in the research area; however, four of them were out of the LiDAR data coverage area ( Figure 1c). Therefore, we reported the result on the 31 plots while conducting a comparison with the LiDAR results.
In the field-collected LAI dataset, the LAI, mean tilt angle, and diffuse non-interceptance were reported for every measurement. Though the gap fraction is one of the key variables in the raw LAI-2000 output file, unfortunately, it is not included in the published dataset. Thus, we reproduced the field gap fraction following the below steps.
Firstly, the field mean tilt angle was used to calculate the leaf angle distribution (LAD) parameter. According to the generalized LAD ellipsoidal model [28], the leaf mean tilt angle can be calculated using the following equation: From the inverted form of the Equation (1), LAD parameter χ can be derived as where χ is the LAD model parameter and is defined as the ratio of vertical to horizontal projections of canopy elements, and α is the leaf mean tilt angle which is reported in the LAI-2000 output file. Secondly, the extinction coefficient can be approximated as [21] k(θ) = [χ 2 + tan 2 θ] where k(θ) is the extinction coefficient at the zenith angle θ. Then, based on the light extinction model of the Beer-Lambert law, in the case of a random spatial distribution of infinitely small leaves and assuming leaf random azimuth distribution, the gap fraction can be related to the LAI and extinction coefficient using the Poisson Model [12], where P gap (θ) is the gap fraction at the scan zenith angle θ. Thus, the gap fraction value corresponding to the measured LAI is reproduced. It should be noted that we use the concept of effective LAI here without the consideration of clumping correction for the nonrandomized forest leaves [29], so the clumping index is not included in the equation. The statistics of measurements taken at individual sample points within a plot were taken as the indicator of the spatial variance of the forest structure. At the plot level, we take the averaged LAI and gap fraction from the sample point data within the plot as the reference data while evaluating the LiDAR result.

LiDAR Data
The LiDAR survey was performed with an Optech Gemini Airborne Laser Terrain Mapper (ALTM) which operated at a wavelength of 1064 nm. The LiDAR survey took place over 2 days, starting on August 5 and finishing on 6 August 2014. The 3-D point cloud in LAS format (Version 1.2) was classified as ground or non-ground in 1 km square tiles. Survey parameters for the LiDAR portion are provided in Table 1. The LiDAR elevation data were normalized to the height above ground according to the classification of ground and vegetation and the algorithm of Parkan [30]. Two levels of LiDAR gap fractions were generated: plot and tile level. For the field plot level (red square in Figure 2), to match the field LAI-2000 measurement, the LiDAR data were clipped as 40 × 40 m 2 plots, where the centroids were the field plot locations. The size is larger than the field plot diameter, which ensures that the LiDAR data encompasses the ground measurements considering the geo-registration error. Figure 2 visualizes the difficulty of obtaining a sufficient range of angular gap fractions for plot-level LiDAR data due to the limited spatial size of the plots. Thus, at the plot level, only an averaged gap fraction was calculated. On the other hand, at the tile level, larger areas were available, with more variation in the scan zenith angles. At the tile level, all points were grouped into different sliced bins as a function of the scan zenith angle. Then, the points in each bin were used to calculate the angular gap fraction at that specific angle. Here, the range of scan zenith angles within each bin was set as 3 degrees.

Airborne Hyperspectral Data
We extracted the soil-vegetation reflectivity ratio from the airborne hyperspectral survey acquired on 30 July 2014. The hyperspectral-imaging mission was acquired by the CASI 1500 (ITRES, Calgary, Alberta, Canada) sensor, which was flown at an altitude of 2150 m above ground level. The central wavelengths of the hyperspectral image bands ranged from 368 to 1041 nm, and the Full Width Half Maximum (FWHM) of bands were approximately 7.2 nm. The spatial resolution was one meter. The hyperspectral data were delivered as digital numbers (DN) that were directly used to identify the soil and vegetation pixels and to calculate the spectral ratio of the targets [31].

Overview
The focus of this study is to analyze the LiDAR and hyperspectral data to extract the forest gap fraction and extinction coefficient independently to ultimately estimate the LAI. Gap fractions are obtained by utilizing a combination of the LiDAR penetration ratio and soil-vegetation reflectivity ratio obtained from hyperspectral data. The forest extinction coefficient is obtained from tiled LiDAR data through a nonlinear optimization procedure, which utilizes the large angular variation present in tiled LiDAR data. The flowchart of the methodology is outlined in Figure 3 and is explained in detail below.

Spectral Correction Coefficient
When the hyperspectral image is available, the calculation of the spectral ratio coefficient is reduced to the extraction of the target optical properties in the given image. We begin this procedure by classifying the image pixels as vegetation and soil. The classification is implemented by the identification of soil and vegetation lines.
The soil and vegetation lines represent a linear relationship for bare soil and fully vegetated pixels, respectively [32]. As indicated in Figure 4, in the two-dimensional space of the red and NIR band, the soil or vegetation line can be characterized by slope and intercept parameters: where r nir and r red are the DN value red and near infrared bands, respectively, a is the slope of the soil or vegetation line and b is the intercept of the soil or vegetation line. The automatic extraction method estimates the soil and vegetation line parameters by deriving a set of maximum and minimum red bands across the range of near-infrared band values [33].
Assuming that all the pixels located around the soil line or the vegetation line are pure pixels, we can calculate the averaged digital number of soil and vegetation in the near-infrared band; i.e., DN grd and DN veg , respectively. Since we focus on the relative value of the soil and vegetation reflectance rather than the absolute value, we take the ratio of DN value as the approximation to the ratio of reflectance of soil and vegetation (r grd , r veg ); i.e., r nir = a × r red + b. We use an iterative procedure for soil and vegetation line extraction in the two-dimensional space of the red and near-infrared band [33]. The wavelength range of the hyperspectral image was 368-1041 nm. Here, we select the 24th band (697.9 nm) as the red band, and the 48th band (1041 nm) as the NIR band. The wavelength of the NIR band differs from the LiDAR sensor that operates at 1064 nm, but the difference in reflectance caused by the bias of the wavelength can be ignored, according to the statistic on the reflectance data from the spectral library of SPECCHIO (Spectral Input/Output) [34].

Gap Fraction from Spectrally Corrected LiDAR Penetration Ratio
The foundation of calculating the gap fraction from the LiDAR penetration ratio is based on the assumption that the number of returns is related to the gap fraction; i.e., a bigger forest gap results in more returns from the ground and fewer returns from leaves, and vice versa [35]. The above concept can be formulated as where N grd and N veg are the return numbers of the ground and vegetation, respectively, and P gap is the gap fraction. However, the above concept formulation neglects the spectral difference of the vegetation and ground. We can imagine an extreme case where the ground reflectance is zero and will produce a N grd value of zero; i.e., all the signals will be absorbed by ground surface even though they can penetrate to pass leaves through canopy gap. Apparently, this result contradicts the concept of the equation.
To account for the target optical difference while calculating the LiDAR gap fraction, we revised the above conceptual formulation equation as where ρ veg is the leaf volume backscattering coefficient, which is the function of the LAD, the phase function of leaf scattering and the optical properties [18], and ρ grd is the backscattering coefficient of the ground. It is worth pointing out that the number of returned points not only relates to the gap fraction and the optical property of the target, but also to other parameters, such as the flight height [36] and the LiDAR sensor characteristics. However, the influence of these parameters can be neglected because in a certain survey campaign they will be kept stable or less varied across different fly strips. The above formulation can be further simplified as where P lidar is the penetration ratio, which can be directly derived from the return numbers of ground and vegetation according to the following expression: and γ is the spectral correction coefficient as the rate of the soil and vegetation backscattering coefficient; i.e., γ = ρ grd /ρ veg . The spectral correction coefficient can be estimated from hyperspectral images as with the assumption of Lambertian ground and randomly oriented Lambertian leaves [18,37], where R is the ratio of target reflectance.
The equation provides an opportunity to combine the LiDAR and hyperspectral data in this study. It should be noted that the equation is similar to the gap fraction calculation equation of Ni-Meister et al. [18], except that we come to this equation using the discrete LiDAR return count rather than the return power.
To calculate N grd and N veg , a weight was calculated for each return as 1/n, where n is the number of returns detected for the given pulse [38]: where v 1 , v 2 , v 3 , v 4 are the vegetation return numbers with 1, 2, 3, . . . n returns, and g 1 , g 2 , g 3 are the ground return numbers with 1, 2, 3, . . . n returns. One special case can be inferred from the equation: when λ = 1, i.e., ρ g = ρ v , then P gap = P lidar . This deduction implies that the LiDAR penetration ratio can be taken as a gap fraction only when the vegetation and soil have an identical back scatter coefficient [39].

LiDAR Extinction Coefficient
When the sensor scan zenith angle is known, the extinction coefficient is a function of the canopy leaf angle distribution [21], which is mainly controlled by plant species [23]. With the assumption that the forest in a certain area is comprised of a constant combination of tree species, the extinction coefficient for this forest is relatively stable compared with the variation of LAI [35]. Thus, we hypothesize that the forests in the research area share a stable but unknown leaf angle distribution but with varied LAI. Therefore, from a statistic point of view, the extinction coefficient in the tile level is one of the samples from the population of the whole research area, and so the averaged value of all the samples will be the best estimation of the extinction coefficient in this research area.
When the forest multi-directional gap fraction is calculated from LiDAR and spectral data in the equations, then the leaf area index and the parameter (χ) of the leaf angle distribution can be retrieved through a nonlinear optimization of the equations, with y minimizing the following cost function where θ i is the scan zenith angle at the ith bin, as illuminated in Figure 2, and P gap (θ i ) and P gap (θ i ) BLL are the gap fraction derived from LiDAR data and the Beer-Lambert law based model in the equation, respectively.
To reduce the search space and improve the probability of finding the optimized solution, we constrain the search space within the boundary [0.5, 2.5] for χ, which corresponds to a mean leaf angle from 30 to 70 degrees, and this leaf angle range can cover most of naturally developed tree types [38]. With a similar method, the boundary of LAI is defined as [0. 5, 9.0], and this range covers 85% of the LAI data in the global LAI dataset [40].
The above procedure was implemented in Matlab R2017 (The Mathworks, Massachusetts, USA) and the constrained nonlinear multivariable function with the sequential quadratic programming algorithm was selected to find the minimum of cost function (Equation (13)). The initial values for χ and LAI were set to the mid-points of the boundary constraints, i.e., 1.25 and 4.75 respectively, and we found the result was insensitive to the selection of initial values.
After the LiDAR extinction coefficient is estimated at tiled level, then the mean tilt angles for all the tiles can be calculated from the equation. The averaged leaf angle was taken as the result of the mean tilt angle in the research area.

Plot Level LiDAR LAI
The plot-level LAI is estimated using the LiDAR data corresponding to the ground plot location using the Equation (14): where P gap (θ) is calculated using the Equation (9), and θ is the averaged scan zenith angle in the plot level LiDAR data. The extinction coefficient is calculated using the mean value of the leaf angle distribution parameter χ which is derived from nonlinear optimization.

Evaluation
We calculate the root mean square error (RMSE) and relative RMSE (RRMSE) as a measure of the differences between values estimated by the Beer-Lambert law and the values actually observed by LAI-2000: where n is number of data points,ŷ i and y i are the estimated and ground-measured LAI, respectively, and y is the mean value of reference LAI. We also studied the correlations between errors from the fitted model using simple linear models. We calculate the R 2 as a statistical measure of how close the LiDAR-estimated LAIs are to the fitted line regression.

Statistical Features of Field Measurement
Among all the field point measurements, there were six outliers that had invalid mean tilt angle, so we did not include those points in the statistical results. As a result, there were 204 field LAI values available. The shape of the LAI histogram (Figure 5a) approximates to a normal distribution but with slight positive skewness of the mean value of 3.45, indicating that the forests in the research area are in the near-normal natural growing status, and there is no extreme external disturbance. The mean tilt angle and the vertical extinction coefficient share common features that are concentrated to the mean values (Figure 5b,c). These indicate that, for naturally growing trees, there is a predominant leaf angle distribution, and this feature results in a more concentrated distribution of extinction coefficient due to the nonlinear relationship between the mean tilt angle and leaf angle distribution. This finding is also demonstrated by the decreasing coefficient of variation (CV) of the LAI, mean tilt angle and extinction coefficient from 0.27 to 0.18 and 0.16. The result of the field LAI and mean tilt angle supports our assumption that there is stable leaf angle distribution with a less varied extinction coefficient and much more varied LAI in the study area.

Spectral Ratio of Soil and Vegetation
The pixels of hyperspectral data corresponding to the field plots were extracted and the soil and vegetation lines were identified on hyperspectral images. The calculation of the spectral ratio of soil and vegetation (r grd /r veg ) is illustrated in Figure 6. The value of r grd /r veg lies in a narrow range from 0.35 to 0.68, but the highest probability occurs at 0.53 (Figure 6a). The 25th and 75th percentiles are 0.50 and 0.63, respectively (Figure 6b). We calculated the mean value of the spectral ratio as 0.55.

LiDAR Gap Fraction
The nearly uniform LiDAR scan zenith angles lie in the range from -20 to 24 • (Figure 7a), but the frequency of the bins in the larger zenith angles was lower than the frequency of the bins in the smaller zenith angles. The total number of points in the zenith angles larger than 16 degrees is approximately one-quarter of all the points. The gap fraction is nonlinear to the forest LAI, so we calculated the statistic for the negative logarithm of the gap fraction (-lnP gap ) as the proxy variable to indicate the distribution of the forest LAI. As shown in Figure 7b, the histogram shape of -lnP gap is normally distributed with a mean of 1.55 and a CV of 0.29. The similar CVs of -lnP gap and field LAI (0.27 in Figure 5a) suggest that the data measured by LiDAR has an ability to describe the variance of forest LAI on a regional scale.

Tile-Level LiDAR Extinction Coefficient and LAI
The histogram of the LiDAR extinction coefficient in Figure 8a suggests that, in most cases, the LiDAR extinction coefficients fall into the range of field measured values, but have a slightly lower mean compared with the result in Figure 5c. The estimated LAI shows a normal distribution with a mean value of 3.36 and a CV of 0.32 (Figure 8b). A t-test demonstrated that there was no significant difference in mean LAI estimates using LiDAR (3.36, σ = 1.08) and LAI-2000 (3.45, σ = 0.94).

Plot-Level LiDAR LAI and Gap Fraction
The statistical results on the input parameters of calculating the plot-level LAI are shown in Table 2 and the plot-level LAI is shown in Figure 9.  At the plot level, the scan zenith angle of LiDAR data is limited to a small dynamic range, as shown in Table 2, which means that, at the plot level, it is difficult to invert the physical model to derive the plot-wise extinction coefficient. This is the main reason we inverted the model at the tile level and used the mean extinction coefficient obtained from tile-level LiDAR data. We compared the plot-level LiDAR LAI with the field measurement in Figure 9. The scatter plot (Figure 9a) shows that there is a significant positive correlation between LiDAR and ground results with a coefficient of determination (R 2 ) of 0.66 and an RMSE of 0.60.
In this study, the regression line displays good adherence to the 1:1 line, which indicates that there is no apparent bias between the LiDAR-based LAI and the field LAI, and it can also be manifested by the Q-Q plot in Figure 9b. A Q-Q plot is a probability plot, and if the two distributions are similar, the points will lie on or near the line y = x. As suggested by the Q-Q plot, we can conclude that the LAI datasets have similar distributions.

Performance of LiDAR LAI Estimation
In this study, we proposed a simple but valid method to estimate forest LAI from airborne discrete return LiDAR data. As described in the introduction, the estimation of LAI from LiDAR data can be based on empirical and physical methods. In the same research area, Majasalmi et al. [41] estimated LAI using an empirical method and the same dataset as in our work. They used the LiDAR penetration ratio without spectral correction, which was calculated using the method of Korhonen et al. [15] and Solberg et al. [42]. Their result showed that the empirical method produced LiDAR LAI estimates with a lower coefficient of determination of 0.56.
In our study, we explored the application of a physically-based method. For comparison purposes, we collected previously published results estimated by a physically-based method using discrete or waveform LiDAR data in other research areas and then compared the published work with our results (Table 3). Though there may be incompatible factors that will influence the algorithm, the comparison shown in Table 3 provides an overview of the performance of different algorithms. Among all the collected datasets, the LAI range in this study can represent the common situation in a natural forest, except for the two cases that have an apparently larger maximum LAI (9.60 in recorder 1# and 8.77 in recorder 5#). As shown in Table 3, our results have a relatively high determination coefficient R 2 of 0.66 and the highest accuracy, with a lower RMSE of 0.60 and the lowest RRMSE of 0.15.

Impact of Target Optical Property on LiDAR-Derived LAI
The variable describing the spectral properties of our method is the ratio of the vegetation and soil reflectivity ρ g /ρ v . This parameter directly determines the gap fraction in the equation and further influences LAI. For a randomly oriented leaf canopy, the ratio is a function of leaf reflectance and soil albedo, as indicated in the Equation (9) [37]. In the context where there is no prior knowledge of soil and leaf reflectance measurements, Ni-Meister et al. [37] suggested that a unity value should be used for ρ g /ρ v . Results from Armston et al. [17] showed that the stability of ρ g /ρ v depends on the footprint size of the LiDAR sensor and flight altitude. Chen et al. [19] reached a similar conclusion, but their investigation on the aggregated 5 m footprint size showed that there is a relatively stable value of ρ g /ρ v .
For the natural forest, on the fine scale, there may be a varied backscatter spectral property of forest and background. However, in a specific research area and in the short flight duration time, it is hypothesized that there is high probability that the spectral property retains stability and results in a fixed ρ g /ρ v value. The statistical results support our assumption that there is a relatively stable ρ g /ρ v value ( Figure 6).

On the Constraints on Nonlinear Optimization
We set the prior knowledge with the bounds to constrain the target parameters (χ and LAI) for the nonlinear optimization. A priori knowledge is useful when noise and poor angular sampling reduce the accuracy of model inversion given a limited number of observations [47]. Prior knowledge may be difficult to obtain for the specific research area. Fortunately, the retrieved result is still reasonable even with a wider range of prior knowledge. This implies that, in the context where there is no exact information on the target parameters, a relaxed constraint can produce acceptable results. In this study, we provide loosely constrained knowledge. From the point of view of information accumulation [47], a long-term field experiment and remote sensing observations may provide a useful source of the prior knowledge. There are large volumes of datasets available, such as the historical leaf angle distribution recorder [38], global leaf area index data from the field [40], the library of forest structure [48] and the recently reported canopy structure in Finland [49]. All the above observation data can be used to narrow the boundary of the knowledge and thus improve the probability of reaching an optimized solution for the targeted parameters.

Conclusions
Our research is one of the few attempts to derive LAI using discrete LiDAR data based on a physically-based model rather than empirical methods. This study has shown that forest LAI can be directly estimated from discrete LiDAR data using the spectral corrected LiDAR penetration ratio and a Beer-Lambert law-based physical model. The contribution of this work can be summarized as follows: (1) It is feasible to retrieve forest LAI from discrete LiDAR data without field data when the forest gap fraction and extinction coefficient can be appropriately calculated. (2) Angular gap fractions can be obtained from large tiles of LiDAR data, which usually have a much larger range of scan zenith angles than plot-level data. (3) For tiled data, the inversion of the Beer-Lambert law-based model provides a feasible method to retrieve tile-level LAI and extinction coefficients.
(4) Statistics for the tile-level extinction coefficient are valid for the plot-level LiDAR to estimate LAI corresponding to field measurements.
The proposed method is independent of the reference data; i.e., the method does not need a field-based training dataset of the LAI or gap fraction. This feature makes it feasible to be developed into an automatic method, which may be a promising solution to produce forest LAI from historical, archived discrete LiDAR datasets, because there was scarce forest LAI measurement when LiDAR data were sensed in the early period of discrete LiDAR development. For example, the Sustainable Landscapes Project has documented a huge volume of discrete LiDAR data over 10 years in the Americas [50]; however, only few ground LAI data are available. We believe that, with the help of the direct estimation of forest LAI, the historic remotely sensed LiDAR data will play an import role in supporting forest sustainable development.