Retrieval of Vertical Foliage Proﬁle and Leaf Area Index Using Transmitted Energy Information Derived from ICESat GLAS Data

: The vertical foliage proﬁle (VFP) and leaf area index (LAI) are critical descriptors in terrestrial ecosystem modeling. Although light detection and ranging (lidar) observations have been proven to have potential for deriving the VFP and LAI, existing methods depend only on the received waveform information and are sensitive to additional input parameters, such as the ratio of canopy to ground reﬂectance. In this study, we proposed a new method for retrieving forest VFP and LAI from Ice, Cloud and land Elevation Satellite (ICESat) Geoscience Laser Altimeter System (GLAS) data over two sites similar in their biophysical parameters. Our method utilized the information from not only the interaction between the laser and the forest but also the sensor conﬁguration, which brought the beneﬁt that our method was free from an empirical input parameter. Speciﬁcally, we ﬁrst derived the transmitted energy proﬁle (TEP) through the lidar 1-D radiative transfer model. Then, the obtained TEP was utilized to calculate the vertical gap distribution. Finally, the vertical gap distribution was taken as input to derive the VFP based on the Beer–Lambert law, and the LAI was calculated by integrating the VFP. Extensive validations of our method were carried out based on the discrete anisotropic radiative transfer (DART) simulation data, ground-based measurements, and the Moderate Resolution Imaging Spectroradiometer (MODIS) LAI product. The validation based on the DART simulation data showed that our method could e ﬀ ectively characterize the VFP and LAI under various canopy architecture scenarios, including homogeneous turbid and discrete individual-tree scenes. The ground-based validation also proved the feasibility of our method: the VFP retrieved from the GLAS data showed a similar trend with the foliage distribution density in the GLAS footprints; the GLAS LAI had a high correlation with the ﬁeld measurements, with a determination coe ﬃ cient


Introduction
Forest canopies play a vital role in determining the energy cycle, regulating the climate, and providing habitats for various species; most of these functions are affected by foliage respiration and photosynthesis [1][2][3]. Therefore, it is essential to quantitatively describe the foliage in the canopy. The leaf area index (LAI) and vertical foliage profile (VFP) are two key canopy structure parameters that quantify the distribution of foliage in the horizontal and vertical directions, respectively. Accurate estimation of the LAI and VFP is essential for promoting understanding of the role of forest canopies in terrestrial ecosystems [4][5][6].
Due to the importance of the LAI and VFP, various efforts have been made to estimate them from a variety of data sources. Passive remote-sensing data are some of the broadest data sources for LAI estimation because of their easy availability and global observation capabilities [7][8][9][10][11]. However, the LAI from passive remote-sensing data suffers from some limitations. For example, the LAI tends to saturate in dense forests because of the limited penetrability of sunlight [12][13][14]. In addition, passive remote-sensing data do not include detailed vertical measurement information, which limits its application areas, such as VFP inversion. In contrast, measurements from light detection and ranging systems (lidar) have good penetration ability and can provide detailed vertical observations from the top of the canopy to the ground. The above reasons promote the application of lidar observations to forest structure parameters [15][16][17], especially the LAI and VFP estimations [18][19][20][21][22][23][24][25][26][27][28]. In general, lidar systems can be classified into three types depending on the platform, including terrestrial, airborne, and spaceborne lidar. Due to the limited spatial coverage, terrestrial and airborne lidar are only suitable for small-or middle-scale research. The spaceborne lidar has global observation capabilities, which supports its utility as a data source for global-scale research.
Existing spaceborne lidar systems include the Geoscience Laser Altimeter System (GLAS) on board the Ice, Cloud and land Elevation Satellite (ICESat), ICESat-2, and Global Ecosystems Dynamics Investigation (GEDI), of which the ICESat GLAS is the most widely used spaceborne lidar system because it is the first spaceborne lidar for Earth observation. More importantly, the office of ICESat GLAS can provide free and rich datasets corresponding to the lidar observation, which makes it an ideal data source for estimating the LAI and VFP at the global scale [23][24][25][26][27][28]. Meanwhile, methods that serve to estimate the LAI or VFP using GLAS data come into rescue. Currently, two types of method have been developed in accordance with GLAS data to estimate the LAI or VFP: empirical methods and physical methods. The empirical methods were developed based on a statistical regression relationship between the GLAS-derived metrics and field-measured LAI [23,24]. Although these empirical methods are easy to implement and understand, they cannot be free from the need for the associated field data. Thus, this type of method is not adapted for application in large-scale LAI mapping because the empirical relationships usually vary with locations, species, canopy structures, etc.
Physical methods, which are constructed around the gap fraction obtained based on radiative transfer theory [25][26][27][28], open the door to large-scale LAI or VFP mapping using lidar data. For example, the LAI and VFP have been successfully retrieved over the entire state of California based on a geometric optical and radiative transfer (GORT) model [25]. By summarizing the current physical methods [25][26][27][28], we find that all methods are constructed depending on the received waveform information. However, this type of method has a limitation, that is, a ratio of canopy and ground reflectance needs to be set; the usual practice is to set a constant value for all GLAS measurements, such as 2.0 [27,28]. Although it is simple to operate, this strategy can introduce error because the ratio value often varies for different sites or even within sites due to different environmental conditions. In addition, a previous study concluded that the variation in this ratio impacts model performance [18]. At present, no one has tried to develop methods from other perspectives, such as using the transmitted Remote Sens. 2020, 12, 2457 3 of 23 energy profile (TEP). The TEP is an intuitive representation of the distribution density of forest canopy components along the lidar observation direction, which may be a new perspective for estimating the LAI and VFP.
The interaction between the laser and vegetation structure should be understood for better use of lidar measurements in estimating forest structure parameters. This reason drives the development of lidar modules in radiative transfer models (RTMs). Examples include FLIGHT [29], Digital Imaging and Remote Sensing Image Generation (DIRSIG) [30], RAYTRAN [31], LIBRAT [32], Forest Light Environmental Simulator (FLiES) [33], and the Discrete Anisotropic Radiative Transfer (DART) [34,35]. The lidar module in the DART model has now been extended to enable simulations of lidar sensitivity to different scene elements (e.g., vegetation, urban, and atmosphere) and sensor configuration (e.g., GLAS). Its performance in modeling lidar has been proven in the radiation transfer model intercomparison (RAMI) project [36]. Therefore, the lidar module in the DART model is an ideal tool for evaluating the constructed inversion methods based on the lidar data.
In this study, we proposed a new method for retrieving VFP and LAI free from the need for the ratio of canopy and ground reflectance-a major limitation with the existing methods that had not been addressed so far. First, we derived the TEP from the GLAS data based on the radiative transfer theory, which in turn was used to calculate the vertical gap distribution. Then, the obtained vertical gap distribution was utilized as an input to calculate the VFP based on the Beer-Lambert law, and the LAI could be obtained by integrating the VFP information. Next, we comprehensively validated the results with various realistic structure scenes (i.e., DART model) and field measurements (i.e., Tracing Radiation and Architecture of Canopies (TRAC) measurements), and the accuracies and sensitivities of our method were analyzed and discussed. Finally, we implemented our method to create footprint-level LAI over Heilongjiang Province, and comparative analysis was conducted between the GLAS LAI and Moderate Resolution Imaging Spectroradiometer (MODIS) LAI products.

Study Area
Our study was conducted in two study areas, that was sites A and B as shown in Figure 1. Site A is Heilongjiang Province (43 • 25 −53 • 33 N, 121 • 11 −135 • 05 E), which is located in northern China ( Figure 1a). Heilongjiang is a land of varied topography, with the average elevation varying from 200 m to 350 m. The forest survey in 2005 reported that Heilongjiang Province with forestry land area was 23.86 million hm 2 , and the forest volume was 1.432 billion m 3 , the province's forest coverage was as high as 41%; more than 60.4% of the forests were middle and young forests with an average diameter at breast height (DBH) of 12.8 cm [37]. The dominant forest types include larch (Larix gmelinii (Rupr.) Rupr.), birch (Betula platyphylla Suk), and scotch pine (Pinus sylvestris var. mongolica); the proportion of these three forest types is about 50% [38]. GLAS campaigns of forest over this study area were used to conduct LAI inversion; here, we termed it province-level LAI inversion.
Site B is the Saihanba National Forest Park (42 • 24 N, 117 • 15 E) located in Heibei Province, China (Figure 1b), which is state-owned and managed by the government. This park is the largest artificial forest in the world and serves as an important forest service research site in China. The climate is a temperate continental climate, with an average annual temperature of -1.2°C and an average annual rainfall of 452 mm. The terrain is complicated, with an elevation varying from 1500 m to 2067 m. Forest types are mainly composed of larch (artificial forest) and birch (natural forest), among which larch has an absolute advantage; larch area accounts for more than half of the total forest area [39]. At present, the forests are dominated by middle-aged and near-mature forests: the average age of larch is 29 with an average DBH of 19.04 cm and an average height of 12.71 m; the average age of birch is 30 with an average DBH of 16.84 cm and an average height of 11.61 m [40]. Validation work based on ground-based measurements was implemented in this area.
. Figure 1. Land cover distribution and locations of Geoscience Laser Altimeter System (GLAS) footprints over study areas A and B. (a) A is Heilongjiang Province, China. We conducted provincelevel leaf area index (LAI) inversion in this study area. (b) B is Saihanba National Forest Park in Heibei Province, China, which is our ground investigation area. ). GLA01 is a Level-1A product and provides waveform information [i.e., received echo waveform (r_rng_wf) and emitted transmitted waveform (r_tx_wf)] and time information [i.e., starting address of the transmit pulse sample (i_TxWfStart) and ending address of the range response (i_RespEndTime)]. GLA05 is a Level-1B product and provides the gain values of the received waveform and transmitted waveform (i_gval_rcv and i_gval_tx). GLA14 is a Level-2 product that contains information about the properties of the GLAS measurements, including the footprint centroid coordinate (i_lat and i_lon), max amplitude of the received echoes (i_maxRecAmp), standard deviation of the background noise (i_sDevNsObl), reflectivity (d_reflctUC), and reflectivity correction factor for atmospheric effects (d_reflCor_atm).

GlobeLand30 Landcover Data
The GlobeLand30 landcover images from paths 112-123 rows 23-30 were used to mask out the non forest-covered GLAS data in Heilongjiang Province so that we could focus on the forest. This product was produced by using multisource high-resolution remote sensing data, which include Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) multispectral images and Chinese Environmental Disaster Alleviation Satellite (HJ-1) multispectral images [41]. Ten land cover types were classified in GlobeLand30, as shown in Figure 1, and the accuracy of GlobeLand30 was approximately 83% at the global scale [42].

Shuttle Radar Topography Mission (SRTM) Data
The version 4 Shuttle Radar Topography Mission (SRTM) digital elevation dataset was used to produce the slope maps of the study areas with a resolution of 90 m. 5° × 5° SRTM tiles in GeoTIFF format corresponding to our study areas were downloaded from http://srtm.csi.cgiar.org/srtmdata/, and the downloaded data have been processed to fill data voids.  ). GLA01 is a Level-1A product and provides waveform information [i.e., received echo waveform (r_rng_wf) and emitted transmitted waveform (r_tx_wf)] and time information [i.e., starting address of the transmit pulse sample (i_TxWfStart) and ending address of the range response (i_RespEndTime)]. GLA05 is a Level-1B product and provides the gain values of the received waveform and transmitted waveform (i_gval_rcv and i_gval_tx). GLA14 is a Level-2 product that contains information about the properties of the GLAS measurements, including the footprint centroid coordinate (i_lat and i_lon), max amplitude of the received echoes (i_maxRecAmp), standard deviation of the background noise (i_sDevNsObl), reflectivity (d_reflctUC), and reflectivity correction factor for atmospheric effects (d_reflCor_atm).

GlobeLand30 Landcover Data
The GlobeLand30 landcover images from paths 112-123 rows 23-30 were used to mask out the non forest-covered GLAS data in Heilongjiang Province so that we could focus on the forest. This product was produced by using multisource high-resolution remote sensing data, which include Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) multispectral images and Chinese Environmental Disaster Alleviation Satellite (HJ-1) multispectral images [41]. Ten land cover types were classified in GlobeLand30, as shown in Figure 1, and the accuracy of GlobeLand30 was approximately 83% at the global scale [42].

Shuttle Radar Topography Mission (SRTM) Data
The version 4 Shuttle Radar Topography Mission (SRTM) digital elevation dataset was used to produce the slope maps of the study areas with a resolution of 90 m. 5 • × 5 • SRTM tiles in GeoTIFF format corresponding to our study areas were downloaded from http://srtm.csi.cgiar.org/srtmdata/, and the downloaded data have been processed to fill data voids.

Moderate Resolution Imaging Spectroradiometer (MODIS) LAI Product
We used the latest eight-day version of the MODIS LAI product (MCD15A2H Collection 6), which temporal phase is consistent with GLAS product, as the comparative dataset. Here, five MODIS tiles (i.e., H25V03, H25V04, H26V03, H26V04, H27V04) were used to mosaic the LAI map for Heilongjiang Province. Compared with the previous version (MCD15A2H Collection 5), the Collection 6 product was generated by using the improved reflectance and land cover product, which made it more accurate than the Collection 5 product [43]. Moreover, the Collection 6 product was generated at a native resolution of 500 m rather than the 1000 m resolution of Collection 5.

Ground-Based Data from Tracing Radiation and Architecture of Canopies (TRAC)
We collected 50 field-measured LAIs corresponding to nearly mature forest and mature forest across the Saihanba National Forest Park in August 2018 to validate our method (Figure 1b). Field measurement work was carried out with TRAC equipment, which could measure sunflecks for a whole canopy by walking along a transect that was in the range of tens of meters. Since the GLAS sensor emits energy from the zenith direction, therefore, we chose a time window with a small solar zenith angle (i.e., 10:00 a.m. to 3:00 p.m.) to conduct field work, to ensure the consistency of GLAS and TRAC observations. For each sampling, the Trimble GEO7X handheld GPS was first used to determine the position of GLAS footprint. Then the sunflecks within the GLAS footprint was measured by using the TRAC through walking along a transect with 30 m, which in turn was used to convert to LAI by TRACWin software installed on Microsoft Windows [44]. Here, the LAI from TRAC was an effective LAI because we did not consider the clumpiness effect. The absence of the clumpiness information did not affect the validation results because the clumpiness effect has the same contribution to both GLAS retrieval and TRAC measurement. In addition, photos of the field investigation footprints were obtained as evidence to describe the vertical foliage profile within footprints, which were used to evaluate the accuracy of the retrieved VFP.

Models and Methods
Our proposed method is a physically-based method based on the 1-D radiative transfer modeling of GLAS waveforms, which can fully utilize the forest canopy structure information contained in GLAS lidar data compared with commonly used empirical models. To assess how the new method can potentially improve on existing VFP and LAI products, we apply three strategies for algorithm assessments. The workflow of our method is shown in Figure 2. The details of the used model and method are introduced in Sections 3.1 and 3.2.

Physical Model of Light Detection and Ranging (Lidar) Echoes for Forest
We developed a physical model to derive the TEP for forest canopies, which linked lidar echoes to forest scenes based on the 1-D radiative transfer theory. In this model, foliage was thought to be the main cover to intercept solar radiation, and only the first collision between the laser and canopy was considered. The mathematical expression of this model is shown as follows: where R n is the echo energy of the recorded layer n, k is the extinction coefficient, LAI n describes the foliage between two successive recorded layers, ρ is the reflectance of the observed materials (e.g., forest canopy), and I n is the lidar transmitted energy received by recorded layer n and is defined as follows:

Physical Model of Light Detection and Ranging (Lidar) Echoes for Forest
We developed a physical model to derive the TEP for forest canopies, which linked lidar echoes to forest scenes based on the 1-D radiative transfer theory. In this model, foliage was thought to be the main cover to intercept solar radiation, and only the first collision between the laser and canopy was considered. The mathematical expression of this model is shown as follows: where Rn is the echo energy of the recorded layer n, k is the extinction coefficient, LAIn describes the foliage between two successive recorded layers, ρ is the reflectance of the observed materials (e.g., forest canopy), and In is the lidar transmitted energy received by recorded layer n and is defined as follows: ( )

Discrete Anisotropic Radiative Transfer (DART) Model
The DART model is a 3-D physically-based RTM that has been under development since 1992 [45]. After years of development, it is now a comprehensive model with multiple functions that can

Discrete Anisotropic Radiative Transfer (DART) Model
The DART model is a 3-D physically-based RTM that has been under development since 1992 [45]. After years of development, it is now a comprehensive model with multiple functions that can simulate not only radiative budgets but also remote-sensing images from different sensors (e.g., passive and active sensors), different platforms (e.g., airborne and spaceborne), and different bands for any Earth scene (e.g., vegetation and urban landscapes). The lidar module in the DART model was implemented based on the so-called "box method" and "Ray Carlo method" [35], which can provide accurate lidar simulations for any landscape and lidar sensor configuration (e.g., view direction, footprint size, and pulse characteristics). In this study, the lidar module of the DART model was used to validate our method because the simulated data usually had corresponding ground truth values that have an absolute consistent temporal phase with lidar observation. These so-called ground truth values were an ideal dataset for verifying the method. To ensure that the lidar data were simulated more realistically, we used the same parameter configuration as the GLAS device as the input parameters of the DART model (Table 1). In addition, the forest structure parameters set in the model were directly based on the forest structure characteristics in our study areas. The two following realistic canopy architecture scenarios were constructed in the DART model to validate our method: (1) homogeneous turbid scene ( Figure 3a): the scene only included the gap within the canopy; this scene was composed of a horizontal ground and a cube plant canopy with randomly distributed foliage. The thickness of the canopy was 9 m, and its bottom was 4 m above the ground. In addition, to verify the performance of our method under different vegetation abundance, we set the LAI values to three levels in this scene: high (i.e., LAI = 8), medium (i.e., LAI = 6), and low (i.e., LAI = 4). (2) Discrete individual-tree scene ( Figure 3b): the scene included large gaps between tree crowns; this scene consisted of trees with a canopy height of 9 m. The 3-D tree model is shown in Figure 3c. In this scene, the LAI of the GLAS footprint was 0.98 [i.e., the area covered by the red circle, as shown in Figure 3b].

Vertical Foliage Profile (VFP) and LAI Calculation
The basic principle for calculating the VFP and LAI from GLAS lidar data is based on the Beer-Lambert law. According to the Beer-Lambert law, we can establish a relationship between the leaf area density (LAD) and the degree of penetration of light within unit space, which is the space between two successive recorded layers as follows: where T(α) is the gap fraction of the unit space when the observation zenith angle is α; ∆l(α) is the

Vertical Foliage Profile (VFP) and LAI Calculation
The basic principle for calculating the VFP and LAI from GLAS lidar data is based on the Beer-Lambert law. According to the Beer-Lambert law, we can establish a relationship between the Remote Sens. 2020, 12, 2457 8 of 23 leaf area density (LAD) and the degree of penetration of light within unit space, which is the space between two successive recorded layers as follows: where T(α) is the gap fraction of the unit space when the observation zenith angle is α; ∆l(α) is the distance of the light traveling through the unit space along the zenith angle α; and G(α) is the leaf projection coefficient. For GLAS lidar, the vertical resolution is equal to 0.15 m; this means the thickness of the unit space is 0.15 m (i.e., ∆l(α) = 0.15). Therefore, we only needed the data on the leaf projection coefficient and gap fraction to obtain the LAD. The leaf projection coefficient G(α) is defined as the mean projection of a unit foliage area in the direction of the solar beam. For a canopy with a spherical leaf distribution, that is, the leaf inclination angle is 57.3 • , G is equal to 0.5 and independent of the solar zenith angle [46]. Due to the difficulty of estimating real leaf inclination angles, a spherical leaf inclination angle distribution is commonly assumed for a vegetative canopy. This assumption has been followed in studies on the inversion of LAI based on the GLAS data and performs well [25,27]. Therefore, the projection coefficient G(α) in this study was set to 0.5. The gap fraction T(α) of the unit space can be calculated by the TEP (i.e., I 0 , I 1 , . . . , I n ), and the calculation equation is shown as follows: where I n is the incident energy of the unit space of layer n and I n+1 is the energy passing through this space. The calculation of the TEP was introduced in Section 3.2.2. When the LAD is obtained, the corresponding LAI can be calculated by the relationship between the LAD and LAI as follows: where m and n are the layers of the LAD. The total LAI requires the integration of the LAD from the canopy top to ground. In addition, the cumulative LAI profiles can be obtained based on Equation (5) by calculating the increment in LAD at each height interval. In this study, we took the LAD profile and cumulative LAI profile as the VFP.

Transmitted Energy Profile (TEP) Calculation
The TEP was derived based on material reflectance characteristics. For lidar observation, the observed material reflectance characteristics usually depend on the sensor configuration and the relationship between the detected object and the laser. The sensor configuration can be expressed by the lidar equation [47], which is shown as follows: where R return is the received return energy, I trans is the emitted energy, τ opt is the optics transmission, A scope is the telescope area, D is the range between lidar sensor and ground, ρ is the reflectance of the target, τ atm is the roundtrip atmosphere transmission, and Ω is the scattering solid angle. The relationship between the observed object and laser can be given by the physical model of lidar echoes for forest (i.e., Equation (1)). By combining the lidar equation [i.e., Equation (6)] and the physical model of lidar echoes for a forest [Equation (1)], we can obtain an equation to express the observed material reflectance characteristics, which is shown as follows: Remote Sens. 2020, 12, 2457 9 of 23 Here, we assumed that the observed surfaces were Lambertian surfaces, so the Ω was equal to π. According to the relationship of the transmitted energy of two successive recorded layers (i.e., Equation (2)), we rewrote Equation (7); the result is shown as follows: In our application scene, two types of materials were involved: ground and vegetation. Based on Equation (8), the corresponding mathematical expression including these two materials can be given as follows: where I 0 is the total energy emitted by the lidar, R n−1 is the lidar received energy reflected from the vegetation, ρ vegetation is the vegetation reflectance, R n is the lidar received energy reflected from the ground surface, and ρ ground is the ground reflectance. The values of τ opt , A scope , τ atm , D, I 0 , and R i (I = 1,2, . . . ,n) are known, as shown in Table 2. Therefore, if we know the ground position in the recorded received waveform and the reflectance of the observed materials (i.e., ρ vegetation and ρ ground ), the TEP can be derived from Equation (8).
where C is the light transmission speed.
GLAS sensor emitted total pulse energy I 0 J I 0 is calculated from the GLAS transmitted waveform recorded in the GLAS product (i.e., r_tx_wf), and the calculation details are shown in the Appendix A.
GLAS sensor received echo energy at each recorded layer R i J R i is calculated from the GLAS received waveform recorded in the GLAS product (i.e., r_rng_wf), and the calculation details are shown in the Appendix A.
The ground position can be identified from the GLAS received waveform through the Gaussian decomposition method ( Figure 2); the process includes filtering, Gaussian decomposition, and determining the ground position. First, the raw GLAS echo waveform was denoised by Gaussian filtering. Then, the denoised GLAS echo waveform was decomposed into multiple Gaussian waveforms using the Gaussian decomposition method. Finally, the peak position of the last decomposed Gaussian waveform was assumed to correspond to the ground surface.
Equation (9) shows that if the ground reflectance is known, the canopy reflectance can be derived. Here, we used the average reflectance results (i.e., GLAS reflectance = d_reflctUC * d_reflCor_atm) from 12 non forest-covered GLAS footprints in Saihanba National Forest Park as ρ ground . The location of non forest-covered GLAS footprints is shown in Figure 2, the reflectance corresponding to these GLAS footprints is shown in Figure 4, which has an average result of 0.21 (i.e., ρ ground = 0.21). This strategy can be close to the actual reflectance of the GLAS laser for the ground.
Equation (9) shows that if the ground reflectance is known, the canopy reflectance can be derived. Here, we used the average reflectance results (i.e., GLAS reflectance = d_reflctUC * d_reflCor_atm) from 12 non forest-covered GLAS footprints in Saihanba National Forest Park as ρground.The location of non forest-covered GLAS footprints is shown in Figure 2, the reflectance corresponding to these GLAS footprints is shown in Figure 4, which has an average result of 0.21 (i.e., ρground = 0.21). This strategy can be close to the actual reflectance of the GLAS laser for the ground.

LAI Retrieval in Virtual Scenes
The VFP and LAI were first retrieved in the homogeneous turbid scene (Figure 3a), which had LAI values of 4, 6, and 8, respectively. Here, we used the cumulative LAI profile to represent the VFP. The vertical resolution of the cumulative LAI profile from the waveform data simulated by DART model was 0.9 m, matching that of the canopy in the virtual scene. Cumulative LAI profiles from our method showed the same trend as true values given from the DART model for all levels of LAI distributions [i.e., in Figure 5a-c: LAI ≈ 4, LAI ≈ 6, LAI ≈ 8]. Then, we applied our method in the discrete individual-tree scene (Figure 3b), and the result is shown in Figure 5d The cumulative LAI profile derived from our method still had strong consistency with the true value given by the DART model (i.e., LAI ≈ 0.98).
Remote Sens. 2020, 07, x FOR PEER REVIEW 11 of 23 The VFP and LAI were first retrieved in the homogeneous turbid scene [ Figure 3a], which had LAI values of 4, 6, and 8, respectively. Here, we used the cumulative LAI profile to represent the VFP. The vertical resolution of the cumulative LAI profile from the waveform data simulated by DART model was 0.9 m, matching that of the canopy in the virtual scene. Cumulative LAI profiles from our method showed the same trend as true values given from the DART model for all levels of LAI distributions [i.e., in Figure 5a-c: LAI ≈ 4, LAI ≈ 6, LAI ≈ 8]. Then, we applied our method in the discrete individual-tree scene [ Figure 3b], and the result is shown in Figure 5d The cumulative LAI profile derived from our method still had strong consistency with the true value given by the DART model (i.e., LAI ≈ 0.98). The above validation theoretically proves the reliability of our method because the scene constructed in the DART model was generated for an ideal situation, such as randomly distributed leaves and a spherical leaf inclination angle distribution. However, there are many uncertainties in the real environment. For example, the forest canopy is usually a mixture of branches and leaves, and it is difficult to separate them. Furthermore, the leaf inclination angle is sometimes not spherically distributed in an actual environment. Therefore, in addition to using virtual scenes for validation, ground-based measurements and other LAI products were used for cross-validation.  The above validation theoretically proves the reliability of our method because the scene constructed in the DART model was generated for an ideal situation, such as randomly distributed leaves and a spherical leaf inclination angle distribution. However, there are many uncertainties in the real environment. For example, the forest canopy is usually a mixture of branches and leaves, and it is difficult to separate them. Furthermore, the leaf inclination angle is sometimes not spherically distributed in an actual environment. Therefore, in addition to using virtual scenes for validation, ground-based measurements and other LAI products were used for cross-validation. Figure 6a presents the comparison between the GLAS LAI (i.e., the result of LAD integration from the canopy top to the ground) and the TRAC LAI. There was moderate consistency between them, with a determination coefficient (R 2 ) = 0.79 and root mean square error (RMSE) = 0.83. However, the GLAS LAI was higher than the TRAC LAI (bias = 0.68), which may be because TRAC works. In the measurement process, the TRAC is on a transect approximately 1 m above the ground, and the sensor on the TRAC is upward; therefore, the vegetation near the ground (i.e., understory vegetation) may be ignored. To prove our conjecture, we chose a set of comparative data, and the corresponding plot photos are shown in Figure 6a. Compared to the LAI derived from the footprint with almost no understory vegetation, the LAI from the footprint with dense understory vegetation has a considerable bias with the TRAC LAI. Therefore, the difference between ground observations and satellite observations was the main reason for inconsistencies in the comparison. To ensure consistency between the observations, we recalculated the GLAS LAI by integrating the LAD from the canopy top to the layer 1 m above the ground. We then obtained improved validation results (Figure 6b From Figure 6b, we can find that there were some visually obvious outliers in the validation results. It has been well documented by previous studies that slope and SNR (signal-to-noise ratio = i_maxRecAmp/i_sDevNsObl) have a great influence on the accuracy of lidar measurements [48][49][50][51][52][53]. Some studies concluded that GLAS observations with slopes of less than 15 and SNR values greater than 60 could accurately represent the vegetation structure [27,54]. Therefore, we analyzed the role of slope and SNR in our method to determine whether they caused these outliers. The distribution of the slope and SNR values of the GLAS footprints in Saihanba National Forest Park, which provide by SRTM data and GLAS product, are shown in Figure 7a and Figure 7b. We found four footprints with slopes greater than 15 and five footprints with SNR values less than 60. Among these potentially problematic GLAS observations, three observations not only had large slopes (i.e., slope > 15) but also low SNR values (i.e., SNR < 60); that is, six of 50 GLAS observations might have been problematic. After removing these six potentially problematic GLAS observations, the validation accuracy was improved, with R 2 = 0.85, RMSE = 0.35, and bias = 0.10 [ Figure 7c]. The above result indicated that the From Figure 6b, we can find that there were some visually obvious outliers in the validation results. It has been well documented by previous studies that slope and SNR (signal-to-noise ratio = i_maxRecAmp/i_sDevNsObl) have a great influence on the accuracy of lidar measurements [48][49][50][51][52][53]. Some studies concluded that GLAS observations with slopes of less than 15 and SNR values greater than 60 could accurately represent the vegetation structure [27,54]. Therefore, we analyzed the role of slope and SNR in our method to determine whether they caused these outliers. The distribution of the slope and SNR values of the GLAS footprints in Saihanba National Forest Park, which provide by SRTM data and GLAS product, are shown in Figure 7a,b. We found four footprints with slopes greater than 15 and five footprints with SNR values less than 60. Among these potentially problematic GLAS observations, three observations not only had large slopes (i.e., slope > 15) but also low SNR values (i.e., SNR < 60); that is, six of 50 GLAS observations might have been problematic. After removing these six potentially problematic GLAS observations, the validation accuracy was improved, with R 2 = 0.85, RMSE = 0.35, and bias = 0.10 ( Figure 7c). The above result indicated that the uncertainty caused by the slope and SNR were probably two primary error sources in the LAI estimation using the method proposed in this study.

GLAS LAI vs. TRAC LAI
by SRTM data and GLAS product, are shown in Figure 7a and Figure 7b. We found four footprints with slopes greater than 15 and five footprints with SNR values less than 60. Among these potentially problematic GLAS observations, three observations not only had large slopes (i.e., slope > 15) but also low SNR values (i.e., SNR < 60); that is, six of 50 GLAS observations might have been problematic. After removing these six potentially problematic GLAS observations, the validation accuracy was improved, with R 2 = 0.85, RMSE = 0.35, and bias = 0.10 [ Figure 7c]. The above result indicated that the uncertainty caused by the slope and SNR were probably two primary error sources in the LAI estimation using the method proposed in this study.  We further explored how these two factors affected our method. Here, we selected typical GLAS observations with large slopes and low SNRs for further analysis. Figure 8 presents an example of GLAS observations with complex terrain (i.e., slope = 21 • ). Notably, the echo energy from the canopy was mixed with the echo energy from the bare ground on the top of the mountain, which resulted in considerable canopy echo energy (Figure 8b). In addition, the obvious scattering phenomenon on the slope area led to the echo energy of the last Gaussian component being very small (Figure 8b). The above two reasons caused the illusion that the observation footprint had dense canopy (i.e., large LAI); however, the reality was that the observation footprint had a sparse forest (i.e., small LAI), as shown in Figure 8a. The ground validation in this footprint was very poor; the difference between the GLAS LAI and the TRAC LAI was equal to 1.84, which was the largest of all validation footprints. Figure 9 shows an example of a comparison of GLAS observations with low and high SNR, respectively (i.e., SNR = 21 vs. 119). Low SNR means poor data quality, and it usually increases the difficulty in preprocessing. From Figure 9, we can note that the GLAS observation with low SNR has difficulty recognizing useful information from the received waveform compared to the GLAS observation with a high SNR value. Therefore, uncertainty may be introduced during the preprocessing process, such as determining the ground position, thereby affecting the accuracy of LAI derived from our method. The validation based on the TRAC measurement also proves the impact of SNR on our method, that is, the GLAS observation with high SNR value performed much better than the one with low SNR value; the difference between the GLAS LAI and TRAC LAI was 0.17 vs. 1.26. was mixed with the echo energy from the bare ground on the top of the mountain, which resulted in considerable canopy echo energy [ Figure 8b]. In addition, the obvious scattering phenomenon on the slope area led to the echo energy of the last Gaussian component being very small [ Figure 8b]. The above two reasons caused the illusion that the observation footprint had dense canopy (i.e., large LAI); however, the reality was that the observation footprint had a sparse forest (i.e., small LAI), as shown in Figure 8a. The ground validation in this footprint was very poor; the difference between the GLAS LAI and the TRAC LAI was equal to 1.84, which was the largest of all validation footprints.  Figure 9 shows an example of a comparison of GLAS observations with low and high SNR, respectively (i.e., SNR = 21 vs. 119). Low SNR means poor data quality, and it usually increases the difficulty in preprocessing. From Figure 9, we can note that the GLAS observation with low SNR has difficulty recognizing useful information from the received waveform compared to the GLAS observation with a high SNR value. Therefore, uncertainty may be introduced during the preprocessing process, such as determining the ground position, thereby affecting the accuracy of LAI derived from our method. The validation based on the TRAC measurement also proves the impact of SNR on our method, that is, the GLAS observation with high SNR value performed much better than the one with low SNR value; the difference between the GLAS LAI and TRAC LAI was 0.17 vs. 1.26.   It can be noted that the ground echo signal (Rg) was mixed with the canopy echo signal (Rv) because of the complex terrain. In addition, the obvious scattering phenomenon (S) on the slope area caused the echo energy of the last Gaussian component to be very small. Figure 9 shows an example of a comparison of GLAS observations with low and high SNR, respectively (i.e., SNR = 21 vs. 119). Low SNR means poor data quality, and it usually increases the difficulty in preprocessing. From Figure 9, we can note that the GLAS observation with low SNR has difficulty recognizing useful information from the received waveform compared to the GLAS observation with a high SNR value. Therefore, uncertainty may be introduced during the preprocessing process, such as determining the ground position, thereby affecting the accuracy of LAI derived from our method. The validation based on the TRAC measurement also proves the impact of SNR on our method, that is, the GLAS observation with high SNR value performed much better than the one with low SNR value; the difference between the GLAS LAI and TRAC LAI was 0.17 vs. 1.26.

GLAS VFP Inversion
We selected two typical inversion results of the VFP for display and analysis to verify our method, and the results are shown in Figure 10, including those for site A and site B. Site A corresponds to a birch forest, in which the foliage is sparse, and there is almost no understory vegetation ( Figure 10). Site B corresponds to a dense larch forest with lush understory vegetation ( Figure 10). Here, we used the LAD profile to represent the VFP. The vertical resolution of the VFP was 0.15 m, which was consistent with the resolution of the echoes recorded by GLAS lidar.
At site A, the LAD gradually increases as the canopy height decreased and exhibited peak leaf density at a height of approximately 7 m; then, the LAD shows a downward trend and gradually increases when approaching the ground. The LAD distribution in the vertical direction was nearly consistent with the foliage distribution density shown in the photo. From the ground measurement side, the integration of LAD from the canopy top to the ground is almost consistent with TRAC measurement, which is 3.36 vs. 3.09.
At site B, the LAD distribution exhibited the same trend as site A except for the part close to the ground; this is because site B had denser understory vegetation than site A (Figure 10), which caused the LAD for site B to be larger than that for site A when approaching the ground. In addition, we can note that the LAD value of site B performed larger than that of site A in the middle parts of the canopy (e.g., at canopy heights equal to 8 m); this result was roughly consistent with the density of foliage shown in the photos. The LAI from the TRAC measurement for site B shows a bias with the integration of LAD from the canopy top to the ground (i.e., 3.28 vs. 4.50). Recalling the ground validation analysis in Section 4.2.1, the analysis result explains the understory vegetation caused the positive bias between the integration of LAD from the canopy top to the ground and ground measurement. Here, our validation result fits this reason, which also indirectly proves the rationality of LAD inversed based on our method. measurement, which is 3.36 vs. 3.09.
At site B, the LAD distribution exhibited the same trend as site A except for the part close to the ground; this is because site B had denser understory vegetation than site A (Figure 10), which caused the LAD for site B to be larger than that for site A when approaching the ground. In addition, we can note that the LAD value of site B performed larger than that of site A in the middle parts of the canopy (e.g., at canopy heights equal to 8 m); this result was roughly consistent with the density of foliage shown in the photos. The LAI from the TRAC measurement for site B shows a bias with the integration of LAD from the canopy top to the ground (i.e., 3.28 vs. 4.50). Recalling the ground validation analysis in Section 4.2.1, the analysis result explains the understory vegetation caused the positive bias between the integration of LAD from the canopy top to the ground and ground measurement. Here, our validation result fits this reason, which also indirectly proves the rationality of LAD inversed based on our method.

Uncertainty Analysis for Ground Validation
In the ground validation, there is a 10-year temporal difference between the ground measurement and the GLAS data, which will introduce uncertainties in the validation process. Despite this weakness, our validation results are still useful and informative for two major reasons. First, the main purpose of our ground-based validation is to assess how well GLAS can capture spatial patterns and heterogeneity in LAI. This spatial variability contributed to the majority of the variation in both the in situ and GLAS LAI data-a signature that should be independent of time. In other words, our validation should provide valid quantification of the GLAS accuracies over space. Second, to overcome the limitation associated with the temporal discrepancy, we used the MODIS LAI product (i.e., MCD15A2H.A2018169 vs. MCD15A2H.A2006169) as an intermediary to evaluate the variation of LAI in 10 years for GLAS footprints located in Saihanba National Forest Park. The result shows that the LAI of 2018 has increased compared with the LAI of 2006, but the overall change

Uncertainty Analysis for Ground Validation
In the ground validation, there is a 10-year temporal difference between the ground measurement and the GLAS data, which will introduce uncertainties in the validation process. Despite this weakness, our validation results are still useful and informative for two major reasons. First, the main purpose of our ground-based validation is to assess how well GLAS can capture spatial patterns and heterogeneity in LAI. This spatial variability contributed to the majority of the variation in both the in situ and GLAS LAI data-a signature that should be independent of time. In other words, our validation should provide valid quantification of the GLAS accuracies over space. Second, to overcome the limitation associated with the temporal discrepancy, we used the MODIS LAI product (i.e., MCD15A2H.A2018169 vs. MCD15A2H.A2006169) as an intermediary to evaluate the variation of LAI in 10 years for GLAS footprints located in Saihanba National Forest Park. The result shows that the LAI of 2018 has increased compared with the LAI of 2006, but the overall change is not very large (Figure 11a). Besides, it can be noted that the LAI of 2018 has a moderate correlation with the LAI of 2006 (i.e., R 2 = 0.41) ( Figure 11b); thus, the error caused by the time discrepancy in the validation could be taken as a systematic error. Therefore, we can reasonably consider the impact of this error on all results is almost the same. This does not necessarily affect the overall validity of the correlation between in situ and GLAS data.

LAI Retrieval in Heilongjiang Province
We derived a total of 4071 GLAS LAI over Heilongjiang Province. Among these GLAS LAI values, approximately 35% were derived from the GLAS observations with large slope (i.e., slope > 15 • ) or low SNR (e.g., SNR < 60). According to the analysis results in Section 4.2.1, we thought that 35% of the inversion results might have uncertainty. To ensure the accuracy of the analysis, we only used the remaining 65% of inversion results for analysis. Figure 12 shows the distribution of GLAS LAI after removing these the 35% of values that were potentially problematic results, which totaled 2628 GLAS LAI values, with a mean value of 2.64 and a mean square error of 1.15.
The vertical LAI integrations from 0-4 m, 4-8 m, and 8-18 m were also mapped at the scale of GLAS footprints (~65 m) ( Figure 13). We noted that the LAIs at different height levels performed differently; in particular, the LAI integrated from 0-4 m was relatively large, which indicated that there was understory vegetation in the forest of Heilongjiang Province. In addition, we also found that in some cases, although the total LAI (e.g., LAI ≈ 4) was the same, the LAIs at different height levels had different contributions to the total, as seen in Figure 14.
is not very large [ Figure 11a]. Besides, it can be noted that the LAI of 2018 has a moderate correlation with the LAI of 2006 (i.e., R 2 = 0.41) [ Figure 11b]; thus, the error caused by the time discrepancy in the validation could be taken as a systematic error. Therefore, we can reasonably consider the impact of this error on all results is almost the same. This does not necessarily affect the overall validity of the correlation between in situ and GLAS data.

LAI Retrieval in Heilongjiang Province
We derived a total of 4071 GLAS LAI over Heilongjiang Province. Among these GLAS LAI values, approximately 35% were derived from the GLAS observations with large slope (i.e., slope >15°) or low SNR (e.g., SNR < 60). According to the analysis results in Section 4.2.1, we thought that 35% of the inversion results might have uncertainty. To ensure the accuracy of the analysis, we only used the remaining 65% of inversion results for analysis. Figure 12 shows the distribution of GLAS LAI after removing these the 35% of values that were potentially problematic results, which totaled 2628 GLAS LAI values, with a mean value of 2.64 and a mean square error of 1.15.

LAI Retrieval in Heilongjiang Province
We derived a total of 4071 GLAS LAI over Heilongjiang Province. Among these GLAS LAI values, approximately 35% were derived from the GLAS observations with large slope (i.e., slope >15°) or low SNR (e.g., SNR < 60). According to the analysis results in Section 4.2.1, we thought that 35% of the inversion results might have uncertainty. To ensure the accuracy of the analysis, we only used the remaining 65% of inversion results for analysis. Figure 12 shows the distribution of GLAS LAI after removing these the 35% of values that were potentially problematic results, which totaled 2628 GLAS LAI values, with a mean value of 2.64 and a mean square error of 1.15.  The vertical LAI integrations from 0-4 m, 4-8 m, and 8-18 m were also mapped at the scale of GLAS footprints (~65 m) ( Figure 13). We noted that the LAIs at different height levels performed differently; in particular, the LAI integrated from 0-4 m was relatively large, which indicated that there was understory vegetation in the forest of Heilongjiang Province. In addition, we also found that in some cases, although the total LAI (e.g., LAI ≈ 4) was the same, the LAIs at different height levels had different contributions to the total, as seen in Figure 14.

GLAS LAI vs. MODIS LAI
The GLAS and MODIS LAIs were found to correlate poorly ( Figure 15). The overall correlation coefficient was only sqrt (0.003), although the comparison of the two gave relatively low RMSE (1.69) and bias (0.89). The differences between the GLAS and MODIS LAI values were highly skewed because of large positive differences at low and middle LAIs [ Figure 15b]. The main reason for the poor correlation was likely attributed to the differing characteristics of the two LAI products. In particular, the MODIS LAI retrieval algorithm accounted for vegetation clumping through 3-D radiative transfer modeling [55,56] and was supposed to provide the true LAI. In contrast, our GLAS method did not consider the clumping effect and derives only the effective LAI. Therefore, the MODIS LAI tended to be higher than the GLAS LAI. Interestingly, for the LAIs higher than 4, the pattern was reversed with higher values for GLAS [ Figure 15a]. This reversal was due to the saturation of MODIS LAI at approximately 5.5 [ Figure 15a]. However, the GLAS LAI estimates showed no signs of saturation, even up to 8 [ Figure 15a]. Overall, the comparison clearly highlighted the unique value or even potential superiority of lidar data for measuring LAI for dense forests-an expected finding consistent with the majority of lidar studies over forests of medium to high biomass [57].

GLAS LAI vs. MODIS LAI
The GLAS and MODIS LAIs were found to correlate poorly ( Figure 15). The overall correlation coefficient was only sqrt (0.003), although the comparison of the two gave relatively low RMSE (1.69) and bias (0.89). The differences between the GLAS and MODIS LAI values were highly skewed because of large positive differences at low and middle LAIs (Figure 15b). The main reason for the poor correlation was likely attributed to the differing characteristics of the two LAI products. In particular, the MODIS LAI retrieval algorithm accounted for vegetation clumping through 3-D radiative transfer modeling [55,56] and was supposed to provide the true LAI. In contrast, our GLAS method did not consider the clumping effect and derives only the effective LAI. Therefore, the MODIS LAI tended to be higher than the GLAS LAI. Interestingly, for the LAIs higher than 4, the pattern was reversed with higher values for GLAS (Figure 15a). This reversal was due to the saturation of MODIS LAI at approximately 5.5 (Figure 15a). However, the GLAS LAI estimates showed no signs of saturation, even up to 8 (Figure 15a). Overall, the comparison clearly highlighted the unique value or even potential superiority of lidar data for measuring LAI for dense forests-an expected finding consistent with the majority of lidar studies over forests of medium to high biomass [57].

Model Uncertainty Analysis
The ground reflectance (ρg) is a necessary input parameter in our method; therefore, incorrect input value may lead to errors in LAI estimation. Ground reflectance often varies among different sites or even within sites due to different environmental conditions, such as water content and ground cover. We used a mean value from the non forest-covered GLAS footprints as the input for the whole study area, which may have introduced errors to the results. To see what role the ground reflectance played in our method, we analyzed the sensitivity of our method to variations in input ground reflectance. We varied the ground reflectance from 1.5 to 2.9 to evaluate its effect on LAI and

Model Uncertainty Analysis
The ground reflectance (ρ g ) is a necessary input parameter in our method; therefore, incorrect input value may lead to errors in LAI estimation. Ground reflectance often varies among different sites or even within sites due to different environmental conditions, such as water content and ground cover.
We used a mean value from the non forest-covered GLAS footprints as the input for the whole study area, which may have introduced errors to the results. To see what role the ground reflectance played in our method, we analyzed the sensitivity of our method to variations in input ground reflectance. We varied the ground reflectance from 1.5 to 2.9 to evaluate its effect on LAI and found that its effect on the inversion results was minimal and nearly negligible ( Figure 16).
legend explains the number of data points that fell into certain square areas. (b) Histogram of LAI difference between GLAS and MODIS; the red line indicates the bias (0.89).

Model Uncertainty Analysis
The ground reflectance (ρg) is a necessary input parameter in our method; therefore, incorrect input value may lead to errors in LAI estimation. Ground reflectance often varies among different sites or even within sites due to different environmental conditions, such as water content and ground cover. We used a mean value from the non forest-covered GLAS footprints as the input for the whole study area, which may have introduced errors to the results. To see what role the ground reflectance played in our method, we analyzed the sensitivity of our method to variations in input ground reflectance. We varied the ground reflectance from 1.5 to 2.9 to evaluate its effect on LAI and found that its effect on the inversion results was minimal and nearly negligible ( Figure 16).

Discussion
Our proposed method for retrieving forest VFP and LAI is based on physical principle, rather than based on a statistical regression, which is currently a preferential method by modelers. In addition, our method used more information provided by GLAS products compared with the previous methods only using received waveform information and solved a problem of existing methods, which will likely improve the accuracy of VFP and LAI estimation. Although there is a limitation that we cannot provide accurate ground reflectance as input, the analysis in Section 4.4 has shown that it has a small effect on the result. Overall, our research provided a new perspective for retrieving the VFP and LAI.
To test the performance of our method, we used ground-based data to validate our method. However, there were some limitations in our validation process. First, the different times of GLAS data acquisition and TRAC measurements will introduce errors. As shown in the comparison scattering plot [ Figure 7c], the comparison points were roughly distributed on both sides of the 1:1 Figure 16. An example of a LAD profile derived from GLAS data with different values of ρ g .

Discussion
Our proposed method for retrieving forest VFP and LAI is based on physical principle, rather than based on a statistical regression, which is currently a preferential method by modelers. In addition, our method used more information provided by GLAS products compared with the previous methods only using received waveform information and solved a problem of existing methods, which will likely improve the accuracy of VFP and LAI estimation. Although there is a limitation that we cannot provide accurate ground reflectance as input, the analysis in Section 4.4 has shown that it has a small effect on the result. Overall, our research provided a new perspective for retrieving the VFP and LAI.
To test the performance of our method, we used ground-based data to validate our method. However, there were some limitations in our validation process. First, the different times of GLAS data acquisition and TRAC measurements will introduce errors. As shown in the comparison scattering plot (Figure 7c), the comparison points were roughly distributed on both sides of the 1:1 line, which may have been caused by the above reason. Second, due to the limitations of the field measurement equipment, we were unable to obtain the LAD profile data through measurement; therefore, we used photos of the sites along the observation footprints as a reference to evaluate the rationality and accuracy of the VFP but did not conduct a quantitative evaluation. Regardless, the quantitative evaluation of the VFP has been addressed by a variety of virtual scenes. In view of the limitation in the validation of this work, future studies should place an emphasis on field data collection to ensure that the measurement time is consistent with the satellite LiDAR data.
We also extended our method to Heilongjiang Province. For this province-level LAI inversion, there was an error source, which may come from landcover information provided by the GlobeLand30 landcover product. Our research focused on forests, so the non forest-covered GLAS data had to be filtered out. In Saihanba National Forest Park, we selected forest-covered GLAS data by visual interpretation of high-resolution images and on site investigation. However, this strategy is challenging to conduct in large study areas, such as those at the province level, because of time and labor costs. Here, we adopted an effective strategy to select the forest-cover GLAS data that used the high-resolution land-cover product to filter out the non forest cover GLAS data. The GlobeLand30 product has an accuracy of 83% at a global scale [42]; therefore, the misclassification of landcover types in the GlobeLand30 product may lead to errors in the LAI inversion in Heilongjiang Province. For example, the low LAI values (e.g., LAI < 1) shown in Figure 11 may not have resulted from the GLAS observations of the forest because of the misclassification of land-cover types.
Due to the lack of available validation datasets at the province level, we therefore limited the assessment of the GLAS LAI differences in Heilongjiang Province by using the MODIS LAI product. The result showed a weak correlation between the two; the reason is that our method did not consider canopy clumpiness. For any given effective LAI value, the error in the true LAI varied linearly with the clumping index. Previous studies have used a strategy that builds a clumping index look-up table by assigning each land cover class an average clumping index from coarse resolution clumping index products such as MODIS or Polarization and Directionality of the Earth's Reflectance (POLDER) to revise this error. However, this simple and crude method does not take into account the difference in the level of clumping of the same landcover type, so it may introduce error. In addition, the scale effect, which is caused by the direct application of coarse-resolution products to high-resolution products, needs to be considered [58]. Currently, a major challenge is that no high-resolution clumping index products are available, especially those products that match the GLAS footprint level, although several moderate-resolution clumping index products have been generated and released to users [59,60], based on the hotspot-corrected bidirectional reflectance distribution function (BRDF) model [61,62] and a linear model [63]. Considering the potential problems of using coarse resolution data, we focus on the effective LAI at present. The problem of the clumpiness effect will be solved in the future, for example, to develop a method that can estimate the clumping index with GLAS data.
One of the major limitations of our method was that assuming all the reflecting components within the GLAS footprint behaved as Lambertian scatterers because of the lack of knowledge of the directional reflectance properties of various components at the GLAS resolution. This assumption has been widely used in the calculation of surface albedo [64,65]. Regardless of this, the assumption of Lambertian scattering is a practical approximation of the potentially complex character of land footprints.
The advantage of our work is that we can obtain the vertical foliage profiles at the landscape scale, which allows us to better understand vertical structural parameters in the forest. For example, our method can describe the understory vegetation very well, as shown in Figure 12 (i.e., LAI integration from 0-4 m). Understory vegetation is usually a forest ecosystem driver, and its productivity is probably comparable to that of trees [66]. Previous studies have suggested that the lack of vertical foliage distribution information caused approximately 50% underestimation of GPP [67,68]. Our method provides an opportunity to obtain foliage profile information on the global scale, which should be helpful for improving the accuracy of ecological process simulations.
Recently two new lidar observation platforms GEDI and ICESat-2 were successfully launched and operated. They have improved instrument technologies and spatial coverage compared with ICESat GLAS. In the next step, our method will combine with these two new platforms, which will allow for enhanced coverage and improved accuracy by leveraging the advantages of the individual instrument technologies and spatial coverage.

Conclusions
The advances in complexity and resolution in land-surface modeling lead to increasing demand for more detailed descriptions of canopy architecture, such as the VFP. Traditional methods for estimating LAI only focus on horizon information and do not provide vertical structure information throughout the canopy. In this work, we developed a new method to retrieve the VFP and LAI using the transmitted energy information derived from the ICESat GLAS data. The advantage of our method is that it is physically based rather than an empirical method; it follows the radiative transfer theory. We have also listed the limitations with our method, namely, the terrain factor and signal noise greatly affected our results. How these two factors affect the inversion results has been analyzed and discussed; the results can guide future efforts in this research direction. Comprehensive validations are conducted for our method from virtual and real scenes, and the results indicate that our method can be successfully used to retrieve VFP and total LAI. Our method can be used as the basis for obtaining continental-and global-scale VFP and total LAI datasets from ICESat GLAS observations, or can even be extended to new satellite lidar platforms such as GEDI and ICESat-2.
Funding: This work was supported by the National Natural Science Foundation of China (Nos. 41971288 and 41801237).

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
The GLAS emitted and received pulse energy can be computed from the digitized pulse waveforms recorded in the GLAS product, and the mathematical formula is as follows: where E is the pulse energy, α is the GLAS receiver calibration constant, η e is the electronic throughput, η is the optical throughput, T Det is the detector responsibility, G is the variable gain amplifier (VGA) gain, W(i) is the pulse waveform recorded in volts, N is the waveform sampling number, and ∆t is the sampling time of the GLAS receiver. The above parameters are different for transmitted and received pulse energy; The corresponding parameters are given in Table A1.