Effectiveness of the Reconstructed MODIS Typical-Angle Reflectances on Forest Biomass Estimation

Abstract

Multi-angle optical reflectance measurements such as those from the NASA moderate resolution imaging spectroradiometer (MODIS) are sensitive to forest 3D structures, potentially serving as a useful proxy to estimate forest structural variables such as aboveground biomass (AGB)—a potential theoretically recognized but rarely explored. In this paper, we examined the effectiveness of the reconstructed MODIS typical-angle reflectances—reflectances observed from the hotspot, darkspot, and nadir directions—for estimating forest AGB from both theoretical and practical perspectives. To gain theoretical insights, we first tested the sensitivities of typical-angle reflectances to forest AGB through simulations using the 4-scale bidirectional reflectance distribution function (BRDF) model. We then built statistical models to fit the relationship between MODIS multi-angle observations and field-measured deciduous-broadleaf/mixed-temperate forest AGB at five sites in the eastern USA, assisted by a semivariogram analysis to determine the effect of pixel heterogeneity on the MODIS–AGB relationship. We also determined the effects of terrain and season on the predictive relationships. Our results indicated that multi-angle reflectances with fewer visible shadows yielded better AGB estimates (hotspot: R² = 0.63, RMSE = 54.28 Mg/ha; nadir: R² = 0.55, RMSE = 59.95 Mg/ha; darkspot: R² = 0.46, RMSE = 65.66 Mg/ha) after filtering out the effects of complex terrain and pixel heterogeneity; the MODIS typical-angle reflectances in the NIR band were the most sensitive to forest AGB. We also found strong sensitivities of estimated accuracies to MODIS image acquisition dates or season. Overall, our results suggest that the current practice of leveraging only single-angle MODIS data can be a suboptimal strategy for AGB estimation. We advocate the use of MODIS multi-angle reflectances for optical remote sensing of forest AGB or potentially other ecological applications requiring forest structure information.

1. Introduction

Forest aboveground biomass (AGB) is a key vegetation feature for quantifying terrestrial carbon dynamics, and has been widely used for explaining climate change, and describing habitats and biodiversity [1,2,3]. Continuous forest AGB maps are effective means to monitor vegetation dynamics across a large scale, such as a mainland or global scale [4]. The ecosystem community has suggested that continuous AGB maps with sufficient accuracy are desirable, which if deployed can avoid representation issues that often arise when using discontinuous field-based observations [5]. With the development of remote sensing technology, a variety of obtained remote sensing data provided opportunities for mapping large-scale forest AGB, such as optical remote sensing data, radio detection and ranging (radar) data, and light detection and ranging (lidar) data.

Radar and lidar, as active remote sensing, have been widely used for large-scale forest AGB estimation because of their ability to penetrate the canopy. In terms of radar, the most frequently used strategies for estimating forest AGB are based on the backscatter values and interferometry technique [6]. For example, radar backscatters in the P and L bands are strongly correlated with AGB and have been successfully used to estimate forest AGB [7,8,9]; adding C-band HV or HH polarization data to the regression equations significantly improved the performance of forest AGB estimation [10,11]. Spaceborne lidar is also a suitable means for large-scale forest AGB estimation—using spaceborne lidar observation combined with continuous remote sensing data, such as radar and multispectral data, for mapping forest AGB [12,13,14,15]. Recently, some researchers have started to test the performance of combining the newly developed space-borne lidar system (i.e., global ecosystem dynamics investigation) with InSAR observations to map forest AGB [16].

Although lidar and radar instruments are promising approaches for forest AGB estimation—data alone or combined with other ancillary data for AGB estimation, these active instruments are limited in temporal and geographic coverage—a problem that cannot be solved in the short term. Multi-angle remote sensing with long temporal and global scales coverage provides an alternative for forest AGB estimation because the multi-angular remotely sensed signals are theoretically sensitive to vegetation structures [17,18,19,20]. Multi-angle remote sensing refers to the observation of surface reflectance from multiple view angles together with an account of the solar position as well; this kind of technology provides an opportunity to measure the anisotropic reflectance of the land surface. The observed anisotropic reflectances, especially from the typical viewing directions of hotspot (observation with no visible shadows), darkspot (observation with a maximum of visible shadows), and nadir, contain important information about the forest structure for driving the physical models. These physical-based canopy geometric-optical models were developed to estimate the bidirectional reflectance distribution function (BRDF) analytically or empirically [21,22,23], as a function of the canopy shape and size, canopy density, tree height, and spatial distribution of trees. Therefore, it can reasonably be assumed that these multi-angle anisotropic reflectances have the ability to retrieve forest structures, such as biomass.

Until now, multi-angle reflectances obtained from airborne and spaceborne sensors have been widely used for forestry applications, such as leaf area index (LAI), clumping index (CI), and canopy height mapping [24,25,26]. Due to canopy height, LAI, and CI, these canopy structure parameters are related to AGB estimation in some respect [27,28,29]. Accordingly, some researchers have attempted to use BRDF anisotropic information to estimate forest AGB by using multi-angle surface reflectance data from a multi-angle imaging spectro radiometer (MISR), successfully mapping forest AGB for the United States [30]. The moderate resolution imaging spectroradiometer (MODIS), one of the important components of the global optical multi-angle observation system, provides massive multi-angle observation information on global-scale long-term series; it should be a good choice for independently, or as auxiliary data combined with radar or lidar observation, mapping forest AGB. Some researchers attempt to use vegetation indices derived from MODIS, such as a normalized difference vegetation index (NDVI), an enhanced vegetation index (EVI), a soil-adjusted vegetation index (SAVI), or an optimized soil-adjusted vegetation index (OSAVI) to estimate AGB [31,32]; among these indices, the NDVI usually has the best performance [33]. However, some researchers considered that indices such as the NDVI produced less spectral information; they therefore introduced more bands of MODIS spectral information to improve the accuracy of the AGB inversion [34]. As research progressed, landcover types and meteorological data were used as auxiliary variables to combine the MODIS vegetation indices to model AGB [35,36]. It should be noted that all the above methods based on MODIS observations for estimating the AGB focus on the nadir spectral information. The multiangle anisotropic information (e.g., hotspot and darkspot directions anisotropic reflectance), which contains important data about the forest structure information, has rarely been explored for forest AGB estimation.

In this study, we designed an experiment that exploits the MODIS multi-angle reflectances observed from hotspot, nadir, and darkspot directions for forest AGB estimation. Specifically, we first used the 4-scale model simulations and performed a global sensitivity analysis based on the extended Fourier amplitude sensitivity test (EFAST) to understand the sensitivity of the multi-angle reflectances to the variation of forest AGB in theory. We subsequently took the field-based forest AGB data at the hectare scale as the true biomass value for MODIS pixels. We used a semivariance function to analyze the heterogeneity of MODIS pixels, which in turn helps to determine whether the field-based biomass data can represent the true biomass value of MODIS pixels. A hotspot-revised kernel-driven BRDF model was used to reconstruct the MODIS multi-angle reflectances; the model was operated based on a version 6.1 MCD43A1 MODIS parameter product as input. According to the reconstructed multi-angle reflectances and BRDF shape indicators constructed based on these reflectances, we applied a linear regression strategy to examine the relationship between BRDF information and forest AGB. Furthermore, we analyzed the influence of factors such as terrain slope and tree-cover percentage on this relationship. Finally, we discussed the seasonal effects on the BRDF, while these effects on the relationship between the multi-angle reflectances and forest AGB were also discussed.

2. Materials

2.1. Study Area

Our study areas are located on the eastern coast of the United States and include a total of 55 field-measurement sites distributed in five forests (Figure 1), namely Howland Forest (lat: 42.20°, lon: −68.73°), Penobscot Forest (lat: 44.87°, lon: −68.65°), Bartlett Forest (lat: 44.04°, lon: −71.16°), Hubbard Brook Forest (lat: 43.56°, lon: −71.45°), and Harvard Forest (lat: 42.53°, lon: −72.17°) [37]. Howland Forest is located in central Maine and was designated as a research forest in 1986. The forest stands located in the boreal-northern hardwood transitional zone; the dominant species were red spruce (Picea rubens), eastern hemlock (Tsuga canadensis), balsam fir (Abies balsamea), white pine (Pinus strobus), and northern white cedar (Thuja occidentalis). The terrain at this site was flat to gently rolling with a maximum elevation of less than 68 m. Penobscot Forest is also located in Maine, across the Penobscot River from Orono and the University of Maine. State Highway 178 parallels the eastern bank of the river and provides access to the forest. The dominant forest type is mixed northern conifers, and the average growing season is 156 days. Bartlett Forest was established in 1931 as a research site within the U.S. Forest Service, which was located at White Mountain National Forest in New Hampshire— it extends from the village of Bartlett in the Saco River valley at 210 m to about 915 m at its upper reaches. The primary forest cover of this site is the sugar maple-beech-yellow birch association. Hubbard Brook Forest lies in the southern part of the White Mountain National Forest in central New Hampshire; forest types are mainly deciduous northern hardwoods, such as sugar maple (Acer saccharum), beech (Fagus grandifolia), and yellow birch (Betula alleghaniensis). Harvard Forest is an ecological research site owned and managed by Harvard University and located in Petersham, Massachusetts. The forest stands are located in the hardwoods-white pine-hemlock transition zone; forest types mainly include red oak (Quercus rubra), white birch (B. papyrifera), white pine (Pinus strobus), and beech (Fagus grandifolia). Field measurement plots in this forest are located at the Prospect Hill tract; the terrain is relatively complex.

2.2. MODIS MCD43A1 BRDF Product

The MODIS MCD43A1 Version 6.1 daily BRDF/Albedo model parameters dataset was produced using 16 days of Terra and Aqua MODIS data at 500 m resolution, which provided the three model weighting parameters (isotropic, volumetric, and geometric) for the kernel-driven BRDF model [38]. Here, we used these three model weighting parameters to reconstruct the reflectances in the hotspot, nadir, and darkspot directions in the red and NIR bands based on a hotspot-revised kernel-driven BRDF model. Compared to the previous version of the MODIS products, the version 6.1 product was improved in several aspects through various calibration changes. In particular, a new version of the response-versus-scan angle (RVS) method was employed to minimize optical crosstalk in the Terra MODIS infrared (IR) bands. In addition, the Terra MODIS forward look-up table (LUT) was updated for the 2012–2017 period.

The MODIS MCD43A1 Version 6.1 product includes a quality assessment parameter to help users to perform quality filtering. We only used the data flagged as high-quality for our research, that is, where the “QA” flag was equal to 0. Due to cloud contamination, high-quality data were unavailable at a given time. Here, we followed a strategy that fills in these missing data with high-quality full inversion data (QA = 0) from a month close to the field-measurement time. The time range of the used MODIS data was from June to August and was located in the growing season. Due to the fact that the structure of the forest does not change much in a month, our strategy will in theory have limited impact on study results.

2.3. MODIS MOD44B Vegetation Continuous Field (VCF) Data

The MOD44B vegetation continuous fields (VCF) annual product layers include percentage tree cover, percentage non-tree cover, and percentage non-vegetation, providing a continuous, quantitative representation of land surface cover at 250 m resolution in bands 1 and 2 [39]. Here, we used the layers of percentage tree cover from the year 2009, which were time-matched with field measurements, to characterize the forest coverage situation of the corresponding MODIS pixel. The purpose of characterizing the forest cover information at the pixel level is to analyze the influence of forest coverage on variations in surface BRDFs, which, in turn, helps us to understand its effect on forest AGB estimation based on BRDF data.

2.4. MODIS Land Cover Data

The MODIS land cover type (MCD12Q1) Version 6 data product provides global land cover types at annual intervals with 500 m pixel resolution [40]. We used the International Geosphere-Biosphere Programme (IGBP) classification system of the MODIS MCD12Q1 product to filter out the non-tree-covered study sites and help to analyze the influence of forest types on the BRDF surface-reflectance patterns, which further analyzes if these influences will pass to forest biomass estimation. The acquired time of the used MCD12Q1 data (i.e., year) was consistent with field-based data.

2.5. SRTM Data

The shuttle radar topography mission (SRTM) collected elevation data on a near-global scale to generate the most complete high-resolution digital topographic database on Earth [41]. Previous studies have shown that BRDF effects would be weak when topographic effects were significant [42,43,44]. Therefore, using BRDF information to estimate forest structure in complex terrain areas, such as study areas with a large slope, will have certain limitations. Here, we used SRTM 1 Arc-Second Global elevation data to calculate the slope information of the MODIS pixel, which in turn was used to analyze whether terrain effects influence using BRDF observations to estimate forest biomass.

2.6. Landsat Surface Reflectance Data

The Landsat clear sky surface reflectance data describe the percentage information about the incoming solar radiation reflected from the Earth’s surface [45,46]. These data contain information describing the differences in the material on the Earth’s surface; thus, they have been widely used as a reference to conduct pixel heterogeneity analysis [47,48]. In this study, we used the field-measured forest biomass as the true value for the MODIS pixel; the problem here is therefore whether the data observed on the ground can represent the real situation of the MODIS pixel. For homogeneous pixels, the ground observations can represent the real situation of the pixel scale to a large extent. To ensure that the ground observation data we used can reasonably indicate the real situation of the MODIS pixel scale, we used the 30 m clear sky Landsat surface reflectance data, which had the same acquired time as field measurements, as input to the semivariogram function to perform the pixel heterogeneity analysis.

2.7. Field-Based Forest Biomass Data

The field-based forest AGB data used in this study were provided by the North American Carbon Program (NACP). The NACP is a multi-disciplinary research program designed to obtain a scientific understanding of North America’s carbon sources and sinks [37]. The forest AGB datasets in this study were collected in 2009 from Bartlett Experimental Forest, Harvard Forest, Howland Research Forest, Hubbard Brook Experimental Forest, and Penobscot Experimental Forest, including 55 field measurement sites (partial collection data shown in Appendix A 

Table A1. The field-measured sites selected in the Howland Forest (lat: 42.20°, lon: −68.73°), Penobscot Forest (lat: 44.87°, lon: −68.65°), Bartlett Forest (lat: 44.04°, lon: −71.16°), Hubbard Brook Forest (lat: 43.56°, lon: −71.45°), and Harvard Forest (lat: 42.53°, lon: −72.17°). Here, we provided some property information for these sites, including site name, location information, measured biomass, measured time, and the composition of measured biomass (i.e., dead tree counts (D) and live tree counts (L)).

Site	Plot_ID	Latitude (°)	Longitude (°)	Biomass (Mg ha−1)	Time (yyyy-mm-dd)	D/L
Bartlett	13L	44.053185	−71.310543	200.59	2009.7.11	75/648
Bartlett	21L	44.054021	−71.300303	229.44	2009.7.13	17/411
Bartlett	30N	44.054162	−71.289812	255.46	2009.7.12	82/800
Harvard	PH1	42.534074	−72.182013	305.44	2009.6.25	179/936
Harvard	PH10	42.536547	−72.175841	139.03	2009.7.24	55/785
Harvard	PH2	42.538096	−72.177597	208.88	2009.7.27	63/753
Harvard	PH3	42.536557	−72.172724	256.43	2009.7.23	24/475
Harvard	PH4	42.536516	−72.179817	271.89	2009.7.14	27/551
Harvard	PH5	42.540983	−72.170486	219.34	2009.7.28	66/702
Harvard	PH6	42.540467	−72.183034	127.38	2009.7.16	61/817
Harvard	PH7	42.539223	−72.187066	282.08	2009.7.17	70/834
Harvard	PH8	42.551416	−72.176897	145.36	2009.7.26	22/359
Harvard	SC1	42.480697	−72.174601	206.87	2009.7.25	83/763
Harvard	SF2	42.508234	−72.250973	309.28	2009.7.27	77/612
Harvard	TS2	42.512857	−72.205741	236.60	2009.7.25	20/520
Howland	H2	45.22755	−68.725911	26.83	2009.8.20	12/242
Howland	H3	45.225188	−68.724381	34.39	2009.8.24	26/629
Howland	H5	45.222658	−68.716496	91.87	2009.8.25	42/571
Howland	H6	45.214881	−68.735791	57.80	2009.8.24	18/432
Howland	H8	45.214646	−68.709366	18.65	2009.8.26	0/148
Howland	H9	45.210844	−68.737554	105.80	2009.8.19	14/541
Howland	H12	45.203327	−68.741371	167.57	2009.8.19	76/1212
Howland	H17	45.152076	−68.735178	131.73	2009.8.27	35/687
Howland	H18	45.147732	−68.718229	122.98	2009.8.27	30/677
Hubbard Brook	1	43.936143	−71.741518	267.26	2009.7.22	52/614
Hubbard Brook	200	43.940344	−71.778636	261.25	2009.7.20	57/833
Hubbard Brook	339	43.945148	−71.709622	257.82	2009.7.27	97/850
Hubbard Brook	354	43.941246	−71.703841	246.54	2009.7.18	63/628
Hubbard Brook	349–350	43.947527	−71.704189	213.38	2009.7.24	60/618
Penobscot	P1	44.871236	−68.626076	233.43	2009.8.25	97/687
Penobscot	P4	44.858001	−68.620421	44.76	2009.8.24	12/886
Penobscot	P5	44.851611	−68.618074	124.65	2009.8.18	17/484
Penobscot	P6	44.850592	−68.613788	51.60	2009.8.19	29/604
Penobscot	P7	44.848417	−68.615501	122.27	2009.8.19	13/491
Penobscot	P10	44.84406	−68.619475	120.84	2009.8.20	19/672
Penobscot	P11	44.844779	−68.614519	93.37	2009.8.20	10/549
Penobscot	P13	44.835663	−68.599269	199.65	2009.8.26	94/994

). At each site, live trees and standing dead trees with a diameter at breast height (DBH) greater than or equal to 10 cm within a 50 m by 200 m plot were counted to calculate the AGB. The AGB for each stem was calculated in megagrams (Mg) using the general equations from Young et al. (1980) [49].

3. Models and Methods

We used two strategies to assess the ability of the MODIS reflectances observed from the hotspot, nadir, and darkspot directions on forest AGB estimation. Firstly, we used the simulations from the 4-scale model and based on the EFAST global sensitivity analysis method to explore the relationship between multi-angle reflectances and forest AGB. Then, according to the field-measured forest AGB data, we used the regression method to analyze the correlation between MODIS multi-angle reflectances and the forest AGB; the MODIS multi-angle reflectances were reconstructed based on a hotspot-revised kernel-driven BRDF model. The performance of the BRDF shape indicators constructed from these multi-angle reflectances in the forest AGB estimation was also analyzed. Furthermore, we analyzed the influence of complex terrain and pixel heterogeneity on their relationship. Finally, we analyzed the seasonal effects of MODIS data on the estimation of the forest AGB. The details of the models and methods used in this study are as follows.

3.1. Models

3.1.1. The RTCLSR Kernel-Based BRDF Model

We employed a hotspot-revised kernel-based BRDF model, named RTCLSR [50], to reconstruct MODIS multi-angle reflectances. The kernel-based BRDF model usually consists of three components: isotropic scattering, volumetric scattering, and geometric-optical surface scattering [51,52]. The equation is shown as follows:

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

(1)

where R(θ, ϑ, φ, Λ) is the BRDF in waveband Λ; $f_{iso}$, $f_{vol}$, and $f_{geo}$ are the weight parameters that determine the shape of the BRDF, that is, the proportion of the above three parts of scattering; K_geo and K_vol are geometric-optical and volume-scattering functions, respectively; all are dependent on the geometric configurations, including the view zenith (ϑ), illumination zenith (θ), and relative azimuth (φ). These geometric configurations can be used to define the concerned sun-sensor geometries, such as the hotspot (ϑ = 35°, θ = 35°, φ = 0°), nadir (ϑ = 35°, θ = 0°, φ = 0°) and darkspot (ϑ = 35°, θ = 35°, φ = 180°).

For the RTCLSR model, K_geo adopts a kernel function derived from a sparsely vegetated canopy surface, named the LiSparseReciprocal kernel (K_LSR) [23], which is shown 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 represents the overlap area between the view and solar shadows, and $\xi$ is the phase angle between the incident and observe direction. In terms of K_vol, a hotspot-revised volume-scattering kernel named the RossThickChen kernel (K_RTC) was adopted; its mathematical function was as follows:

$$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)

This function was developed from the RossThick kernel, which introduced a hotspot factor, $1+C_1{e}^{-\frac{\xi }{C_2}}$, compared with its original term. C₁ and C₂ are two adjustable parameters related to the height and width of the hotspot effect. Here, the suitable input values of C₁ and C₂ for red and NIR bands reference previous studies [25,53].

Subsequently, the reconstructed multi-angle reflectances were used to calculate the MODIS BRDF shape indicators, which were further used for forest AGB estimation. Here, five kinds of BRDF shape indicators were selected for our study (

Table 1. Formulating MODIS BRDF shape indicators based on the multi-angle reflectances observed from the hotspot, nadir, and darkspot directions of the red and NIR bands.

Names	Abbr.	Formulas	Application Areas
Normalized difference between hotspot and darkspot index	NDHD	$NDHD=\frac{Hotspo{t}_{\mathrm{NIR}/\mathrm{red}}-Darkspo{t}_{\mathrm{NIR}/\mathrm{red}}}{Hotspo{t}_{\mathrm{NIR}/\mathrm{red}}+Darkspo{t}_{\mathrm{NIR}/\mathrm{red}}}$	Clumping index
Hotspot darkspot index	HDS	$HDS=\frac{Hotspo{t}_{\mathrm{red}}-Darkspo{t}_{\mathrm{red}}}{Darkspo{t}_{\mathrm{red}}}$	Clumping index
Anisotropic factor	ANIX	$ANIX=\frac{Hotspo{t}_{\mathrm{NIR}/\mathrm{red}}}{Darkspo{t}_{\mathrm{NIR}/\mathrm{red}}}$	Landcover types
Normalized difference anisotropic index	NDAX	$NDAX=\frac{ANI{X}_{\mathrm{red}}-ANI{X}_{\mathrm{NIR}}}{ANI{X}_{\mathrm{red}}+ANI{X}_{\mathrm{NIR}}}$	Landcover types
Hotspot darkspot NDVI	NDVI-HD	$NDVI\_HD=\frac{Hotspo{t}_{\mathrm{NIR}}-Darkspo{t}_{\mathrm{red}}}{Hotspo{t}_{\mathrm{NIR}}+Darkspo{t}_{\mathrm{red}}}$	Leaf area index
Hotspot incorporated NDVI	NDVI-HS	$NDVI\_HS=NDVI\_HD\times \left(1-Hotspo{t}_{\mathrm{red}}\right)$	Leaf area index

), which have been successfully used for retrieving CI [54], LAI [55], and landcover types [56].

3.1.2. 4-Scale Model

The 4-scale model is a geometric optical (GO) BRDF model, that is widely used to simulate the multi-angle reflectances of vegetation surfaces observed from remote sensing sensors [21]. We used the 4-scale model to perform the simulation of measuring the BRDF behavior of forest canopies; the output will be used as input for the EFAST method to assess the effectiveness of the BRDF on forest biomass estimation from a theoretical perspective. There are four types of input data that need to be inserted to realize the simulation of the four-scale model at the stand scale, crown scale, branch scale, and leaf scale, such as site parameters, forest canopy parameters, optical parameters for canopy and background, and observation geometry parameters (

Table 2. Input parameters for 4-scale model, including parameters value and ranges.

Input Parameter	Symbol	Unit	Values and Ranges
Site parameters
Stand density	SD	Trees/ha	500–5000
Canopy parameters
Leaf area index	LAI	m2	0–8
Clumping index	CI	dimensionless	0.33–1
Canopy height	HC	m	5–60
Crown based height	HB	m	1–10
Crown radius	RC	m	0.5–5
Newman 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 transitivity in red	REDTT	dimensionless	0.05
Leaf transitivity in NIR	NIRTT	dimensionless	0.35
Background reflectance in red band	REDG	dimensionless	0.1
Background reflectance in NIR band	NIRG	dimensionless	0.25
Observation geometry parameter
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

). The setting of these input parameters in our study references our study areas’ field measurements and previous studies [21,54], which can make these input parameters reasonable and suitable for our study areas to a certain extent.

As we can see from Table 2, there is no biomass of this option as input parameters for the 4-scale model. Existing studies have shown that canopy height and crown diameter are highly correlated with the variation in the forest biomass [57,58]; we can note that these two parameters are in the input table of the 4-scale model. Therefore, we adopted an indirect strategy to analyze the sensitivity of the BRDF information on the variation of forest biomass—we analyzed the relationship of the BRDF information with the canopy height and crown diameter. Specifically, we used the EFAST method to quantitatively evaluate the sensitivity of the BRDF information to the variation in canopy height and crown diameter.

3.2. Methods

3.2.1. EFAST Global Sensitivity Analysis Method

We employed the EFAST global sensitivity analysis method to evaluate the importance and apportioned the uncertainty of the input factors to the output factors of the model. The EFAST method was developed from the classical Fourier amplitude sensitivity test (FAST) method [59]. Compared with FAST, the new feature of EFAST is that it can compute not only the first-order sensitivity index (first-order) but also the mutual sensitivity index (high-order) and the total sensitivity index.

The sensitivity index of the EFAST method is calculated based on the predictions variances of the model, which are shown as follows:

$$D={\displaystyle \sum _{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, D_ij and D_1,2,…n are the variances in the output caused by the different dimensions of the input parameter; D is the total variance of the output. For the main sensitivity index (S_i), also known as the first-order sensitivity index, characterizes a single input parameter to the sensitivity of the output; its mathematical expression is as follows:

$$S_i=D_i/D$$

(5)

The mutual sensitivity index, which describes the influence of the interactions between one specific parameter and all other parameters in the output, is calculated as follows:

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

(6)

Conducting a sum process for the main sensitivity index and mutual sensitivity index can produce the total sensitivity index:

$$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)

3.2.2. Assessment of Pixel Homogeneity

The field-based measurements of whether we can represent the real situation at pixel scale is an unavoidable problem in remote sensing study fields, as well as in our study, as we mentioned above. In relatively homogeneous pixels, field-based measurements can usually represent the real state of the landcover of the pixel [47]. Therefore, evaluating the homogeneity of pixels is key to carrying out remote sensing authenticity evaluation experiments. Here, we employed a strategy based on a semivariogram function [48] to assess the homogeneity of MODIS pixels corresponding to field observations. Specifically, using Landsat-TM finer resolution images as an input for semivariogram functions, if the output yields a sill value of less than 5.0 × 10⁻⁴, we considered the landcover of MODIS pixels to be homogeneous; this conclusion was reached by the previous study [60]. We give an example that shows our employed strategy to assess the homogeneity of MODIS pixels (Figure 2): using a cloud-free Landsat surface reflectance product centered on a MODIS pixel as input for the semivariogram function. We obtained three different scales (i.e., 1 km, 1.5 km, and 2 km) of Landsat-TM subsets’ variogram estimators; the sill values of these three spatial elements are 0.95 × 10⁻⁴, 2.25 × 10⁻⁴, and 2.95 × 10⁻⁴, respectively, which are all less than 5.0 × 10⁻⁴. Therefore, we can reasonably consider that the observed targets within this MODIS pixel are homogeneously distributed.

3.3. Accuracy Validation

The ability of MODIS BRDF observations to retrieve forest AGB was evaluated by the reference datasets. Specifically, we employed ordinary least squares regression (OLS) and regression through the origin (RTO) analysis [61] to validate the accuracy; we additionally utilized three statistical indices; they were the mean bias, the coefficient of determination, and the root mean square error:

$$R^2=\frac{{\displaystyle \sum _{i=1}^n{\left(V_i-\overline{V}\right)}^2-{\displaystyle \sum _{i=1}^n{\left(V_i-\overline{\overline{V_i}}\right)}^2}}}{{\displaystyle \sum _1^n{\left(V_i-\overline{V}\right)}^2}}$$

(8)

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(V_i-\overline{\overline{V_i}}\right)}^2}}{n}}$$

(9)

$$nRMSE=\frac{RMSE}{\overline{V}}$$

(10)

where V_i is the estimated value, $\overline{\overline{V_i}}$ is the reference value, $\overline{V}$ is the average result of reference values, and n is the used data amount.

4. Results and Analysis

4.1. Sensitivity Analysis of the BRDF Information to Forest Biomass

Here, we used the EFAST method to assess the sensitivity of the BRDF information on the variation of the canopy height and crown diameter based on the simulations of the 4-scale model. Since the canopy height and crown diameter are closely related to AGB, this section can indirectly help us to understand the feasibility of using BRDF data to estimate the AGB.

4.1.1. Sensitivity Analysis of Typical-Angle Reflectances to Forest Biomass from Canopy Height and Crown Diameter

Our results show that the reflectances obtained in the typically observed direction in the red and NIR bands show a high sensitivity to the variation in crown diameter and canopy height (Figure 3). For example, the total sensitivity index ($S_i^T$) of the hotspot from the NIR band and darkspot from the red band both exceed 0.6 for the canopy height; in terms of the crown diameter, every observed typical-angle reflectances’ total sensitivity index ($S_i^T$) exceeds 0.7. Furthermore, it can be noted that this sensitivity degree performs differently for red and NIR bands. The above results may relate to the formation mechanism of the anisotropic pattern of BRDF—BRDF anisotropic patterns that are usually not only determined by vegetation structure but also by spectral characteristics [56,62]. For a red band, leaf transmittance activity is usually weak, which in turn leads to weak multiple scattering effects within the canopy; in other words, the anisotropic pattern of the BRDF performs strongly in the red band. For the NIR band, the results are exactly the opposite—the spectral properties of the NIR band with a strong transmittance characteristic usually caused strong multiple scattering effects, and thus exhibited relatively weak anisotropic patterns. However, the main sensitivity index does not show much difference for the darkspot in the red and NIR bands; this result occurs because the dependence of shading on the spectrum is not very strong. In addition, it can be noted that the reflectances observed from the nadir direction performed differently to the variation of crown diameter and canopy height; although these reflectances showed sensitivity to both parameters, they are more sensitive to a variation in the crown diameter.

Furthermore, our results show that the total sensitivity is not only affected by the main sensitivity but also by the mutual sensitivity (Figure 3). For example, the mutual sensitivity index of a hotspot in the NIR band contributes around 80% of the total sensitivity to the variation in crown diameter and canopy height. Multiple scattering effects in the NIR band should be the main reason leading to the mutual sensitivity—the multiple scattering effects within the canopy describe the interactions in the forest canopy structure from the spectral aspects; these interactions finally contributed to the mutual sensitivity index. Therefore, we can reasonably consider that the BRDF anisotropic patterns for the NIR band are mostly affected by common effects of various structural parameters.

Overall, the reflectances obtained from the typically observed directions are sensitive to the variation in forest structure variables closely related to the forest biomass from the perspective of sensitivity studies. Therefore, the multi-angle reflectances obtained from typically observed directions should have the potential ability to estimate the forest AGB.

4.1.2. Sensitivity Analysis of BRDF Shape Indicators to Forest Biomass from Canopy Height and Crown Diameter

Our results show that all of these selected BRDF shape indicators display high sensitivity to variations in the canopy height and crown diameter (Figure 4), although these indexes were not previously used for forest AGB estimation. Compared to the typical-angle reflectances, the sensitivity of these indexes to variations in the forest canopy structure mostly comes from mutual sensitivity. As shown in Table 1, all of these BRDF shape indicators are derived from multi-angle reflectances in red and NIR bands. These multi-angle reflectances result from the combination of geometrical optics scattering and volume scattering. For our studies, volume scattering is in general multiple scattering caused by randomly distributed components that form gaps within the canopy; the geometrical optics scattering effects can be understood as a single scattering caused by the gap between different canopies. In other words, these BRDF shape indicators are constructed based on the multi-angle reflectances resulting from the coupling of various forest canopy structures, such as canopy height, LAI distribution, and crown diameter. In general, the geometrical optical effect in the red band is more pronounced than the volume scattering effects [23,42], while in the NIR band, the multiple scattering effects are enhanced [63]. Previous studies documented that the vegetation usually exhibits different anisotropic scattering patterns due to differences in the canopy structure and main scattering types [64,65]. Therefore, the forest canopy structure presented by the BRDF anisotropic pattern will behave differently due to the variations in scattering type; this theory is consistent with our results—the main and mutual sensitivities of the different BRDF shape indicators to the canopy structure parameters behave differently. Overall, the above results indicate that these BRDF shape indicators constructed from typical-angle reflectances have the potential to estimate the forest AGB.

4.2. Performance of the MODIS BRDF Information on Forest Biomass Estimation

4.2.1. Performance of the MODIS Typical-Angle Reflectances on Forest Biomass Estimation

Linear regression results for estimating the forest AGB using MODIS typical-angle reflectances based on 55 sets of field-based biomass data are summarized in

Table 3. Regression results for estimating forest AGB using MODIS typical-angle reflectances based on field-based forest biomass data.

Multi-Angle Reflectances	Using All Field-Based Biomass Data			Using Selected Field-Based Biomass Data
	R2	RMSE (Mg/ha)	nRMSE	R2	RMSE (Mg/ha)	nRMSE
Hotspot-NIR	0.20	77.18	0.41	0.63	54.28	0.31
Nadir-NIR	0.21	76.54	0.41	0.55	59.95	0.34
Darkspot-NIR	0.12	80.92	0.43	0.46	65.66	0.38
Hotspot-red	0.07	83.27	0.45	0.25	77.65	0.45
Nadir-red	0.02	85.59	0.46	0.21	79.44	0.46
Darkspot-red	0.001	86.33	0.46	0.06	86.82	0.50

. It is obvious that there is an un-ideal correlation between MODIS typical-angle reflectances and the forest AGB. For example, the hotspot reflectance of the NIR band can only explain 20% of the forest AGB variation, when RMSE = 77 Mg/ha and nRMSE = 0.41, which is the best performance of typical-angle reflectances. The above results indicated that the reflectances of MODIS under typical observation directions were worthless for forest AGB estimation, as they were completely inconsistent with our sensitivity analysis results.

As shown in Figure 5a, some of the sites measured in the field are in some heterogeneous MODIS pixels due to water bodies, bare land, and building land. Furthermore, we counted the terrain slope information of the MODIS pixels where the ground measurement sites were located; a total of 10 sites were located in rugged areas, that is, the terrain slope was greater than 8 degrees (Figure 5b). The above evidence indicates two noteworthy issues in our linear regression analysis: first is the representation issue—whether the field-based forest AGB can represent the three-dimensional surface that characterizes the forest AGB across the MODIS pixel; second is the rugged terrain’s influence on BRDF anisotropic reflectances. These issues’ influence on remote sensing inversion and authenticity verification experiments has been widely reported [66,67,68]. We have a total of 19 sites with the above issues. Once these sites were excluded from the regression analysis, the results were significantly improved (Table 3). For example, the explanatory ability of the hotspot of the NIR band for the forest AGB increased from 20% to 63%, with an improved RMSE from 77.18 Mg/ha to 54.28 Mg/ha; furthermore, the normalized RMSE improved from 0.41 to 0.31, indicating a better fitting model. The above results are also consistent with our sensitivity analysis based on 4-scale model simulation data; that is, the hotspot reflectance of the NIR band has the best performance among typical-angle reflectances for explaining forest AGB variation. Although reflectances in the NIR band observed from the nadir and darkspot direction are not as good as the hotspot for estimating AGB, they are still effective.

The X-Y scatter plot of the linear regression shows that the MODIS typical-angle reflectances in the NIR band yielded a positive correlation with the forest AGB variation, while the parameters in the red band showed a negative correlation (Figure 6). Furthermore, the MODIS multi-angle reflectances in the NIR band can explain more forest AGB variation than those in the red band. For the NIR band, it can be noted that multi-angle reflectances with fewer shadowing parts (i.e., shadowed crown and shadowed ground) can explain more forest AGB variation. Reflectances observed from the hotspot are better than those observed from the nadir, and the reflectances observed from the nadir are better than those observed from the darkspot (i.e., hotspot: R² = 0.63, RMSE = 54.28 Mg/ha, nRMSE = 0.31; nadir: R² = 0.55, RMSE = 59.95 Mg/ha, nRMSE = 0.34; darkspot: R² = 0.46, RMSE = 65.66 Mg/ha, nRMSE = 0.38). The effect of shadow on the spectral signal has drawn much attention in previous studies because of the links between shadow and forest structure parameters [22,69,70]. Therefore, the forest AGB estimation using BRDF spectral information seemed to be influenced by the shadowing effect.

Although the X-Y scatter plot of the linear regression shows a strong correlation between the hotspot reflectance of the NIR band and forest AGB, most of the scatters are distributed on both sides of the fitted line. This result may be due to several reasons. First, although we conducted a pixel heterogeneity analysis, all MODIS pixels’ tree cover percentage we used was around 60% (Figure 7a); this means that the reflectance from non-forest parts within the pixel, such as grass, will be contained in the total reflectance of the MODIS pixel. Therefore, some biases will be involved and unavoidable when using these reflectances to explain the forest structure. Second, the field-based forest AGB data, which we used as the actual biomass value for the MODIS pixel, are collected from both the live tree and the dead tree (Figure 7b); most of the biomass is contributed to by the live tree. The spectral characteristics of live and dead trees are quite different during the growing season, as red and NIR bands are sensitive to variations in green leaves within the canopy. In other words, the BRDF anisotropy pattern in red and NIR bands performs differently for the canopy with or without green leaves, which in turn will have a certain influence on forest AGB estimation based on BRDF information.

4.2.2. Performance of the MODIS BRDF Shape Indicators on Forest Biomass Estimation

We also evaluated the performance of some classical BRDF shape indicators constructed from typical-angle reflectances on the forest AGB estimation. We used the selected field-based biomass data (Appendix A Table A1) as the reference forest biomass and these study sites without the influence of terrain slope and pixel heterogeneity. From assessment results, it is evident that the BRDF shape indicators NDVI-HD and NDVI-HS can capture the variation in the forest AGB (Figure 8). These two parameters originally designed to estimate the LAI achieved a good performance as the density of leaf elements of the canopy is related to the AGB estimation to a certain extent. Among the two, the NDHD-HS performs better than the NDVI-HD, which can explain 52% of the AGB variation, with RMSE = 60.55 Mg/ha and nRMSE = 0.35. The NDHD-HS is an improved version of the parameter NDVI-HD, where a new term (1-hotspot_red) was introduced. This new term represents the forest coverage to some extent because the hotspot reflectance in the red band is sensitive to the hidden area of the ground which is greatly affected by the tall and dense canopy. Therefore, introducing this new term can reduce the interference of non-forest elements within the MODIS pixel to the forest AGB estimation, further improving the accuracy of the forest AGB estimation.

4.3. Seasonal Effects of Using MODIS BRDF Information to Estimate Forest Biomass

The BRDF anisotropy effect usually varies over the seasons due to the solar angle variation and land cover changes [65,71]; using BRDF information to derive the forest AGB should therefore be influenced by this variation (i.e., seasonal effects). Here, we selected the MODIS BRDF data observed from the summer (August) and winter (December) to summarize BRDFs variation between seasons. It can be noted that the BRDFs show differences in both the NIR band and the red band, of which the hotspot reflectance in the NIR band, which is most sensitive to the forest AGB variation, performs very differently in summer and winter (Figure 9). Furthermore, we quantitatively assessed this seasonal effect based on BRDF data observed in winter and biomass data measured in summer (Figure 10). The result indicated that the BRDF information showed limitations to reflecting the forest AGB variation when the experimental results from ground and satellite observations were not time-matched (a hotspot of the NIR band: RMSE = 54.77 Mg/ha vs. 76.84 Mg/ha, nRMSE = 0.31 vs. 0.44; NDVI-HS: RMSE = 60.55 Mg/ha vs. 72.92 Mg/ha, nRMSE = 0.35 vs. 0.42), although the correlation still existed (a hotspot of the NIR band: R² = 0.61 vs. 0.23; NDVI-HS: R² = 0.52 vs. 0.3). Therefore, this seasonal effect needs to be taken into account when estimating the forest AGB based on BRDF information.

5. Discussion

The performance of MODIS multi-angle observations in estimating the forest AGB has rarely been studied, in particular the typical directional reflectances in the hotspot, nadir, and darkspot directions. However, BRDF anisotropic information to estimate forest structure parameters has been widely employed, such as the CI, the canopy height, and the LAI. For the above reason, we explored the feasibility of BRDF anisotropy information in retrieving the forest AGB. The overall results show that BRDF anisotropic information can estimate the forest AGB, especially the hotspot reflectance in the NIR band; however, some potential questions of uncertainty require further discussion, as shown below.

5.1. Uncertainty Factors of This Study

Due to the limited ability of MODIS sensors sampling, some typical directional reflectances (e.g., hotspot and darkspot) are somewhat challenging to capture most of the time. However, our results show that the accuracy of these reflectances, especially the hotspot reflectances, is critical for estimating the forest AGB accurately. To address this issue, we used a hotspot-revised kernel-driven BRDF model, RTCLSR [50], to reconstruct these typical directional reflectances. Although the multi-angle reflectances reconstructed based on RTCLSR have been successfully used for the CI [25] and canopy height mapping [26], there are still some potential uncertainties. For example, the RTCLSR model used two adjustable hotspot parameters, C₁ and C₂, to accurately capture hotspot signatures; the optimized values of the C₁ and C₂ for the red and NIR bands that we used were provided by the previous study, which was derived from the polarization and directionality of the Earth’s reflectances (POLDER) observations with wide-angle sampling capabilities, including the hotspot direction. Clearly, the resolution is very different between the 6 × 7 km POLDER and the 500 × 500 m MODIS, with the result that there is a scale effect when applying prior knowledge derived from the POLDER directly to MODIS observations [72], which in turn brings uncertainty to subsequent applications, such as the biomass estimation.

Another issue that needs attention is whether the field-based biomass data can represent the real situation of the MODIS pixel. Our ground surveys were carried out on a 1-hectare scale, which was smaller than the 500 m spatial resolution MODIS pixels. In this study, we consider that the field-based data can represent the real situation of the MODIS pixel when the MODIS pixel is homogeneous, which is consistent with the previous conclusions [47,73,74]. The homogeneity of the MODIS pixel is assessed by 30 m Landsat surface reflectance data using the semivariance function. Although the strategy we adopted has been widely used in pixel homogeneity evaluation, this strategy is largely based on the spectral information provided by remote sensing images; some non-forest elements within a pixel usually have similar spectral properties as a canopy, such as grassland, which may lead to uncertainties for our biomass estimation results. In further research, we will consider replacing the field-based measurements with almost continuous distributed lidar data as the “true value” for MODIS pixels, which may avoid such representative problems to a certain extent.

5.2. Prospects of Using BRDFs for Large-Scale AGB Mapping

The results of this study represent an important advance by providing an option for mapping the forest AGB over large areas and at high temporal frequencies from space that can complement the data analysis of orbital active remote sensing data. However, several issues still need attention before large-scale mapping applications. First, our method does not consider the effects of differences between tree species. Different tree species usually have different canopy structures, which means that the biomass is closely related to tree species. For example, allometric functions for biomass estimation are frequently given separately for different tree species [75]. In this study, we hoped to obtain a universal relationship between the forest biomass and BRDF, which means that we hoped the differences in canopy structure could be explained directly from BRDF information. The theoretical basis for this claim is that BRDF observations contain not only spectral information, but also information about the canopy structure (i.e., BRDF anisotropic information) [23]. Although our results prove that using BRDF information can directly carry out the biomass inversion, in view of the limited observation capability of the canopy structure by passive remote sensing, we think the effect of tree species differences in the inversion still exists; how this affects the performance of estimating biomass using BRDFs needs to be further explored.

The second point is that our study results only focus on temperate forests, with the consequence that our method’s performance in retrieving other forest types’ AGB, such as tropical forests, needs to be explored. MODIS-derived anisotropy information has been proven to have a strong correlation with tropical forest canopy roughness, even in high biomass regions [18]. The canopy roughness usually has a close relationship with variability in canopy heights; this variability leads to a variability in geometric scattering components especially of NIR reflectance [63]—a key basis for our method to estimate the AGB. The above evidence gives us confidence that our method has the potential to estimate the AGB for tropical forests. One further point is that our method still requires more evaluation before application because our method was developed in a specific land cover environment. In other words, the generality of our method requires further verification, such as conducting a comparison with widely used biomass estimation approaches. In the next stage, we will investigate the performance of our method in mapping the AGB of tropical forests and conduct cross-validation experiments; we will then combine them with the results of this study to carry out a large-scale forest AGB mapping, such as on a global scale.

6. Conclusions

MODIS observed anisotropy information has been acknowledged as an effective means for the inversion of vegetation structure parameters in passive optical remote sensing areas; however, the estimation of forest AGB is rarely reported. In this work, we explored the value of MODIS BRDF anisotropy information for forest AGB estimation, especially from the typical directional reflectances (i.e., reflectances observed from the hotspot, darkspot, and nadir directions). Overall, the BRDF anisotropy information obtained from MODIS observations has the potential to retrieve the forest biomass, and the obtained results are as follows:

(i)

The typical directional reflectances in the red and NIR bands and the constructed BRDF shape indicators show sensitivity to capture the variation in biomass-related canopy structure parameters (i.e., canopy height and crown diameter) in terms of the sensitivity analysis using the 4-scale model simulations.

(ii)

The MODIS typical directional reflectances in the NIR band show a good linear relationship with the field-based forest AGB after filtering the influence of terrain slope and pixel heterogeneity; in particular, the hotspot reflectance from the NIR band can explain up to 62% of the biomass variations. It is also worth noting that the BRDF shape indicators (i.e., NDVI-HD and NDVI-HS) that are constructed from MODIS multi-angle observations and were originally designed for the inversion of LAI have a good linear relationship with field-based forest AGB. In particular, NDVI-HD with a ground vegetation cover index term yielded a better performance than NDVI-HS and can explain up to 52% of the biomass variation.

(iii)

Seasonal effects on biomass estimation using BRDF data are noteworthy. Seasonal changes lead to changes in the spectral characteristics of the land cover and the observed geometry at the corresponding locations; the above changes will eventually influence the observed multi-angle anisotropic information, further affecting the estimation accuracy of forest AGB. Therefore, it is necessary to ensure the consistency of the ground observation and satellite observation time as fully as possible when constructing a biomass inversion model based on BRDF data.