Estimating Forest Canopy Height Using MODIS BRDF Data Emphasizing Typical-Angle Reﬂectances

: Forest-canopy height is an important parameter for the estimation of forest biomass and terrestrial carbon ﬂux and climate-change 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.

. Scheme of the hotspot, nadir, and darkspot directions along the principal plane (PP) and the corresponding directional reflectances. RHS refers to hotspot reflectance, and RDS refers to darkspot reflectance. Negative values along the abscissa axis refer to the backward scattering direction, and positive values refer to the forward scattering direction.
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 6, respectively.

Study Area
The study areas are located on the east coast of the United States, including three forest regions ( Figure 2). The Howland Forest (  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 Sections 5 and 6, respectively.

Study Area
The study areas are located on the east coast of the United States, including three forest regions ( Figure 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) 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

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) 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 shows the LVIS-derived information for each research forest region. 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 shows the LVIS-derived information for each research forest region.  Figure 3. A schematic diagram of canopy height (RH100) derived from LVIS data. RH100 is the difference between the elevation of the highest-detectable return and the ground elevation.

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, Figure 3. A schematic diagram of canopy height (RH100) derived from LVIS data. RH100 is the difference between the elevation of the highest-detectable return and the ground elevation. 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

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.

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.

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.

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 Remote Sens. 2019, 11, 2239 7 of 21 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.

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-samplingstrategy parameters. Table 3 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]. The kernel-driven model has a general form [66,67] for a soil-vegetation system as follows: 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: 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]: where 1 + C 1 e − ξ 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 1 and C 2 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].

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: 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 T i ) 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.
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.
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 forestcanopy height, as shown in the processing workflow given in Figure 4.

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.  Figure 4. Workflow chart of the study methods, indicating that a necessary sensitivity analysis was performed prior to the wall-to-wall canopy-height estimation. . Workflow chart of the study methods, indicating that a necessary sensitivity analysis was performed prior to the wall-to-wall canopy-height estimation.

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 1 and C 2 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.
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.

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 T i > 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 multiplescattering 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 ( 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., T i S for hotspot and darkspot reflectances both exceed 0.4, Figure 5). 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 i···k ( j, k i) = 0.71), which in turn causes the hotspot reflectance in the NIR band to have the largest total sensitivity index value (S T i = 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 T i 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.

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 (Figures 6-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 2 . 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. 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.

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 (Figures 6-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 6(a) 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 2 . 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 6(b)). 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 6(a) reveals some differences between the two maps of the LVIS-derived and modelestimated 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 6 (a)), 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 6(b) are marked by gray points. Once such outliers are removed in regions A and B (gray spots shown in Figure 6(b)), 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 2 (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. We then analyzed another large area in the Harvard Forest. This region covers approximately 39 × 9 km 2 (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 2 (Figure 8(a)). 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 9(a)), 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 8(b)); 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 Finally, the proposed method was applied to the mountainous Bartlett Forest with an area of approximately 27 × 53 km 2 (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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 21 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 9(b)). 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.   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 9(b)). 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.

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 (C1 and C2). The optimized values of the C1 and C2 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 Figures 6 and 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

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 1 and C 2 ). The optimized values of the C 1 and C 2 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 Figures 6 and 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.

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.