Assessment of Active LiDAR Data and Passive Optical Imagery for Double-Layered Mangrove Leaf Area Index Estimation: A Case Study in Mai Po, Hong Kong

Abstract

Remote sensing technology is a timely and cost-efficient method for leaf area index (LAI) estimation, especially for less accessible areas such as mangrove forests. Confounded by the poor penetrability of optical images, most previous studies focused on estimating the LAI of the main canopy, ignoring the understory. This study investigated the capability of multispectral Sentinel-2 (S2) imagery, airborne hyperspectral imagery (HSI), and airborne LiDAR data for overstory (OLe) and understory (ULe) LAI estimation of a multi-layered mangrove stand in Mai Po, Hong Kong, China. LiDAR data were employed to stratify the overstory and understory. Vegetation indices (VIs) and LiDAR metrics were generated as predictors to build regression models against the OLe and ULe with multiple parametric and non-parametric methods. The OLe model fitting results were typically better than ULe because of the dominant contribution of the overstory to the remotely sensed signal. A single red-edge VI derived from HSI data delivered the lowest RMSE of 0.12 and the highest R²_adj of 0.79 for OLe model fitting. The synergetic use of LiDAR metrics and S2 VIs performed best for ULe model fitting with RMSE = 0.33, R²_adj = 0.84. OLe estimation benefited from the high spatial and spectral resolution HSI that was found less confounded by the understory. In addition to their penetration attributes, LiDAR data could separately describe the upper and lower canopy, which reduced the noise from other components, thereby improving the ULe estimation.

1. Introduction

Leaf area index (LAI), defined as one-half the total leaf area per unit ground surface area [1,2,3], can reflect foliage density and the health condition of plants. LAI is a critical variable for reflecting canopy structure, estimating photosynthetic primary production, and measuring vegetation responses to global environmental change in process-based models [4]. LAI can be most accurately measured by direct methods such as harvesting, but it is extremely time-consuming, labor-intensive, and operation-constrained [5]. In recent decades, indirect non-contact LAI measurement techniques (based on the measurement of light transmission through canopies) have offered an alternative for LAI measurement. Coupled with remote sensing techniques, they have provided a timely and labor-efficient solution for LAI estimation over large spatial and long temporal scales [6,7,8,9,10]. Applying remote sensing approaches to estimate LAI can be categorized into physically based models and empirical statistical regression models. The former method involves the inversion of a canopy radiative transfer model (RTM) with remote sensing data for LAI retrieval. The latter method establishes a regression model between remote sensing variables and LAI. Empirical statistical regression models are the most commonly used and simplest method for LAI estimation, and they are highly flexible (which is critical when combining multi-source data) [11,12]. To apply a statistical approach, a regression method is selected to relate the observed LAI to several predictor variables. Parametric regression methods, which are defined by mathematically linear or non-linear functions, have been frequently employed [13,14,15]. These methods can describe the apparent relationship between LAI and predictors, but may also suffer from overfitting and poor accuracy when dealing with multiple predictors [12]. Meanwhile, non-parametric regression methods have been recommended for model fitting for various LiDAR metrics [16,17]. Non-parametric regression methods such as partial least squares (PLS), random forests (RFs), and artificial neural networks (ANNs) have been found to produce accurate and stable results for biophysical parameter retrieval over crop and forest areas [18,19,20,21,22]. Since the regression models tend to be sensor- and site-specific and the regression results also might change in time and space, the best regression method could vary case by case.

Vegetation indices (VIs) derived from multispectral and hyperspectral images are commonly used as predictor variables to build regression models against LAI. VIs are calculated from two or more spectral bands in order to enhance the sensitivity of spectral features to a vegetation property and reduce disturbance from the background or atmosphere [23]. Conventional broadband indices are mainly developed from the ratio, slope, or distance between the red reflectance valley and near-infrared (NIR) reflectance peak. They include the Normalized Difference Vegetation Index (NDVI), Simple Ratio (SR), Soil Adjusted Vegetation Index (SAVI), Enhanced Vegetation Index (EVI), and Atmospherically Resistant Vegetation Index (ARVI). Such indices have been found useful for estimating LAI in the literature [24,25,26]. Since novel sensors such as WorldView-2/3 and Sentinel-2 have been providing red-edge bands, a series of modified or novel red-edge indices have been proposed, including the Modified Red-Edge Normalized Difference (MNDre), Modified Simple Ratio Red Edge (MSRre), and Inverted Red-Edge Chlorophyll Index (IRECI). These indices have proven to be more useful than traditional VIs in canopy property monitoring [27,28,29,30]. VIs derived using red-edge bands have been found to alleviate the saturation effect and improve LAI estimation accuracy when compared with VIs constructed from traditional NIR bands [20]. Sentinel-2 data have been widely used to estimate local and regional scale LAI [31,32,33,34].

Narrowband indices are derived from hyperspectral data, including those derived from field spectroradiometers, spaceborne (EO-1 Hyperion), and airborne instruments (APEX, CASI, Headwall, etc.). They can be designed to emphasize subtle features along the spectral curve such as the reflectance slope in the red-edge region. Some narrowband indices have been found to perform consistently and robustly for LAI or chlorophyll content estimation. For example, the Non-Linear Vegetation Index (NLI) and Modified NLI, which exploit NIR and shortwave-infrared (SWIR) wavelengths [14,35], and the Transformed Chlorophyll Absorption in Reflectance (TCARI) and Modified Chlorophyll Absorption in Reflectance Index (MCARI), which are related to chlorophyll absorption [36,37,38]. When these narrowband indices were further modified and adopted for multispectral data, strong correlations with LAI were still obtained [20,29,39]. However, the limitation of VIs (such as the saturation effect) has been noted in other studies [40,41].

Most previous studies estimated the LAI of forests’ main canopy but neglected the understory. Photosynthesis mostly occurs in overstories that are closely related to biomass and productivity [42]. Moreover, understory LAI has remained inherently difficult to predict with passive optical images, as vegetation reflectance is mainly dominated by signals from the upper crowns in optical images. Understories aid soil nutrient cycling [42], provide habitat for wildlife [43], and contribute to the diversity in both species and structure [42]. Importantly, understories also make a contribution to biomass, net primary production and carbon stocks [44]. Despite the importance of understories, the exploration and quantification of their distribution has been limited in previous research. The estimation of LAI for both the understory and the main canopy would provide complete knowledge, enabling us to better understand forest communities, model biogeochemical cycles, and evaluate ecosystem function [8,9]. A few studies demonstrated that the presence of the understory can affect observed reflectance and LAI estimation [24,45,46], but the effect of background reflectance (understory vegetation, senescent leaves, and soil) is difficult to measure [47,48].

LiDAR data provide 3D information that provides opportunities to describe forest structure and improve LAI estimation. It is common to use geometric and radiometric metrics such as tree height, return elevation distribution statistics, and mean return intensity to estimate LAI [49,50]. Both canopy volume and foliage density metrics have been found useful for LAI estimation because LAI can be considered as a proportion of the crown volume containing foliage, estimated from the penetration length of lasers into the top canopy surface [51,52]. Penetration metrics indicate the canopy gap probability, which is non-linearly related to LAI according to the Beer–Lambert law [53]. Other useful metrics include the Last Echo Cover Index (LCI) and Solberg’s Cover Index (SCI), which quantify the rate of pulse penetration through the canopy [17,49,54]. Some studies attempted to combine LiDAR data and spectral images for biophysical parameter retrieval, in order to improve LAI estimation. For example, ref. [55] reported that the combination of VIs and LiDAR metrics could potentially better estimate biophysical parameters (e.g., canopy cover, water content, biomass) in tropical riparian zones, and [46,56,57] demonstrated that the combination of these two data sources improves LAI estimation accuracy for crops and forests. Meanwhile, ref. [58] suggested that LiDAR data can provide satisfactory estimates of LAI, but no significant improvement occurs when VIs are included for model fitting.

In addition to its utility for estimating the LAI of the main canopy, the penetrability of LiDAR data also enables the exploration of the understory layer. LiDAR data assisted in detecting the distribution of understory vegetation with accuracies of between 83% and 88% [43,59,60]. Furthermore, studies demonstrate the accurate stratification of the overstory and understory with LiDAR data based on canopy vertical profiles [61,62,63], which facilitate the separate estimation of overstory and understory LAI. Recently, ref. [64] improved overstory and understory LAI estimation with modified LiDAR indices in loblolly pine forest. Similarly, ref. [65] quantified understory LAI by combining full-waveform and point cloud LiDAR in a temperate forest. Although LiDAR data demonstrate the capacity to predict LAI in each forest component, there are limited studies that investigate the role of both optical and LiDAR data in a separated estimation of overstory and understory LAI. In boreal, deciduous, and planted forests, recent research has demonstrated that overstory and understory LAI/fractional cover are strongly related, mainly due to the reliance of the understory on light transmission through the overstory [66,67,68,69]. However, studies have rarely explored the overstory and understory structure and their relationship in coastal mangrove areas, which can be affected by varying tides, mudflat soils, disturbances, etc.

A double-layered structure was reported in the mangroves of Mai Po, Hong Kong [70], but no research has investigated the overstory LAI without neglecting the understory LAI. Based on the aforementioned research background, the overall objective of this study is to explore the feasible methods for LAI estimation in double-layered mangrove areas, so as to obtain a comprehensive understanding of mangrove growth. More specifically, this study attempts to investigate:

(1)

The ability of hyperspectral imagery, multispectral imagery, and LiDAR data to estimate the LAI of the overstory and understory, respectively.

(2)

The important features derived from remote sensing data to estimate each LAI category.

(3)

Possible factors such as remote sensing data, canopy structure, field measurement might affect the LAI estimation, especially for understories.

2. Materials and Methods

2.1. Study Area

The study area was the mangrove conservation core zone of Mai Po Nature Reserve, Ramsar wetland, (22°29′N–23°31′N, 113°59′E–114°03′E) in Hong Kong, China (Figure 1). A double-layered canopy structure was commonly found in the study area (Figure 1a). Overstories are commonly composed of mature arbor mangroves including Kandelia obovata (KAOB), Avicennia marina (AVMA), Aegiceras corniculatum (AECO), and Sonneratia (SO). The average canopy height of overstory mangroves varies from 5.5 m to 8.5 m for different species [61]. Understories were mainly found as shrub mangroves Acanthus ilicifolius (ACIL), which are commonly observed with a canopy height of between 1 m and 1.5 m [61]. A single-layered canopy structure was observed as two types. The first type consists of short mangroves such as shrub and small mangroves (Figure 1c) and another type is comprised of tall mangroves (most found as KAOB) without an understory beneath (Figure 1b).

2.2. Remote Sensing Data and Processing

Multispectral Sentinel-2 imagery, airborne hyperspectral imagery, and discrete return LiDAR data of the study area were used as the main remote sensing data.

The Sentinel-2 (S2) imagery includes 13 spectral bands in 3 spatial resolutions: 10 m in the blue, green, red, and near-infrared (NIR) bands, 20 m in 3 vegetation red-edge, a narrow NIR, and 2 SWIR bands, and 60 m in coastal aerosol and water vapor bands. An S2 Level-2A image captured on 26 December 2018 was downloaded from the Sentinel Scientific Data Hub (https://scihub.copernicus.eu/ (accessed on 15 January 2019)), which was the clearest and most concurrent image with the in-situ LAI measurement period. Atmospheric correction was conducted using Sen2Cor to obtain surface reflectance. The S2 image was reprojected into the Hong Kong 1980 Grid. Furthermore, 10 m was the finest spatial resolution of the S2 captures, and was also the same size as the LAI field-measured plots [30]. Therefore, the S2 image was resampled to a common 10 m spatial resolution using nearest neighbor resampling to preserve the original spectral values in ENVI 5.3.

Airborne hyperspectral imagery (HSI) was acquired using a headwall Nano-Hyperspec push-broom sensor on 13 March 2019. The HSI consisted of 271 spectral channels from 400 to 1000 nm with a full-width at a half maximum (FWHM) of 6 nm. The flying height was less than 300 m above the ground level, resulting in a spatial resolution of between 0.27 and 0.38 m. A total of 21 hyperspectral strips were acquired for the study area. The raw HSI was radiometrically corrected and then orthorectified in the WGS 1984 geographic coordinate system with inertial motion data of the flight campaign and a digital elevation model, using the SpectralView software (Headwall Photonics Inc., Bolton, MA, USA, 2018). After that, the HSI was registered to the Hong Kong 1980 Grid using digital aerial orthophotos (Lands Department, Hong Kong SAR Government) with manually selected tie points, obtaining an of accuracy less than 1 m in ArcMap 10.6. The spatial resolution of all HSI was homogenized to 0.5 m with nearest neighbor resampling. Atmospheric conditions were corrected to derive surface reflectance using the empirical line method (1) coupled with field reflectance measurements in ENVI.

Surface Reflectance = Gain × Radiance + Offset

(1)

Reflectance calibration tarpaulins in a white color (~70% reflectance) and a in black color (~5% reflectance) with a size of 3.5 m × 3 m that were placed within the flight area were measured by a spectroradiometer (ASD FieldSpec 4, ASD Inc., Boulder, CO, USA) throughout the flight. The spectroradiometer was calibrated using a Spectralon panel (SRT-99-100, Labsphere, North Sutton, NH, USA) before each recording. A total of 34 and 24 measured points were evenly distributed over the white and black tarpaulins, respectively. With the known reflectance of the white and black tarpaulins (as measured with the spectroradiometer) and their radiance extracted from the HSI, the gain and offset could be computed (1) [71].

Airborne LiDAR data were acquired on 21 March 2018 with a Riegl VUX-1LR laser scanner (RIEGL Inc., Horn, Austria) emitting a NIR 905 nm beam with a laser pulse rate of 600 kHz. The system can collect up to five discrete returns. The LiDAR data were acquired at an altitude between 155 and 270 m above the ground level, resulting in an average point density of 20 points per m². The LiDAR data were projected on the Hong Kong 1980 Grid and were further aligned with the digital aerial orthophotos. Subsequently, ground returns were classified with the Progressive Triangular Irregular Network (TIN) Densification filtering method [72], and then the return height was normalized by subtracting the TIN derived from ground returns in LAStools (https://rapidlasso.com/lastools/ (accessed on 15 January 2019)).

An object-based vertical stratification was applied to classify overstory and understory returns. The multiresolution segmentation algorithm in eCongition developer 9 (Trimble, Westminster, CA, USA) was applied to segment the 2 m spatial resolution canopy height model (CHM) into objects. After multiple trials and visual interpretation, the shape and compactness parameters were set to 0.1 and 0.8, respectively, and a scale parameter of 5 was determined. As a result, the CHM was segmented into 3088 objects. The frequency histogram of the return height was plotted in 0.2 m intervals for each object. A kernel regression was fit to smooth the distribution of histogram, where peaks would represent the canopy layers while valleys would represent the threshold for stratifying different layers [61,64]. The local minimum was searched as the stratified threshold in the bimodal distribution (Figure 1a), after which returns below the threshold were classified as understory and those beyond the threshold were regarded as overstory. Histograms showing unimodal distribution indicated that areas were comprised of a single-layer canopy, such as shrub (Figure 1c).

It is worth noting that although the airborne hyperspectral, LiDAR, and satellite-borne data were collected at different times, most of the mangroves within the study area are mature, and their canopies have not undergone drastic changes during recent years. Additionally, all data were collected during the dry season, reducing differences due to phenology. The airborne data represent the latest hyperspectral and LiDAR data available for Mai Po at the time of the study. A comparison of three data can be found in Supplementary File Figure S1.

2.3. Field Data Collection and Processing

A digital hemispherical photography (DHP) system comprised of a Nikon Coolpix P5000 digital camera and a 180° field of view fisheye lens (FC-E8) was used to measure LAI. Fieldwork was undertaken between 14 and 30 January in 2019 under overcast weather conditions. Mangrove fieldwork is particularly challenging due to the harsh environment (i.e., swampy mudflats and dense and prickly understories). As a result, most plots were located around the boardwalk. Sample plots were designed in 10 × 10 m quadrants at least 10 m away from the boardwalk and at least 10–20 m between each other (Figure 1). The understory vegetation had denser and larger leaves than the overstory, with understory canopy heights often taller than 1 m. Therefore, DHP images of both the overstory and understory (if present) were taken facing upward at five locations in each plot: one at the center of the plot and four at the centers of four sub-grids [9,73].

Plots were positioned using GPS (Garmin, Montana 600) with an average positioning uncertainty of +/− 3.65 m, and nearby landmarks were also recorded to help visually identify the plot location in the high-resolution images. Tree species in each plot were recorded, and tree heights of the understory at each sampling point were measured. Since AECO only occasionally appeared as one or two individual crowns within the accessible area, and SOAP only presented in the inaccessible seaward edge, there was no pure species plot for AECO and SOAP. Accessible single-layered shrub areas were not large enough to contain a 10 m × 10 m plot or Sentinel-2 pixel. As a result, a total of 51 plots including 34 KAOB plots and 17 AVMA plots were measured. Sample plots represented various canopy types, including mature and young mangroves with sparse to lush shrub beneath. One plot had no understory and seven plots had sparse and low (0.52–1.06 m) shrub ACIL beneath. The remaining 43 plots had a dense and relatively tall (1.1–1.65 m) thicket underneath.

The CAN-EYE software (https://www6.paca.inra.fr/can-eye/ (accessed on 15 January 2019)) was used to post-process the DHP images for effective LAI calculation. Undesirable objects in understory DHPs such as operator or obvious overstory components in the DHPs of the understory were masked out using the masking tool in CAN-EYE. Due to the large occurrence of mixed pixels at higher zenith angles, all classified binary images were limited to 0°–60°. Woody areas in the DHP images (such as trunks, branches, and twigs) were difficult to clearly distinguish from leaves, especially for AM, which was characterized by crooked trunks. Additionally, since is not possible to recognize if some leaves are present behind these woody areas, only classifying visible leaves could lead to a large underestimation of the actual LAI value [74]. Therefore, the woody areas were included in the vegetation class. With the inclusion of wood material, the plant area index (PAI) was derived, which was used as a proxy for LAI in this study. The statistical summary of measured LAI can be found in Supplementary File Table S1.

2.4. Generation of Vegetation Indices and LiDAR Metrics

A large set of different vegetation indices that were reported as useful for LAI estimation in previous studies were calculated as predictors [11,14,30,75,76]. VIs were categorized into the following groups: conventional NIR, atmospheric, red-edge, and SWIR indices (see

Table A1. A list of vegetation indices calculated for LAI estimation.

Index	Description	Formula	Reference
Conventional NIR indices
NDVI	Normalized Difference Vegetation Index	(NIR − Red)/(NIR + Red)	[88]
DVI	Difference Vegetation Index	NIR − Red	[89]
CIg	Green Chlorophyll Index	(NIR/Green) − 1	[90]
WDRVI	Wide Dynamic Range Vegetation Index	(a × NIR − Red)/(a × NIR + Red)	[91]
NR	Normalized Red Band Index	Red/(NIR + Red + Green)	[92]
NNIR	Normalized NIR Band Index	NIR/(NIR + Red + Green)	[92]
NLI	Non-Linear Index	(NIR2 − Red)/(NIR2 + Red)	[93]
RDVI	Renormalized Difference Vegetation Index	$(NIR-Red)/\surd (NIR+Red)$	[94]
SPVI	Spectral Polygon Vegetation Index	$0.4\times 3.7\times (\rho 800-\rho 680)-1.2\times \sqrt{{(\rho 530-\rho 670)}^2}$	[95]
Atmospheric indices
ARVI	Atmospherically Resistant Vegetation Index	$\frac{NIR-Red-(Blue-Red)}{NIR+Red-(Blue-Red)}$	[96]
EVI	Enhanced Vegetation Index	$\frac{2.5\times (NIR-Red)}{NIR+6\times Red-7.5\times Blue+1}$	[97]
GARI	Green Atmospherically Resistant Index	$\frac{NIR-(Green-\gamma (Blue-Red\left)\right)}{NIR+(Green-\gamma (Blue-Red\left)\right)}$	[41]
VARIg	Visible Atmospherically Resistant Index	$\frac{(Green-Red)}{(Green+Red-Blue)}$	[98]
Red edge indices
CIre	Red-Edge Chlorophyll Index	(B7/B5) − 1	[90]
WDRVI-re	Red-Edge Wide Dynamic Range Vegetation Index	(a × Red Edge − Red)/(a × Red Edge + Red)	[99]
PSRI	Plant senescence Reflectance Index	$(\rho 650-\rho 550)/\rho 750$	[100]
MTCI	MERIS Terrestrial Chlorophyll Index	$(\rho 750-\rho 708)/(\rho 708-\rho 680)$	[101]
MCARI	Modified Chlorophyll Absorption Ratio Index	$\frac{\rho 700–(\rho 1640-\rho 2150)}{\rho 860+\rho 1640-\rho 2130}$	[102]
MCARI2	Modified Chlorophyll Absorption Ratio Index Improved	$\frac{1.5\times \left(1.2\times \rho 800-\rho 670\right)-2.5\times (\rho 800-\rho 550))}{\sqrt{{(2\times \rho 800+1)}^2}-\left(6\times \rho 800-5\times \sqrt{\rho 670}\right)-0.5}$	[75]
TCARI	Transformed Chlorophyll Absorption Reflectance Index	$3\times \left(\rho 700-\rho 670\right)-0.2\times (\rho 700-\rho 550)\times \frac{\rho 700}{\rho 670}$	[103]
TCARI2	Transformed Chlorophyll Absorption Reflectance Index 2	$3\times \left(\rho 705-\rho 770\right)-0.2\times (\rho 705-\rho 550)\times \frac{\rho 705}{\rho 770}$	[37]
TVI	Triangular Vegetation Index	$\frac{120\times \left(\rho 750-\rho 550\right)-200\times (\rho 670-\rho 550)}{2}$	[77]
MTVI2	Modified Triangular Vegetation Index—Improved	$\frac{1.5\times \left(1.2\times (\rho 800-\rho 550\right)-2.5\times (\rho 670-\rho 550))}{\sqrt{{\left(2\times \rho 800+1\right)}^2-(6\times \rho 800-5\times \sqrt{\rho 670})-0.5}}$	[75]
IRECI	Inverted Red-Edge Chlorophyll Index	$\frac{\rho 783-\rho 665}{\rho 705/\rho 740}$	[29]
S2REP	Sentinel-2 Red-Edge Position Index	$705+35\times \frac{\left(\frac{\rho 783+\rho 665}{2}\right)-\rho 705}{\rho 740-\rho 705}$	[29]
NDre1m	Modified Red-Edge Normalized Difference 1	$\frac{\rho 750-\rho 705}{(\rho 750+\rho 705-2\times \rho 445)}$	[28]
NDre2m	Modified Red-Edge Normalized Difference 2	$\frac{\rho 783-\rho 705}{(\rho 783+\rho 705-2\times \rho 445)}$	[28]
SRre1	Red-Edge Simple Ration 1	$\frac{\rho 750–\rho 445}{\rho 700-\rho 445}$	[28]
SRre2	Red-Edge Simple Ration 1	$\frac{\rho 783–\rho 445}{\rho 700-\rho 445}$	[28]
MSRren	Modified Simple Ratio Red-Edge Narrow	$\frac{\left(\frac{\rho 865}{\rho 700}\right)-1}{\sqrt{\left(\frac{\rho 865}{\rho 700}\right)+1}}$	[39]
BSItian	Tian’s three-band spectral index	$\frac{\rho 605-\rho 521-\rho 682}{\rho 605+\rho 521+\rho 682}$	[104]
SWIR indices
NDII	Normalized Difference Infrared Index	$\frac{\rho 819–\rho 1649}{\rho 819+\rho 1649}$	[105]
NBR	Normalized Burn Ratio	(NIR − SWIR)/(NIR + SWIR)	[106]
NMDI	Normalized Multi-Band Drought Index	$\frac{\rho 860–(\rho 1640-\rho 2150)}{\rho 860+\rho 1640-\rho 2130}$	[107]

). The individual S2 pixel that intersected with the central point of the plot was selected to calculate VIs, as both the pixel size of the resampled S2 data and the plot size were equal to 10 m. To match this extent, for the 0.5 m HSI, a 21 × 21 pixel window was used to represent the plot. The average reflectance of the window was used to calculate the VIs.

LiDAR returns within a 10 × 10 m square around the plot center were extracted to derive LiDAR metrics. Widely used standard metrics that statistically describe return height were generated (

Table A2. A list of LiDAR metrics calculated for leaf area index estimation.

Abbreviation of LiDAR Metric	Metric Details	Reference
Statistical metrics
zmax	Maximum height
zmean	Mean height
zsd	Standard deviation of height distribution
zskew	Skewness of height distribution
zkurt	Kurtosis of height distribution
pzabovezmean	Percentage of returns above mean height
pzabove 2	Percentage of returns above 2 m
zqx (x = 5, 10, 15… 95)	Xth percentile (quantile) of the height distribution
zpcumx (x = 1, 2, 3… 9)	For each plot, the LiDAR height range was divided into 10 equal intervals. Cumulative percentage of return in the Xth layer	[108]
den	Point density of all returns
den_o	Point density of overstory returns
den_u	Point density of understory returns
Penetration related metrics
CRR	Canopy relief ratio,	[109]
cover_o	Overstory canopy cover is the proportion of overstory returns over vegetation returns
cover_u	Understory canopy cover is the proportion of understory returns over vegetation returns
LCI	Last Echo Cover Index	[82]
SCI	Solberg’s Cover Index	[49]
SCI_O	Modified Solberg’s Cover Index using overstory returns
SCI_U	Modified Solberg’s Cover Index using understory returns
LPIFR	${LPI}_{FR}=\frac{\sum R_{O\left(first\right)}}{\sum R_{frist}}$	[110]
ABRI	$ABRI=\frac{\sum R_O}{\sum R_U}$	[111]
Pgap_weight	$P_{gap\_weight}=1-\frac{\sum 1/NoR}{N_{total}}$	[112]
pground	$P_{ground}=\frac{\sum R_{ground}}{\sum R_{total}}$
BL	$BL=\frac{\left(\frac{\sum I_{<T\left(single\right)}}{\sum I_{\left(total\right)}}\right)+\sqrt{\frac{\sum I_{<T\left(last\right)}}{\sum I_{\left(total\right)}}}}{\left(\frac{\sum I_{\left(frist\right)}+\sum I_{\left(single\right)}}{\sum I_{\left(total\right)}}\right)+\sqrt{\frac{\sum I_{\left(intermediate\right)}+\sum I_{\left(last\right)}}{\sum I_{\left(total\right)}}}}$	[110]
BLcanopy	$BLcanopy=\frac{\left(\frac{\sum I_{O\left(frist\right)}+\sum I_{O\left(single\right)}}{\sum I_{\left(total\right)}}\right)+\sqrt{\frac{\sum I_{O\left(intermediate\right)}+\sum I_{O\left(last\right)}}{\sum I_{\left(total\right)}}}}{\left(\frac{\sum I_{\left(frist\right)}+\sum I_{\left(single\right)}}{\sum I_{\left(total\right)}}\right)+\sqrt{\frac{\sum I_{\left(intermediate\right)}+\sum I_{\left(last\right)}}{\sum I_{\left(total\right)}}}}$	[64]
BLunder	$BLunderstand=\frac{\left(\frac{\sum I_{U\left(Frist\right)}+\sum I_{U\left(Single\right)}}{\sum I_{\left(Total\right)}}\right)+\sqrt{\frac{\sum I_{U\left(Intermediate\right)}+\sum I_{U\left(Last\right)}}{\sum I_{\left(Total\right)}}}}{\left(\frac{\sum I_{\left(Frist\right)}+\sum I_{\left(Single\right)}}{\sum I_{\left(Total\right)}}\right)+\sqrt{\frac{\sum I_{\left(Intermediate\right)}+\sum I_{\left(Last\right)}}{\sum I_{\left(Total\right)}}}}$	[64]

$R_{frist}$: first returns; $R_O$: overstory returns, $R_U$: understory returns; $R_{ground}$: ground return; $N_{total}$: total number of returns; $NoR$: number of returns; I: intensity

). Additionally, various metrics describing the penetration rate (i.e., canopy cover) were calculated and some of them were modified with consideration of the overstory and understory.

2.5. Statistical Methods to Estimate LAI and Validation

Both univariate and multivariate regression methods were examined for overstory effective leaf area index (OLe) and understory effective leaf area index (ULe) fitting. Univariate methods including linear and non-linear (exponential and logarithm) regression were applied so that the useful variables (VIs or LiDAR metrics) and their relationship with overstory and understory LAI could be explored. Multivariate regression methods including parametric multiple linear regression (MLR), non-parametric method partial least squares regression (PLSR), and random forest regression (RFR) were employed to understand whether the combination of data could improve the LAI model fitting.

Every additional predictor in a multivariate regression model might decrease the RMSE and increase the R², and also increase the risk of overfitting [30]. Therefore, several rules were applied to select variables for multivariate regression:

(1)

The maximum number of predictor variables from single-source data (i.e., VIs from S2 or HSI and LiDAR metrics) was restricted to four, because the RMSE was observed to be stable with more than four predictors. Therefore, the maximum number of combination predictors (i.e., four VIs + four LiDAR metrics) would be eight.

(2)

To include an additional predictor variable, the RMSE of the model had to be decreased by at least 2% when compared with the original value.

Predictors were then selected by the recursive feature elimination (RFE) feature selection algorithm based on RF. RFE-RF iteratively eliminates redundant variables to fit the RFR model until the best root-mean-square error (RMSE) is obtained. Feature selection was performed with S2 VIs, HSI VIs, and LiDAR metrics separately, to identify important predictors. Subsequently, the selected VIs and LiDAR metrics were input into the REF-RF again to determine the optimal combination of predictors.

Regression models were fitted and assessed through five-fold cross-validation (CV) due to the limited number of field samples. The coefficient of determination (R²) and the RMSE of cross-validation were used to assess the model fit. Because adding additional predictors tends to increase the R², the adjusted R² (R²_adj) was also calculated to better evaluate models with a different number of input variables. The regression fitting was implemented in R (https://www.r-project.org/ (accessed on 15 January 2019)) with additional packages: “pls”, “randomForest”, and “caret”. Figure 2 summarizes the workflow of above-mentioned data process and analysis.

3. Results

Table 1. OLe and ULe regression model fitting results. The best model fitting result of each regression method is highlighted in bold. For multiple regression, variables are listed in descending order of importance, indicated by feature selection processing.

	Data	Method	Variables	RMSE	R2	R2adj
OLe
Univariate regression	HSI	Exp	IRECI a	0.12	0.79	0.79
	S2	Exp	NDII	0.15	0.68	0.68
	LiDAR	Exp	zsd	0.17	0.53	0.52
Multivariate regression	HSI	MLR	IRECI, RDVI, GARI	0.16	0.79	0.77
	S2	RFR	NBR, TCARI2, Cig	0.20	0.70	0.68
	LiDAR	RFR	zsd, LCI, zkurt	0.21	0.63	0.61
	HSI + LiDAR	PLSR	IRECI, RDVI, GARI, zsd, LCI b	0.16	0.81	0.78
	S2 + LiDAR	PLSR	NBR, TCARI2, Cig, zsd, LCI, zkurt	0.16	0.78	0.76
ULe
Univariate regression	HSI	Exp	MSRren	0.37	0.41	0.40
	S2	Exp	NLI	0.30	0.57	0.56
	LiDAR	Exp	Pground c	0.24	0.57	0.56
Multivariate regression	HSI	MLR	MSRren, MTCI, SRre2,	0.66	0.44	0.40
	S2	RFR	BSItian, MTVI2, MCARI2	0.47	0.64	0.61
	LiDAR	RFR	den_u, zpcum6, zq40, zpcum7, zkurt	0.40	0.71	0.67
	HSI + LiDAR	RFR	zpcum7, zpcum6, zq40, den_u, zkurt, MTCI, SRre2, MSRren	0.40	0.70	0.65
	S2 + LiDAR	RFR	zpcum7, zpcum6, MTVI2, MCARI2, den_u d	0.33	0.84	0.82

Abbreviations: HSI, hyperspectral imagery; S2, Sentinel-2; Exp, exponential regression; MLR, multiple linear regression; RFM, random forest regression; PLSR, partial least squares regression. See Table A1 and Table A2 for explanation of variables. The superscript a, b, c, and d suggest the predictors and methods used to generate LAI map in Figure 3.

presents the models achieving the highest R² and lowest RMSE by using each remote sensing data source, univariate regression, and multivariate regression, respectively. Important features determined by the feature selection process for multivariable regressions are listed by decreasing order of importance in Table 1. For OLe estimation, the top three VIs come from different groups, which indicates that the decorrelation among variables may benefit the OLe model fitting. In our case, the optimal combination of IRECI (Inverted Red-Edge Chlorophyll Index, red-edge VI), RDVI (Renormalized Difference Vegetation Index, NIR indices), and GARI (Green Atmospherically Resistant Index, atmospheric indices), derived from HSI and NBR (Normalized Burn Ratio, SWIR indices); and TCARI2(Transformed Chlorophyll Absorption Reflectance Index 2, red-edge VI) and Cig (Green Chlorophyll Index, NIR indices), derived from S2 data, were selected as input. LiDAR metrics reflecting the LiDAR return distribution of canopy structure, including zsd (standard deviation of height distribution), LCI (Last Echo Cover Index), and zkurt (kurtosis of height distribution), were selected. Both HSI VIs and S2 VIs were considered more important than LiDAR metrics in OLe estimation if combining both spectral and spatial features together. For ULe estimation, all selected VIs were from the red-edge VI group. This result may reflect that the synergy of multiple red-edge spectral features works better for the understory with weak reflectance signals. LiDAR metrics that reflected the distribution of LiDAR returns in lower canopy layered such as zpcum7 (cumulative percentage of return in the seventh layer), zq40 (40th percentile of the height distribution), den_u (point density of understory returns) were selected. When combing VIs and LiDAR metrics to estimate ULe, LiDAR metrics were considered more important than VIs as the penetration property of LiDAR data is able to capture more signals from the understory.

For OLe estimation, using only the HSI red-edge VI IRECI produced the best model fitting result, with the highest R²_adj of 0.79 and the lowest RMSE of 0.12. IRECI alone performed better than using multiple HSI VIs or combining spectral VIs and LiDAR metrics. The individual S2 VI NDII (Normalized Difference Infrared Index), which incorporates the NIR and SWIR bands, could also produce a good model fitting result, with an R²_adj of 0.68 and RMSE of 0.15. The combination of S2 VIs and LiDAR metrics significantly improved the R²_adj to 0.76, which was better than using either S2 VIs (R²_adj = 0.68) or LiDAR metrics (R²_adj = 0.61) alone.

Overall, ULe models were observed with larger RMSEs than OLe models. Applying LiDAR metric P_ground (percentage of ground return) and S2 NLI (Non-Linear Index of Red and NIR bands) alone generated a low RMSE of 0.24 and 0.3, respectively, with an R²_adj = 0.56. It was observed that although increasing the number of predictors (such as using multiple S2 VI or LiDAR metrics) could improve R²_adj to 0.61–0.67, the RMSE increased to 0.4–0.47. The combination of S2 red-edge VIs and LiDAR metrics produced the best R²_adj of 0.82 with an RMSE of 0.33. Although HSI VIs performed best in OLe estimation, they seem less useful in ULe estimation compared with S2 VIs and LiDAR metrics.

As shown in Figure 3, the best univariate and multivariate regression models were applied to generate maps of Ole and Ule, respectively. Maps generated by the multivariate regression model (Figure 3b,d) are more homogenous than the univariate ones for both OLe and ULe (Figure 3a,c). A similarity in spatial distribution of both OLe and ULe values was observed, although different predictors and methods were employed. For double-layered mangroves in region 1 and 2 (region are shown in Figure 1), OLe values are moderate (between 2 and 4), where OLe values of AMVA are slightly higher than KAOB. Meanwhile, the understories beneath them have relatively high ULe (ranging from 2.1 to 2.9 or 2.1 to 4.2). The OLe distribution in region 3 is complex. Young mangroves near the seaside and along the fringe have the highest OLe (3–5), while the adjacent KAOB region 2 is found to have lower OLe (less than 2). The ULe in region 3 is relatively low, showing a gradual decrease toward the seaside. Prediction bias can be indicated by comparing the measured and predicted LAI values. The OLe map produced by the combination of predictors has a slightly higher OLe compared to the map predicted by the single HSI VI, especially for double-layered areas. Areas with low OLe (below 1.2) may be overestimated while areas with OLe higher than 2 may be underestimated. The estimated LAI of areas where ULe was less than 2 has a larger bias than where ULe was 3 or higher. For high ULe (2.5–3.5), the bias between observation and estimation could be even larger if using the single LiDAR metric P_ground.

4. Discussion

4.1. Estimation of OLe and ULe

OLe and ULe were predicted with varying degrees of accuracy by using different data sources. OLe exhibited better relationships with predictors because both spectral reflectance and LiDAR pulses were dominated by the upper canopy. HSI was found as the most useful data to estimate OLe in the experiments, which may benefit from its high spatial and spectral resolution. The high spatial resolution HSI provided more detail of the mangrove stand. For instance, the overstory canopy, gaps, and shadow could be distinguished (Figure S1). Additionally, the narrow bands of the HSI are thought to be less sensitive to the background when compared to broad bands [77]. As a result of these factors, in the HSI, the spectral properties of the overstory might be less confounded by background information.

Table 2. A comparison of the correlation coefficients between field-measured leaf area indices (OLe: overstory leaf area index; ULe: understory leaf area index) and Sentinal-2 (S2) spectral bands/selected hyperspectral (HSI) bands.

S2 Bands	B1 443 nm	B2 490 nm	B3 560 nm	B4 665 nm	B5 705 nm	B6 740 nm	B7 783 nm	B8 842 nm	B8A 865 nm	B9 945 nm	B11 1610 nm	B12 2190 nm
OLe	−0.154	−0.551 ***	−0.528 ***	−0.606 ***	−0.347 *	0.683 ***	0.730 ***	0.736 ***	0.713 ***	0.314 *	−0.563 ***	−0.678 ***
ULe	−0.495 ***	−0.645 ***	−0.592 ***	−0.609 ***	−0.399 **	0.230	0.338 *	0.474 ***	0.400 **	0.385 **	0.037	−0.160
HSI bands	444 nm	491 nm	560 nm	666 nm	706 nm	740 nm	782 nm	842 nm	864 nm	945 nm
OLe	−0.564 ***	−0.556 ***	0.030	−0.527 ***	0.168	0.773 ***	0.837 ***	0.840 ***	0.838 ***	0.847 ***
ULe	−0.334 *	−0.307 *	−0.286 *	−0.226	−0.255 .	0.057	0.168	0.186	0.192	0.206

The significance of correlation is indicated by * symbol, where p-value ranges to symbols is listed as: 0–0.001 to ***; 0.001–0.01 to **; 0.01–0.05 to *; 0.05–0.1 to .; 0.1–1.0 no symbol.

compares the correlation coefficient between field-measured LAI (OLe and ULe) and S2 spectral bands and selected HSI spectral bands. For easy comparison, hyperspectral bands with wavelengths closest to the central wavelength of each S2 band were selected. As the wavelength range of HSI is 400 nm–1000 nm, SWIR bands are not available in HSI. In the visible wavelengths range from 443 nm to 705 nm, S2 and HSI spectral bands show similar correction levels to OLe. However, in the red-edge to NIR bands (740 nm–945 nm), HSI spectral bands produced better correlation coefficients (0.773–0.847) to OLe when compared with S2 spectral bands.

The ULe model fitting obtained less accuracy compared with OLe, which may result from inherent factors in both remote sensing data and environment. Understory information could be partially blocked by the upper canopy or affected by shadow, leading to weaker relationships between ULe and predictors. Therefore, it is particularly useful to combine spectral and spatial features to provide as much information as possible for the understory. LiDAR data are more adaptable for ULe prediction as a result of their penetration characteristics. S2 VIs appeared to be more helpful than HSI VI for ULe estimation. As shown in Table 2, visible bands (443 nm–705 nm) and red-edge to NIR bands (783 nm–945 nm) of S2 show significant correlations to Ule, with correlation coefficients between 0.338 and 0.645. Although the blue, green, and red-edge bands of HSI are negatively correlated to ULe, the correlation and significance are weaker than the corresponding S2 spectral bands. The broad spectral bands may have received more spectral signal from the background than narrow spectral regions [77]. The bandwidth of the hyperspectral data is much narrower than the corresponding S2 spectral band shown in Table 2. As a result, S2 bands commonly showed better correlations with ULe, and similarly, S2 VIs were more capable of reflecting ULe than narrowband VIs. Future studies may attempt to apply multiple contiguous hyperspectral bands to better extract understory signals. Meanwhile, as parts of the understory can be seen from the gaps of canopy, coarse pixels in S2 imagery are more likely to capture a certain amount of understory signal. The understory in HSI was often affected by the shadow (Figure S1), which lowers the sensitivity to the understory. This might explain the relatively poor regression results between HSI VIs and ULe. Additionally, the canopy gaps and shadow captured in HSI could introduce noise and led to a heterogeneous OLe map (Figure 3a).

Unlike previous studies mainly composed of low herbs, ferns, grasses, and moss [24,45], understories in the study area were composed of shrub mangroves and poaceae, with substantial variation in canopy density and height. Previous studies found that ULe was negatively correlated with OLe in both deciduous and coniferous forests because the overstory gap fraction was directly related to light availability, which was a key factor restricting the growth of the understory [24,45]. However, OLe and ULe showed no significant correlation in our case, which may relate to factors such as species and historical disturbances apart from light availability. The average OLe of AVMA was observed to be higher than KAOB, but the understory LAI under the AVMA is not necessarily lower than under the KAOB (Table S1). For example, only the AVMA species defoliated in late spring due to moth pests. The understory of defoliated areas was able to receive abundant sunlight and grow rapidly. After that, AVMA leaves grew back into a healthy and lush canopy in summer [78,79].

4.2. Important Features and Regression Methods

According to the regression results in Table 1, VIs derived from red-edge, NIR, and SWIR that are highly correlated to OLe (Table 2) are most useful for OLe estimation. Specifically, the contiguous spectral sampling provided by the HSI allows for subtle spectral details to be resolved, especially for red-edge slope/ratio characteristics, which might make the HSI red-edge indices less liable to saturation at high LAI levels [14,20,80]. For instance, IRECI derived from HSI outperformed other variables as it characterizes the red-edge slope by using two red-edge bands while also making use of the minimum and maximum vegetation reflectance found in red and NIR bands, respectively [29]. For S2 VIs, compared with the NIR band (842 nm) that can significantly correlate to both OLe and Ule, SWIR bands (1610 nm and 2190 nm) only show negative correlations with OLe. This allows SWIR bands to reflect the canopy closure and meanwhile eliminate the effect from the background [81], specific to the understory in this study. As a result, SWIR VIs such as NDII and NBR were found to be more useful than other red-edge or NIR VIs. Although LiDAR metrics were less important than VIs for estimating OLe, a few useful LiDAR metrics were identified. Specifically, zsd (the standard deviation of return height) reflected the roughness of the canopy surface and zkurt (kurtosis of return height) indicated the number of returns in the main canopy. We calculated LCI (last echo cover index) by using the proportion of last and single returns in overstory layers, which could more accurately describe both between-crown and within-crown canopy gaps [82].

The penetration attribute of LiDAR data played an important role in estimating ULe. Some LiDAR metrics could represent the understory or low layered returns, for example, den_u reflected the return density of the understory and P_ground and zpcum6 emphasized the percentage of return in the low layer. These LiDAR metrics using returns within a height range corresponding to the understory could help to eliminate the noise from other components and improve the ULe estimation accuracy. S2 VIs such as MTVI2 and MCARI2 have similar mathematical forms that calculate the difference between the red-edge (783 nm) and the red (665 nm) or green (560 nm) band. Previous studies indicated that MTVI2 and MCARI2 were less sensitive to chlorophyll concentration variations and exhibited the best linearity with green LAI. Consequently, VIs that were robustly coupled with LAI were more likely to be selected to present the weak signals from the understory.

It was noted that exponential regression consistently produced the best results between LAI and individual VIs or LiDAR metrics. The results demonstrated a robust non-linear relationship between predictors derived from remote sensing data and field-measured LAI in our case. Similarly, non-parametric regression methods such as RFR and PLSR commonly performed better when handling multiple variables [83,84,85]. RFR is well-suited to handling collinear features by searching for the best split for nodes of the decision tree [86]. In PLSR, it is possible to interpret the underlying physical variables to latent variables [87]. However, the explanatory power of the regression model may be limited by training data for a relatively good fitting result. For instance, it is difficult to furthermore improve the OLe regression model, even by combining the HSI VIs and LiDAR metrics. With limited field observations for this study, we applied feature selection to reduce input variables and five-fold cross-validation to overcome overfitting problem. While in some cases, especially when combining the VIs and LiDAR metrics, the overfitting problem may happen with the small dataset. More field data will be collected in the future to improve the explanatory power and prevent overfitting.

4.3. Technical Challenges and Outlooks

Actual LAI might be overestimated by PAI but removing woody areas from DHP images was difficult. Nevertheless, it is important to recognize that satellite and airborne data always contain a mixed signal of leaf and woody components, and this might balance the relationship between predictors and observations. The direct methods would acquire the most accurate LAI. However, destructive sampling is commonly prohibited in mangrove reserves; therefore, accurate overstory, understory, and total LAI could not be obtained in this study. Fieldwork is particularly challenging in mangrove areas due to the swampy mudflats and tidal fluctuation. Limited by the field environment and manpower, we could not measure understory LAI in a more ideal condition at the time—for instance, raising the lens at least 2 m above the understory or either acquired more observations. It would be ideal to design a solution to better measure the understory LAI and collect more samples for building a robust estimation model. It is worth estimating the total LAI in the multi-layered forest in future studies, so that the contribution of each component to the total LAI could be understood [42]. Rigorous approaches to measure total LAI should be designed to collect the total LAI, which will need to consider the saturation in LAI of the equipment and the measurement distance for both the overstory and understory. The above-mentioned factors, such as species and disturbances, that might affect the relationship between the upper and lower layers should also be considered when building the total LAI estimation model.

5. Conclusions

Our results demonstrated the potential of multi-source remote sensing data and regression methods in OLe and ULe estimation over double-layered mangrove areas. Both spatial and spectral resolution could affect the interpretation of overstory and understory reflectance signals from optical images. High spatial and spectral resolution HSI could better eliminate the influence from the background and performed best in OLe modeling. Meanwhile, it was important to be aware of the background signal captured by multispectral fine resolution S2 images, which was useful for ULe estimation. Novel VIs constructed by red-edge, NIR, and SWIR bands have the potential to improve LAI estimation. The penetrability of LiDAR data alleviated the challenges in ULe estimation. LiDAR data stratified the overstory and understory, which helped to eliminate noise from other forest components. The synergetic use of LiDAR metrics and VIs demonstrates their complementary nature for overcoming limitations such as weak signals and shadow effects for the understory. The results indicated the feasibility of statistical regression methods for both OLe and ULe estimation. Non-linear and non-parametric regression methods frequently produced the best model fitting in this study.