LAI estimation based on physical model combining airborne LiDAR waveform and Sentinel-2 imagery

Leaf area index (LAI) is an important biophysical parameter of vegetation and serves as a significant indicator for assessing forest ecosystems. Multi-source remote sensing data enables large-scale and dynamic surface observations, providing effective data for quantifying various indices in forest and evaluating ecosystem changes. However, employing single-source remote sensing spectral or LiDAR waveform data poses limitations for LAI inversion, making the integration of multi-source remote sensing data a trend. Currently, the fusion of active and passive remote sensing data for LAI inversion primarily relies on empirical models, which are mainly constructed based on field measurements and do not provide a good explanation of the fusion mechanism. In this study, we aimed to estimate LAI based on physical model using both spectral imagery and LiDAR waveform, exploring whether data fusion improved the accuracy of LAI inversion. Specifically, based on the physical model geometric-optical and radiative transfer (GORT), a fusion strategy was designed for LAI inversion. To ensure inversion accuracy, we enhanced the data processing by introducing a constraint-based EM waveform decomposition method. Considering the spatial heterogeneity of canopy/ground reflectivity ratio in regional forests, calculation strategy was proposed to improve this parameter in inversion model. The results showed that the constraint-based EM waveform decomposition method improved the decomposition accuracy with an average 12% reduction in RMSE, yielding more accurate waveform energy parameters. The proposed calculation strategy for the canopy/ground reflectivity ratio, considering dynamic variation of parameter, effectively enhanced previous research that relied on a fixed value, thereby improving the inversion accuracy that increasing on the correlation by 5% to 10% and on R2 by 62.5% to 132.1%. Based on the inversion strategy we proposed, data fusion could effectively be used for LAI inversion. The inversion accuracy achieved using both spectral and LiDAR data (correlation=0.81, R2 = 0.65, RMSE=1.01) surpassed that of using spectral data or LiDAR alone. This study provides a new inversion strategy for large-scale and high-precision LAI inversion, supporting the field of LAI research.

Leaf area index (LAI) is an important biophysical parameter of vegetation and serves as a significant indicator for assessing forest ecosystems. Multi-source remote sensing data enables large-scale and dynamic surface observations, providing effective data for quantifying various indices in forest and evaluating ecosystem changes. However, employing single-source remote sensing spectral or LiDAR waveform data poses limitations for LAI inversion, making the integration of multi-source remote sensing data a trend. Currently, the fusion of active and passive remote sensing data for LAI inversion primarily relies on empirical models, which are mainly constructed based on field measurements and do not provide a good explanation of the fusion mechanism. In this study, we aimed to estimate LAI based on physical model using both spectral imagery and LiDAR waveform, exploring whether data fusion improved the accuracy of LAI inversion. Specifically, based on the physical model geometric-optical and radiative transfer (GORT), a fusion strategy was designed for LAI inversion. To ensure inversion accuracy, we enhanced the data processing by introducing a constraint-based EM waveform decomposition method. Considering the spatial heterogeneity of canopy/ground reflectivity ratio in regional forests, calculation strategy was proposed to improve this parameter in inversion model. The results showed that the constraint-based EM waveform decomposition method improved the decomposition accuracy with an average 12% reduction in RMSE, yielding more accurate waveform energy parameters. The proposed calculation strategy for the canopy/ground reflectivity ratio, considering dynamic variation of parameter, effectively enhanced previous research that relied on a fixed value, thereby improving the inversion accuracy that increasing on the correlation by 5% to 10% and on R 2 by 62.5% to 132.1%. Based on the inversion strategy we proposed, data fusion could effectively be used for LAI

Introduction
Leaf area index (LAI) is one of the prime determinants of photosynthesis, which makes it an important quantity controlling physical and biological processes of plant canopies and assessing forest growth potential (Chen and Black, 1992; Barclay and Goodman, 2000). As a fundamental attribute of global vegetation, LAI has been listed as an essential climate variable by the global climate change research community (Fang et al., 2019). The ability to accurately and rapidly acquire LAI is an indispensable component of process-based ecological research facilitating the understanding of gas-vegetation exchange phenomenon at an array of spatial scales from the leaf to the landscape (Zheng and Moskal, 2009).
Remote sensing data provides large-scale, systematic land surface observations consistently over the globe (Pan et al., 2008). With remote sensing technology, LAI can be mainly derived from a variety of sensors including passive optical sensors, active light detection and ranging (LiDAR) instrument, and microwave sensors (Fang et al., 2019). In passive optical sensors, multispectral and hyperspectral sensors provide spectral measurements across the electromagnetic spectrum, which are sensitive to subtle variations in reflected energy and, therefore, have a giant potential for detecting differences in vegetation (Mananze et al., 2018). However, LAI retrieval using multispectral and hyperspectral data has potential problems, such as the low signal-to-noise ratios (SNRs) of some remote sensing data, the "curse of dimensionality" and problems of saturation (Liu et al., 2016). In addition, optical remote sensing is only capable of capturing information from the horizontal canopy, resulting in a lack of information pertaining to vertical canopy (Xu et al., 2022). In active LiDAR instrument, full-waveform LiDAR systems can digitize the entire reflected energy, resulting in complete waveforms from the top of the canopy to the ground which reflect vertical profiles (Lefsky et al., 1999;Mallet and Bretar, 2009). It has been used to estimate LAI based on canopy structure and radiation transfer principles, especially by means of the correlation with the gap fraction (Wang and Fang, 2020). The primary advantage of LAI estimation using full-waveform LiDAR lies in its ability to capture detailed structural information beneath the canopy through the complete energy waveform, thereby mitigating estimation errors stemming from leaf aggregation. However, the high-density data obtained from airborne LiDAR is limited to the measurement range, whereas spaceborne full-waveform LiDAR data is characterized by its substantial volume. Multi-source remote sensing data have their own advantages and disadvantages in LAI estimation. On this basis, researches on using multi-source remote sensing data fusion to estimate LAI is becoming a research hotspot for which it can give full play to the advantages of different remote sensing data (Clevers and vanLeeuwen, 1996;Yang et al., 2011;Ma et al., 2014;Qu et al., 2015).
Existing methods for fusing spectral and LiDAR data are mostly based on empirical models. The empirical model directly relates inputs to outputs by pure statistical means, the advantage of which lies in its simplicity (Weiss et al., 2020). Thomas et al. (2011) estimated LAI with LiDAR and multispectral data by constructing fused LiDAR-optical indices. Ma et al. (2014) combined LiDAR data with the MODIS and MISR products to retrieve canopy height and LAI by multivariate linear regression model and geometricoptical mutual-shadowing (GOMS) model. Luo et al. (2019) estimated maize LAI using the combined hyperspectral imagery and LiDAR pseudo-waveforms by random forest (RF) regression algorithm. Zhang et al. (2022) estimated the LAI of a short-crop using UAS-based SfM and LiDAR point clouds, as well as the spectral information from multispectral imagery. Zhang et al. (2023) used six typical machine learning algorithms to construct prediction models of LAI, among which the XGBoost model showed the best performance. It also showed that the fusion of data could significantly improve the predictive ability of the models. However, the empirical model is strongly grounded with a large amount of statistical data and can only be applied within a relatively localized area because their performance is highly dependent on vegetation types, canopy structures, sensors, and temporal change . Empirical models, even superior deep learning models, can also suffer from statistical problems such as overfitting (Verrelst et al., 2015;Neinavaz et al., 2016).
In contrast to empirical methods, the physical model can be better generalized and its physical principles is helpful for the analysis of fusion mechanism (Myneni et al., 1997;Verrelst et al., 2019;Kennedy et al., 2020). Based on LiDAR waveform, commonly used models include gap fraction models (homogeneous canopy or clumping-aware canopy) and three-dimensional (3-D) radiative transfer models. Within multiple gap fraction models, gap fraction can be performed by directly computing its aggregate value (Luo et al., 2013;Fieber et al., 2014;Tseng et al., 2016;Jiang et al., 2022) or by considering the vertical accumulation of gaps within the canopy layers (Ni et al., 1999;Yang et al., 2019). Gap fraction models, with minimal input parameters, streamline LAI computation through forward modeling . The geometric-optical and radiative transfer (GORT) model, one of the gap fraction models, capitalizes on the capability of lidar waveform data to characterize the underlying canopy structure (Ni-Meister et al., 2001). Accurate vertical profiles of LAI can be derived by GORT (Tang et al., 2012). This methodology has been employed in the derivation of products for the spaceborne lidar GEDI and exhibits a strong applicability of large-scale LAI estimation (Tang et al., 2014;Wang et al., 2023). Regarding the 3-D radiative transfer models such as Discrete Anisotropic Radiative Transfer (DART), they are proficient in simulating lidar waveforms effectively (Gastellu-Etchegorry et al., 2016;Gastellu-Etchegorry et al., 2017). However, due to its multitude of model parameters and complex scenarios, time-consuming aspects arise during the inversion process. Taking the above reasons into consideration, we consider conducting a data fusion inversion study based on the GORT model. It is still a challenge to combine spectral imagery and LiDAR waveform for LAI retrieval based on physical models. In LAI estimation, spectral features including sensitive bands' reflectance and spectral indices, play a pivotal role in accuracy (Potithep et al., 2013;Liang et al., 2015). Waveform information such as height, echo energy ratio and leaf coverage are key parameters of LAI estimation using full-waveform LiDAR (Pope and Treitz, 2013;Ma et al., 2015). Improving the physical model to use the above parameters so as to make use of the respective advantages of multi-source data is the key for LAI estimation using both spectral and LiDAR data. In addition, the accuracy of the physical model is susceptible to the initial assignment value of model parameters and the quality of the input data (Houborg et al., 2007). Adjusting the input parameters of the model based on the study area and source data is also an important measure to ensure the accuracy of LAI estimation.
In response to the above problems, the main objectives of this paper are: (1) Developing an estimation strategy based on GORT model to achieve data fusion. (2) Extracting accurate parameters for LAI estimation from both spectral and LiDAR waveform data. (3) Assessing the performance of the joint data fusion for LAI estimation. Specifically, for the first time, we attempt to fuse spectral and waveform data within the GORT model, achieving a joint LAI estimation. We enhance the waveform decomposition method to extract more accurate waveform parameters. Furthermore, on the basis of the existing retrieval, we improve the model input parameter canopy/ground reflectivity ratio (r v =r g ) to make it more suitable for the large-scale forest with the heterogeneity of the spectrum. The significance of this study lies in providing novel insights into the fusion of active and passive remote sensing data, thereby contributing to the enhancement of accurate large-scale LAI estimation.

Study area and data
The study area is located in Harvard Forest, a 4000 acres forest in Petersham, Mass., which is now among the most studied forests in the world. As a critical node in USA national ecological network (LTER, NEON and ForestGEO), Harvard Forest department gathers and produces various datasets from its ecological scientific researches. At the same time, many remote sensing ecosystem projects choose Harvard Forest as study area, which provide multi-source remote sensing data. For these reasons, we chose the Harvard Forest to carry out remote sensing research of multisource data.

LVIS airborne LiDAR data
NASA's Land, Vegetation, and Ice Sensor (LVIS), is an airborne, wide-swath, full waveform imaging laser altimeter system, which emits 1064nm wavelength laser pulses to collect data on surface topography and 3-d structure with medium 25m footprint.
During summer 2021, LVIS operated as a NASA Facility to calibrate and validate the space-based LiDAR sensor GEDI (Global Ecosystem Dynamics Investigation) by conducting overflow ground tracks over the Eastern United States and French Guiana. The LVIS Classic instrument was flown on Gulfstream V at a flight altitude of 41,000', covering Harvard Forest completely on August 6, 2021. The data products of LVIS include Level 1B Geolocated LVIS Waveforms (HDF format) and Level 2 Geolocated Surface Elevation and Height Product (ASCII Text format), from which ecosystem structure parameters can be derived (Blair et al., 1999;Blair and Hofton., 2020).
In this research, we used data products of the LVIS flight on August 6, 2021, obtained from 'https://nsidc.org/data/LVISC1B/ versions/1'. From these, we extracted multiple parameters for each pulse into a comprehensive dataset (.csv). The parameters included: laser shot (shotnumber), longitude (lon), latitude (lat), elevation of the highest detected signal (zt), elevation of the lowest detected mode within the waveform (zg), return waveform (rxwave), signal mean noise level (sigmean).

Sentinel-2 multispectral images
The Copernicus Sentinel-2 mission comprises a constellation of two polar-orbiting satellites. It offers free multi-spectral images with high spatial resolution (four bands at 10 m, six bands at 20 m and three bands at 60 m spatial resolution). The orbital swath width is 290 km with high revisit time (5 days with 2 satellites under cloudfree conditions which results in 2-3 days at mid-latitudes), which support to accurately monitor land surface changes especially vegetation changes and are beneficial to biophysical indicators estimation (Drusch et al., 2012). The Sentinel-2 data has good temporal and spatial resolution with high-quality. Numerous studies have shown that it provides accurate retrieval of LAI. (Korhonen et al., 2017;Hu et al., 2020;Zhou et al., 2020;Sun et al., 2022).
Sentinel Applications Platform (SNAP), released by European Space Agency (ESA), has been accelerating Earth observation innovation since 2014. It can help to process and analyze Sentinel-2 imagery. In this experiment, the Sentinel-2 image was processed by the Sentinel-2 Toolbox (S2TBX) in SNAP. We acquired an L2A image of Sentinel-2 covering the Harvard Forest on July 31, 2021 (https://doi.org/10.5270/S2_-znk9xsj). The base map of Figure 1 is the Sentinel-2 RGB image of the study area processed by SNAP. Based on SNAP, we obtained canopy cover map from Sentinel-2 L2A image. Subsequently, using the latitude and longitude information of the laser pulse from the LVIS, we extracted parameters at each pulse point, including spectral reflectance and canopy cover (CC) values.

Ground based LAI in Harvard Forest
The Harvard Forest Data Archive contains various datasets from scientific research at the Harvard Forest, in which HF150 dataset collects Leaf Area Index at Harvard Forest HEM and LPH Towers since 1998 (https://harvardforest.fas.harvard.edu/harvardforest-data-archive). Leaf area index is measured with the LAI-2000 canopy analyzer with one LAI sensor made at multiple plots within each forest type -usually 12 plots within the old-growth hemlock forest, and 36 plots on Little Prospect Hill (Orwig and Hadley, 2022). The time, distance from tower, compass direction from tower from geographic north, LAI value of each plot is given in the dataset, which help to correlate exact coordinates and values of plots. The spatial distribution of ground plots is shown in Figure 1.

Methods
To better utilize the advantages of spectral imagery and LiDAR waveform, and considering the scale difference between satellite and airborne data, we first designed a LAI estimation strategy based on a physical model GORT, using airborne LiDAR waveform as the main body and satellite spectral data as support. In order to improve the LAI inversion accuracy, we optimized the model input data and parameters: optimizing waveform decomposition algorithm for more precise waveform energy data, improving the method for obtaining the canopy/ground reflectivity ratio as model parameter to obtain values more consistent with the actual research area. Figure 2 shows an overview of the methods used in this paper.

A fusion strategy proposal based on GORT model deconstruction
In order to perform joint estimation of LAI using both spaceborne multispectral images and airborne LiDAR waveforms, we decomposed the GORT model and developed a data fusion strategy.
The geometric-optical and radiative transfer (GORT) model is for the bidirectional reflectance of a vegetation cover combines principles of geometric optics and radiative transfer (Li et al., 1995). Ni-Meister et al. (2001) developed a method based on the modified The spatial distribution of ground plots in Harvard Forest (The base map is Sentinel-2 RGB image.). Shi et al. 10.3389/fpls.2023.1237988 Frontiers in Plant Science frontiersin.org GORT model to derive gap probability and canopy cover from LiDAR waveforms. In the modified GORT model, the cumulative canopy height profile (CHP) was calculated by using a logarithmic transformation of (1-canopy cover). The canopy cover and gap probability could be calculated as follows (Ni-Meister et al., 2001): In formula (1), P gap (h) and fcoVer(h) represented the gap probability and canopy cover percentage above a particular height h within canopy respectively. The terms R V (h), R V (0) and R g were the integrated laser energy returns from the canopy top to height h, from canopy top to canopy bottom, and from the ground return individually. The canopy and ground reflectance were r v and r g respectively.
Canopy cover (CC) is typically defined as the extent of ground area covered by the foliage of trees or other vegetation, as projected from a vertical viewpoint onto a horizontal plane (Fiala et al., 2006). From the perspective of remote sensing approaches, two types of canopy cover estimates are commonly derived: metrics describing the 2D horizontal extent of canopy, which is often expressed for a given cover type as a percentage of pixels (Silvan-Cardenas and Wang, 2010); or as 3D LiDAR metrics that represent the transmission of light through the canopy (Morsdorf et al., 2006;Korhonen et al., 2011;Moran et al., 2020). Researches have demonstrated a strong correlation between canopy cover extracted from spectral data and LiDAR-derived inversion values, with the latter frequently exhibiting higher accuracy (Smith et al., 2009;Ma et al., 2017). Figure 3 shows a schematic diagram of the principles of acquiring forest canopy information using active and passive remote sensing sensors.
Based on the above researches, we assumed that the canopy cover obtained by spectral imagery and LiDAR waveform are equal in this study: where fcoVer(z) represented canopy cover percentage above height z. We set the accumulated canopy cover percentage above height z as equal to the estimated canopy cover (FVC) in the pixel at the position of the LiDAR pulse using spectral imagery. z was set to the height corresponding to 80% of the canopy energy in the return. Therefore, the estimate of r v =r g could be expressed as: Based on the above formulas, the canopy/ground reflectivity ratio value of each pulse could be calculated according to LiDAR waveform and spectral imagery. Tang et al. (2012) deduced the formula of LAI derivation based on the GORT model. The effectiveness of the method has been proved by experiments. Total LAI can be calculated as (Tang et al., 2012): The overview of the methods used in this paper. Shi et al. 10.3389/fpls.2023.1237988 Frontiers in Plant Science frontiersin.org where C represented the clumping index which adjusted the linear relationship between effective LAI and true LAI (Chen, 1996). We chose the clumping index value of 1.58 for the in Harvard Forest (Tang et al., 2012). G was the projection coefficient and was set to be 0.5 assuming a random foliage distribution within the canopy (Ni-Meister et al., 2001). R V (0) and R g were the integrated laser energy returns from canopy top to canopy bottom, and from the ground return individually, which were obtained by waveform processing. r v =r g was calculated as formula (3) using both spectral data and LiDAR waveforms. Based on formula 3 and 4, LAI was estimated by combining spaceborne and airborne data. The strength of this strategy is that we used the physical model to realize the data fusion, and improve the parameter r v =r g by using the spaceborne spectral reflectance.

Optimization of waveform decomposition method for accurate waveform energy parameters extraction
Waveform energy parameters are the main input data in LAI inversion based on the GORT model, including canopy backward energy, ground backward energy and waveform energy integral returned at different altitudes, which are useful for segmentation, classification and inversion purposes, in both forested and urban areas (Mallet and Bretar, 2009). Selecting an appropriate processing method to "purify" the original waveform data is vital to extract structure parameters of forest, so as to accurately invert physical parameters of vegetation. In this paper, waveform processing procedure was designed for the processing of the Level 1B Geolocated LVIS Waveforms and Level 2 Geolocated Surface Elevation and Height Product. Specifically, to extract energy parameters more accurately from waveforms, the waveform decomposition algorithm was improved. The processing flow chart is shown in Figure 4.
Firstly, the background noise of the echo was removed based on the average noise parameter "sigmean" of each echo calculated in flight provided by LVIS Level 1B product. Values less than the noise threshold was eliminated.
Secondly, the Gaussian filtering algorithm was used to remove other types of noise and smoothed the waveform. The methods of waveform denoising mainly include Gaussian filtering, mean filtering, Fourier low-pass filtering, etc. . Gaussian filtering has small time-frequency window area and a simple design, which makes it widely used in the field of signal processing. By measuring the echo denoising effect and adjusting the parameters, a Gaussian filter with better denoising effect was finally selected for echo denoising.
Finally, waveform decomposition method was applied to decompose the echo and extracted effective waveform parameters. Since the backscattered echo signal can be considered as the superposition of multiple Gaussian signals, the Gaussian decomposition method was used to fit the original signal to the superposition of multiple Gaussian function curves (Zhou et al., 2022). The backscattered echo can be expressed as: In formula (5), W(t) is the amplitude of the waveform at time t; e is the bias of the Gaussian waveform; N p is the number of Gaussian components; A m , t m , s m are the amplitude, peak position and waveform width of the waveform of the mth Gaussian component respectively.
There are two main steps in the waveform decomposition: 1) estimation of the initial parameters; 2) optimization of the parameters and fitting the waveform (Zhou et al., 2021). After extracting accurate initial values of parameters, the commonly used waveform fitting methods include LM (Levenberg-Marquardt optimization algorithm) method (Wagner et al., 2006) and EM (Expectation-Maximization algorithm) method (Persson et al., 2005), the accuracy of which has no significant difference (Zhou et al., 2022). In the global optimization algorithm, waveform components that are too close in distance are prone to being merged during optimization, leading to significant fitting errors. In this paper, in order to address recognition errors caused by some echo components being close to the target, constraints were placed The schematic diagram illustrating the principle of acquiring forest canopy information through active and passive remote sensing methods. Shi et al. 10.3389/fpls.2023.1237988 Frontiers in Plant Science frontiersin.org on the values of each peak to improve the accuracy of the fitting based on EM algorithm. The specific approach of the constraintbased EM was designed as follows: E-step: Compute the posterior probabilities for each component given the data points using Bayes' rule. The posterior probability for the jth component of the Gaussian mixture model for the ith data point was given by: where g ij represented the probability that the ith data point belonged to the jth component, p j (x i ) was the probability density function of a Gaussian distribution at x i of the j th component. M-step: Update the parameters of the Gaussian mixture model using the posterior probabilities computed in the E-step. The update equations for the means (m j ), and standard deviations (s j ) of the jth component were given by: where P ½a,b (x) was the projection operator that maps x onto the interval ½a, b, m j 0 , s j 0 were the initial values of the mean and standard deviation for the jth component.
Repeat E-step, M-step until convergence was achieved. The constrained EM algorithm for LiDAR waveform decomposition imposed constraints on the parameter range within the optimization problem, leading to enhanced stability and accuracy of the algorithm.
After waveform processing, we identified the last waveform component as the ground component, and the rest as the canopy components. R g was the area enclosed by the amplitude of the last waveform component and the coordinate axis within its start-stop The flow chart of waveform processing (R v (z), R v (0) and R g are the integrated laser energy returns from the canopy top to height z, from canopy top to canopy bottom, and from the ground return individually.).
range. R V (0) was the integral value of the other waveform components within their start-stop range. R V (z) was the integral value of the waveform from the initial position of canopy component to the height of z.
In this section, we implemented the processing of LiDAR waveforms, especially by adding constraints on the peak positions in the waveform decomposition algorithm to obtain better waveform decomposition results and calculate more accurate waveform energy parameters.

Optimization of model parameter r v =r g for large scale forest by gridding study area
When using the strategy proposed in Section 3.1 to calculate the LAI of each pulse position, both spectral and LiDAR data are used to calculate the model parameter r v =r g at this point, that is, r v =r g varies with the input data in each position. Although the parameter setting strategy is more accurate than taking a fixed value for the entire study area, it causes large data uncertainty and increases the computation of inversion. For example, abnormal waveform of a pulse will result in abnormal calculated value of r v =r g , thus leading to abnormal LAI inversion results at this point.
To solve this problem, we proposed a method to optimize the model parameter r v =r g , that was, gridded the study area and calculated the mean value of r v =r g in each grid. Then the LAI estimation model was constructed by each grid. In a large area, the canopy/ground reflectivity ratio varies with forest environmental conditions. Gridding the study area not only accounts for the spatial heterogeneity within Harvard Forest but also reduces the computational complexity of LAI modeling, thereby mitigating uncertainty arising from anomalous input data. Theoretically, the optimization method is helpful to improve the inversion accuracy.
Specifically, we divided the study area into 2331 rectangular grids (63*37), the coordinate size of each grid was 0.002°* 0.002°( about 36118m 2 ). In order to reduce the uncertainty, the statistical method of histogram distribution was used to eliminate extreme abnormal ratio values that without the range of the mean plus or minus three standard deviations in each grid. Then an average reflectivity ratio was assigned for all the pulses in grid to reduce the uncertainty caused by pulse quality. The reflectivity ratio of each grid was the trimmed mean of all pulses in the grid.

Waveform energy parameters based on optimized waveform decomposition method
For the purpose of obtaining accurate canopy and ground energy parameters, we processed a total of 187,302 LiDAR waveforms in the study area through denoising, smoothing, and waveform decomposition based on least squares optimization. During this process, the waveform energy parameters Rg and RV obtained using different least squares optimization methods were compared.
The correspondence between the original waveforms and Gaussian model estimates of Rg and RV using different optimization algorithms is shown in Table 1. Overall the correspondence is high for both RV and Rg, the R 2 of which are all above 0.99. Further comparative analysis between the canopy and ground components reveal that the overall canopy R 2 is lower than that of the ground. The RMSE for the canopy is significantly higher than that of the ground, with a maximum difference of 35.656. It can be attributed to two main factors. Firstly, the canopy component often has one or more echo components in LiDAR waveforms, whereas the ground component typically only has one. Multiple echo components result in larger errors in waveform decomposition, leading to poorer correspondence of RV. Secondly, due to the significant differences in canopy coverage, the waveform energy fluctuation of canopy is also much larger, resulting in a much higher RMSE compared to the ground.
Comparing the results obtained from two different optimization methods, it is found that correspondence is higher for both Rg and RV when the bound of peaks are constrained, compared to the result without parameter constraints (higher R 2 and lower RMSE). To investigate the reasons for the result, we compare the waveform fitting results of two different optimization methods. The results of waveform decomposition using unconstrained optimization algorithm and optimization algorithm with peak boundary constraints are shown in the Table 2. We conducted an average statistical analysis of the fitting results for all waveforms in the study area and found that the constrained optimization algorithm yielded a fitting waveform with R 2 of 0.979, MAE of 6.907, MSE of 103.016, and RMSE of 8.596. The fact that the average R 2 is almost close to 1  Waveform fitting and waveform energy parameters results using two different optimization methods. (A) with two canopy echoes and one ground echo; (B) with two canopy echoes and one ground echo; (C) with three canopy echoes and one ground echo (The left column: the original Expectation-Maximization algorithm, the right column: the Expectation-Maximization algorithm with peaks boundary constraints).
R 2 is not significant. This suggests that the overall fitting accuracy of the unconstrained method is lower than that of the constrained method.
To further analyze the waveform decomposition result using different methods and its effects on the value of Rg and RV, we extracted some waveforms with significant differences in inversion accuracy using two different optimization methods ( Figure 5). The left column shows results obtained by EM method, while the right column shows results obtained by constraint-based EM method. The waveform energy parameters result of the three waveforms are shown in Table 3. Waveform (a) contains two canopy echoes and one ground echo. The waveform decomposition using EM method identifies only one canopy echo, and the calculated energy of the canopy echo is 2179.58, which differs significantly from the actual value of 2166.87. The constraintbased EM method identifies two canopy echoes, and the calculated value of 2166.21 is almost the same as the actual value. Waveform (b) contains two canopy echoes and one ground echo, and the second canopy echo and ground echo are combined into one waveform component using the EM method, resulting in a large difference between the canopy and ground energy calculation results of 1510.62 and 805.17 and the actual 1562.44 and 775.86, while the added constraint method does not show this phenomenon. Waveform (c) contains three canopy echoes and one ground echo, which are partially combined by the EM method due to the close distance of each canopy echo, resulting in a large difference between the calculated canopy energy (1939.33) and the actual value (1945.25). The results reveal that the waveform energy parameters calculated by constraint-based EM method are more precise than the unconstrained EM method. The waveform diagram shows that the main reason for the difference in accuracy between the two waveform decomposition methods is that after constraining the peak position of each waveform component, the merge of the closed waveform echoes can be avoided.
These results provide evidence that compared to optimization algorithms without parameter boundary constraints, the proposed The heat map of gridded canopy/ground reflectivity ratios in the study area (left) (The right image shows the Sentinel-2 RGB image of the study area).

The gridded r v =r g result
After gridding the study area, the r v =r g value of each grid is shown in Figure 6, and the statistical result is shown in the Table 4. Comparing the value heatmap ( Figure 6 left) with the RGB images ( Figure 6 right) of the study area, it can be found that the r v =r g is correlated with the degree of vegetation coverage. The value of r v =r g in non-vegetated and sparsely vegetated areas is generally lower than that of areas with higher vegetation coverage. The spatial distribution of r v =r g coincided with the actual vegetation distribution, indicating the accuracy of the calculation method we proposed. The statistical results show that the final ratio of 2331 grids in the study area is within the range of [0, 3.68], the average value is 2.45, and the root mean square is 0.80.
In previous studies, Lefsky et al. (1999) suggested using a constant (r v =r g = 2) for 1064 nm. Tang et al. (2012) obtained the r v =r g value of 2.5 at 1064 nm and used it as the mean value for the whole study area. In addition to determining the ratio by empirical field measurements, extracting from LiDAR waveforms using statistical methods (Ni-Meister et al., 2010;Armston et al., 2013) is also a way to calculate this ratio. The value of canopy/ground reflectivity ratio is basically between [0, 3]. Due to the lack of measured data, and no research has used the method of combining spectral and LiDAR data to calculate r v =r g , it is currently impossible to accurately demonstrate the accuracy of the calculated ratios. However, the average and mean square deviation results show that most of the reflectance ratios are [0, 3] with only a few abnormal values, it can be proved that accurate gridded r v =r g of the study area can be obtained by this strategy. These provide a new idea for the calculation of r v =r g .

Comparison of LAI inversion results based on different data and inversion methods
In this section, we explored the LAI retrieval results based on different datasets and different methods for calculating r v =r g . The accuracy is compared with the existing field measurement data.
Using the fusion strategy we proposed, we obtained the LAI map using LiDAR waveform and multispectral image (Figure 7). The result shows that directly using the data of entire study area without land cover classification can estimate LAI well. For example, LAI values of the road in the lower left of this research area are 0. We can clearly distinguish the vegetation and nonvegetation areas from the LAI map. In the vegetated area, the LAI values are generally around 3-7, and reach above 7 in a few dense areas, indicating that the area is relatively heavily forested, which is basically consistent with the actual situation. Comparing the inversion results with the true LAI value provided by Harvard Forest HEM plots, the correlation, R 2 and RMSE are 0.81, 0.65, 1.01 respectively (Figure 8), which shows that the LAI map obtained have high accuracy. The method we proposed to invert LAI by fusing LiDAR and spectral data is feasible.
We conducted LAI estimation based on four different strategies for comparative analysis. These strategies encompassed two that exclusively utilized LiDAR data, one that solely relied on spectral data, and one that integrated both LiDAR and spectral data. Using only LiDAR data, we reproduced the methods used by Tang et al. (2012) and Armston et al. (2013) respectively. They both performed the inversion based on the GORT model, only some of the parameters in the model were determined in different ways. The LAI estimation maps ( Figures 9A, B) were performed according to the parameters they set. Using only spectral data, we adopted the most traditional empirical model to construct the linear LAI map of the study area using both spectral and LiDAR data. Shi et al. 10.3389/fpls.2023.1237988 relationship between NDVI-LAI for LAI inversion using the LAI true value of the Harvard Forest LPH flux tower's 36 plots ( Figure 9C). Using LiDAR waveform and spectral data, we employed the method proposed by Yang et al. (2019) and the obtained result was depicted in Figure 9D. Based on the true LAI value of Harvard Forest, the accuracy evaluation of the estimation results obtained by various methods was carried out. The overall accuracy results are shown in Table 5. Figures 9A, B indicates the LAI inversion results using only LiDAR data. Similar to LAI map obtained by the method proposed in this study, the results show obvious differentiation between vegetation and non-vegetation areas. However, the difference lies in that, LAI values of invention are too high in places with dense vegetation only using LiDAR data. By Tang et al. (2012) 's method, the LAI inversion result shows the correlation of 0.63, R 2 of 0.40 and RMSE of 2.01 with ground plot LAI, indicating a moderate correlation between the two. By Armston et al. (2013) 's method, the inverted LAI is weakly correlated with the true value, with a correlation coefficient of 0.53, R 2 of 0.28 and RMSE of 2.85. It shows that the accuracy of the two LiDAR-only inversion methods is lower than that of the fusion of LiDAR and spectral data. The LAI map inverted by the empirical model is shown in Figure 9C. The results show that there is a large difference between the LAI inversion results obtained by using only spectral data and those obtained by fusing the two data. Based only on spectral data, the inverted LAI is an underestimate with the highest value being only 4.8. It shows a correlation of 0.48, R 2 of 0.25 and RMSE of 2.72 with ground plot LAI, which is poor compared to the inversion result that combines the two data. Figure 9D illustrates a low accuracy in LAI estimation (the correlation of 0.51, R 2 of 0.30 and RMSE of 2.76), with LAI values consistently underestimated. This discrepancy may be attributed to the unsuitability of the model for the tree species and data sources in this study area. The original LiDAR waveform used by Yang et al. (2019) is acquired from a large-footprint LiDAR system (70m), which significantly differs from the experimental data used in this study (25m). Comparing various estimation strategies, it is evident that the fusion of both active and passive remote sensing data contributes to improved LAI estimation accuracy. The enhancements we have proposed for the GORT model further enhance LAI estimation precision.
We successfully estimate LAI based on the GORT model combining LiDAR and spectral data with a correlation of 0.81 and R 2 of 0.65, which shows a large accuracy improvement compared to both LiDAR data alone and spectral data alone. These improvements can be attributed to the addition of spectral data to improve the parameters of the model. Comparing the three inversion methods based on the GORT model, the main difference between them is the r v =r g value of the model, which may be the main reason for the difference in inversion accuracy. The r v =r g of the study area is determined as a constant value of 2.5 by experience, which helps to reduce the amount of computation. But this empirical value may not necessarily be applicable to Harvard Forest and a fixed value cannot adequately represent the forest conditions of the entire study area, which may be the prime causes of the low accuracy. The use of least squares methods provides a new approach for the calculation of the ratio, which does not rely on manual measurements, but rather on the energy returned by the LiDAR. Based on this approach, we obtain r v =r g of 1.17 for the study area. The main reasons for the low accuracy can be attributed to two factors. Firstly, similar to using experienced value, calculating a single value for the entire Harvard Forest will cause abnormal results due to the complexity of the forest canopy. Secondly, the quality of the LiDAR can greatly affect the results. Such as areas in low point density in canopy that do not reflect enough energy (Chauve et al., 2009). Combining the spaceborne spectral data and airborne LiDAR data to calculate the reflectance ratio can use highquality spectral data to a certain extent to eliminate the abnormal phenomenon of reflectance ratio caused by the abnormal collection of some LiDAR footprints. Compared with the above two methods, the gridded r v =r g calculation method we proposed considers the influence of these factors. By combining spaceborne spectral data and airborne LiDAR waveform to calculate the r v =r g , the influence of abnormal waveforms on the value of r v =r g can be eliminated to some extent. At the same time, laborious ground truth measurements of reflectance are no longer needed (Yang et al., 2006). Furthermore, dividing the study area into grids and calculating the average r v =r g in each grid can not only further eliminate abnormal values through statistical method, but also calculate different r v =r g in view of the canopy heterogeneity in large-scale complex forest. It can be found from the LAI maps of the three methods that the method of taking a constant value would lead to higher LAI values in the areas of dense vegetation, which is due to the fact that the actual r v =r g in these areas are higher than the determined value of the model. LAI results obtained by the method we proposed are basically within 7, with only a few outliers. The results also confirm the validity of the proposed method.
In addition to the inversion based on the GORT physical model, the LAI inversion based on the empirical NDVI-LAI relationship (Turner et al., 1999) is also carried out based on the spectral data. The Scatterplots of field-observed LAI against estimated LAI using both spectral and LiDAR data.
results show that the inversion accuracy of the empirical model using spectral data (correlation of 0.48 and R 2 of 0.25) is much lower than that of the physical model using LiDAR waveforms. The empirical model needs a certain amount of truth values to ensure the accuracy of the inversion equation. However, there are only 48 small plots in the study area, which in theory can cause large errors when used for empirical model construction and verification. Also, the existing plots of Harvard Forest HEM and LPH Towers are concentrated in a small area, which cannot well represent the NDVI-LAI relationship of the entire study area. It is the main reason for the lower precision., Our

Limitations
While this study has successfully achieved the joint inversion of spectral and LiDAR data for LAI estimation based on a physical model, there are still some limitations. Firstly, due to the limited number of ground measurement points in our study area, we are unable to achieve the fusion of LiDAR and spectral data for LAI retrieval based on empirical models, and compare it with the results from physical models. Nevertheless, this also partly demonstrates the advantages of developing data fusion inversion based on physical models, which helps us reduce reliance on ground measurements, lower manual labor costs, and facilitates the widespread application of large-scale regions. Secondly, while both LiDAR data and spectral data are commonly employed for retrieving canopy cover, there remains a disparity between the values obtained through these two methods (Smith et al., 2009;Li et al., 2023). This discrepancy, though overlooked in this experiment, may cause errors in inversion. It might be one of the contributing factors to the limited precision in LAI inversion. Research is warranted in future experiments to address this issue and enhance the inversion accuracy. Additionally, due to time and resource constraints, we do not validate the effectiveness of this method in regional scale. In future research, we will further explore the contribution of data fusion to LAI based on theoretical analysis.

Conclusion
Spectral imagery and full-waveform LiDAR data can provide reflectance information and echo energy information reflecting the vertical structure of the forest canopy respectively. Joint active and passive remote sensing data has great potential for accurate inversion of forest canopy LAI. Our research is one of the few attempts to derive LAI using both spectral imagery and LiDAR waveform based on physical model retrieval rather than through empirical methods. We proposed a useful data-joint LAI inversion strategy based on the GORT model using LiDAR waveform and spectral data. For the large-scale heterogeneous forest, we further accurately extracted the waveform energy parameters as the model input data and optimized the model input parameter canopy/ ground reflectivity ratio to improve the inversion accuracy. The results show that comparing with only using LiDAR or spectral imagery, the LAI calculated by the proposed strategy using both LiDAR waveform and spectral imagery has a higher accuracy, indicating the effectiveness of the proposed strategy. Overall, our study confirms that optimizing the input parameter and data of the model for the study area can help improve the inversion accuracy, and the combined LiDAR waveform and multispectral imagery have potential for improving prediction accuracies of LAI.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.