Retrieval and validation of vertical LAI profile derived from airborne and spaceborne LiDAR data at a deciduous needleleaf forest site

ABSTRACT

Leaf area index (LAI) is defined as one half of the total green leaf area per unit ground surface area. Its vertical profile is critical for understanding the remote sensing radiative transfer processes. LAI profile has been derived from airborne and spaceborne LiDAR data, such as the Global Ecosystem Dynamics Investigation (GEDI) installed on the International Space Station. However, the capability of various algorithms for the LAI profile estimation with airborne LiDAR is not clearly evaluated, and the estimated LAI profiles, including the GEDI LAI products, are not been fully validated. This study conducted a quantitative retrieval and validation of the LAI profiles using terrestrial and airborne laser scanning (TLS and ALS) and spaceborne GEDI data over a deciduous needleleaf forest site in northern China. The vertical LAI profile was estimated in the field using an upward digital hemispherical photography (DHP) attached to a portable measurement system in 2020 and 2021. A suite of new LiDAR indices combining both LiDAR return number and return intensity was explored for the LAI profile estimation. All LAI profiles obtained from the DHP, TLS, ALS, and GEDI during the leaf-on season and leaf-off season were compared. The DHP shows a good agreement with the TLS LAI profiles (R² = 0.97). The LAI profile derived from the ALS data using the combined light penetration index (LPI_RI) agrees well (R² ≥0.86) with the DHP, TLS, and GEDI estimates. In general, the LPI_RI is advantageous for regional LAI profile mapping from ALS. The GEDI cumulative LAI corresponds well with the DHP during the leaf-on season (R² = 0.90, RMSE = 0.23), but underestimates during the leaf-off season (R² = 0.70, RMSE = 0.14, bias=−0.13). The underestimation is attributed to the higher canopy and ground reflectance ratio (ρ_v/ρ_g) assigned in the algorithm and the height discrepancy between the GEDI and field measurements. For the GEDI LAI profile product, further validation and improvement are necessary for other biome types and landscape conditions, especially during the leaf-off season.

KEYWORDS:

• Leaf area index (LAI)
• vertical profile
• digital hemispherical photography (DHP)
• airborne laser scanning (ALS)
• Global Ecosystem Dynamics Investigation (GEDI)
• terrestrial laser scanning (TLS)

1. Introduction

The leaf area index (LAI) is an important vegetation biophysical parameter and a critical variable in photosynthesis and respiration processes (Alton 2016; Asner et al. 1998; Chen and Black 1992). Forests display a vertical stratification of foliage, which affects the light used in photosynthesis, and water interception (Goetz et al. 2007; Hopkinson et al. 2013). Knowledge regarding the forest LAI profile is significant for understanding remote sensing radiative transfer processes (Béland and Baldocchi 2021; Gastellu-Etchegorry, Martin, and Gascon 2004; Parker, Lefsky, and Harding 2001; Qi et al. 2022; Yang et al. 2020).

The LAI can be measured by direct sampling methods and indirect optical methods (Fang et al. 2019; Jonckheere et al. 2004; Weiss et al. 2004; Zheng and Moskal 2009). Both methods can be used to obtain the forest LAI profile by utilizing observation towers and mobile elevators (Clark et al. 2008; Wang et al. 1992). For example, Clark et al. (2008) harvested foliage at 1.86 m height intervals via a modular tower in a tropical rain forest to obtain the landscape LAI profile. Ryu et al. (2014) conducted a tower-based measurement using LAI-2200 Plant Canopy Analyzer (PCA) and digital hemispherical photography (DHP) to monitor the variation of forest vertical LAI. However, ground measurement of the LAI profile is labor-intensive and large area deployment is difficult.

On the other hand, remote sensing provides an optimal tool for LAI estimation on large scales (Fang et al. 2019). LAI has been derived from passive remote sensing data and active light detection and ranging (LiDAR) data (Baret et al. 2013; Zheng and Moskal 2009). However, passive remote sensing systems do not adequately capture vertical variation in LAI. In contrast, LiDAR provides unique advantages in capturing the vertical dimension information (Guo et al. 2021; Wang et al. 2019; Yan et al. 2019). The terrestrial laser scanning (TLS) provides a 3D point cloud from which LAI can be extracted at the individual tree or forest stand levels (Beland et al. 2011; Guo et al. 2021; Zheng et al. 2017; Zheng and Moskal 2012). However, TLS beams sampled in the upper canopy are easily obstructed by the middle and lower canopy components, which limits its broad application in forests.

The airborne laser scanning (ALS) can better characterize the upper layer and is commonly acquired at a large scale, which facilitates the LAI estimation at the regional level (Wang and Fang 2020). Various LiDAR indices, such as the light penetration index (LPI), canopy fractional cover index (FCI), and the above and below ratio index (ABRI), have been proposed to estimate the canopy LAI and LAI profile using various regression models (Hopkinson and Chasmer 2009; Sumnall et al. 2016, 2016, 2021; Zhao and Popescu 2009). However, the statistical regression models are specific to the area in which they are developed. The LiDAR indices have shown different performances in various forests (Alonzo et al. 2015; Sumnall et al. 2021; Zhao and Popescu 2009). Recently, Sumnall et al. (2021) evaluated the capacity of eight LiDAR indices to estimate the forest overstory and understory LAIs. In their study, the forest was divided into overstory and understory with a threshold height corresponding to the height to live crown (HTLC). LiDAR indices derived from each layer were compared with field LAI, and the highest correlations are produced for overstory (R² = 0.71) and understory (R² = 0.49), respectively. However, the methods for the two-layer LAI may not be applicable for estimating the multi-layer LAI profiles. Moreover, it is still unclear which index is optimal for the LAI profile estimation. More extensive validation studies are necessary for the LiDAR LAI profile data estimated from various methods.

Compared to the TLS and ALS data, the spaceborne LiDAR is advantageous for forest measurement on a global scale. The Geoscience Laser Altimeter System (GLAS) was the first spaceborne full waveform LiDAR. GLAS provided the capability to map LAI and LAI profiles over large areas (Tang et al. 2014, 2016), but the process was affected by the slope and the mission was terminated in 2009. The Global Ecosystem Dynamics Investigation (GEDI) installed on the International Space Station is a spaceborne full waveform LiDAR with a 25 m footprint size and a 60 m along-track sampling interval (Dubayah et al. 2020; Roy, Kashongwe, and Armston 2021). GEDI provides high-quality measurements of vertical forest structure in temperate and tropical forests between 51.6° N and 51.6° S. A variety of footprint and gridded products (Level 1 to Level 4) were derived from GEDI observations, including canopy height, canopy cover and vertical profile, canopy LAI and LAI profile, topography, and biomass (Dubayah et al. 2020). GEDI Level 2B product provides the first global footprint level LAI profile data at a vertical resolution of 5 m since April 2019 (Dubayah et al. 2020). The LAI profile is derived from the modified Beer-Lambert method, which employs a global constant canopy and ground reflectance ratio (ρ_v/ρ_g = 1.5) at the footprint level (Dhargay et al. 2022; Tang and Armston 2020). However, whether this ratio is applicable for different biome types and temporal periods is questionable. Moreover, validation of the new GEDI LAI profile is lacking.

GEDI data are acquired throughout the year, and LAI data are available during both the leaf-on and leaf-off seasons. The LAI acquired during the leaf-off season is actually the woody area index (WAI) (Fang et al. 2019). WAI is an important vegetation biophysical parameter in canopy reflectance and land surface models (Bonan et al. 2002; Lawrence and Chase 2007). WAI accounts for an important component of light interception in forests (Sánchez-Azofeifa et al. 2009), and its presence significantly affects the snow surface albedo because of the absorption of nonphotosynthetic vegetation, the decrease of gaps in illumination, and the increase in shadows (Tian et al. 2004). However, optical remote sensing retrieval of WAI is difficult because of the mixed leaf and woody components in forest canopies. GEDI data products obtained during the leaf-off season provide a new data source for WAI. Therefore, it is necessary to assess the performance of the new GEDI LAI profile products during the leaf-off seasons.

As mentioned above, two questions are crucial: (1) What is the capability of various LiDAR indices and algorithms for the LAI profile estimation? (2) How does the new GEDI LAI profile product perform during both the leaf-on and leaf-off seasons? This study conducted a quantitative retrieval and validation of the LAI profiles derived from TLS, ALS, and GEDI. We aim to (1) develop a portable field measurement system for the LAI profile measurement, (2) obtain an optimal LiDAR index for the LAI profile estimation from ALS, and (3) validate the GEDI LAI profile product during both the leaf-on and leaf-off seasons. The study enhances the understanding of the capability of LAI profile estimation from different observation platforms and provides insight into the accuracy of the GEDI LAI profile product.

2. Materials and methods

The framework of this study is illustrated in Figure 1, including field data collection, evaluation of LiDAR indices, and the performance of GEDI LAI profile. The vertical LAI profile was first estimated in the field using an upward DHP attached to a portable measurement system in 2020 and 2021. Then, a suite of new LiDAR indices combining both the return number and return intensity were explored for the estimation of the LAI profile from ALS. As an illustration, the optimal ALS model was used to map the LAI profile over the study area. Finally, the GEDI LAI profiles were validated by the DHP, TLS, and ALS during the leaf-on and leaf-off seasons.

Figure 1. The framework for retrieval and validation of the LAI profiles derived from DHP, TLS, ALS, and GEDI.

2.1. Study area

The location of the study area (Figure 2) is within the Saihanba National Forest Park (SNFP) in the Hebei Province, North China (42°24′N, 117°18′E). The dominant species of SNFP are larch (Larix principis-rupprechtii), birch (Betula platyphylla) and Mongolian pine (Pinus sylvestris var. mongolica) (Zeng and Wang 2015). The deciduous needleleaf forest site with larch (Larix principis-rupprechtii) was chosen as an area of interest (AOI) and for field data collection (Figure 2). The terrain of the needleleaf forest site is generally flat (Figure S1).

Figure 2. The location of the study area. The left panel (a) shows the locations of all measurements and the GEDI footprints (Background Google Earth Image obtained by CNES/Airbus on June 29, 2018). The right panel shows two field pictures taken near a measurement tower (42°24′N, 117°18′E, yellow triangle in (a) during (b) leaf-on (September 6, 2020) and (c) leaf-off (April 24, 2021) seasons, respectively.

2.2. Field data collection

A total of 46 plots (Table 1) were collected in September 2020 (leaf-on, 23 plots) and April 2021 (leaf-off, 23 plots), corresponding to the leaf-off GEDI footprints (gray circles in Figure 2). Twelve plots (Table 1) were collected in September 2021 (leaf-on), corresponding to the leaf-on GEDI footprints (green circles in Figure 2) and the TLS site (orange rectangle in Figure 2). The plot is located at the center of the GEDI footprints and the TLS site, with a radius of 15 m (Figure 3a). The GEDI version 1 (V1) footprint locations were considered as the version 2 (V2) data were not available at the experiment design.

Figure 3. (a) The sampling strategy within a plot with five elementary vertical sampling points (VSPs) located at the center and four cardinal points, respectively; (b) a diagram of the portable measurement system at different heights (1.5–13.5 m) with 2 m interval; (c) field pictures were taken at 1.5–13.5 m with 2 m interval, respectively, on September 9, 2020. The upward DHP is within the red circle.

Table 1. DHP, TLS, ALS, and GEDI data used in the study. N is the number of plots for DHP and TLS and the number of footprints for GEDI, respectively.

	Acquisition date
	2018 (leaf-on)	2019 (leaf-off)	2020 (leaf-on)	2021 (leaf-off)	(leaf-on)
DHP	-	-	September 4–10 (N=23)	April 21–24 (N=23)	September 3–9 (N=12)*
TLS	-	-	September 5 (N=1)	-	-
ALS	September 5–17	-	-	-	-
GEDI	-	May 21 (N=23)	August 24 (N=11)	-	-

Note: *There are 11 plots in the GEDI footprints, one in the TLS site.

An extensible lifting mast was used for field vertical measurement. A Nikon D5100 camera and a fisheye converter were attached to the top of the mast. The calibration of the camera is carried out following the CAN-EYE software manual (Weiss and Baret 2017). The field measurements were made in near dusk or cloudy conditions (Fang et al. 2014; Zhang et al. 2020).

Within each plot, five elementary vertical sampling points (VSPs) were selected. The VSPs are located at the center and four cardinal points, respectively, and a brick was used to mark the VSP position for consecutive measurements (Figure 3a). The measurement heights range from 1.5 m to 13.5 m in 2 m intervals (Figure 3b–c). Four upward-facing DHP images were taken at each height in VSP to ensure the measurement and environment stability. Sample upward DHP images acquired from a VSP during different seasons are shown in Figure 4.

Figure 4. Upward DHP images acquired at different heights (1.5–13.5 m) with 2 m interval from the same vertical sampling point (VSP) in (a) leaf-on season (September 8, 2020) and (b) leaf-off season (April 25, 2021), respectively.

All DHP images at the same height of the five VSPs (i.e., 4 pictures/height × 5 VSPs = 20 pictures) were processed by the CAN-EYE software (version 6.49) to calculate the plot’s LAI. For the DHP images, only zenith angles varying from 0–60° were utilized to avoid the distortions of the edge. The zenith resolution and azimuth resolution were set as 2.5°, respectively. The LAI was estimated using a lookup-table technique in the CAN-EYE software (Weiss et al. 2000).

2.3. Derivation of LAI profile from TLS

A Riegl VZ-400 LiDAR system was used for TLS measurement at a 30 m × 30 m plot on 5 September 2020 (leaf-on, orange rectangle in Figure 2). The system employs a 1550 nm laser with a beam divergence of 0.35 mrad, and maximum range of 350 m (Zhu et al. 2018). The zenith angles ranged from 30° to 130°, and the azimuth angles ranged from 0° to 360°. A central scan position and four cardinal scan positions were used to minimize the occlusion effect. Eight cylindrical reflectance targets were marked as control points to register the TLS data acquired from multiple scans.

The LAI profile was derived from the TLS data using a voxel-based method which has been frequently used in similar studies (Beland et al. 2014; Beland, Kobayashi, and Goslee 2021; Hosoi and Omasa 2006; Zheng and Moskal 2012). The leaf angle distribution is assumed to be spherical and the voxel size is set to 0.05 m (Wei et al. 2020).

2.4. LAI profile estimation from ALS

2.4.1. ALS data and processing

A Riegl LMS-Q680i installed on the CAF-LiCHy platform was used to obtain the ALS data on September 5–17, 2018 (Pang et al. 2016). The detailed information about the LiDAR system is shown in Table 2.

Table 2. The parameters of airborne LiDAR system and flight.

Acquisition date	September 5–17, 2018
Sensor	Riegl LMS-Q680i
Beam wavelength (nm)	1550
Beam divergence (mrad)	0.5
Footprint (m) @1000 m altitude	0.5
Altitude (m)	1000
Pulse repetition frequency (kHz)	240
Point density (m^−2) @1000 m altitude	4
Scan angle	≤30°

The LiDAR return height was normalized by the digital elevation model (DEM) extracted from the ALS data. Returns from large off-nadir angles (>23°) were excluded to avoid potential uncertainties caused by the underestimated gap fraction (Liu et al. 2018). To reduce the impact of laser power, target area, atmosphere, and the range, the LiDAR intensity value was corrected to a user-defined standard range following (García et al. 2010; Gatziolis 2011; Hofle and Pfeifer 2007):

(1) $I_{corrected}=I_{raw}\times \frac{{R_0}^2}{R_s^2}$(1)

where $I_{corrected}$ and $I_{raw}$ are the corrected intensities and raw intensities, respectively. $R_0$ is the range between the ALS instrument and the return, and $R_s$ is the standard range, i.e. the average flying altitude ( = 1000 m in this study).

For the estimation of the LAI profile, the canopy and ground reflectance ratio (ρ_v/ρ_g) was first derived from the return energies of two adjacent footprints (Armston et al. 2013; Ni-Meister et al. 2010):

(2) $\frac{{\rho }_v}{{\rho }_g}=\frac{R_{v1}-R_{v2}}{R_{g2}-R_{g1}}$(2)

where $R_{v1}$, $R_{v2}$, $R_{g1}$, and $R_{g2}$ are the canopy and ground return energies from two adjacent vegetation and ground footprints, respectively.

In this study, ten pairs of adjacent footprints (red adjacent forest and ground paired circles in Figure 2), each containing a forest and ground return, were randomly selected. The ALS data were clipped with a diameter size of 25 m, similar to the GEDI. $R_v$ and $R_g$ are the mean intensity of the canopy and ground returns, respectively. To properly calculate $R_v$ for the forest footprints (marked by a number in Figure 2), the 0.5 m height threshold was used to mask the ground grass contribution (Armston et al. 2013; Chen et al. 2014). The average canopy and ground reflectance ratio (ρ_v/ρ_g) at 1550 nm was calculated as 1.04 ± 0.12, which was used in the subsequent analysis. More details about the calculation of the canopy and ground reflectance ratio from ALS is shown in Appendix A.

2.4.2. LAI estimation from LPI with the Beer-Lambert model

The theoretical basis of indirect LAI estimation is the Beer-Lambert law (Nilson 1971):

(3) $P\left(\theta \right)=e^{-G\left(\theta \right)\cdot LAI/cos\theta }$(3)

where $P\left(\theta \right)$ is the gap fraction and $G\left(\theta \right)$ is the leaf projection function.

The cumulative LAI is generally estimated from the LiDAR data following a modified Beer-Lambert method (Tang and Armston 2020; Tang et al. 2012):

(4) $LA{I}_{cum}\left(z\right)={\int }_{z_0}^z\frac{1}{G}\ast \frac{\mathrm{d}lnP\left(z\right)}{\mathrm{d}z}\mathrm{d}z$(4)

(5) $P\left(z\right)=1-\frac{R_v\left(z\right)}{R_v\left(0\right)}\frac{1}{1+\frac{{\rho }_v}{{\rho }_g}\hspace{0.17em}\frac{R_g}{R_v\left(0\right)}}$(5)

where $LA{I}_{cum}\left(z\right)$ is the cumulative LAI, $P\left(z\right)$ is the gap fraction, and z₀ is the height of the canopy bottom. $R_v\left(z\right)$, $R_v\left(0\right)$ and $R_g$ are the energy from the top of the canopy to the height z, from the canopy to bottom, and from the ground, respectively.

Based on Eqs. (4) and (5), the $LA{I}_i$ for layer $i$ is calculated as:

(6) $LA{I}_i=\frac{1}{G}\ast ln\left(1+\frac{{\rho }_g}{{\rho }_v}\ast \frac{1-LPI\left(i\right)}{LPI\left(i\right)}\right)$(6)

where the light penetration index (LPI) for layer $i$ is calculated as the ratio of the LiDAR return number (LPI_R) or intensity (LPI_Int) going through a specific layer (dz) to those reaching the top of the layer (z + dz) (Table 3). The ground and canopy reflectance ratio $({\rho }_g/{\rho }_v)$ was derived during the ALS data processing procedure (see 2.4.1).

Table 3. Various LiDAR indices derived from the ALS return number and intensity*.

Name	Equation	Description	References
Light penetration index from returns (LPIR)	$LP{I}_R\left(i\right)=\frac{{\sum }^{R_{\le z}}}{{\sum }^{R_{\le z+dz}}}$	The ratio of returns going through the layer dz to those reaching the height z + dz.	(Bouvier et al. 2015)
Light penetration index from return intensity (LPIInt)	$LP{I}_{Int}\left(i\right)=\frac{{\sum }^{I}n{t}_{\le z}}{{\sum }^{I}n{t}_{\le z+dz}}$	The ratio of return intensity going through the layer dz to those reaching the z + dz.	(Luo et al. 2018)
Light penetration index combining returns and return intensity (LPIRI)	$LP{I}_{RI}\left(i\right)=\frac{1}{2}\left(1+\frac{1}{LP{I}_R\left(i\right)}+\frac{{\rho }_g}{{\rho }_v}\ast \frac{1-LP{I}_{Int}\left(i\right)}{LP{I}_{Int}\left(i\right)}\right)$	Combine the LPIR and LPIInt at each height interval.	This study
Fractional cover index from returns (FCIR)	$FC{I}_R\left(i\right)=\frac{{\sum }^{R_{>z}}}{{\sum }^{R_{all}}}$	The ratio of returns above height z to all returns.	(Riano et al. 2004)
Fractional cover index from return intensity (FCIInt)	$FC{I}_{Int}\left(i\right)=\frac{{\sum }^{\hspace{0.17em}}In{t}_{>z}}{{\sum }^{\hspace{0.17em}}In{t}_{all}}$	The ratio of return intensity above height z to all return intensity.	(Hopkinson and Chasmer 2009)
Fractional cover index combining returns and return intensity (FCIRI)	$FC{I}_{RI}\left(i\right)=\frac{1}{2}\left(\frac{{\sum }^{R_{>z}}}{{\sum }^{R_{all}}}+\frac{{\sum }^{I}n{t}_{>z}}{{\sum }^{I}n{t}_{all}}\right)$	Combine the FCIR and FCIInt at each height interval.	This study
Above and below ratio index from returns (ABRIR)	$ABR{I}_R\left(i\right)=\frac{{\sum }^{R_{>z}}}{{\sum }^{R_{\le z}}}$	The ratio of returns above height z to those below height z.	(Sumnall et al. 2021)
Above and below ratio index from return intensity (ABRIInt)	$ABR{I}_{Int}\left(i\right)=\frac{{\sum }^{I}n{t}_{>z}}{{\sum }^{I}n{t}_{\le z}}$	The ratio of return intensity above height z to those below height z.	(Sumnall et al. 2016)
Above and below ratio index combining returns and return intensity (ABRIRI)	$ABR{I}_{RI}\left(i\right)=\frac{1}{2}\left(\frac{{\sum }^{R_{>z}}}{{\sum }^{R_{\le z}}}+\frac{{\sum }^{I}n{t}_{>z}}{{\sum }^{I}n{t}_{\le z}}\right)$	Combine the ABRIR and ABRIInt at each height interval.	This study

Note: * Subscripts R and Int denote ALS return and intensity, respectively, and RI indicates the new index combining both ALS return and intensity. $z$ represents the height, and $dz$ is the height interval.

To harness the advantages of LiDAR return number and intensity, a suite of new LiDAR indices calculated as a weighted linear average of the original indices, was explored for the LAI profile estimation. The combined light penetration index (LPI_RI) for layer $i$ is expressed as:

(7) $LP{I}_{RI}\left(i\right)=\frac{1}{2}\left(1+\frac{1}{LP{I}_R\left(i\right)}+\frac{{\rho }_g}{{\rho }_v}\ast \frac{1-LP{I}_{Int}\left(i\right)}{LP{I}_{Int}\left(i\right)}\right)$(7)

where the subscripts R and Int denote the ALS return and intensity, respectively, and RI indicates the new index combining both ALS return and intensity.

Accordingly, the $LA{I}_i$ for layer $i$ was calculated as:

(8) $LA{I}_i=\frac{1}{G}\ast ln\left(LP{I}_{RI}\left(i\right)\right)$(8)

2.4.3. LAI estimation from linear regression model

LAI was also estimated from the FCI and ABRI using the linear regression method (Luo et al. 2015; Richardson, Moskal, and Kim 2009; Zhao and Popescu 2009). In this case, the $LA{I}_i$ for layer $i$ was calculated as:

(9) $LA{I}_i=a_1\ast FCI\left(i\right)+b_1$(9)

(10) $LA{I}_i=a_2\ast ABRI\left(i\right)+b_2$(10)

where $a_1$, $a_2$, $b_1$ and $b_2$ are the regression coefficients. The FCI and ABRI are derived from the ALS return number, intensity, and both (Table 3).

Likewise, the combined fractional cover index (FCI_RI) for layer $i$ is expressed as:

(11) $FC{I}_{RI}\left(i\right)=\frac{1}{2}\left(\frac{\sum R_{>z}}{\sum R_{all}}+\frac{\sum In{t}_{>z}}{\sum In{t}_{all}}\right)$(11)

and the combined above and below ratio index (ABRI_RI) for layer $i$ is expressed as:

(12) $ABR{I}_{RI}\left(i\right)=\frac{1}{2}\left(\frac{\sum R_{>z}}{\sum R_{\le z}}+\frac{\sum In{t}_{>z}}{\sum In{t}_{\le z}}\right)$(12)

where R and Int denote the ALS return and intensity, respectively. The subscript all refer to all returns, and $z$ represents a height threshold. For example, $\sum R_{>z}$ is the number of returns above the height $z$.

For the linear regression method, all leaf-on DHP data (34 plots in September 2020 and 2021, 238 heights) were divided into two parts: one half (119 heights) for training the ALS model to estimate LAI and the other half (119 heights) for the validation purpose.

2.4.4. Mapping the vertical LAI profile

Nine LiDAR indices (Table 3) were investigated to estimate the LAI profile from 1.5 m to the maximum return height. All LiDAR returns were used to calculate the indices as previous studies have indicated that it would increase the point density throughout the canopy and may improve the estimation of LAI (Riano et al. 2004; Sumnall et al. 2016). Subsequently, the estimated ALS profile was validated using both the DHP and TLS estimates and compared with the GEDI product.

To illustrate the vertical LAI profile, the optimal ALS model was used to map the LAI profile for the 1 km × 1.5 km region (red box in Figure 2). The grid resolution and height intervals are 25 m and 5 m, respectively, similar to the GEDI footprint size and vertical resolution.

2.5. Validation of the GEDI LAI profile product

The cumulative vertical LAI profile was distributed in the GEDI Level 2 product (L2B) at 5 m vertical resolution (Dubayah et al. 2020). The GEDI L2B data product explored in the study were acquired on 21 May 2019 and 24 August 2020, respectively (gray and green circles in Figure 2). A total of 34 sample footprints were selected by the quality flag parameter (quality_flag = 1) for subsequent analyses (Table 1). The slopes of the GEDI footprints are less than 7° (Table S1). The GEDI V2 LAI profiles were explored to examine their performance (Dubayah et al. 2021). The GEDI V2 LAI profile was compared to those derived from the DHP, TLS, and ALS during the leaf-on season. During the leaf-off season, because of the lack of the TLS and ALS data, the GEDI V2 profile was compared to the DHP profile only.

2.6. Statistical analysis

All LAI profiles obtained from the DHP, TLS, ALS, and GEDI were inter-compared. The R², RMSE, and bias metrics were used to evaluate the LAI estimation.

(13) $R^2=1-\frac{{\sum }_{i=1}^n{\left({\hat{y}}_i-y_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\bar{y}\right)}^2}$(13)

(14) $RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left({\hat{y}}_i-y_i\right)}^2}{n}}$(14)

(15) $bias=\frac{1}{n}\sum _{i=1}^n\left({\hat{y}}_i-y_i\right)$(15)

where $y_i$ and ${\hat{y}}_i$ are the reference LAI and estimated LAI, respectively, and $\bar{y}$ is the mean reference LAI value.

3. Results

3.1. Comparison of ALS indices for the estimation of LAI profile

This section first presents the relationships between the DHP LAI and various ALS indices based on the training data, and then compares the LAI estimated from ALS with the DHP data. Finally, the LAI derived from different ALS indices is compared with the TLS measurements.

Figure 5 shows the linear relationships constructed between the DHP LAI profile and the FCIs and ABRIs from the training data. For the layered LAI (Figures 5a–f), the FCI_R shows a good correspondence with the DHP LAI (R² = 0.88, RMSE = 0.10, Figure 5a), while the ABRI indices show a weaker relationship (R² <0.83, Figures 5d–f). For the cumulative LAI (Figures 5g–I), all of the correlation values are similar to those for the layered LAI, with R² and RMSE values ranging from 0.75–0.89 and 0.19–0.29, respectively. The FCI_R and FCI_RI are similar, and both show a significant correspondence with the DHP measurements (Figures 5g–I), whereas the ABRI indices show a slightly weaker relationship, especially for the ABRI_Int (Figure 5k).

Figure 5. Relationships between the DHP layered and cumulative LAI (September 2020 and 2021) and various ALS (September 2018) indices from the training data. (a-c) and (g-i) the canopy fractional cover indices from return (FCI_R), intensity (FCI_Int), and both (FCI_RI), respectively; (d-f) and (j-l) the above and below ratio indices from return (ABRI_R), intensity (ABRI_Int), and both (ABRI_RI), respectively. The vertical color bar in (a) and (g) indicates the different heights for layered LAI and cumulative LAI, respectively. The index values are normalized to 0–1. The black line represents the regression line.

The FCI and ABRI models (Figure 5) were used to derive the vertical LAI profile, which was then compared with the layered DHP LAI. Figure 6 compares the layered LAI estimated from the ALS indices with those derived from DHP. The layered LAI from ALS is generally consistent with the DHP LAI (R² >0.80, RMSE <0.15, Table 4a). The LPI_RI model shows good correspondence with the DHP measurements (R² = 0.86, RMSE = 0.13) with a slope close to the 1:1 line (Figure 6c), while the LPI_R model slightly overestimates (bias = 0.06) and the LPI_Int slightly underestimates (bias = −0.07). Similar results were obtained from the FCI_R, FCI_Int, and FCI_RI models (R² ≥0.87, Figures 6d–f). The ABRI_Int model shows the lowest correspondence (R² = 0.81, RMSE = 0.14, Figure 6h), and the results are more scatterly distributed than those of the other models, especially at >11.5 m. The ABRI_R and ABRI_RI models perform similarly and show better results than the ABRI_Int (Figure 6g–i).

Figure 6. Comparison of the layered LAI estimated from ALS (September 2018) with the validation DHP data (September 2020 and 2021) during the leaf-on season. (a-c) 34 plots with the Beer-Lambert method; (d-i) 17 plots with the linear regression method.

Table 4. Comparison of statistics between the LAI profiles derived from various ALS indices (Table 3) and those from DHP (34 plots, 238 heights) during the leaf-on season.

Indices	Model	Modeling			Validation
		N	R^2	RMSE	N	R^2	RMSE	Bias
(a) Layered LAI (Figures 5 & 6)
LPIR	equation (6)(6) $LA{I}_i=\frac{1}{G}\ast ln\left(1+\frac{{\rho }_g}{{\rho }_v}\ast \frac{1-LPI\left(i\right)}{LPI\left(i\right)}\right)$(6)	-	-	-	238	0.86	0.15	0.06
LPIInt	equation (6)(6) $LA{I}_i=\frac{1}{G}\ast ln\left(1+\frac{{\rho }_g}{{\rho }_v}\ast \frac{1-LPI\left(i\right)}{LPI\left(i\right)}\right)$(6)	-	-	-	238	0.86	0.14	−0.07
LPIRI	equation (8)(8) $LA{I}_i=\frac{1}{G}\ast ln\left(LP{I}_{RI}\left(i\right)\right)$(8)	-	-	-	238	0.86	0.13	0.00
FCIR	$LAI=1.94\ast FCI+0.08$	119	0.88	0.10	119	0.88	0.11	−0.01
FCIInt	$LAI=2.37\ast FCI+0.09$	119	0.86	0.11	119	0.87	0.11	0.00
FCIRI	$LAI=2.13\ast FCI+0.08$	119	0.87	0.11	119	0.88	0.11	−0.01
ABRIR	$LAI=0.52\ast ABRI+0.09$	119	0.82	0.13	119	0.86	0.12	0.00
ABRIInt	$LAI=0.95\ast ABRI+0.10$	119	0.77	0.15	119	0.81	0.14	−0.01
ABRIRI	$LAI=0.66\ast ABRI+0.09$	119	0.81	0.13	119	0.83	0.13	−0.01
(b) Cumulative LAI (Figures 5 & 7)
LPIR	equation (6)(6) $LA{I}_i=\frac{1}{G}\ast ln\left(1+\frac{{\rho }_g}{{\rho }_v}\ast \frac{1-LPI\left(i\right)}{LPI\left(i\right)}\right)$(6)	-	-	-	238	0.89	0.27	0.06
LPIInt	equation (6)(6) $LA{I}_i=\frac{1}{G}\ast ln\left(1+\frac{{\rho }_g}{{\rho }_v}\ast \frac{1-LPI\left(i\right)}{LPI\left(i\right)}\right)$(6)	-	-	-	238	0.86	0.19	−0.28
LPIRI	equation (8)(8) $LA{I}_i=\frac{1}{G}\ast ln\left(LP{I}_{RI}\left(i\right)\right)$(8)	-	-	-	238	0.89	0.23	−0.10
FCIR	$LAI=2.36\ast FCI+0.24$	119	0.89	0.19	119	0.92	0.17	−0.03
FCIInt	$LAI=2.95\ast FCI+0.27$	119	0.86	0.22	119	0.91	0.19	−0.04
FCIRI	$LAI=2.64\ast FCI+0.26$	119	0.88	0.20	119	0.91	0.19	−0.05
ABRIR	$LAI=0.61\ast ABRI+0.28$	119	0.83	0.25	119	0.89	0.21	−0.05
ABRIInt	$LAI=1.06\ast ABRI+0.32$	119	0.75	0.29	119	0.84	0.26	−0.06
ABRIRI	$LAI=0.78\ast ABRI+0.30$	119	0.81	0.26	119	0.87	0.23	−0.06

Figure 7 shows the cumulative LAI derived from various ALS indices compared with the DHP measurements. The ALS outputs generally show good consistency with the DHP measurements (R² ≥0.84, RMSE ≤0.26, Table 4b). The LPI_RI and FCI models show similar performance (R² ≥0.89, Figures 7c–f), whereas the ABRI_Int model results in the lowest correspondence (Figure 7h). The LPI_R shows a small overestimation, especially at >7.5 m (Figure 7a). The LPI_Int systematically underestimates (bias = −0.28, Figure 7b), while the FCI and ABRI underestimate for LAI >2.

Figure 7. Comparison of the cumulative LAI estimated from ALS (September 2018) with the validation DHP data (September 2020 and 2021) during the leaf-on season. (a-c) 34 plots for Beer-Lambert; (d-i) 17 plots for linear regression method.

The statistics for the different ALS indices are summarized in Table 4. For the layered LAI (Table 4a), the FCI_RI displays the highest agreement with DHP (R² = 0.88), while the ABRI_Int produces the lowest agreement. For the cumulative LAI (Table 4b), the LPI_RI and FCIs perform similarly and show the best correspondence with the DHP (R² ≈0.90). The ABRI_Int model shows the lowest correspondence. The performance of the return intensity models is similar to that of the return number models, especially for the LPI and FCI models, whereas the ABRI_Int model performance is slightly inferior to the ABRI_R and ABRI_RI models.

Figure 8 compares the ALS LAI profiles derived from different LiDAR indices with the TLS measurements at 2 m height interval. In general, the LAI profiles estimated from ALS agree well with those from the TLS (R² ≥0.78). The LAI profile from the LPI_RI model shows the best correspondence with those from the TLS (R² = 0.95, RMSE = 0.05, Figure 8c), while the LPI_R shows a slightly weaker performance (R² = 0.88, RMSE = 0.07, Figure 8a). The ABRI_Int model shows the weakest correlation and slightly overestimates the TLS (R² = 0.78, RMSE = 0.09, Figure 8h). All the other models show moderate correlations (R² >0.80, RMSE ≤0.08).

Figure 8. Comparison of the layered LAI derived from ALS (September 2018) and TLS (September 2020) during the leaf-on season.

3.2. Vertical LAI distribution maps from ALS

The vertical LAI distribution maps of the study area were generated at 25 m resolution from ALS. The LAI profile maps derived from the LPI_RI method are presented in Figure 9. The LAI is close to 0 in the bare land and road areas, and the highest total LAI value reaches about 3.07. The LAI of each layer shows a geographic pattern generally similar to the total LAI map. However, different patterns exist for different layers. The highest LAI is distributed at 10–15 m (0.73 ± 0.35), while the lowest LAI is observed at 0–5 m (0.11 ± 0.11). The 5–10 m (0.38 ± 0.26) and >15 m (0.36 ± 0.31) layers show moderate LAI values. The vertical LAI maps were also generated using the other LiDAR indices (Figure S2). The results are generally similar to those of the LPI_RI. However, the FCI and ABRI methods show a slight underestimation of the total LAI compared to those of the LPI models.

Figure 9. Vertical LAI maps (1 km × 1.5 km) Figure 2derived from the ALS using LPI_RI (September 2018) over the study area (red box in Figure 2).The color bar indicates the LAI values.

3.3. Comparison of LAI profiles derived from different systems

The LAI profiles derived from the GEDI V2, ALS, TLS, and DHP are compared in Figure 10 and Table 5. The LAI values derived from ALS (1.98 ± 0.42) and DHP (1.92 ± 0.35) are generally similar during the leaf-on season (Figure 10a and Table 5). The LAI values derived from DHP (1.92 ± 0.35, Table 5) and GEDI (1.90 ± 0.29, Table 5) during the leaf-on season are clearly larger than those from the leaf-off season (DHP: 1.12 ± 0.17, GEDI: 0.68 ± 0.16, Table 5), especially for >5.5 m. At <5.5 m, the LAI profiles derived from all systems and seasons are nearly the same (~0.10–0.20, Figure 10 and Table 5) because of the dominant tree trunks at this height. For TLS, the maximum LAI profile appears at 13.5–15.5 m (Figure 10b). As shown in Appendix B, the TLS and DHP LAI profiles are highly correlated (R² = 0.97 and RMSE = 0.06, Figure B1).

Figure 10. The vertical LAI profiles derived from (a) DHP (September 2020 and 2021) and ALS (September 2018) during the leaf-on (34 plots) and leaf-off (23 plots in April 2021) seasons, and (b) TLS and GEDI V2 during the leaf-on (1 plot for TLS in September 2020, and 11 plots for GEDI in August 2020) and leaf-off (23 plots for GEDI in May 2019) seasons. The ALS LAIs were derived from LPI_RI. The line and shades represent the mean and standard deviation values, respectively. Subscripts “on” and “off” indicate the leaf-on and leaf-off seasons, respectively.

Table 5. The mean LAI at different heights and the total LAI obtained from DHP, TLS, ALS, and GEDI V2 at different seasons, respectively. The ALS LAIs were derived from LPI_RI. The values of the standard deviations are shown in the brackets. All LAI profile data are from Figure 10.

System	Acquisition date	Height (m)							Total LAI
		1.5–3.5	3.5–5.5	5.5–7.5	7.5–9.5	9.5–11.5	11.5–13.5	>13.5
DHP	September 2020 and 2021 (leaf-on)	0.12 (0.05)	0.12 (0.04)	0.11 (0.05)	0.14 (0.06)	0.20 (0.10)	0.32 (0.14)	0.92 (0.34)	1.92 (0.35)
	April 2021 (leaf-off)	0.11 (0.03)	0.08 (0.03)	0.08 (0.03)	0.07 (0.05)	0.09 (0.04)	0.09 (0.04)	0.59 (0.14)	1.12 (0.17)
TLS	September 2020 (leaf-on)	0.03	0.01	0.02	0.04	0.22	0.43	1.12	1.87
ALS	September 2018 (leaf-on)	0.04 (0.05)	0.04 (0.03)	0.08 (0.06)	0.15 (0.09)	0.26 (0.12)	0.37 (0.14)	0.99 (0.35)	1.98 (0.42)
GEDI	-	0–5	5–10			10–15		>15	-
	21 May 2019 (leaf-off)	0.07 (0.06)	0.27 (0.08)			0.18 (0.07)		0.16 (0.05)	0.68 (0.16)
	24 August 2020 (leaf-on)	0.14 (0.06)	0.33 (0.23)			0.72 (0.12)		0.71 (0.29)	1.90 (0.29)

3.4. The performance of the GEDI LAI profile data products

Figure 11 compares the LAI profiles derived from the GEDI V2 and DHP during the leaf-on season. The GEDI profile generally shows moderate correspondence with the DHP LAI (R² = 0.49, RMSE = 0.21, Figure 11a). At <15 m, the GEDI profile correlates well with the DHP value (R² = 0.77, RMSE = 0.14, Figure 11b). At >15 m, the GEDI values are generally lower than the DHP values (bias = −0.24, Figure 11c). The cumulative LAI from the GEDI presents good correspondence with the DHP LAI (R² = 0.90, RMSE = 0.23, Figure 11d). The LAI profile derived from GEDI also agrees very well with the TLS LAI during the leaf-on season (Figure B2).

Figure 11. Comparison of the layered LAI at (a) all heights, (b) 0–15 m, and (c) > 15 m and (d) the cumulative LAI derived from GEDI V2 (August 2020) with the DHP measurements (September 2021) during the leaf-on season (11 plots).

The LAI profiles derived from the GEDI V2 and DHP are also compared during the leaf-off season (Figure 12). The GEDI profile shows poor correspondence with the DHP LAI (R² = 0.04, RMSE = 0.10, Figure 12a). At <15 m, the GEDI values correlate moderately with the DHP (R² = 0.40, RMSE = 0.08, Figure 12b). However, at >15 m, the GEDI LAI significantly underestimates (~0.43) the DHP LAI (Figure 12c). For the cumulative LAI values, the correlation coefficient is better than that for the layered LAI (R² = 0.70), although the GEDI still shows a slight underestimation (bias = −0.13, Figure 12d).

Figure 12. Comparison of the layered LAI at (a) all heights, (b) 0–15 m, and (c) > 15 m and (d) the cumulative LAI derived from GEDI V2 (May 2019) and DHP (April 2021) during the leaf-off season (23 plots).

Figure 13 compares the cumulative LAI of ALS and GEDI V2 during the leaf-on season. In general, the LAI profiles estimated from the various LiDAR indices agree well with those from the GEDI (R² >0.85). The LPI_RI model presents the best results among all models (R² = 0.93, RMSE = 0.20, Figure 13c). The LPI_R model presents a systemic overestimation (bias = 0.22), whereas the LPI_Int presents a systemic underestimation (bias = −0.25). Similar performance is observed for the other models (Figures 13d–i). Compared with the LPI_RI model, the FCI and ABRI models show a slight underestimation of the cumulative LAI.

Figure 13. Comparison of the cumulative LAI derived from ALS (September 2018) and GEDI V2 (August 2020) during the leaf-on season (11 footprints).

4. Discussion

4.1. LAI profile measurement with DHP

In this study, a portable field measurement system was developed to estimate the vertical LAI profile using upward-facing DHP images at different heights (Figure 3). This system mitigates the challenges of tower-based measurements in representing large sampling areas and provides a promising method for validating the vertical LAI profile derived from ALS. The maximum height of the DHP system in this study is 13.5 m. A system utilizing a greater height would improve the LAI measurement in the upper canopy, but the system is prone to leaning at a higher level. Therefore, the maximum measurement height is a tradeoff between the canopy height and system stability.

4.2. Estimation of the LAI profile from LiDAR indices

Various LiDAR indices based on the return number and intensity were assessed in the LAI profile estimation. The LAI profiles estimated from ALS agree well with those from the DHP (R² ≥0.81) (Table 4). The ALS performance is comparable with those reported in other similar studies (Sumnall et al. 2016, 2021; Zhao and Popescu 2009). In this study, a new LiDAR index combining both the return number and return intensity was proposed to estimate the LAI profile. Other researchers have used only return number or intensity in similar studies (Arnqvist, Freier, and Dellwik 2020; Hopkinson and Chasmer 2009; Solberg et al. 2009; Sumnall et al. 2021). Compared with the original indices based only on the return number (R² = 0.86–0.88, RMSE = 0.11–0.15) or return intensity (R² = 0.81–0.87, RMSE = 0.11–0.14), the new combined index improves the LAI estimation (R² = 0.83–0.88, RMSE = 0.11–0.13, Figure 6). Among the combined indices, the performance of LPI_RI is superior to those of FCI_RI and ABRI_RI, indicating that the LPI_RI with the Beer-Lambert law may be more suitable for LAI estimation than the linear regression model with FCI_RI and ABRI_RI (Figures 6–7). The LPI_RI model is physically-based and less affected by site conditions and generally obtains the best results. The model may be useful for LAI mapping from ALS on a larger scale.

With the FCI models, the LAI values are slightly underestimated for LAI >2 (Figures 7d and 13d–f). This may possibly be attributed to the limited number of samples for LAI >2 in this study. Similar findings with FCI have been reported in other studies (Sumnall et al. 2016). The results from ABRI are similar to those from FCI, possibly because of FCI and ABRI are arithmetically exchangeable. Nevertheless, the FCI and ABRI regression models require sufficient amount of field data during the training process and are site-specific.

The LAI profiles estimated from ALS also show good agreement with TLS (Figure 8). The LPI_RI is superior to other indices (R² = 0.95, RMSE = 0.05). The small differences between ALS and TLS are mainly attributed to the different observation systems. The ALS on fixed wing aircraft may better characterize the upper layer when looking downward, while the lower canopy is easily obstructed by the upper canopy. Using both TLS and ALS may help alleviate the occlusion effects and improve the LAI estimation (Hosoi, Nakai, and Omasa 2010). In this study, the TLS data were available in only one site, and more TLS data are necessary to evaluate the performance of different LiDAR indices in the future.

All of the ALS first, intermediate, and last returns were used to calculate the LiDAR indices. Using all returns not only provides increased point sampling density throughout the canopy to the ground, but also minimizes the LAI bias usually induced by using only the first or last return (Hopkinson and Chasmer 2009; Lovell et al. 2003). Table 4 and Figure 8 show that slightly better results can be obtained using the return number than using the return intensity. In theory, the return intensity derived from radiometric information is more suitable for LAI estimation than the return number derived from point density (Yin, Cook, and Morton 2022; Yin et al. 2020). In practice, however, the intensity data are easily affected by atmospheric and aircraft conditions which greatly hampers the LAI estimation (García et al. 2010; Luo et al. 2019). Nonetheless, our results show that the LiDAR intensity data can be used to estimate the forest LAI profile satisfactorily (R² >0.75).

4.3. Performance of the GEDI LAI profile

The GEDI LAI product was compared to the DHP LAI at different seasons. During the leaf-on season, the GEDI layered LAI generally showed moderate correspondence with the DHP LAI profile (R² = 0.49, RMSE = 0.21, Figure 11a). However, the GEDI layered LAI values were smaller than the DHP values at >15 m (bias = −0.24, Figure 11c). On the other hand, the cumulative LAI from GEDI showed good correspondence with those of the DHP (R² = 0.90, RMSE = 0.23, Figure 11d). The main reason for the negative GEDI layered LAI (>15 m) discrepancy is attributed to the height mismatch between GEDI and DHP. The maximum height of the DHP system in this study is 13.5 m, and the layered LAI >15 m derived from DHP actually represents the LAI from 13.5 to the canopy top. Moreover, as shown in Appendix C, the layered LAI usually peaks at 13.5 − 15.5 m in the study area (0.39 ± 0.09, Figure C1 and Table C1). Another reason may be because of the increasing LAI spatial variation at higher canopy. The LAI spatial variation generally increases with canopy height with the increasing of large gaps in tall trees (Wang et al. 1992).

During the leaf-off season, the GEDI products underestimated the layered and cumulative LAI values, especially at the upper canopy layers (>15 m) (Figure 12). This underestimation is mainly attributed to the global static canopy and ground reflectance ratio (ρ_v/ρ_g = 1.5) adopted in the GEDI algorithm (Tang and Armston 2020). In the study area, the typical ρ_v/ρ_g value is around 1.04 ± 0.12 during the leaf-on season (Table A1). Other studies, such as the RAdiation transfer Model Intercomparison IV (RAMI-IV) project has indicated even smaller ρ_v/ρ_g values (0.71 ± 0.09, Table D1) during the leaf-off season (Appendix D). These values are significantly smaller than the GEDI assignment. Like the leaf-on season, the underestimation is also caused by the height mismatch between GEDI and DHP. Furthermore, the stronger ground return in the absence of leaves and the weak return energy from the small twigs and branches in the upper canopy may also bring retrieval uncertainties during the leaf-off season.

4.4. Woody area measurement during the leaf-off season

In the study, the LAI profile was determined from DHP and LiDAR at different seasons. During the leaf-on season, the LAI is actually the plant area index (PAI) since it consists of both leaf and woody components, while the LAI obtained during the leaf-off season is the woody area index (WAI) (Fang et al. 2019). As shown in this study, the GEDI data obtained during the leaf-off season can be used to represent the WAI. However, estimation of the WAI from LiDAR data during the leaf-off season is usually difficult because of the relatively lower WAI value and weaker canopy return. The DHP system experimented in this study provides a promising way to determining the WAI profile in the leaf-off season (Figure 4b). For the study area, the WAI values obtained from DHP are around 1.12 ± 0.17 (Table 5). The value is slightly lower than the WAI value (1.4) reported for a similar temperate larch forest in Japan (Hirata et al. 2007). Indeed, more studies are necessary to obtain the WAI profile which is required in various canopy reflectance and land surface models.

4.5. Limitations

In this study, the ALS data were acquired in 2018, and the DHP data in 2020–2021. The two-year difference between the ALS and DHP data may introduce errors in the evaluation of the ALS indices (Shao et al. 2019; Sumnall et al. 2021). However, this impact is small in this study because of the short growing season (from May to September) and the slow forest growth rate in the SNFP (Li et al. 2020). Moreover, there is no deforestation in the protected study area during the study period. In this study, the TLS data were available in only one site. For a robust analysis, more TLS data are necessary to demonstrate the suitability of the DHP measurement system.

5. Conclusions

In this study, we introduced a quantitative retrieval and validation of the LAI profile estimation from TLS, ALS, and GEDI over a deciduous needleleaf forest site in northern China. The LAI profile data were obtained in the field with a portable and extensible DHP measurement system at different seasons. With enhanced heights and stability, similar measurements can be conducted in more areas to obtain field reference data for the validation of the emerging LAI profile product. The LAI profile estimation from ALS can be improved with the proposed LiDAR index, which combines both the return number and intensity. The Beer-Lambert model with the combined light penetration index (LPI_RI) performs best when compared with the DHP, TLS, and GEDI data (R² ≥0.86), and is recommended for the LAI profile mapping over larger areas.

The GEDI cumulative LAI product performs well with the field measurements during the leaf-on season, but shows a negative bias during the leaf-off season, especially for the upper canopy layer. The negative bias is attributed to the global constant canopy and ground reflectance ratio (ρ_v/ρ_g = 1.5) and the height discrepancies between the GEDI and field measurements. Further improvement of the GEDI LAI algorithm can be attained with locally adjusted ρ_v/ρ_g value and ancillary data. Further validation studies are necessary for the GEDI LAI product across various biome types and landscape conditions, especially during the leaf-off season.

Highlights

• A portable field measurement system was developed to measure the LAI profile.

• The combined light penetration index (LPI_RI) is recommended for the LAI profile estimation from ALS.

• GEDI cumulative LAI performs well during the leaf-on season.

• GEDI cumulative LAI shows an underestimation during the leaf-off season.

Nomenclature

θ	=	Zenith angle
z	=	Height in the canopy
$G\left(\theta \right)$	=	Leaf projection function
J0	=	Transmitted laser pulse energy
$P\left(\theta \right)$	=	Gap fraction at zenith angle θ
$P\left(z\right)$	=	Gap fraction at height z
Iraw	=	The raw return intensity
Icorrected	=	The corrected return intensity
R0	=	The range between the ALS instrument and the return
Rs	=	The standard range (the average flying altitude)
ρg	=	The ground reflectance
ρv	=	The canopy reflectance
Rg	=	The laser energies from the ground return
Rv(0)	=	The laser energies from the canopy top to bottom
Rv(z)	=	The laser energies from the canopy top to height z
Rv1 and Rv2	=	The canopy return energies from two adjacent footprints
Rg1 and Rg2	=	The ground return energies from two adjacent footprints
R2	=	Coefficient of determination
ALS	=	Airborne Laser Scanning
DEM	=	Digital Elevation Model
DHP	=	Digital Hemispherical Photography
GEDI	=	Global Ecosystem Dynamics Investigation
GLAS	=	Geoscience Laser Altimeter System
LAI	=	Leaf Area Index
LAIcum	=	Cumulative Leaf Area Index
LAIi	=	Leaf Area Index for layer i
LiDAR	=	Light Detection and Ranging
PAI	=	Plant Area Index
RMSE	=	Root Mean Square Error
SNFP	=	Saihanba National Forest Park
TLS	=	Terrestrial Laser Scanning
VSP	=	Vertical Sampling Point
WAI	=	Woody Area Index
ABRI	=	Above and Below Ratio Index
ABRIR	=	Above and Below Ratio Index from Returns
ABRIInt	=	Above and Below Ratio Index from return Intensity
ABRIRI	=	Above and Below Ratio Index combining Returns and Intensity
FCI	=	canopy Fractional Cover Index
FCIR	=	canopy Fractional Cover Index from Returns
FCIInt	=	canopy Fractional Cover Index from return Intensity
FCIRI	=	canopy Fractional Cover Index combining Returns and Intensity
LPI	=	Light Penetration Index
LPIR	=	Light Penetration Index from Returns
LPIInt	=	Light Penetration Index from return Intensity
LPIRI	=	Light Penetration Index combining Returns and Intensity

Acknowledgments

This study was mainly supported by the National Natural Science Foundation of China (42171358) and the National Key Research and Development Program of China (2016YFA0600201) to H.F. The TLS data were provided by Dr. Jie Zou, Fuzhou University. The GEDI data are available from the NASA’s Land Processes Distributed Active Archive Center (LP DAAC). Drs. Lu Xu and Lixia Ma helped with the field data collection and the TLS data, respectively. The insightful comments provided by the anonymous reviewers helped improve the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The field measurement data in this study are available from the corresponding author, Y.W., upon reasonable request.

Supplementary Material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2023.2214987

Additional information

Funding

The work was supported by the National Natural Science Foundation of China [42171358]; National Key Research and Development Program of China [2016YFA0600201].

Appendix A.

Calculation of the canopy and ground reflectance ratio (ρ_v/ρ_g, 1550 nm) from ALS

The vegetation and ground return energies (R_v and R_g) at each footprint are calculated as (Ni-Meister, Jupp, and Dubayah 2001):

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

(A2) $R_g=J_0{\rho }_{g}P$(A2)

where $J_0$ is the laser pulse energy and $P$ is the canopy gap fraction.

Set $R_v^{\prime }=R_v/J_0$ and $R_g^{\prime }=R_g/J_0$, the above equations become:

(A3) $\frac{R_v^{\prime }}{{\rho }_v}+\frac{R_g^{\prime }}{{\rho }_g}=1$(A3)

(A4) $\frac{P}{1-P}=\frac{R_g^{\prime }{\rho }_v}{R_v^{\prime }{\rho }_g}$(A4)

Assuming constant ${\rho }_v$ and ${\rho }_g$ for two adjacent footprints, ${\rho }_v/{\rho }_g$ can be estimated from the $R_v$ and $R_g$ of two adjacent footprints (Ni-Meister et al. 2010):

(A5) $\frac{R_{g1}^{\prime }}{{\rho }_g}+\frac{R_{v1}^{\prime }}{{\rho }_v}=1$(A5)

(A6) $\frac{R_{g2}^{\prime }}{{\rho }_g}+\frac{R_{v2}^{\prime }}{{\rho }_v}=1$(A6)

(A7) $\frac{{\rho }_v}{{\rho }_g}=\frac{R_{v1}^{\prime }-R_{v2}^{\prime }}{R_{g2}^{\prime }-R_{g1}^{\prime }}=\frac{R_{v1}-R_{v2}}{R_{g2}-R_{g1}}$(A7)

In this study, ten pairs of footprints, each consisting of two adjacent footprints, were randomly selected (red circles in Figure 2). Figure A1 shows an example intensity image of three pairs of adjacent footprints (#1, #9 and #10). The return intensities and reflectance ratio (ρ_v/ρ_g) for the 10 pairs are shown in Table A1.

Figure A1. An example intensity image (in raw digital number) for 3 pairs of adjacent footprints (#1, #9 and #10) selected from Figure 2. Left: bare ground returns; right: forest returns. The average return intensity of each footprint is shown in the lower right corner.

Table A1. The average return intensities (in raw digital number) for ρ_v/ρ_g calculated from 10 pairs of adjacent footprints (Figure 2). The values of the standard deviations are shown in the brackets.

	#1	#2	#3	#4	#5	#6	#7	#8	#9	#10	Average
Rv1	0	0	0	0	0	0	0	0	0	0	0 (0)
Rv2	37.56	35.30	32.34	30.23	34.81	33.14	31.07	57.33	33.26	33.94	35.90 (7.42)
Rg1	127.86	99.61	108.12	119.27	118.15	105.73	117.65	195.61	102.10	155.07	124.90 (28.03)
Rg2	87.02	68.70	79.51	95.15	78.08	71.91	88.87	136.61	67.49	125.22	89.86 (22.34)
ρv/ρg	0.92	1.14	1.13	1.25	0.87	0.98	1.08	0.97	0.96	1.14	1.04 (0.12)

Appendix B.

Comparison of the DHP, TLS, and GEDI LAI profiles

  Figure B1 shows the LAI profiles derived from the TLS and DHP from 1.5 m to the canopy top. Both TLS and DHP profiles are nearly the same, while the TLS LAI are slightly lower than the DHP LAI at <5.5 m. In general, the TLS and DHP LAI profiles are highly correlated (R² = 0.97 and RMSE = 0.06). Figure B2 shows that the LAI profile derived from GEDI V2 agrees well with the TLS LAI for the layered (R² = 0.87, RMSE = 0.13) and cumulative (R² = 0.95, RMSE = 0.16) LAI values, respectively.

Figure B1. (a) The vertical LAI profiles derived from TLS (September 2020) and DHP (September 2021) during the leaf-on season. (b) The scatterplot between the TLS and DHP. The subscript “on” indicates that the data were obtained in the leaf-on season.

Figure B2. Comparison of the layered (a) and cumulative (b) LAI profiles derived from GEDI V2 (August 2020) and TLS (September 2020) during the leaf-on season.

Appendix C.

The 2 m LAI profiles derived from the GEDI waveform data

We derived the 2 m LAI profiles from the GEDI V2 waveform data using the same method as the GEDI standard product (Tang et al. 2014, 2012). The layered LAI was estimated as (Boucher et al. 2020):

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

(A9) $LAI\left(z\right)=\frac{ln\left(P\left(z+dz\right)\right)-ln\left(P\left(z\right)\right)}{G}$(A9)

where dz is the height interval for LAI calculation, dz = 2. The constants were set the same as in the GEDI standard product: G = 0.5, and ρ_v/ρ_g = 1.50.

The resultant LAI profiles show that the layered LAI peaks at 13.5 − 15.5 m in the study area during the different seasons (Figure C1 and Table C1).

Figure C1. The 2 m LAI profiles derived from the GEDI V2 waveform data during the leaf-on (11 footprints) and leaf-off (23 footprints) seasons.

Table C1. The mean and total LAI values obtained from the 2 m GEDI LAI profiles (Figure C1). The standard deviation values are shown in the brackets.

GEDI	Height (m)										Total LAI
	1.5–3.5	3.5–5.5	5.5–7.5	7.5–9.5	9.5–11.5	11.5–13.5	13.5–15.5	15.5–17.5	17.5–19.5	19.5–21.5
Leaf-on	0.07 (0.04)	0.05 (0.05)	0.08 (0.12)	0.15 (0.12)	0.21 (0.07)	0.32 (0.07)	0.39 (0.09)	0.32 (0.09)	0.21 (0.11)	0.12 (0.13)	1.94 (0.36)
Leaf-off	0.04 (0.02)	0.02 (0.01)	0.05 (0.01)	0.06 (0.02)	0.08 (0.02)	0.09 (0.03)	0.14 (0.04)	0.09 (0.04)	0.07 (0.06)	0.04 (0.04)	0.65 (0.16)

Appendix D.

Typical field measured canopy and ground reflectance ratio used in the RAdiation transfer Model Intercomparison IV (RAMI-IV)

Table D1. Typical field measured canopy and ground reflectance ratio (ρ_v/ρ_g) for different species used in the RAMI-IV during the leaf-off season* (Widlowski et al. 2015). The last row shows the mean (standard deviation) of the values.

Species	ρv/ρg
Maple	0.69
Birch	0.90
Alder	0.69
Ash	0.68
Linden	0.68
Poplar	0.74
Norway spruce	0.59
All	0.71 (0.09)

Note: *Data source: http://rami-benchmark.jrc.ec.europa.eu/HTML/RAMI-IV/RAMI-IV.php.