Evaluation of the MODIS LAI/FPAR Algorithm Based on 3D-RTM Simulations: A Case Study of Grassland

Abstract

Uncertainty assessment of the moderate resolution imaging spectroradiometer (MODIS) leaf area index (LAI) and the fraction of photosynthetically active radiation absorbed by vegetation (FPAR) retrieval algorithm can provide a scientific basis for the usage and improvement of this widely-used product. Previous evaluations generally depended on the intercomparison with other datasets as well as direct validation using ground measurements, which mix the uncertainties from the model, inputs, and assessment method. In this study, we adopted the evaluation method based on three-dimensional radiative transfer model (3D RTM) simulations, which helps to separate model uncertainty and other factors. We used the well-validated 3D RTM LESS (large-scale remote sensing data and image simulation framework) for a grassland scene simulation and calculated bidirectional reflectance factors (BRFs) as inputs for the LAI/FPAR retrieval. The dependency between LAI/FPAR truth and model estimation serves as the algorithm uncertainty indicator. This paper analyzed the LAI/FPAR uncertainty caused by inherent model uncertainty, input uncertainty (BRF and biome classification), clumping effect, and scale dependency. We found that the uncertainties of different algorithm paths vary greatly (−6.61% and +84.85% bias for main and backup algorithm, respectively) and the “hotspot” geometry results in greatest retrieval uncertainty. For the input uncertainty, the BRF of the near-infrared (NIR) band has greater impacts than that of the red band, and the biome misclassification also leads to nonnegligible LAI/FPAR bias. Moreover, the clumping effect leads to a significant LAI underestimation (−0.846 and −0.525 LAI difference for two clumping types), but the scale dependency (pixel size ranges from 100 m to 1000 m) has little impact on LAI/FPAR uncertainty. Overall, this study provides a new perspective on the evaluation of LAI/FPAR retrieval algorithms.

1. Introduction

Leaf area index (LAI), defined as half of the total green leaf area of per unit horizontal ground area, is a basic parameter for measuring the vegetation canopies [1,2]. This variable plays a key roles in hydrology, biogeochemistry, and ecosystem models that connect vegetation to the climate observing system through the carbon, water cycles, and radiation [3]. Fraction of photosynthetically active radiation (0.4–0.7 μm) absorbed by vegetation (FPAR) measures the proportion of the solar radiation entering at the top of the plant canopy that contributes to the photosynthetic activity [3,4,5,6]. LAI/FPAR retrieved from remote sensing observations in the reflective solar domain, are used as input parameters for models monitoring the Earth’s surface continuously and are key parameters recognized by the global climate observing system (GCOS) to describe climatic characteristics [3,7]. LAI/FPAR products, derived from atmospherically corrected surface reflectances, have entered a new era since the moderate resolution imaging spectroradiometer (MODIS) became operational in 1999 [8,9,10]. The MODIS LAI/FPAR products (MOD15), based on the radiative transfer (RT) model [11], have been widely used to corroborate global climate change [12], to serve as key inputs for terrestrial carbon cycle models [13], and to support the research of both phenomena and possible reasons of large scale vegetation dynamics [14,15,16]. Moreover, the generation of MODIS LAI/FPAR products does not depend on other LAI/FPAR datasets and they are commonly used as input and reference data for the generation and intercomparison of other products [17,18].

Intensive evaluation and validation efforts have been carried out to examine the uncertainty of MODIS LAI/FPAR products and the corresponding retrieval algorithm. These works mainly included: (1) theoretical derivation based on model mechanisms and error propagation [19]; (2) intercomparison with other LAI/FPAR products or related variables (e.g., GLASS, CYCLOPES, VIIRS) [20,21,22,23,24,25]; (3) direct validation using ground LAI/FPAR measurements [25,26]. The theoretical derivation has an explicit mathematical basis and does not require other datasets; however, this approach is highly correlated with the algorithm itself and is easily affected by model limitations and uncertainties [27,28]. Intercomparison with other LAI/FPAR products is an approach that can effectively analyze the spatio-temporal consistency of long-term LAI/FPAR, but the results cannot meet the requirement of product usage and algorithm refinement. Ground-based validations are essential as the basis of all validations, but the accuracy of this validation method includes the uncertainty of the ground measurements, the spatial heterogeneity-caused uncertainty [29] in the upscaling process from the point measurement to the pixel scale, and the product uncertainty. Above all, the previous studies mainly focused on the evaluation of product uncertainty, which introduces the coupled uncertainties from the model, inputs, and assessment method. Therefore, it would hinder the process of evaluating the uncertainty of the algorithm itself and understanding the deficiencies of the algorithm, thus hampering future improvements to the algorithm.

In the above context, real scene computer simulations provide a new approach for remote sensing evaluation and validation [30]. As computing power improves, several 3D RT models have been developed for scene simulation [31,32], such as DART (discrete anisotropic radiative transfer) [33], RAPID (radiosity applicable to porous individual objects) [34], and LESS (large-scale remote sensing data and image simulation framework) [35]. These models have become an important tool in the field of quantitative remote sensing, particularly for studying the radiometric properties of the Earth’s surface [31,36]. 3D RT models can analyze the detailed interactions between solar radiation and vegetation canopies [37], analyze the radiative properties of specific biome types [38], and help the science team define the characteristics of optical sensors through model simulation [39]. Data from simulations based on 3D RT models are widely used for model validation and evaluation. The DART model has been used in studies on the surface energy budget [40], the impact of canopy structure on satellite image texture [41], the 3D distribution of photosynthesis and primary production rates of vegetation canopies [42], and forest biophysical parameter retrieval [43,44]. The LESS model can synergistically use spectral and angular information to simulate the radiation properties of complex realistic landscapes, which can be used for simulating datasets of 3D landscapes [45]. The outputs of LESS can serve as benchmarks for retrieval algorithm evaluation since it has a solid theoretical foundation and its accuracy has already been well-assessed by comparison with other models of radiation transfer model intercomparison (RAMI) [35] and field measurements [46].

This study aimed to provide a new perspective on the evaluation of MODIS LAI/FPAR retrieval algorithms, which differs from previous research by evaluating the algorithm itself rather than the product. In this paper, a computer simulation of a real grassland scene is performed using the ray-tracing LESS model to analyze the uncertainty of the MODIS LAI/FPAR retrieval algorithm. The advantage of simulation-based model evaluation is that the uncertainty caused by a single variable can be analyzed to avoid the effects caused by the mixing of multiple factors. The uncertainty of the MODIS LAI/FPAR algorithm was evaluated by separating the model and input uncertainties. In addition, further analysis was conducted to understand the impact of scale dependency and clumping. The results can serve as guidance for improving this algorithm continuously.

The structure of this paper is organized as follows. Section 2 briefly describes the MODIS LAI/FPAR retrieval algorithm, how we use LESS to analyze the retrieval algorithm for uncertainty, and the methodologies for uncertainty evaluation. Section 3 details the results of LAI/FPAR uncertainty caused by inherent model, reflectance, and biome type uncertainties as well as the clumping effect, and scale dependency. The discussions, including the analysis of the experiment results, are detailed in Section 4. Finally, Section 5 provides some concluding remarks.

2. Materials and Methods

2.1. MODIS LAI/FPAR Retrieval Algorithm

The MODIS LAI/FPAR retrieval algorithm consists of a main algorithm based on the radiative transfer equation (RTE) and a backup algorithm using the relationship between vegetation index and LAI/FPAR. The retrieval algorithm exploits the spectral information content of MODIS surface reflectances at up to 7 spectral bands (band 1: 620–670 nm; band 2: 841–876 nm; band 3: 459–479 nm; band 4: 545–565 nm; band 5: 1230–1250 nm; band 6: 1628–1652 nm; band 7: 2105–2155 nm) [4,8]. Inputs of this algorithm include BRFs at red and near-infrared (NIR) bands (band 1 and 2), their uncertainties, sun–sensor geometry (SZA: solar zenith angle, SAA: solar azimuth angle, VZA: view zenith angle, VAA: view azimuth angle), and a biome classification map. Note that in the current algorithm version, different biome types use different RT models. Herbaceous biomes (B1: grasses and cereal crops; B2: shrubs; B3: broadleaf crops;) were modelled using 1D RT due to the good continuity of the grass distribution and in consideration of the computational efficiency. Savannas (B4) were modelled by a stationary Poisson germ-grain stochastic process (so called stochastic radiative transfer (SRT) model) [47,48]. Forest biomes (B5: evergreen broadleaf forests; B6: deciduous broadleaf forests; B7: evergreen needleleaf forests; and B8: deciduous needleleaf forests) were based on a 3D RTM (3D structures were represented by columns uniformly (deterministically) spaced on the ground). With these RTMs, the science team constructed an LAI/FPAR main algorithm based on angular information, biome type, and spectral information in which the mean and standard deviation values of the LAI and FPAR selected in the spectral retrieval space are reported for retrieval value and its uncertainty. The main look up table (LUT)-based algorithm was designed as follows. Firstly, the main algorithm evaluates a weight coefficient as a function of sun–sensor geometry, wavelength, and LAI by using a field-tested canopy reflectance model. Then it calculates the BRFs by using the weight coefficient and the same model [4,8]. The algorithm tests the eligibility of a canopy radiation model to generate the LUT file where a subset of coefficients is satisfied within a given accuracy [9]. The given atmosphere-corrected BRFs are then compared with the modeled BRFs, which are stored in the biome-specific LUT files. Finally, all candidates of LAI/FPAR are used to calculate the mean values and uncertainty of the retrieval [9]. In the case of highly dense canopies, reflectance will be saturated and insensitive to changes in canopy properties. Therefore, LAI and FPAR values acquired under saturated conditions are less reliable than those generated by unsaturated BRFs. When the main algorithm fails to localize a solution, the backup algorithm is used to retrieve values through an empirical relationship between the normalized difference vegetation index (NDVI) and the canopy LAI/FPAR [11,21]. Such retrievals are flagged in the algorithm path quality assessment (QA) variable [8], which consists of two values for the main algorithm and two values for the backup algorithm (from high quality to low): the main algorithm without saturation (QA = 0), the main algorithm with saturation (QA = 1), the backup algorithm due to sun–sensor geometry (QA = 3), and the backup algorithm due to other reasons (QA = 4) [9,14,49].

2.2. Three-Dimensional Grassland Scene Simulation

We used the newly proposed but well validated 3D RT model LESS to simulate the interaction between the solar radiation and landscape elements based on the spectral response functions (SRFs) (from ENVI software) of MODIS and calculated the scene BRFs [35,45]. LESS simulates BRFs by a weighted forward photon tracing method as well as simulated energy transfer and generates images by a backward path tracing method [35]. Qi et al. [35] described the comparison between BRFs simulated by LESS and average BRF results from other models (e.g., SPRINT3, RAYTRAN, and RAYSPREAD) over several different homogeneous and heterogeneous canopies from the RAMI website to evaluate the accuracy of LESS.

The input parameters of LESS include 3D landscape elements, optical properties, and sun–sensor geometries. The simulated scenes are covered by grass (Johnson grass) and its component spectra were obtained from the LOPEX93 dataset on the OPTICLEAF website [50]. The soil (grayish brown loam) spectra were selected from the soil spectral library in ENVI software, and the transmittance of the soil is 0 (Figure 1). Then we calculated the two MODIS bands (red: band 1 and NIR: band 2) reflectance and transmittance (see

Table 1. Broad-band reflectance and transmittance of grass and soil used in this study. R and T are abbreviations for broad-band reflectance and transmittance, respectively.

	R (Red)	T (Red)	R (NIR)	T (NIR)
Johnson grass	0.0738	0.0577	0.4276	0.4607
Grayish brown loam	0.1755	0	0.3021	0

) using SRFs (the shaded part of Figure 1) of the MODIS sensor by the following equation [51]:

$$\{\begin{array}{c}R={\displaystyle {\displaystyle \sum }_{\lambda = {\lambda }_{min}}^{{\lambda }_{max}}} S_{\lambda }{R}_{\lambda }/ {\displaystyle {\displaystyle \sum }_{\lambda = {\lambda }_{min}}^{{\lambda }_{max}}} S_{\lambda }\\ T={\displaystyle {\displaystyle \sum }_{\lambda = {\lambda }_{min}}^{{\lambda }_{max}}} S_{\lambda }{T}_{\lambda }/{\displaystyle {\displaystyle \sum }_{\lambda = {\lambda }_{min}}^{{\lambda }_{max}}} S_{\lambda }\end{array}$$

(1)

where, R and T are MODIS band reflectance and transmittance, respectively. The $R_{\lambda }$ and $T_{\lambda }$ are mean narrow-band reflectance and transmittance derived from the spectral curves. The $S_{\lambda }$ is the SRF value of the MODIS sensor. $\lambda$ is the value of wavelength, which has a specific upper (${\lambda }_{max}$: red = 670 $\mu m$, NIR = 876 $\mu m$) and lower (${\lambda }_{min}$: red = 620 $\mu m$, NIR = 841 $\mu m$) limit for each band.

The 3D landscape elements were created with the third-party software OnyxTree, which uses the calculated reflectance and transmittance (Table 1) to make a grass 3D model (obj format file). As shown in Figure 2, we created nine randomly distributed grasslands with different LAIs (0.25, 0.50, 0.75, 1.0, 1.25, 1.5, 2.5, 3.5, and 4.5) using the LESS and grass 3D model. Moreover, LESS calculates FPAR by performing a band integration of the PAR between 380 nm and 710 nm and dividing by the incident radiation (slightly different from MODIS for which the wavelength interval is 400–800 nm) based on the LESS simulation of the collision of photons and the transfer of energy. In addition, to match the canopy structure of grasses in the MODIS LAI/FPAR retrieval algorithm (all organs other than leaves are ignored), only foliage is present in the scene. There is also only direct radiation in these scenes. The size of these scenes is 500 m × 500 m, which matches the spatial resolution of the MODIS LAI/FPAR products.

2.3. Experimental Design

We utilized the standard deviation of all LAI/FPAR candidates (StdLAI and StdFPAR), the retrieval index (RI), and the relative and absolute LAI/FPAR differences as the indicators of LAI/FPAR uncertainty. According to the uncertainty theory, StdLAI and StdFPAR are the standard deviations of all acceptable LAI/FPAR solutions in the LUT, which are the function of both the input uncertainty (biome type and BRF uncertainty) and model uncertainty [4,8]. StdLAI and StdFPAR have been proven and evaluated as quality metrics for MODIS LAI/FPAR products [28,52]. However, these two metrics have limitations due to the regularization introduced by the LUT algorithm and are artificially lowered at large LAIs [8]. Therefore, in this paper, we have also selected the RI (see Equation (2)) as an uncertainty metric, which is defined as the percentage of pixels for which the main RTE-based algorithm generates retrieval results. We note that the RI is used to characterize the overall uncertainty of all pixels [21,25,53], while StdLAI and StdFPAR are used to characterize individual pixel uncertainty.

$$\mathrm{RI}=\frac{Number of pixels retrieved by the main algorithm }{Total number of processed pixels }$$

(2)

To evaluate the consistency between true LAI/FPAR and MODIS retrievals, the difference between the simulation results of LAI (input to the LESS)/FPAR (output from the LESS) and the LAI/FPAR retrieved by the MODIS retrieval algorithm were used. The relative difference (RD, see Equation (3)) and absolute difference (AD, see Equation (4)), were utilized to quantify any differences.

$$\mathrm{RD}=\left(\mathrm{Retrieval}-\mathrm{Truth}\right)/\mathrm{Truth}$$

(3)

$$\mathrm{AD}=\mathrm{Retrieval}-\mathrm{Truth}$$

(4)

Based on the uncertainty theory, the retrieval uncertainty is a function of both model and input uncertainty and is embedded in the MODIS algorithm. In this study, we explored the relationship between retrieval uncertainty and the retrieval space, sun–sensor geometry, surface reflectance uncertainty, and biome type uncertainty using the variable-controlling approach (see

Table 2. Parameter configuration for designed experiments. SZA, SAA, VZA, VAA means solar zenith angle, solar azimuth angle, view zenith angle, and view azimuth angle, respectively.

Experiment	LAI	SZA	SAA	VZA	VAA	Uncertainty Metrics
Retrieval Space	/	0°	0°	0°	0°	StdLAI, StdFPAR
Sun–Sensor Geometry	0.50, 1.5, 3.5	30°/ 0°:10°:60°	90°/ 0°:30°:330°	0°:10°:60°/30°	0°:30°:330°/ 90°	RD, StdLAI, StdFPAR
BRF Uncertainty	1.5	0°	0°	−60°:10°:60°	0°	RD, StdLAI, StdFPAR
Biome Type Uncertainty	0.25, 0.50, 0.75, 1.0, 1.25, 1.5, 2.5, 3.5, 4.5	0°	0°	0°:10°:60°	0°:30°:330°	RI, AD
Clumping and Scale Effect	1.5	30°	0°	0°:30°:60°	0°:60°:300°	RI, StdLAI, StdFPAR

). We analyzed the inherent model uncertainty in two steps: 1) analysis of the retrieval space; 2) uncertainty changed with sun–sensor geometry. We obtained 4000 red-NIR BRF pairs by adding normally distributed errors (errors with 5% and 15% standard deviation) to the LESS simulated red and NIR band BRFs (1000: red without uncertainties and NIR with 5% standard deviation, 1000: red without uncertainties and NIR with 15% standard deviation, 1000: NIR without uncertainties and red with 5% standard deviation, and 1000: NIR without uncertainties and red with 15% standard deviation). Then we analyzed the LAI/FPAR uncertainty caused by BRF uncertainty within the 4 groups of samples. In addition, we analyzed uncertainties due to biome type misclassification, which is one of the main factors affecting the LAI/FPAR retrieval accuracy [4,54]. Each red-NIR BRF pair was sequentially combined with each biome type as the inputs for the MODIS LAI/FPAR algorithm. In this experiment, only B1 (grasses and cereal crops) was correct while the remaining seven combinations represented the biome type misclassification cases. Finally, we analyzed the influence of scale dependency and clumping effect (“tree groups”) [55] on the uncertainty of LAI/FPAR retrievals. We simulated a randomly distributed 1 km × 1 km scene (Figure 9a-1) and two clumping 1 km × 1 km scenes. Clumping type 1 (CT1, Figure 9a-2) had random clumping and Clumping type 2 (CT2, Figure 9a-3) was half bare ground and half grass. The LAI of the three scenes remained constant and these scenes were downscaled into four 500 m × 500 m scenes, sixteen 250 m × 250 m scenes, and one hundred 100 m × 100 m scenes for the discussion of scale dependency and clumping effect.

3. Results

3.1. Inherent Model Uncertainty

To evaluate the inherent model uncertainty of the MODIS LAI/FPAR retrieval algorithm, we analyzed the effect of the retrieval space and sun–sensor uncertainty, separately. In performing the evaluation of the retrieval space, we paid more attention to the changes in the uncertainty of algorithm paths. While the difference between the LAI/FPAR retrieval and LESS simulations were analyzed when evaluating of the sun–sensor geometry.

3.1.1. Analysis of Retrieval Space

Figure 3 indicates the variation of LAI/FPAR and its uncertainty in the retrieval space. As we can see, LAI/FPAR is nonlinearly related to surface reflectance (Figure 3a,b), and FPAR is also nonlinearly related to LAI. Moreover, the relationship between LAI/FPAR and its uncertainty (StdLAI and StdFPAR) is also nonlinear. The StdLAI and StdFPAR are very low for lower LAI/FPAR and then increase to the highest values, and then steadily decrease (from the bottom right to the top left of Figure 3d,e) as the LAI/FPAR gets progressively larger (from the bottom right to the top left of Figure 3a,b). It is also obvious that there is a clear division between the saturated (QA = 1) and unsaturated (QA = 0) parts where the LAI/FPAR values are higher in the saturated part (Figure 3c). Compared to the unsaturated part, the bias of LAI (+4.64) and FPAR (+0.631) are high, but the bias of StdLAI (−0.052) and StdFPAR (−0.169) are low in the saturated part. Figure 3d,e show that StdLAI and StdFPAR are relatively small at the boundaries of the area retrieved by the main algorithm due to the regularization of the algorithm [4,8].

3.1.2. Retrieval Uncertainty as a Function of Sun–Sensor Geometry

The relationship between LAI/FPAR uncertainty and sun–sensor geometry are presented in Figure 4 and Figure 5. In the high LAI scene (Figure 4a, LAI = 3.50), the retrieval results of LAI/FPAR show low consistency with the truth (it yields to an overall uncertainty of 20.01% for RD of LAI and 13.96% for RD of FPAR). The main algorithm shows an averaged 6.61% underestimation of LAI, while the backup algorithm results in an averaged 84.85% overestimation of LAI. In this scene, the backup algorithm appears at the “hotspot” geometry and where the difference between SAA and VAA is large. It can also be seen that the large VZA will lead to saturation. Nevertheless, for the low LAI scene (LAI = 0.50), the retrieved LAI/FPAR showed a significant overestimation (+111.86% RD for LAI, +162.50% RD for FPAR) and large uncertainty (StdLAI = 0.285, and StdFPAR = 0.238). Figure 5 shows the same analysis as above, but we controlled the view position and varied the sun position. Comparing Figure 4 and Figure 5, the distribution of LAI and its uncertainty (Figure 4a,c, and Figure 5a,c) show higher consistency, while FPAR and its uncertainty (Figure 4b,d and Figure 5b,d) are slightly different.

3.2. Input BRF Uncertainty

Here, we calculated the effects of input BRF uncertainty on the LAI/FPAR retrieval. Figure 6 shows that the uncertainty in the LAI/FPAR in the shadow area of the 15% BRF uncertainty is much larger, which means that larger BRF uncertainty will result in larger LAI/FPAR uncertainty. The StdLAI and StdFPAR due to a 5% BRF uncertainty are close to the StdLAI and StdFPAR due to a 15% BRF in the red band. This is because both 5% and 15% BRF uncertainty in the red band will trigger the backup algorithm with no StdLAI and StdFPAR. Comparing the shadow area in panels (a) and (b), we found that the same level of uncertainty in the NIR band BRF has a greater impact on the retrieval than the red band BRF. The main algorithm was not used in the hotspot (VZA = 0) geometry leading to the absence of both StdLAI and StdFPAR in panel (a).

3.3. Input Biome Type Uncertainty

Different biome types have different canopy structures, and the MODIS retrieval algorithm uses photon transport theory and the corresponding RT model for different biome types to parameterize the canopy structures (e.g., reflectance and transmittance of leaves, crown shadowing), which form the LUTs of the MODIS retrieval algorithm. To check the sensitivity of the algorithm to biome type, we modified the input biome type for the retrieval algorithm from the correct type (B1: grasses and cereal crops) to incorrect types (B2: shrubs; B3: broadleaf crops; B4: savannas; B5: evergreen broadleaf forests; B6: deciduous broadleaf forests; B7: evergreen needleleaf forests; and B8: deciduous needleleaf forests). As seen from Figure 7 and Figure 8, the retrieval uncertainty is similar when the input biome types are non-forest biomes (B1-B4) with a greater than 59.5% RI for all four biome types except for B2 in the LAI = 4.5 scene. However, the RI gets much lower when the grassland pixel is misclassified into forest biomes. As shown in Figure 7c and Figure 8c, the uncertainties of the retrieved LAI are high at the scene with high LAI (e.g., LAI = 3.5, 4.5). A significant overestimation (+0.727, +1.434 for AD of LAI) in B2, and a significant underestimation (−1.608, −2.344 for AD of LAI) in B3 is also evident. For B5, the RI is high (>69%) but AD of LAI (>1.656) is also high when LAI is relatively high (e.g., LAI = 2.5, 3.5, 4.5). As shown in Figure 8c, the FPAR calculated from the MODIS algorithm is significantly overestimated for all cases except for B4 high LAI scenes, which appear to be underestimated.

3.4. Impact of Clumping Effect and Scale Dependency

The model scale dependency and clumping effect have attracted much attention from the community in the development of quantitative remote sensing. In this experiment, the model scale dependency refers to the discrepancy between LAI/FPAR uncertainties that are derived from the same algorithm but at different spatial resolutions. The scale dependency determines the adaptive capacity of an algorithm for different pixel size. The clumping effect refers to the discrepancy between retrieved LAI/FPARs with same LAI/FPAR truth but different vegetation spatial distributions. The model nonlinearly and surface heterogeneity together result in the well-known phenomenon called “Inversion first and aggregation later is different from aggregation first and inversion later” [11].

Comparing the algorithm performance at different scales, we found that the MODIS algorithm is nearly scale-invariant from 100 m to 1000 m. Both LAI/FPAR and their uncertainty nearly remain unchanged with increasing pixel size except for the CT2, which shows that the StdLAI and StdFPAR are lower than other scales (

Table 3. The uncertainty metrics of three different clumping scenes and four scales.

	Scene	100 m	250 m	500 m	1000 m
RI (N. of main/N. of all)	Uniform	1700/1800	272/288	68/72	17/18
	CT1	1743/1800	279/288	70/72	18/18
	CT2	1565/1800	251/288	62/72	17/18
StdLAI (mean ± Std)	Uniform	0.147± 0.019	0.148± 0.019	0.149± 0.020	0.150± 0.021
	CT1	0.251± 0.108	0.218± 0.074	0.225± 0.066	0.179± 0.074
	CT2	0.340± 0.160	0.339± 0.160	0.342± 0.159	0.181± 0.027
StdFPAR (mean ± Std)	Uniform	0.088± 0.012	0.088± 0.012	0.089± 0.012	0.089± 0.013
	CT1	0.208± 0.109	0.177± 0.070	0.180± 0.058	0.144± 0.061
	CT2	0.269± 0.188	0.269± 0.188	0.271± 0.187	0.130± 0.022

). The retrieved LAIs for all scenes are less than the LAI truth at 1000 m scale, and the underestimations for Uniform, CT1, and CT2 vegetation distributions are −0.005, −0.846, and −0.525, respectively. Comparing the three clumping scenes, we found that the LAI of a uniform scene is very close to the LAI truth (Figure 9b-1). CT1 shows a significant underestimation, while CT2 shows a significant overestimation except in the 1000 m scale. For FPAR, there is a significant overestimation in all three scenes (Figure 9b-2). At the same spatial resolution, the RI of CT1 is the highest, followed by Uniform, and the lowest is CT2 (Table 3). While the values of StdLAI and StdFPAR are as follows (from small to large): Uniform, CT1, and CT2. The standard deviations of StdLAI and StdFPAR also get larger in this order.

4. Discussion

Because of the different sensitivities of LAI/FPAR to surface reflectances, we note that there would be a gap of uncertainty between the saturated part and the unsaturated part [4,8]. However, Figure 3 indicates that for large LAI/FPAR, their theoretical uncertainty is artificially reduced by the method of regularization, which causes the retrieval to have varying degrees of confidence and leads to a problematic evaluation of high LAI/FPAR scenes using the provided StdLAI and StdFPAR. This also places new requirements on future algorithm refinement that the LUT algorithm should be consistent in the saturated case as in the unsaturated case. The LAI/FPAR values estimated by the backup algorithm and calculated by the main algorithm also show significant discontinuity [4,8,9,10] (Figure 4a). Based on this, we point out that future algorithm refinement should increase the coverage of the main algorithm usage, which will greatly improve the overall accuracy of the product. In addition, according to the 3D RT model, the hotspot means that the radiation field tends to peak around the retro-illumination direction. The results of this study indicate that the uncertainty of the MODIS algorithm in the hotspots is quite large (Figure 4, Figure 5 and Figure 6), due to which the science team decided not to include additional hotspot parameters since their inclusion would make algorithm calibration difficult [56,57]. We note that this will not cause large problems in the MODIS LAI/FPAR production because of the fact that observations near the hotspot are rare for MODIS. However, this points out a new refinement direction of this algorithm to improve the accuracy of hotspot modeling for other sensors.

As is known, the uncertainty of the inputting BRFs has some influence on the uncertainty of the retrieval algorithm. In particular, our results show that the uncertainty of NIR BRFs has a larger effect on LAI/FPAR uncertainty compared to the red BRFs (Figure 6). We know that insufficient input information will lead to the “ill-posed” retrieval problem [11]; however, the inputting BRFs of the MODIS operational algorithm are currently only for the red and NIR bands. Therefore, in the future we may try to make use of BRFs in other bands to improve the retrieval accuracy. The MODIS algorithm depends on a priori information about the land surface given by biome type representing the pattern of the architecture of vegetation, as well as patterns of spectral reflectance and transmittance of vegetation [8]. Figure 7 and Figure 8 confirm that the misclassification of biome types with similar structures will result in smaller LAI/FPAR uncertainty, and vice versa [11,58]. This means that the improvement of biome classification accuracy is an efficient way to improve the LAI/FPAR products. Moreover, different biome types also lead to different clumping types. As Figure 9 shows, the underestimation of LAI is significant for two clumping scenes at 1000 m scale. For the other three scales, however, the overestimation of CT2 is due to the backup algorithm retrievals. As our results show, the algorithm only considers the clumping effect at one scale (e.g., B1 is minimal leaf clumping) [4], which can result in large differences in the retrievals; therefore, we suggest that future algorithms consider the clumping effect at more scales (e.g., leaf, branch, and crown).

We note that there are some problems with the way we use LESS to simulate specific scenes and evaluate the MODIS algorithm. First, according to the algorithm, the retrieved LAI/FPAR is a weighted average of the probability values within the error range. Therefore, the probability distribution of LAI/FPAR within the error range based on a great number of realizations has more statistical significance thus may differ from the specific realization (scene) that was used. Secondly, although the LESS model has been well validated, the confidence of our evaluation results depends on the accuracy of the LESS simulation.

In short, validation in the field of remote sensing utilizing computer simulations has proved feasible. In future studies, we will analyze the other seven biome types, which will provide a more comprehensive evaluation of the MODIS LAI/FPAR retrieval algorithm. In addition, we will change the mode of a single specific scene to obtain retrieval results by simulating multiple scenes. Moreover, evaluation of the algorithm at different levels of vegetation clumping will be the focus of our future research.

5. Conclusions

This paper presents an uncertainty assessment of the MODIS LAI/FPAR retrieval algorithm over B1 (grassland) based on computer simulation. To accomplish this assessment, we first analyzed the theoretical uncertainty caused by inherent model uncertainty, then we calculated the uncertainty caused by input parameters (BRF and biome type) over simulated 3D grass scenes. Finally, we analyzed the effects of vegetation clumping and scale dependency of the MODIS algorithm. The 3D grass scenes were simulated by a well validated 3D RT model (LESS), which helps to separate the model uncertainty and other uncertainties. We found that the uncertainty of the main and backup algorithm varies considerably. In the same scene, there is a −6.61% bias for the main algorithm retrieval, while the backup algorithm retrieval has a +84.85% bias. We noted that the uncertainty of the saturated retrievals is artificially reduced compared with unsaturated retrievals. At the same time, MODIS showed significant overestimation at low LAI scenes, with a maximum bias of +111.86% for LAI and +162.50% for FPAR. In the high LAI scenes, the “hotspot” geometry results in greater retrieval uncertainty from the backup algorithm. Moreover, input uncertainties further increased the uncertainty of LAI/FPAR retrieval. We found that the uncertainties in BRF in the NIR band has a greater impact than in the red band. The biome type uncertainty also leads to great retrieval uncertainty. Large uncertainties occurred when grassland was misclassified into forest biomes, while smaller uncertainties occurred when the misclassification was within the non-forest biomes. In addition, the clumping effect results in underestimation (−0.846 and −0.525 for the two clumping types, respectively) and we found that the MODIS algorithm is nearly scale-invariant from 100 m to 1000 m pixel sizes. Overall, these results, based on novel computer simulation experiments, can guide the future refinements of the MODIS LAI/FPAR algorithm.