Estimating Forest Canopy Height Using MODIS BRDF Data Emphasizing Typical-Angle Reflectances

Abstract

Forest-canopy height is an important parameter for the estimation of forest biomass and terrestrial carbon flux and climate-change research at regional and global scales. Currently, various methods combining Light Detection and Ranging (LiDAR) data with various auxiliary data, particularly satellite remotely sensed reflectances, have been widely used to produce spatially continuous canopy-height products. However, current methods in use for remote sensing reflectances mainly focus on the nadir view direction, while anisotropic reflectances, which are theoretically more sensitive to the forest canopy height in the multiangle remote sensing field, have rarely been explored. Here, we attempted to examine the potential of using modeled multiangle reflectances at three typical viewing angles (i.e., from the hotspot, darkspot, and nadir directions) to estimate forest-canopy height as auxiliary data sources. First, the sensitivities of the typical angular reflectances as a function of forest canopy height were fully examined using the Extended Fourier Amplitude Sensitivity Test (EFAST) method based on the 4-scale Bidirectional Reflectance Distribution Function (BRDF) model simulations. This indicated that reflectances in the off-nadir viewing directions are generally sensitive to canopy-height variations. Then, the canopy heights were extracted from airborne Laser Vegetation Imaging Sensor (LVIS) data, which were further divided into training and validation data. Moderate Resolution Imaging Spectroradiometer (MODIS) multiangle reflectances at typical viewing angles were calculated from the MODIS BRDF parameter product (MCD43A1, version 6) as partial training-input data, based on a hotspot-adjusted, kernel-driven linear BRDF model. Subsequently, the Random Forest (RF) machine learning model was trained to acquire the relationship between the extracted canopy heights and the corresponding MODIS typical viewing reflectances. The trained model was further applied to estimate the canopy height metrics in the study areas of Howland Forest, Harvard Forest, and Bartlett Forest. Finally, the estimated canopy heights were independently validated by canopy heights extracted from the LVIS data. The results indicate that the canopy heights modeled through this method exhibit generally high accordance with the LVIS-derived canopy heights (R = 0.65−0.67; RMSE = 3.63−5.78). The results suggest that the MODIS multiangle reflectance data at typical observation angles contain important information regarding forest canopy height and can, therefore, be used to estimate forest canopy height for various ecological applications.

1. Introduction

Forest-canopy height is a primary vegetation feature in characterizing the roughness of the land surface and plays an important role in the accurate estimation of aboveground biomass [1,2,3]. This, in turn, is critical for acquiring an improved understanding of the earth’s carbon cycle, ecosystem sustainability, biodiversity, and climate change at regional and global scales [4,5,6].

Forest-canopy height, as one of the main structural parameters of forest canopy, is an important input to many physical models in quantitative remote-sensing research fields. Forest-canopy height is particularly important in the geometric-optical (GO) modeling method, which demonstrates that the anisotropic reflectance patterns are particularly sensitive to vertical forest structures [7,8]. Therefore, anisotropic reflectances of the remotely sensed signals acquired by airborne and/or spaceborne sensors provide an opportunity to estimate forest structural parameters, such as the leaf area index (LAI) [9,10,11,12,13] and clumping index [14,15,16,17,18,19,20,21,22,23]. In particular, anisotropic reflectances have potential in estimating forest-canopy height in many applications. For example, Chopping et al. [24] used a simplified GO canopy reflectance model to interpret the multiangle reflectances from a Moderate Resolution Imaging Spectroradiometer (MODIS) to map the mean canopy height. Wang et al. [25] established a multiple regression model between spectrally invariant variables derived from multiangle satellite observations and the canopy heights extracted using the airborne Laser Vegetation Imaging Sensor (LVIS), which was, in turn, used to examine the capability of using the MODIS bidirectional reflectance distribution function (BRDF) parameters to estimate canopy heights. Heiskanen’s study [26] demonstrated the ability of multiangle reflectances to estimate canopy height using remotely sensed multiangle data from the Multiangle Imaging Spectroradiometer (MISR). Although the acquisition of angular measurements from different multiangle sensors (e.g., MISR and MODIS) are different due to the scanning configuration [27,28], these result suggests that directional information acquired by MODIS may also have the potential to improve canopy height estimates, considering the same platform and the similar spectral bands being explored here.

Recently, airborne and spaceborne Light Detection and Ranging (LiDAR), as an important data source of the active remote sensing technique, has exhibited the potential to accurately acquire forest-canopy height over regional, continental, and global scales. Many efforts have demonstrated the ability to estimate vertical forest structural parameters through the use of LiDAR data [29,30]. Specifically, Kimes et al. [31] used airborne LiDAR data (i.e., LVIS) to estimate forest-canopy height in the Howland Forest; these estimates were reported to correlate highly with the field-measured data. Anderson et al. [32] used airborne LiDAR data to estimate the forest-canopy height at the Bartlett Experimental Forest. Brilli et al. [33] and Rahman et al. [34] derived the urban tree canopy height as the difference between the LiDAR obtained Digital Terrain Model (DEM) and Digital Surface Model (DSM). At a larger scale, the spaceborne Geoscience Laser Altimeter System (GLAS) LiDAR data have been reliably used as the main data source to retrieve forest-canopy height [35,36,37,38].

Most recently, forest-canopy height products have been generated at regional and global scales using LiDAR data in conjunction with various auxiliary data (e.g., quantitative remote sensing products such as LAI and vegetation continuous field (VCF), spectral vegetation indices such as the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI), nadir BRDF-adjusted reflectances (NBAR), forest types, and meteorological data) [37,38,39,40,41,42,43,44]. Previously, we explored the capability of the available MODIS BRDF shape indicators in estimating forest-canopy height that was extracted from the airborne LiDAR data, which showed that these BRDF shape indicators could capture the main variances in canopy height [45]. Theoretically, the multiangle reflectances, particularly in the typical viewing directions of the hotspot, darkspot, and nadir, as shown in Figure 1, contain important information about the forest-canopy height for driving the physical models [7,8]; the models are, therefore, reasonably assumed to have the ability to retrieve vertical forest structures. However, such multiangle reflectances, particularly in the typical hotspot and darkspot directions, have been rarely utilized in the investigation of canopy-height retrievals.

In this study, we designed an experiment to explore the potential of using these typical directional reflectances to estimate forest-canopy height. First, the sensitivities of the reflectances in these several directions were completely analyzed as functions of the canopy height using the Extended Fourier Amplitude Sensitivity Test (EFAST) global sensitivity analysis method [46], based on the simulation results of the so-called 4-scale GO BRDF model. Such an analysis has rarely been performed in previous similar studies. Then, the canopy height was extracted from the full waveform LiDAR data from the airborne LVIS as benchmarked data for independent training and validation in the study regions. The MODIS multiangle reflectances in the hotspot, darkspot, and nadir directions were accurately reconstructed from the MODIS BRDF parameter product (MCD43A1, V006), based on the hotspot-revised kernel-driven linear BRDF model. Finally, we used the Random Forest (RF) machine learning method to train the relationship between the canopy heights and the MODIS multiangle reflectances, which was then applied to map the canopy heights in the study regions of the Howland Forest, Harvard Forest, and Bartlett Forest. The validation results using the independently extracted canopy-height values suggest that the MODIS multiangle reflectances in the hotspot, nadir, and darkspot directions contain important information regarding forest-canopy height; therefore, they have potential in estimating the height to a relatively high accuracy in conjunction with the land-cover products for various applications, particularly as important auxiliary data of the LiDAR data.

The organization of this paper is as follows. In Section 2, we introduced the materials of this study, including the study area and data. In addition, we also briefly described the processing of the data used in this study. In Section 3, we introduced the models and methods, and particularly demonstrated the approaches of the sensitivity analysis and wall–to–wall canopy height estimation. Section 4 focuses on the sensitivity analysis result and predicted canopy height validation results, which are appropriately explained and analyzed, based on the theoretical and data evidence. The discussion and conclusions are finally presented in Section 5 and Section 6, respectively.

2. Study Area and Data

2.1. Study Area

The study areas are located on the east coast of the United States, including three forest regions (Figure 2). The Howland Forest (45.2N, 68.74W) is located in central Maine, U.S.A. The dominant land-cover types in the Howland Forest are boreal-northern hardwood transitional forest and palustrine wetlands. The forest species include red spruce, hemlock, and white pine, consisting of approximately 90% conifer and 10% deciduous tree species [47]. The Bartlett Experimental Forest (44.11N, 71.84W) is located within the White Mountain National Forest, which is a heavily forested and mountainous region in north-central New Hampshire, U.S.A. The dominant species in this forest are birch, maple, and pine fir. This forest area has rich precipitation, averaging approximately 50 inches (i.e., 127 cm) annually, distributed throughout the year [48]. The Harvard Forest (42.54N, 72.17W) lies in central Massachusetts, U.S.A. The dominant species are red maple, red oak, birch, and hemlock.

2.2. LVIS LiDAR Data

LVIS is an airborne laser altimeter system that can record full waveform LiDAR data. LVIS was designed and is operated by the Laser Remote Sensing Laboratory at Goddard Space Flight Center [49]. LVIS is also motivated for a spaceborne LiDAR mission: NASA’s Global Ecosystem Dynamics Investigation (GEDI), which deployed on the International Space Station (ISS) in 2018, will provide billions of measurements of ground elevation and vertical forest structure [50]. This study used the LVIS datasets (

Table 1. Datasets used in this study. LVIS, https://lvis.gsfc.nasa.gov/Data/; MODIS (MCD43A1; MCD12Q1; MOD44B), https://search.earthdata.nasa.gov/; SRTM, http://srtm.csi.cgiar.org/srtmdata/.

Dataset	Year	Resolution	Usage	Reference
LVIS	2003 and 2009	20 m	Derive canopy height	Blair et al. [49]
MCD43A1	2003 and 2009	500 m	Derive multi-angle reflectance values	Schaaf et al. [51]; Wang et al. [52]
MOD44B	2003 and 2009	250 m	Describe tree cover (%)	Hansen et al. [53]
MCD12Q1	2003 and 2009	500 m	Distinguish forest types and identify non-forest pixels	Friedl et al. [54]
SRTM	2000	90 m	Produce slope maps	Farr et al. [55]

) that are located in these three forest regions (Figure 2); the datasets were acquired in mid-August 2003 and late August 2009. The RH100 values derived from the LVIS waveform data are shown in Figure 3, which represents the difference between the elevation of the highest detectable return and the ground elevation. The RH100 value was used as the canopy height in this study. The footprint of the LVIS data is approximately 20 m. To match the spatial resolution of the MODIS product, all RH100 values within each MODIS data grid were averaged as the canopy height value at that pixel with a 500 m spatial resolution.

Table 2. The statistics of tree-canopy heights are derived for a 500 m forest-covered pixel scale in the three forest study regions in this study. These canopy heights were first extracted within the footprint of the LVIS data and then were averaged within each MODIS pixel (500 m).

Forest	File Index	Data Acquisition Year	Forest-Covered Pixels (500 m)	Mean Height (500 m)	Min Height (500 m)	Max Height (500 m)	Standard Deviation (500 m)
Howland Forest	HO03	2003	1467	15.98	7.28	25.28	3.96
	HO09	2009	1981	16.13	8.06	24.95	3.28
Bartlett Forest	BA03	2003	4105	22.53	6.17	39.80	5.99
	BA09	2009	2894	24.66	10.68	37.30	4.38
Harvard Forest	HA03	2003	1086	21.46	10.16	36.44	4.34
	HA09	2009	1119	24.44	12.02	34.84	3.65

shows the LVIS-derived information for each research forest region.

2.3. MODIS BRDF Product

The MODIS BRDF parameter product (MCD43A1) (Table 1) provides the weight coefficients of the RossThick-LiSparse Reciprocal (RTLSR) BRDF model, which qualifies the anisotropic reflectance of Earth’s surface at a 500 m spatial resolution [51]. In this study, we used the “MODIS/Terra+Aqua BRDF/Albedo Model Parameters Daily L3 Global-500 m V006” dataset (MCD43A1 Collection 6) to forwardly calculate the reflectance values in several typical observation angles (i.e., hotspot, darkspot, and nadir), which are, in turn, used to train the canopy-height-estimation model. Notably, to provide accurate forward simulations of these reflectances at these typical observation angles, a hotspot-revised kernel-driven model called RossThickChen-LiSparse Reciprocal (RTCLSR), as introduced in Section 3.1.2, was used. Compared to the Collection 5 data, the daily MCD43A1 Collection product has been reported to represent more accurate anisotropic features of Earth’s surface at any solar and viewing geometry through the improved operational RTLSR algorithm [52,56,57]. In spectroscopy, reflectances in red and near-infrared (NIR) bands are usually strongly related to vegetation features [58,59]. Heiskanen’s study explored the sensitivity of the four MISR spectral bands (blue, green, red, and NIR) to forest canopy height; the results indicated that red and NIR bands were the most useful [26]. More importantly, the forest types in Heiskanen’s study region are similar to ours, including coniferous and broadleaved forest. Therefore, the MODIS data in the red (band 1) and NIR (band 2) bands were selected to estimate canopy height in this study.

Quality assessment (QA) information is provided in the Collection 6 MCD43A1 products. Several QA flags are as follows: QA = 0 accounts for the best quality from the full inversion of the main algorithm; QA = 1 refers to the low quality from the magnitude inversion of the backup algorithm; QA = 255 is for the filled values. The number of high-quality full inversions in the MODIS BRDF parameter products is not necessarily sufficient due to cloud contamination, which has a significant impact on the result of this experiment due to the limited availability of multiangular observations. In this study, we adopted an assumption that canopy height does not change much during a two-month period, so that we can fill the gaps caused by filled-value pixels and the low-quality magnitude inversion values only using the high-quality full inversion data (QA = 0), centered at the dates when the LVIS data were obtained.

2.4. MODIS Vegetation Continuous Field (VCF) Data

The MODIS VCF product (MOD44B) (Table 1) is a global representation of Earth’s surface as gradations of land cover, including the percentages of tree cover, nontree cover, and bare area; the product further shows the subpixel mixtures in each 250 m pixel in bands 1 and 2 [27,53]. This product can characterize the forest coverage within a pixel. Previous studies have shown that vegetation coverage has a significant effect on variations in surface BRDFs [7,60,61]. To consider the impact of the vegetation coverage on the estimation of the forest-canopy height based on the reflectance data, MODIS VCF products that are concurrent with the MODIS BRDF parameter product were inputted as the training data in this study.

2.5. MODIS Land Cover Data

We utilized the MODIS International Geosphere Biosphere Programme (IGBP) land-cover product (MCD12Q1) (Table 1) concurrent with the MODIS BRDF parameter product. The MODIS land cover was used mainly to mask out the possible nonforested pixels. It was additionally used as an index to differentiate the forest types because the differences in the vegetation structures between different forest-canopy types not only lead to differences in the LiDAR echo signal [35] but also have a significant impact on the surface-reflectance patterns obtained by remote-sensing sensors [62,63,64]. Therefore, the land-cover type was considered as an input-data source for examining its potential influence on the estimation of forest-canopy height.

2.6. SRTM Digital Elevation Data

The Shuttle Radar Topography Mission (SRTM) digital elevation dataset (Table 1) provides consistent, high-quality elevation data, which covered latitudes from N60° to S56° [55]. In this study, the void-filled 90 m resolution SRTM elevation data for the global Version 4, available from the CGIAR-CSI SRTM 90 m Database (http://srtm.csi.cgiar.org/srtmdata/), was used to produce the slope maps of the study areas. The produced slope maps were used to improve the understanding of the possible influence of the terrain information on the predicted results in the study area. The potential influence of topography on the estimation of forest-canopy height will be further analyzed in this study.

3. Models and Methods

3.1. Models

3.1.1. The 4-Scale GO BRDF Model

The 4-scale GO BRDF model is a geometric optical (GO) model that can be used to simulate the reflectances of vegetation surfaces observed from the viewpoint of remote sensing sensors [8]. The model is designed to consider the influence of four scales of canopy structure on the BRDF pattern of canopy radiation: (1) stand scale, the distribution characteristics of the canopy community; (2) crown scale, the shape of the crown; (3) branch scale, the spatial distribution of branches within the canopy; (4) leaf scale, the angle distribution of leaves within the canopy. The model assumes that the canopy has a gap distribution at various scales to ensure that the model is similar to the natural growth states of trees. In addition, the model considers the overlap and shadow effects between canopies. A multiple scattering scheme that can realize multi-order scatterings among canopy geometrical structures was introduced in the model to obtain reliable simulation results in the NIR band [65].

The 4-scale model has many input parameters, including site parameters, forest-structure parameters, leaf and background optical-property parameters, and hemispheric spatial-sampling-strategy parameters.

Table 3. Values and ranges of input parameters for the 4-scale model used in this study are mainly referred to some field measurements and related studies, as mentioned in the text.

Parameter	Symbol	Unit	Values and Ranges
Site parameters
Stand density	SD	Trees/ha	500–5000
Vegetation parameters
Leaf area index	LAI	m2	0–8
Clumping index	CI	dimensionless	0.33–1
Canopy height	HC	m	5–60
Crown base height	HB	m	1–10
Crown radius	RC	m	0.5–5
Neyman clustering	NC	dimensionless	1–6
Optical property parameters
Leaf reflectance in red	REDT	dimensionless	0.08
Leaf reflectance in NIR	NIRT	dimensionless	0.6
Leaf transmissivity in red	REDTT	dimensionless	0.05
Leaf transmissivity in NIR	NIRTT	dimensionless	0.35
Background reflectance in red	REDG	dimensionless	0.1
Background reflectance in NIR	NIRG	dimensionless	0.25
Hemispheric spatial sampling strategy parameters
Solar zenith angle	SZA	degree	35
Relative azimuth angle	PHI	degree	0, 180
View zenith angle	VZA	degree	0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80

presents the input parameters and their ranges considered in this study. To keep the input parameters within a reasonable range, the setting of the model parameters in the experiment is mainly referred to the field measurements and related studies [8,14].

3.1.2. The Kernel-Driven RTCLSR BRDF Model

The kernel-driven model has a general form [66,67] for a soil-vegetation system as follows:

$$R(\theta ,\vartheta ,\phi ,\mathsf{\Lambda })=f_{iso}(\mathsf{\Lambda })+f_{vol}(\mathsf{\Lambda })K_{vol}(\theta ,\vartheta ,\phi )+f_{geo}(\mathsf{\Lambda })K_{geo}(\theta ,\vartheta ,\phi ),$$

(1)

where R(θ,ϑ,ϕ,Λ) is the BRDF in waveband Λ; K_geo and K_vol are the trigonometric functions of view zenith (ϑ), illumination zenith (θ), and relative azimuth (ϕ), which are also termed the GO kernel and volumetric scattering kernel, respectively. These two functions provide the shapes for GO-scattering and volume-scattering BRDFs. f_iso(Λ) is a constant for isotropic scattering, f_vol(Λ) and f_geo(Λ) are the weight coefficients for K_vol and K_geo, respectively.

The RTCLSR model adopts the LiSparseReciprocal kernel (K_LSR) as K_geo, which is derived from a sparsely vegetated canopy surface [7]. This parameter describes the shadowing and surface scattering from the discrete tree canopy. The mathematical expression of this kernel is given as follows:

$$K_{LiSparseR}=O\left(\theta \prime ,\vartheta \prime ,\phi \right)-\sec \theta \prime -\sec \vartheta \prime +\frac{1}{2}\left(1+\cos \xi \prime \right)\sec \vartheta \prime \sec \theta \prime ,$$

(2)

where O(θ’, ϑ’, ϕ) is the overlap function of view and illumination shadows on the ground and ξ is the phase angle, which represents the angle between the directions of solar radiation and observation.

The RossThickChen kernel (K_RTC) was adopted as K_vol, which is a hotspot-revised volume-scattering kernel. The mathematical expression of this kernel is shown in Equation (3), which is the revised form of the RossThick kernel by considering a hotspot factor for its original term [68]:

$$K_{RTC}=\frac{\left(\frac{\pi }{2}-\xi \right)\cos \xi +\sin \xi }{\cos {\theta }_v+\cos {\theta }_s}\times \left(1+C_1{e}^{-\frac{\xi }{C_2}}\right)-\frac{\pi }{4},$$

(3)

where $1+C_1{e}^{-\frac{\xi }{C_2}}$ is the hotspot function that was originally proposed by Chen and Cihlar [69] and was modified by Jiao et al., [68] and C₁ and C₂ refer to the height and width of the hotspot effect, respectively. Recently, this model has been successfully utilized to estimate the foliage clumping index at a global scale [17,70], and has been further linked with some physical BRDF models (e.g., PROSAIL) for vegetation characterization [71,72]. More recently, the development of this model focuses on accurately modeling the anisotropic reflectance of snow cover [73,74] and further evaluating the possible influence on albedo retrieval for the improved model relative to its original form [75].

3.2. Methods

3.2.1. Sensitivity Analysis

The EFAST method has been developed based on the Fourier Amplitude Sensitivity Test (FAST) [76] and Sobol’s method [77] to estimate the contribution of the input parameters on the output by the model predictions variances, which are calculated from a decomposition of the Fourier-series representation. The total variances of the output (D) are given as follows:

$$D={\displaystyle \sum _{\mathrm{i}=1}^n{D}_i}+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^n{D}_{ij}}}+\dots +D_{1,2,\dots n},$$

(4)

where D_i is the variance of the output caused by input parameter X_i; D_ij is the variance of the interaction between input parameters X_i and X_j; D_1,2…n is the higher-order variance among multiple input parameters.

Three sensitivity indices can be calculated for each input parameter based on the above variances; i.e., the main (first-order) sensitivity index (S_i), which is given in Equation (5), characterizes the contribution of a single input parameter to the output; the mutual (high-order) sensitivity index (S_i…k) given in Equation (6) characterizes the interactions between a specific parameter and all other parameters; the total sensitivity index ($S_i^T$) is the sum of the above two parts (Equation (7)). Thus, this method considers not only individual input parameters but also the effect of the interaction between input parameters on the output.

$$S_{\mathrm{i}}=D_{\mathrm{i}}/D,$$

(5)

$$S_{\mathrm{i}\cdots \mathrm{k}}=\frac{D_{ij}+D_{i\cdots j+1}\dots +D_{i\cdots k}}{D}\left(j,k\ne i\right),$$

(6)

$$S_i^T=S_i+{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^n{S}_{ij}+\cdots +}{S}_{1,2\cdots n}.$$

(7)

In this study, we employed the EFAST method (Simlab2.2) and 4-scale model to analyze the sensitivity of the typical directional reflectances in the red and NIR bands to the variability of forest-canopy height, as shown in the processing workflow given in Figure 4.

3.2.2. Forest-Canopy Height Estimation

We used the RF regression-tree model to train the relationship between the LVIS-derived canopy heights and MODIS typical directional reflectances in the red and NIR bands; the model was then applied to estimate the forest-canopy heights. The processing workflow of canopy-height estimation from the typical directional reflectances is provided in Figure 4.

RF is a machine learning method that was proposed by Breiman [78] based on the classification-tree algorithm. RF runs by constructing multiple decision trees and randomly selecting training samples; it then outputs a mean prediction regression of the individual trees. Each decision tree is generated by the bootstrap-sampling method. The final predictions are determined by a majority vote of the trees. The RF package in the R language was used to implement the RF algorithm in this study.

The reflectances at typical observation angles (hotspot, darkspot, and nadir) in the red and NIR bands were reconstructed using the MODIS BRDF parameter data based on the hotspot-revised kernel-driven linear BRDF model (RTCLSR). The values of C₁ and C₂ used in the RTCLSR model are the optimized values derived by Jiao et al. [68] from the Polarization and Directionality of the Earth’s Reflectances (POLDER) hotspot observations as prior knowledge. This strategy has been successfully used in the generation process of the global MODIS clumping-index product [17].

To ensure the independence of the LVIS canopy-height data used for training the RF prediction model and for validating the estimated results, the LVIS canopy-height validation data were not used in the model-training process. For example, we selected LVIS-derived canopy heights in the Howland Forest as reference data for validation, but we merely used the canopy-height data in the Bartlett and Harvard Forests to train the canopy-height prediction model. The training data and validation data used in this study are shown in

Table 4. The validation and training data for the three study regions. The table shows that the validation data are completely independent of the training data.

Forest	Mapping Region/Validation Data (LVIS File Index in Table 2)	Model Training Data (LVIS File Index in Table 2)
Howland Forest	HO09	BA03, BA09, HA03, HA09
Harvard Forest	HA09	HO03, HO09, BA03, BA09
Bartlett Forest	BA09	HO03, HO09, HA03, HA09

.

In theory, the accuracy of the LiDAR-derived canopy height can be significantly affected by terrain slope [35,79], because the slope broadens the LiDAR waveform and introduces a bias in canopy-height estimation. Therefore, to reduce the slope-induced uncertainty for the LiDAR-derived canopy height as input in the model-training process, we used mainly the LVIS data with a slope < 10° to train the model for the canopy-height estimation, although we appropriately maintained steep-slope LVIS data to test the results and analyze its potential influence on the estimation accuracy in this study.

4. Results and Analysis

4.1. Sensitivity Analysis of the Typical Directional Reflectances to Canopy Height

Within the possible range of the 4-scale model input parameters shown in Table 3, we first used Simlab2.2 to generate 959 pseudo values as each variable input parameter to the 4-scale model (Table 3) and then to simulate typical directional reflectances in the red and NIR bands using the 4-scale model. With these reflectances generated by the 4-scale model as outputs in conjunction with the corresponding sampling values by the Simlab2.2 as inputs, we finally used the EFAST method to analyze the sensitivity of the typical directional reflectances to forest canopy height.

Figure 5 depicts the EFAST sensitivity profiles of the main, mutual, and total sensitivity indices between the canopy height and the typical directional reflectances in the red and NIR bands. In general, the results present that the typical directional reflectances, especially in the direction of the hotspot and darkspot, show a higher sensitivity to canopy-height variations ($S_i^T$ > 0.2) than that in the nadir direction; on the contrary, the sensitivity of the reflectances in the nadir view directions to the canopy-height variation is the lowest among the three viewing angles in both bands.

The sensitivity of these reflectances to canopy height is somewhat different in the red and NIR bands (Figure 5) because the anisotropic reflectance patterns related to vegetation structural characteristics also depend on spectral bands, known as the BRDF spectral effects [60,61], at different spatial resolutions [80,81]. In the red band, the small leaf transmittance results in low multiple-scattering effects within the canopy and generates high reflectance anisotropy for the vegetation canopy. This makes the anisotropic reflectance in the red band more sensitive to the variability in vegetation canopy structures than that in the NIR band. Conversely, the high leaf transmittance in the NIR band leads to strong multiple-scattering effects and decreases the reflectance anisotropy, which therefore makes the reflectance in the NIR band less sensitive to vegetation canopy structure in theory than that in the red band. In conformity with the theoretical analysis, Figure 5 shows that the main sensitivity index of directional reflectances in the red band to the canopy height is generally higher than that in the NIR band. However, the main sensitivity index does not show much difference in darkspot reflectance (Figure 5), which may be attributed to the fact that the information regarding the canopy height in the shadows is not necessarily different between the red and NIR bands.

The mutual-sensitivity index indicates that the interactions among the forest canopy structure parameters (Table 3) on the typical directional reflectances also exhibit significant contributions to canopy height (Figure 5). Theoretically, vegetation canopy reflectances are affected by various canopy structures, e.g., LAI and canopy-coverage fraction. In the NIR band, multiple scattering effects are usually tightly related to the density of green leaves within the canopy that are quantified by the LAI, and more green leaves will generate more multiple scatterings within the canopy. This indicates that vegetation structures are mutually dependent on each other within the canopy, and the anisotropic reflectance patterns of vegetation canopy are usually the result of the common effects of various structural parameters. Figure 5 shows that the value of the mutual sensitivity index of the hotspot reflectance in the NIR band is the largest ($S_{\mathrm{i}\cdots \mathrm{k}}\left(j,k\ne i\right)$ = 0.71), which in turn causes the hotspot reflectance in the NIR band to have the largest total sensitivity index value ($S_i^T$ = 0.82) to the canopy height. This is apparently because of the most pronounced reflectance values in the hotspot direction in the NIR band. In comparison, the mutual-sensitivity index in the red band is substantially smaller than that in the NIR band, most likely due to the overall low reflectance in the red band.

As a summary, following the EFAST sensitivity analysis results, we can conclude that the directional reflectances in the hotspot and darkspot directions are more sensitive to canopy height than in the nadir direction (i.e., $S_i^T$ for hotspot and darkspot reflectances both exceed 0.4, Figure 5). Therefore, the reflectances in the direction of the hotspot and darkspot are probably more helpful in the canopy-height estimation in terms of this sensitivity investigation.

4.2. Canopy-Height Estimation and Validation

Canopy-height maps with a spatial resolution of 500 m over three large areas at the Howland Forest, Bartlett Forest, and Harvard Forest were produced to verify the ability of typical directional reflectances as inputs to estimate canopy height. The training and validation data used in this Section refer to Table 4. In general, the estimated canopy-height maps for these three large areas (Figure 6, Figure 7 and Figure 8) based on this proposed method show visually consistent spatial patterns with the LVIS-derived canopy heights; this verifies the feasibility of this method in canopy-height estimation.

More specifically, Figure 6a presents the comparison of the LVIS-derived canopy heights with the predicted canopy heights at 500 m resolution in the Howland Forest over an area of 64 × 21 km². The correlation coefficient between the model-estimated canopy heights and LVIS-derived canopy heights is ~0.63, indicating a relatively high correlation. The RMSE is 3.93 m, accounting for ~25% relative errors in this study area (Figure 6b). A bias = 0.06 is almost negligible. In general, this result demonstrates a relatively high mapping accuracy based on this proposed method. Close examination of Figure 6a reveals some differences between the two maps of the LVIS-derived and model-estimated canopy heights. Obviously, the model-estimated canopy-height values are somewhat larger than the LVIS extracted canopy-height values, especially at both ends of the areas that are already marked by two black rectangles. Further examination of these two regions, A and B as individual examples (Figure 6a), with Google Maps, shows that these two areas are dominated by some mixed MODIS pixels, including lakes, rivers, and residential areas that are misclassified as forests in the MODIS land-cover product. The corresponding pixels in the scatter plot shown in Figure 6b are marked by gray points. Once such outliers are removed in regions A and B (gray spots shown in Figure 6b), the results have a somewhat improved prediction accuracy with a correlation coefficient of ~0.65 and RMSE of ~3.63 m. This analysis indicates that the uncertainty caused by the land-cover type map, particularly for some mixed pixels, is probably a main error source in the canopy-height estimation using the proposed method in this study.

We then analyzed another large area in the Harvard Forest. This region covers approximately 39 × 9 km² (Figure 7a). The mapping results in Figure 7a show that the model-estimated canopy heights have very consistent spatial patterns with the LVIS-derived canopy heights. Similar to Figure 6, the statistical evaluation results are further provided for this area, which have good agreement between the LVIS-derived and model-estimated canopy heights with a correlation coefficient of ~0.65. The RMSE is 3.89 m, accounting for ~17% relative error (Figure 7b). The bias is 1.58 m. In general, this result further demonstrates the feasibility and reliability of this proposed method. Notably, some prominent outliers still exist in this mapping region, as shown in Figure 7b, and thus it needs a close investigation, as shown in Figure 7a. We found that mixed pixels are again the main error source in this mapping region. As an example, we examined some pixels in region A (Figure 7a), including not only trees but also farmland, buildings, and bare land, which are generally labeled as forests in the MODIS land-cover product. Once we removed these outliers from region A (gray spots in Figure 7b), the results were somewhat improved, with correlation coefficients reaching 0.67, the RMSE reducing to 3.81 m, and the bias decreasing to 1.51 m. Obviously, a better result is expected if such low-quality land-cover pixels are masked out through the use of a high-resolution land-cover map.

Finally, the proposed method was applied to the mountainous Bartlett Forest with an area of approximately 27 × 53 km² (Figure 8a). In this study area, the total correlation coefficient and the RMSE between LVIS-derived and model-estimated canopy heights are ~0.48 and ~6.35 m (Figure 9a), respectively, without taking the effects of terrain into account. Compared with the previous two results of the Howland and Harvard forests, the result accuracy in this region is very challenging. To examine this result, we investigated the slope map in this study area (Figure 8b); it shows that the terrain in this study region is relatively complex with average and maximum slopes of 11.18° and 39.39°, respectively. The terrain effect in this study area is completely different from the first two study areas. Therefore, the influence of topographic factors on the canopy-height estimation in this region must be emphasized and further analyzed. Our specific statistics show that over 50% of the forest-covered pixels in this region are located in areas with slopes greater than 10°.

Figure 10 shows that the LVIS-extracted canopy height has a generally smaller value range (i.e., ~15–31 m) than the model-estimated canopy height (i.e., ~7–35 m), indicating that the LVIS signals are relatively insensitive to canopy-height variations in this mountainous region due to the topographic effects. However, the model-estimated canopy height tends to be overestimated in a larger range of the canopy-height values (e.g., >27 m), but tends to be underestimated in a smaller range of the canopy-height values (e.g., <15 m), compared to the LVIS-extracted canopy height. Once these pixels with steep slopes (slope >10°) were excluded from this map, the validation results greatly improved, with a correlation coefficient and RMSE of ~0.67 and ~5.78 m, respectively (Figure 9b). Obviously, the relatively large RMSE (i.e., ~5.78 m) in this mountainous area indicates that even a relatively small slope (e.g., 5–10°) probably has a relatively large influence on canopy-height estimation. Somewhat surprisingly, the terrain effect does not seem to largely influence the bias, most likely because the bias is offset by the overestimation and underestimation of the canopy heights in this case study.

5. Discussion

For the first time, the typical directional reflectances in the hotspot, nadir, and darkspot directions have been comprehensively tested to estimate forest-canopy height in conjunction with LiDAR signals in several selected study regions. The results indicated that multiangle reflectances in the hotspot, darkspot, and nadir directions contain important information about forest-canopy heights and therefore can be used as an effective data source to estimate forest-canopy height.

In the BRDF sampling capability of the MODIS sensors, typical reflectances are usually difficult to directly observe in most cases, particularly for the hotspot effect. To address this issue, we reconstructed the typical reflectances in the three viewing angles of the hotspot, nadir, and darkspot directions using the operational high-quality MODIS BRDF parameter product (i.e., MCD43A1) based on the RTCLSR model. The RTCLSR model is a hotspot-adjusted version of the kernel-driven RTLSR BRDF model that can accurately capture hotspot signatures by using two adjustable hotspot parameters (C₁ and C₂). The optimized values of the C₁ and C₂ parameters in the RTCLSR model that were derived from the 6 × 7 km POLDER hotspot observations in the previous study [68] were directly used as prior knowledge to reconstruct the 500 m MODIS multiangle reflectances. Obviously, the accuracy of these reconstructed typical reflectances is closely related to the robustness of the kernel-driven model. As is well known, the reconstructed reflectances from MODIS multiangle observations have been widely used in many fields, e.g., derivation of the MODIS global-land-cover products [82], improvement of the accuracy of various land-surface-classification schemes [83,84], and estimation of the structural parameters of vegetation canopies [17,64]; however, the potential uncertainty caused by such typical reflectances through use of the RTCLSR model must be considered and discussed, particularly for the pixel-size difference in hotspot signatures between the MODIS and POLDER sensors at two very different spatial resolutions [72,81].

In addition, in this study, the relatively coarse spatial resolution of MODIS data is definitely an inevitable error source, introducing various uncertainties in the tree-height retrievals, especially by the mixed pixels from the MODIS land-cover product. This study shows that some pixels that are identified as of forest type in the MODIS land-cover product are actually mixed pixels that cover the forest, roads, and water. These mixed pixels tend to have a low prediction accuracy of tree-canopy heights, resulting in a significant bias in the process of tree-height estimations. Indeed, the mixed pixels were previously reported as a strong limitation factor, especially in utilizing coarse-resolution satellite images to monitor vegetation dynamics [85]. Similarly, this study further confirms that the effect of mixed pixels deserves much more attention in the study of canopy-height estimation, as shown in Figure 6 and Figure 7.

Finally, the terrain slope is another important factor determined here to have a large effect on tree-height retrievals. As shown in Figure 9, this study generally performs well in low-to-moderate slope areas (i.e., slope <10°), but the results get worse in mountainous study regions with relatively steep slopes. Indeed, previous studies noted that the topographic factors had a large influence not only on the BRDF patterns [86,87,88,89] but also on the accuracy of the LiDAR-derived canopy height [35,79], which is, in turn, usually used as the true value to train the canopy-height prediction model and/or to validate the predicted results in such studies. Therefore, we must be aware of the potential uncertainty caused by a steep terrain slope not only for using the multiangle reflectances but also for deriving the LVIS-derived canopy height in such studies.

Despite the various uncertainties mentioned above, we must realize that the multiangle reflectances, even in just the three typical directions of the hotspot, nadir, and darkspot being explored here, contain important information regarding canopy heights, and therefore are generally helpful in the retrievals of canopy heights in conjunction with the corresponding land-cover data.

6. Conclusions

Forest-canopy height is an important parameter for the estimations of forest biomass and terrestrial carbon flux that are in turn essential to climate-change research at regional and global scales. In this study, the potential of using three typical directional reflectances (i.e., reflectances obtained from the hotspot, darkspot, and nadir directions) to estimate forest canopy height was comprehensively explored based on the LVIS-extracted canopy-height data in the Howland Forest, Harvard Forest, and mountainous Bartlett Forest. This study highlights several conclusive results as follows:

(1) The reflectances in the off-nadir viewing directions are generally sensitive to canopy height variations in terms of the sensitivity analysis using the 4-scale model simulations.

(2) Generally, a relatively high-accuracy canopy height could be retrieved by using three typical directional reflectances (i.e., hotspot, nadir, and darkspot) in conjunction with corresponding land-cover data in the three study areas. Specifically, the correlation coefficients (R) and RMSE between the modeled canopy height and the LVIS-derived canopy height arrive at 0.63 and 3.93 m for the Howland Forest; 0.65 and 3.89 m for the Harvard Forest; 0.48 and 6.60 m for the mountainous Bartlett Forest. Such results are generally comparable to the previous studies using various auxiliary input data [42,43,44]. The results can be significantly improved if mixed pixels or pixels with a steep slope (e.g., slope >10°) are examined and excluded.

(3) The MODIS multiangle reflectances at the three typical viewing angles show the potential to estimate forest-canopy height in conjunction with corresponding land-cover data. However, the problem of mixed pixels and large terrain slope pixels are still a challenge and need more research, particularly in conjunction with using high-resolution land cover classification maps in the near future.