Analysis of Extracting Prior BRDF from MODIS BRDF Data

Many previous studies have attempted to extract prior reflectance anisotropy knowledge from the historical MODIS Bidirectional Reflectance Distribution Function (BRDF) product based on land cover or Normalized Difference Vegetation Index (NDVI) data. In this study, the feasibility of the method is discussed based on MODIS data and archetypal BRDFs. The BRDF is simplified into six archetypal BRDFs that represent different reflectance anisotropies. Five-year time series of MODIS BRDF data over three tiles are classified into six BRDF archetype classes according to the Anisotropy Flat indeX (AFX). The percentage of each BRDF archetype class in different land cover classes or every 0.1-NDVI interval is determined. Nadir BRDF-Adjusted Reflectances (NBARs) and NDVIs simulated from different archetypal BRDFs and the same multi-angular observations are compared to MODIS results to study the effectiveness of the method. The results show that one land cover type, or every 0.1-NDVI interval, contains all the potential BRDF shapes and that one BRDF archetypal class makes up no more than 40% of all data. Moreover, the differences between the NBARs and NDVIs simulated from different archetypal BRDFs are insignificant. In terms of the archetypal BRDF method and MODIS BRDF product, this study indicates that the land cover or NDVI is not necessarily related to surface reflectance anisotropy.


Introduction
The reflectance of most natural surfaces is anisotropic, resulting in major difficulty regarding the quantitative estimation of vegetation and soil characteristics from ground-based and remotely sensed observations.The Bidirectional Reflectance Distribution Function (BRDF) is used to describe the characteristics of reflectance anisotropy [1].In practice, the Kernel-driven, semi-empirical RossThick-LiSparse-Reciprocal (RTLSR) BRDF model is widely used for BRDF retrieval and reconstruction based on remotely sensed data [2][3][4][5][6].
The BRDF depends on the wavelength and is determined based on the optical properties and structure of the surface [6].It has been widely used for normalizing satellite measurements in the nadir direction and retrieving land surface albedo from sparse angular observations [3,5], performing coupled atmospheric correction [7] and land cover classification [8], and deriving canopy structure and other bio-geophysical parameters [3].Generally, sufficient and well-distributed observations are needed to accurately retrieve BRDF from multi-angular remotely sensed data [9].However, constrained by the observation capacity of sensors, the orbital characteristics of observation platforms [10], clouds, and their shadows, most spacecraft or satellite remote sensors cannot collect sufficient and well-distributed observations over short periods.This constraint is more remarkable for high and medium spatial resolution remotely sensed data because most remote sensors only have near-nadir observation capability.
When observations are sparsely sampled, including the case of measuring only a single observation, we can improve the inversion accuracy of the BRDF and other relevant parameters using an a priori BRDF [11].This is mainly because prior knowledge can limit the variability in the BRDF during the retrieval process.Many investigations have been performed in this field.In the first version of the operational MODIS BRDF/albedo backup algorithm, the prescribed BRDF associated with specific land cover types is used as prior knowledge [6,9].Currently, in the V006 collection, the most recent full inversion BRDF is used as prior knowledge in the BRDF/albedo backup algorithm [12].Moreover, according to Bayesian inference theory, a priori knowledge has been used to improve the retrieval of surface bidirectional reflectance and spectral albedo from satellite observations [11].The BRDF parameters of Multi-angle Imaging SpectroRadiometer (MISR) data were used as prior knowledge to improve the retrieval of the surface BRDF from MODIS observations [2].
Recently, to extract prior BRDF knowledge from the historical BRDF product, many studies have linked surface reflectance anisotropy to land cover or the NDVI [13][14][15]; both types of extracted BRDF knowledge perform well in the retrieval of land surface albedo from nadir reflectance.Land cover is the physical and biological cover on the surface of the Earth, and land surface parameters have been commonly retrieved from remotely sensed data [6,16].The NDVI is a normalized ratio of the near-infrared (NIR) and red bands [17] and is one of the most commonly used vegetation indices [18].However, both land cover and NDVI data are generally calculated from a single directional reflectance; thus, these two data sources contain limited reflectance anisotropy information.Recent studies also indicated a weak correlation between reflectance anisotropy and land cover or the NDVI, e.g., the variability in the estimated POLDER BRDF model parameters for several International Geosphere-Biosphere Program (IGBP) land cover classes were found to be higher within a class than between classes [19].Additionally, the NDVI tends to be saturated in high biomass regions such as dense forest areas [18] that generally have complicated and changing structures, and the NDVI is approximately orthogonal to the Anisotropy Flat indeX (AFX), which can indicate varying reflectance anisotropy [20].
The archetypal BRDF database, which contains six archetypal BRDFs for each spectrum of MODIS, was established based on AFX theory and the MODIS BRDF product [20].The representativeness of these archetypal BRDFs for naturally occurring BRDFs has been proven using the MODIS BRDF product and actual MODIS multi-angular observations [21].These archetypal BRDFs have also been applied in albedo retrieval from small view-angle airborne observations [22].The AFX and archetypal BRDFs offer methods for quantitatively classifying reflectance anisotropy.In this study, archetypal BRDFs were used to represent different levels of surface reflectance anisotropy, and the Version 005 MODIS BRDF/albedo product and actual multi-angular observations were used as study data.The percentage of each BRDF archetype class within several IGBP land cover classes or every 0.1-NDVI interval was used to analyze the feasibility of extracting prior reflectance anisotropy form MODIS BRDF data.An analysis of the difference between the NBARs or NDVIs retrieved from different archetypal BRDFs was performed to assess the effectivity of the approach.This study has important implications for extracting prior BRDF knowledge from the historical MODIS BRDF product.

MODIS Products and Multi-Angular Observations
MODIS BRDF/albedo products (MCD43A) of the global land surface have been routinely available since 2000 [5,23].Time series MODIS BRDF/albedo products over three MODIS tiles from 2008 to 2012 are used in this study.Tile h11v03 is located in northwestern Canada and contains grass, savannah, and forest.Tile h20v11 is located in southern Africa and contains grassland, shrubland, and savannah.Tile h17v06 is located in northwestern Africa, and most of the area of this tile is desert.The MODIS Land Cover Type product (MCD12Q1), which maps global land cover using spectral and temporal information derived from MODIS [24], contains multiple classification schemes.Many studies have confirmed that the accuracy of the land cover product from MODIS is approximately 80% [8,[24][25][26].The NDVI, which can be used to analyze remote sensing measurements and indicate the land cover vegetation condition, is also a routine product of MODIS.The land cover data defined by the IGBP global vegetation classification scheme and the NDVI data are used in this study.
In addition to the abovementioned data, multi-angular data of the cloud-free, atmospherically corrected, high-quality MODIS observations (MCD-obs) [5,9] over tile h20v11 from days 201-216 in 2008 are used for further analysis of the correlation between reflectance anisotropy and the NDVI.The MCD-obs are extracted from the MODIS reflectance products (MOD09GA and MYD09GA) using the operational MODIS BRDF/albedo algorithm [5].Finally, high-quality, snow-free observations of more than four million samples are used in this study.

AFX and Archetypal BRDF Database
The AFX value of each pixel can be easily calculated from the RTLSR BRDF model parameters using Equation (1) [20]: where f iso , f vol , and f geo represent the spectrally dependent BRDF model parameters; the constants 0.189184 and −1.377622 are the bi-hemispherical integrals of the RossThick and LiSparse Reciprocal kernels, respectively; and λ represents the MODIS spectrum.
The AFX provides a feasible method of classifying reflectance anisotropy.As the AFX increase, the land surface scattering pattern changes from geometric optical scattering (AFX < 1) to volumetric scattering (AFX > 1).The shapes of corresponding archetypal BRDFs also change from dome shaped to bowl shaped accordingly.Based on this feature, the archetypal BRDF database [20] was established from the high-quality MODIS BRDF product over the Earth Observing System (EOS) Land Validation Core Sites [27].The BRDF model parameters and corresponding AFX ranges of the BRDF archetype class for the red and near-infrared bands are shown in Table 1 [20].The normalization of BRDF removes the effects of spectral reflectance differences by multiplying the original BRDF by a scale factor K = 0.5/f iso .Both the establishment of the parameters and the normalization method are detailed in [20].The surface reflectance anisotropy changes from strong geometric scattering to strong volumetric scattering as the archetypal BRDFs change from one to six.

Proportion of Each BRDF Archetype Class
The flowchart for calculating the proportion of each BRDF archetype class within different land covers or every 0.1-NDVI interval is shown in Figure 1.The AFX value of each pixel can be calculated from the MODIS BRDF product.The MODIS BRDF product can be divided into six BRDF archetype classes, according to the AFX range of each BRDF archetype class.Based on the number of pixels in each archetype class, the proportion of each BRDF archetype class within different land cover types (P i ) or every 0.1-NDVI interval (P i ) can be calculated.To remove the effect of viewing geometry, the NDVI that is recalculated from simulated NBARs at a certain solar zenith angle (SZA) is used in this process.

Proportion of Each BRDF Archetype Class
The flowchart for calculating the proportion of each BRDF archetype class within different land covers or every 0.1-NDVI interval is shown in Figure 1.The AFX value of each pixel can be calculated from the MODIS BRDF product.The MODIS BRDF product can be divided into six BRDF archetype classes, according to the AFX range of each BRDF archetype class.Based on the number of pixels in each archetype class, the proportion of each BRDF archetype class within different land cover types (Pi) or every 0.1-NDVI interval (P'i) can be calculated.To remove the effect of viewing geometry, the NDVI that is recalculated from simulated NBARs at a certain solar zenith angle (SZA) is used in this process.

Fitting Archetypal BRDFs to Multi-Angular Observations
Using the archetypal BRDFs as prior knowledge, we can predict a probable BRDF for a pixel under observation, and calculate the corresponding albedo.The simulated reflectance (ρ′), which has the same observation geometry as MCD-obs (ρ), can be calculated from the RTLSR BRDF model and the BRDF model parameters of the archetypal BRDF (F*): where θ, υ, and φ denote the illuminating and viewing directions of ρ.Kvol and Kgeo are kernel functions that describe volumetric [28,29] and geometric optical [30,31] scattering, respectively.Based on a multi-angular dataset that has a set of n observations, the factor a, which can adjust and best fit the prior archetypal BRDF (BRDF′) to the observations, can be calculated using the least square method [2,6]: The BRDF that best fits the observations can then be determined as follows: Additionally, the NBARs (ρnadir) associated with prior archetypal BRDFs at a certain SZA θ can be expressed as follows:

Fitting Archetypal BRDFs to Multi-Angular Observations
Using the archetypal BRDFs as prior knowledge, we can predict a probable BRDF for a pixel under observation, and calculate the corresponding albedo.The simulated reflectance (ρ ), which has the same observation geometry as MCD-obs (ρ), can be calculated from the RTLSR BRDF model and the BRDF model parameters of the archetypal BRDF (F * ): where θ, υ, and ϕ denote the illuminating and viewing directions of ρ.K vol and K geo are kernel functions that describe volumetric [28,29] and geometric optical [30,31] scattering, respectively.Based on a multi-angular dataset that has a set of n observations, the factor a, which can adjust and best fit the prior archetypal BRDF (BRDF ) to the observations, can be calculated using the least square method [2,6]: The BRDF that best fits the observations can then be determined as follows: Additionally, the NBARs (ρ nadir ) associated with prior archetypal BRDFs at a certain SZA θ can be expressed as follows: During the fitting process, the reflectance anisotropy characteristics remain constant because the prior BRDF and the adjusted BRDF have the same AFX value (Equations ( 1) and ( 4)).Therefore, the adjusted BRDFs retrieved from different archetypal BRDFs and the same MCD-obs can be regarded as the responses of the observations under different reflectance anisotropies.If we assume that the NDVI can be linked to reflectance anisotropy, then the NBARs and NDVIs that are retrieved from different archetypal BRDFs and the same MCD-obs should have obvious differences.The NDVI refers to the NBARs of the red and NIR bands; thus, to illustrate the relationship more clearly, the MODIS data that belong to a specific BRDF archetype class in the red band are used as the study dataset.During the fitting process, the reflectance anisotropy characteristics remain constant because the prior BRDF and the adjusted BRDF have the same AFX value (Equations ( 1) and ( 4)).Therefore, the adjusted BRDFs retrieved from different archetypal BRDFs and the same MCD-obs can be regarded as the responses of the observations under different reflectance anisotropies.If we assume that the NDVI can be linked to reflectance anisotropy, then the NBARs and NDVIs that are retrieved from different archetypal BRDFs and the same MCD-obs should have obvious differences.The NDVI refers to the NBARs of the red and NIR bands; thus, to illustrate the relationship more clearly, the MODIS data that belong to a specific BRDF archetype class in the red band are used as the study dataset.

Point Assessment
Figure 3 compares the shapes of the MODIS BRDF in the NIR band (a and c) for deciduous broadleaf forest (b) and evergreen needle forest (c) at two locations in tile h11v03.Concurrent high-resolution, truecolor SPOT satellite images of the corresponding area from Google Maps are also shown in Figure 3.All the images were collected in September 2006, and the locations of images are labeled in Figure 3.For the two images of the deciduous broadleaf forest (b), the one on the left is possibly mixed with other land cover types (e.g., shrubland) and is more heterogeneous because its canopy height covers a large dynamic range.For the two images of the evergreen needle forest (d), the image on the left has sparse trees, and the shadows of the trees and the ground can be observed clearly.However, the image on the right appears homogeneous.The two IGBP land cover classes in the two images shown on the left have remarkable geometric scattering and a dome-shaped BRDF, while the two on the right have typical volumetric scattering with a bowl-shaped BRDF.The BRDF shape is determined by the structure of the vegetation canopy [6].Nevertheless, the same IGBP land cover class may have

Point Assessment
Figure 3 compares the shapes of the MODIS BRDF in the NIR band (a and c) for deciduous broadleaf forest (b) and evergreen needle forest (c) at two locations in tile h11v03.Concurrent high-resolution, true-color SPOT satellite images of the corresponding area from Google Maps are also shown in Figure 3.All the images were collected in September 2006, and the locations of images are labeled in Figure 3.For the two images of the deciduous broadleaf forest (b), the one on the left is possibly mixed with other land cover types (e.g., shrubland) and is more heterogeneous because its canopy height covers a large dynamic range.For the two images of the evergreen needle forest (d), the image on the left has sparse trees, and the shadows of the trees and the ground can be observed clearly.However, the image on the right appears homogeneous.The two IGBP land cover classes in the two images shown on the left have remarkable geometric scattering and a dome-shaped BRDF, while the two on the right have typical volumetric scattering with a bowl-shaped BRDF.The BRDF shape is determined by the structure of the vegetation canopy [6].Nevertheless, the same IGBP land cover class may have a completely different structure and contain totally different BRDF shapes (Figure 3), while different land cover classes may have similar structures and reflectance anisotropy.

Spatial Assessment
We analyzed the composition of the BRDF within several IGBP land cover classes using the MODIS BRDF/albedo product from day 201 in 2008.The study area is located in the southwestern portion of tile h11v03.It contains a continuous region of high-quality MODIS data (approximately 1500 × 2000 pixels), and the latitude and longitude of the center of the study area are 54.164583 and −114.223074degrees, respectively.
The left part of Figure 4 shows the IGBP land cover maps of the study area, and the right part shows the normalized BRDF shapes of the red (above) and NIR (below) bands in the principle plane over a large area of evergreen needle leaf forest (ENF), mixed forest (MF), grassland, and cropland.The four different colors represent the four land cover types in the study area.To make the BRDFs comparable, the normalized BRDFs rather than the original BRDFs are shown in Figure 4.Only the normalized BRDFs that belong to BRDF archetype classes No. 1, No. 4, and No. 6 (from left to right) are shown.To generate this figure, the original BRDF is normalized by multiplying by a scale factor of 0.5/fiso [20], and the normalized BRDF shapes of ENF (red lines), MF (green lines), grassland (blue lines) and cropland (cyan lines) are shifted by adding 0, 0.2, 0.4, and 0.6, respectively.
Obviously, each of the four IGBP land cover classes in the study area contains various potential BRDF shapes in both bands, and different IGBP land cover classes generally include similar reflectance anisotropy.The parameters of geometric optical scattering and volumetric scattering are not be strictly independent because the two kernels are somewhat correlated [3].This may explain the discrete distribution of the MODIS BRDF within every BRDF archetype class.The BRDF shapes of cropland are highly consistent with each other, likely due to the homogeneous characteristics of cropland.

Spatial Assessment
We analyzed the composition of the BRDF within several IGBP land cover classes using the MODIS BRDF/albedo product from day 201 in 2008.The study area is located in the southwestern portion of tile h11v03.It contains a continuous region of high-quality MODIS data (approximately 1500 × 2000 pixels), and the latitude and longitude of the center of the study area are 54.164583 and −114.223074degrees, respectively.
The left part of Figure 4 shows the IGBP land cover maps of the study area, and the right part shows the normalized BRDF shapes of the red (above) and NIR (below) bands in the principle plane over a large area of evergreen needle leaf forest (ENF), mixed forest (MF), grassland, and cropland.The four different colors represent the four land cover types in the study area.To make the BRDFs comparable, the normalized BRDFs rather than the original BRDFs are shown in Figure 4.Only the normalized BRDFs that belong to BRDF archetype classes No. 1, No. 4, and No. 6 (from left to right) are shown.To generate this figure, the original BRDF is normalized by multiplying by a scale factor of 0.5/f iso [20], and the normalized BRDF shapes of ENF (red lines), MF (green lines), grassland (blue lines) and cropland (cyan lines) are shifted by adding 0, 0.2, 0.4, and 0.6, respectively.
Obviously, each of the four IGBP land cover classes in the study area contains various potential BRDF shapes in both bands, and different IGBP land cover classes generally include similar reflectance anisotropy.The parameters of geometric optical scattering and volumetric scattering are not be strictly independent because the two kernels are somewhat correlated [3].This may explain the discrete distribution of the MODIS BRDF within every BRDF archetype class.The BRDF shapes of cropland are highly consistent with each other, likely due to the homogeneous characteristics of cropland.5c and 6c).The red band generally has stronger geometric and weaker volumetric scattering than the NIR bands because chlorophyll provides strong absorption in the red band and strong reflection in the NIR band [30].

Temporal and Spatial Assessment
The proportions of different BRDF archetype classes in different IGBP land cover classes over time are used to further analyze the correlation between reflectance anisotropy and land cover.Cumulative histograms regarding the percentage of every BRDF archetype class in primary IGBP land cover classes are shown in Figures 5 and 6 for each tile in the red and NIR bands from 2008 to 2012.
The statistical results show that, except for the desert area, the studied land cover types contain all the potential BRDF shapes in the two bands.For the desert in tile h17v06 (Figure 5e 5c and 6c).The red band generally has stronger geometric and weaker volumetric scattering than the NIR bands because chlorophyll provides strong absorption in the red band and strong reflection in the NIR band [30].The changes in ground vegetation structure over time lead to corresponding changes in the composition of these archetypal BRDFs.For the vegetation area, the distribution of BRDF archetype classes is discrete in both bands because of the large diversity and complexity of the surface structure.The proportion of BRDF archetype class No. 6 in the red band is generally small; however, it is dominant in grassland and cropland in the NIR band during summer (Figure 6a,c).The results for ENF (Figure 5d) are more remarkable because each BRDF archetype class is almost equally weighted throughout most of the year.In addition, BRDF archetype classes No. 1 and No. 6 reach maximums simultaneously in the summer.Cropland (Figure 6c) accounts for the largest percentage of BRDF archetype class No. 6 in summer because crops have high fractional vegetation cover and a uniform spatial distribution.Meanwhile, during the winter, the relatively homogeneous BRDF archetype class No. 4 is dominant (~75%), as there is no vegetation cover.The changes in ground vegetation structure over time lead to corresponding changes in the composition of these archetypal BRDFs.For the vegetation area, the distribution of BRDF archetype classes is discrete in both bands because of the large diversity and complexity of the surface structure.The proportion of BRDF archetype class No. 6 in the red band is generally small; however, it is dominant in grassland and cropland in the NIR band during summer (Figure 6a,c).The results for ENF (Figure 5d) are more remarkable because each BRDF archetype class is almost equally weighted throughout most of the year.In addition, BRDF archetype classes No. 1 and No. 6 reach maximums simultaneously in the summer.Cropland (Figure 6c) accounts for the largest percentage of BRDF archetype class No. 6 in summer because crops have high fractional vegetation cover and a uniform spatial distribution.Meanwhile, during the winter, the relatively homogeneous BRDF archetype class No. 4 is dominant (~75%), as there is no vegetation cover.

Composition of Reflectance Anisotropy in 0.1-NDVI Intervals
The high-quality MODIS BRDF product in the selected tiles is used to analyze the correlation between archetypal BRDFs and the NDVI between 2008 and 2012.The percentage of each BRDF archetype class in the BRDF data of every 0.1-NDVI interval is used to quantitatively explore the associated correlation.All the MODIS datasets, the datasets of different land cover types, and the datasets of different growth periods are used to study the correlation.The statistical results calculated from different datasets are similar; therefore, only the result for all the MODIS datasets is shown as an example.
Figure 7 shows the percentage of each BRDF archetype class for different ranges of the NDVI using all the MODIS data in each tile from 2008 to 2012.There is no obvious relationship between the NDVI values and reflectance anisotropy.In tile h20v11 (Figure 7), most of the pixels have NDVI values between 0.2 and 0.5.In the red band (Figure 7a), the proportion of BRDF archetype classes No. 2 and No. 3 reach maximum levels (~35%) at NDVI values of 0.2-0.3.With the increase in the NDVI, the percentage of BRDF archetype class No. 1 increases first and then decreases, reaching a maximum value at an NDVI of 0.5.Meanwhile, the percentages of BRDF archetype classes No. 2 and No. 3 gradually decrease.When the NDVI is greater than 0.7, the percentages of BRDF archetype classes No. 1 to No. 4 are generally similar (10%~15%), while the proportion of BRDF archetype classes No. 5 and No. 6 increase gradually.In the NIR band (Figure 7b), one of the most remarkable features is that BRDF archetype classes No. 2 and No. 3 account for large proportions when the NDVI is small.When the NDVI is greater than 0.5, BRDF archetype class No. 6 accounts for a large proportion, while the other BRDF archetype classes remain stable.When the NDVI is greater than 0.8, the percentages of the six BRDF archetype classes tend to be similar (10%~20%).
The statistical result in tile h11v03 is similar to that in tile h20v11; however, in the red band (Figure 7c), BRDF archetype classes No. 4 and No. 5 account for large proportions when the NDVI is small, while in the NIR band (Figure 7d), the proportion of BRDF archetype class No. 6 reaches a maximum (~35%) when the NDVI is approximately 0.4.The other five BRDF archetype classes always have equal weights.

Composition of Reflectance Anisotropy in 0.1-NDVI Intervals
The high-quality MODIS BRDF product in the selected tiles is used to analyze the correlation between archetypal BRDFs and the NDVI between 2008 and 2012.The percentage of each BRDF archetype class in the BRDF data of every 0.1-NDVI interval is used to quantitatively explore the associated correlation.All the MODIS datasets, the datasets of different land cover types, and the datasets of different growth periods are used to study the correlation.The statistical results calculated from different datasets are similar; therefore, only the result for all the MODIS datasets is shown as an example.
Figure 7 shows the percentage of each BRDF archetype class for different ranges of the NDVI using all the MODIS data in each tile from 2008 to 2012.There is no obvious relationship between the NDVI values and reflectance anisotropy.In tile h20v11 (Figure 7), most of the pixels have NDVI values between 0.2 and 0.5.In the red band (Figure 7a), the proportion of BRDF archetype classes No. 2 and No. 3 reach maximum levels (~35%) at NDVI values of 0.2-0.3.With the increase in the NDVI, the percentage of BRDF archetype class No. 1 increases first and then decreases, reaching a maximum value at an NDVI of 0.5.Meanwhile, the percentages of BRDF archetype classes No. 2 and No. 3 gradually decrease.When the NDVI is greater than 0.7, the percentages of BRDF archetype classes No. 1 to No. 4 are generally similar (10%~15%), while the proportion of BRDF archetype classes No. 5 and No. 6 increase gradually.In the NIR band (Figure 7b), one of the most remarkable features is that BRDF archetype classes No. 2 and No. 3 account for large proportions when the NDVI is small.When the NDVI is greater than 0.5, BRDF archetype class No. 6 accounts for a large proportion, while the other BRDF archetype classes remain stable.When the NDVI is greater than 0.8, the percentages of the six BRDF archetype classes tend to be similar (10%~20%).
The statistical result in tile h11v03 is similar to that in tile h20v11; however, in the red band (Figure 7c), BRDF archetype classes No. 4 and No. 5 account for large proportions when the NDVI is small, while in the NIR band (Figure 7d), the proportion of BRDF archetype class No. 6 reaches a maximum (~35%) when the NDVI is approximately 0.4.The other five BRDF archetype classes always have equal weights.

Comparison of NBAR and NDVI Values
The MCD-obs used in this study were extracted from the MODIS reflectance product through the operational MODIS BRDF/albedo algorithm.The number of pixels and the percentages of each BRDF archetype class in the red and NIR bands used in this study are shown in Table 2.   2), the NBARs retrieved from one specific archetypal BRDF can still exhibit good consistency with MODIS NBARs.The SZA has a large effect on the retrieval of NBARs.When the SZA is 15 • , the two NBARs are relatively different, and the coefficients of determination (R 2 ) are less than 0.87 and 0.9 in the red and NIR bands, respectively.However, as the SZA increases, the consistency between the two NBARs tends to increase, and the R 2 values are over 0.92 and 0.98.The large difference at a small SZA is mainly caused by the hot spot effect of the BRDF [32,33].Moreover, when the SZA is small, archetypal BRDF No. 1 (No. 6) tends to slightly overestimate (underestimate) the NBARs, compared to the MODIS NBARs; with increasing SZA, this difference tends to disappear.

Comparison of NBAR and NDVI Values
The MCD-obs used in this study were extracted from the MODIS reflectance product through the operational MODIS BRDF/albedo algorithm.The number of pixels and the percentages of each BRDF archetype class in the red and NIR bands used in this study are shown in Table 2.   2), the NBARs retrieved from one specific archetypal BRDF can still exhibit good consistency with MODIS NBARs.The SZA has a large effect on the retrieval of NBARs.When the SZA is 15°, the two NBARs are relatively different, and the coefficients of determination (R 2 ) are less than 0.87 and 0.9 in the red and NIR bands, respectively.However, as the SZA increases, the consistency between the two NBARs tends to increase, and the R 2 values are over 0.92 and 0.98.The large difference at a small SZA is mainly caused by the hot spot effect of the BRDF [32,33].Moreover, when the SZA is small, archetypal BRDF No. 1 (No. 6) tends to slightly overestimate (underestimate) the NBARs, compared to the MODIS NBARs; with increasing SZA, this difference tends to disappear.

Discussion
In this study, the feasibility of extracting the BRDF feature from MODIS BRDF data based on land cover or the NDVI was analyzed using the five-year time series MODIS BRDF/albedo product in three tiles.The AFX value is a new quantitative index that was used to classify reflectance anisotropy.Six archetypal BRDFs established using the MODIS BRDF product and AFX provide different reflectance anisotropy patterns.MODIS BRDF data were classified into six BRDF archetype classes, and

Discussion
In this study, the feasibility of extracting the BRDF feature from MODIS BRDF data based on land cover or the NDVI was analyzed using the five-year time series MODIS BRDF/albedo product in three tiles.The AFX value is a new quantitative index that was used to classify reflectance anisotropy.Six archetypal BRDFs established using the MODIS BRDF product and AFX provide different reflectance anisotropy patterns.MODIS BRDF data were classified into six BRDF archetype classes, and the percentage of each BRDF archetype class in the primary land cover classes and every 0.1-NDVI interval was used to explore the correlation of reflectance anisotropy with land cover and the NDVI.Moreover, each of these archetypal BRDFs was used to fit the same MODIS multi-angular observations, and the simulated NBARs and NDVIs were compared to MODIS results to study this correlation.
The reflectance anisotropy of the land surface is determined based on the surface structure, such as shadow casting and the spatial distribution of components.Land cover and NDVI data, which are generally retrieved from directional reflectance, rarely contain vegetation density and canopy structure information.The density, height, and spatial distribution of the same land cover or similar NDVIs may exhibit large differences.Therefore, they may have completely different BRDF shapes.Our study proved that, except for the desert area, the selected land cover classes or every 0.1-NDVI interval contain various BRDF shapes, and one BRDF archetype class rarely dominates.From the opposite perspective, although the six archetypal BRDFs exhibit obvious differences, the NBARs and NDVIs retrieved from these BRDFs and the same multi-angular observations are similar, especially when the SZA is relatively large.Only the NBARs were used in this comparison, and the reflectance in the other direction may lead to a different outcome.Nevertheless, land cover or the NDVI itself cannot be applied directly to extract prior BRDF knowledge from historical MODIS BRDF data in terms of the archetypal BRDF method.
Considerable effort has been devoted to extracting prior BRDF knowledge from historical BRDF data based on land cover or the NDVI [13,14].Studies have proved that linking the surface reflectance anisotropy to land cover or the NDVI may be useful for extracting prior BRDF knowledge from the historical BRDF product, which can improve the retrieval of land surface albedo.To extract accurate prior knowledge from the historical BRDF product, we must typically filter the existing BRDF data by setting limiting conditions, such as selecting homogenous pixels or pixels that change minimally over time.These additional conditions are difficult to achieve using the MODIS BRDF data at a resolution of 500 m.The improvement in this land cover-or NDVI-based prior BRDF knowledge for albedo retrieval can be explained by the constraints of the prior BRDF knowledge on BRDF model parameters.The retrieved BRDF is limited to a certain range in the entire viewing hemisphere, and the albedo, which can be calculated by integrating the BRDF over both the reflected and incident radiation hemispheres, can be limited to a normal extent.This study shows that an intermediate archetypal BRDF can also be used to correct the reflectance measurements for directional effects.However, when the reflectance anisotropy information in a specific direction must be used, such as the hot spots and dark spots [32,34,35], this prior BRDF knowledge may become unrepresentative.Further investigations are required to extract more accurate prior reflectance anisotropy knowledge from the historical BRDF product.

Summary and Conclusions
MODIS data and archetypal BRDFs were used to study the feasibility of extracting prior reflectance anisotropy knowledge from the BRDF product based on land cover or the NDVI.To account for the complicated anisotropic properties of surface reflectance, the BRDF dimensionality was reduced to six archetypal BRDFs.Their dependency on land cover or the NDVI was investigated by determining the percentage of each BRDF archetype class associated with a land cover type or a certain range of the NDVI.Additionally, this study compared the difference between the simulated NBARs and NDVIs retrieved from different archetypal BRDFs and the same multi-angular observations.The results illustrate only weak relationships between reflectance anisotropy and land cover or the NDVI, indicating that land surface reflectance anisotropy is not necessarily associated with land cover types or NDVIs at the resolution of the MODIS BRDF data, and that they are not reliable standards to accurately distinguish and extract surface reflectance anisotropy from the historical MODIS BRDF product.
A large amount of remote sensing data has been collected during the last two decades.Further work will focus on how to fully utilize these data, mainly concerning how to extract prior BRDF knowledge from historical BRDF data to improve albedo retrieval from insufficient multi-angular observations or even from a single directional observation.

Figure 1 .
Figure 1.Flowchart for calculating the proportion of each BRDF archetype class.

Figure 1 .
Figure 1.Flowchart for calculating the proportion of each BRDF archetype class.

Figure 2
shows the flowchart for the retrieval of the NBARs and NDVIs.The chart contains three parts: the recalculation of the MODIS NBARs from the MODIS BRDF product and RTLSR BRDF model, the retrieval of NBARs from archetypal BRDFs and MCD-obs, and the calculation of NDVIs from MODIS NBARs (red band) and different NBARs (NIR band) calculated from six archetypal BRDFs.Remote Sens. 2016, 8, 1004 5 of 16 Figure 2 shows the flowchart for the retrieval of the NBARs and NDVIs.The chart contains three parts: the recalculation of the MODIS NBARs from the MODIS BRDF product and RTLSR BRDF model, the retrieval of NBARs from archetypal BRDFs and MCD-obs, and the calculation of NDVIs from MODIS NBARs (red band) and different NBARs (NIR band) calculated from six archetypal BRDFs.

Figure 2 .
Figure 2. Flowchart of the retrieval and comparison of the NBARs and NDVIs.The present study first uses the MODIS BRDF/albedo product at specific points and corresponding true-color images of the ground to perform a visual assessment of the composition of surface reflectance anisotropy within different land covers.Then, the composition of the MODIS BRDF within several IGBP land cover types or every 0.1-NDVI interval is analyzed based on a fiveyear time series of MODIS data in three tiles.Finally, the six archetypal BRDFs are taken as prior BRDF knowledge in sequence to fit the MODIS multi-angular observations for tile h20v11 during days 201-216 in 2008.The NBARs and NDVIs that are calculated from different archetypal BRDFs and the same MCD-obs are compared to MODIS NBARs and NDVIs to analyze the correlation.In this study, only the anisotropic characteristics of the MODIS red and NIR bands are analyzed as an example.

Figure 2 .
Figure 2. Flowchart of the retrieval and comparison of the NBARs and NDVIs.

Figure 3 .
Figure 3. MODIS BRDF shape in the principal plane (a) for the deciduous broadleaf forest (b) and the one (c) for evergreen needle forest (d) images at two pixels in the NIR band.The solid lines refer to the BRDF shapes of the images on the left, and the dashed lines refer to the BRDF shapes of the images on the right.

Figure 3 .
Figure 3. MODIS BRDF shape in the principal plane (a) for the deciduous broadleaf forest (b) and the one (c) for evergreen needle forest (d) images at two pixels in the NIR band.The solid lines refer to the BRDF shapes of the images on the left, and the dashed lines refer to the BRDF shapes of the images on the right.

Figure 4 .
Figure 4.The land cover map of the study area (a) and the corresponding normalized BRDF shapes of the red (b) and NIR (c) bands in the principal plane for ENF (red lines), MF (green lines), grassland (blue lines), and cropland (cyan lines).

3. 1 . 3 .Figure 4 .
Figure 4.The land cover map of the study area (a) and the corresponding normalized BRDF shapes of the red (b) and NIR (c) bands in the principal plane for ENF (red lines), MF (green lines), grassland (blue lines), and cropland (cyan lines).

Figure 4 .
Figure 4.The land cover map of the study area (a) and the corresponding normalized BRDF shapes of the red (b) and NIR (c) bands in the principal plane for ENF (red lines), MF (green lines), grassland (blue lines), and cropland (cyan lines).
), BRDF archetype classes No. 4 and No. 5 are dominant (80%~90%), and during the summer, the percentage of BRDF archetype class No. 4 can reach up to 60%.By contrast, the percentages of BRDF archetype classes No. 1 and No. 6 are close to zero.During the winter, most of tile h11v03 is coved by snow, causing a lack of data (white areas in Figures

Figure 5 .
Figure 5. Cumulative histograms regarding the percentage of each BRDF archetype class within grassland (a); shrubland (b); cropland (c); forest (d); and desert (e) in the red band in three tiles from 2008 to 2012.Different grey levels refer to different BRDF archetype classes.The numbers of pixels used in the calculation are shown together (black lines).The white areas in the chart are due to the lack of data.

Figure 5 .
Figure 5. Cumulative histograms regarding the percentage of each BRDF archetype class within grassland (a); shrubland (b); cropland (c); forest (d); and desert (e) in the red band in three tiles from 2008 to 2012.Different grey levels refer to different BRDF archetype classes.The numbers of pixels used in the calculation are shown together (black lines).The white areas in the chart are due to the lack of data.

Figure 6 .
Figure 6.Cumulative histograms regarding the percentage of each BRDF archetype class within grassland (a); shrubland (b); cropland (c); forest (d); and desert (e) in the NIR band in three tiles from 2008 to 2012.

Figure 6 .
Figure 6.Cumulative histograms regarding the percentage of each BRDF archetype class within grassland (a); shrubland (b); cropland (c); forest (d); and desert (e) in the NIR band in three tiles from 2008 to 2012.

Figure 7 .
Figure 7.The average percentage of each BRDF archetype class for different NDVI ranges over five years.(a) and (b) show the results of tile h20v11 in the red and NIR bands.(c) and (d) show the results of tile h11v03 in the red and NIR bands.

Figure 7 .
Figure 7.The average percentage of each BRDF archetype class for different NDVI ranges over five years.(a,b) show the results of tile h20v11 in the red and NIR bands; (c,d) show the results of tile h11v03 in the red and NIR bands.

Figures 8 and 9
Figures 8 and 9 show a comparison of the MODIS NBARs and the NBARs retrieved from archetypal BRDF No. 1, 4, and 6.To emphasize the effect of the SZA on NBARs, the results of three SZAs (15 • , 40 • , and 60 • ) are shown together.In both bands, although the data used in this study includes various potential BRDF shapes (Table2), the NBARs retrieved from one specific archetypal BRDF can still exhibit good consistency with MODIS NBARs.The SZA has a large effect on the retrieval of NBARs.When the SZA is 15 • , the two NBARs are relatively different, and the coefficients of determination (R 2 ) are less than 0.87 and 0.9 in the red and NIR bands, respectively.However, as the SZA increases, the consistency between the two NBARs tends to increase, and the R 2 values are over 0.92 and 0.98.The large difference at a small SZA is mainly caused by the hot spot effect of the BRDF[32,33].Moreover, when the SZA is small, archetypal BRDF No. 1 (No. 6) tends to slightly overestimate (underestimate) the NBARs, compared to the MODIS NBARs; with increasing SZA, this difference tends to disappear.

Figures 8 and 9
Figures 8 and 9 show a comparison of the MODIS NBARs and the NBARs retrieved from archetypal BRDF Nos. 1, 4, and 6.To emphasize the effect of the SZA on NBARs, the results of three SZAs (15°, 40°, and 60°) are shown together.In both bands, although the data used in this study includes various potential BRDF shapes (Table2), the NBARs retrieved from one specific archetypal BRDF can still exhibit good consistency with MODIS NBARs.The SZA has a large effect on the retrieval of NBARs.When the SZA is 15°, the two NBARs are relatively different, and the coefficients of determination (R 2 ) are less than 0.87 and 0.9 in the red and NIR bands, respectively.However, as the SZA increases, the consistency between the two NBARs tends to increase, and the R 2 values are over 0.92 and 0.98.The large difference at a small SZA is mainly caused by the hot spot effect of the BRDF[32,33].Moreover, when the SZA is small, archetypal BRDF No. 1 (No. 6) tends to slightly overestimate (underestimate) the NBARs, compared to the MODIS NBARs; with increasing SZA, this difference tends to disappear.

Figure 8 .Figure 9 .
Figure 8.Comparison of the MODIS NBARs and NBARs retrieved from different archetypal BRDFs in the red band.The black points account for 95% of all selected data.(a-c) refer to the result at solar zenith angle of 15°, 40°, and 60°.Three columns represent the result of archetypal BRDF Nos. 1, 4, and 6.

Figure 8 .Figure 8 .
Figure 8.Comparison of the MODIS NBARs and NBARs retrieved from different archetypal BRDFs in the red band.The black points account for 95% of all selected data.(a-c) refer to the result at solar zenith angle of 15 • , 40 • , and 60 • .Three columns represent the result of archetypal BRDF No. 1, 4, and 6.(c) Figure 8.Comparison of the MODIS NBARs and NBARs retrieved from different archetypal BRDFs in the red band.The black points account for 95% of all selected data.(a-c) refer to the result at solar zenith angle of 15°, 40°, and 60°.Three columns represent the result of archetypal BRDF Nos. 1, 4, and 6.

Figure 9 .
Figure 9.Comparison of MODIS NBARs with NBARs retrieved from different archetypal BRDFs in the NIR band.(a-c) refer to the result at solar zenith angle of 15°, 40°, and 60°.Three columns represent the result of archetypal BRDF Nos. 1, 4, and 6.

Figure 9 .
Figure 9.Comparison of MODIS NBARs with NBARs retrieved from different archetypal BRDFs in the NIR band.(a-c) refer to the result at solar zenith angle of 15 • , 40 • , and 60 • .Three columns represent the result of archetypal BRDF No. 1, 4, and 6.

Table 2 Figure 10 .
Figure 10.Comparison of the MODIS NDVI and NDVIs retrieved from different archetypal BRDFs.(a-c) refer to the result at solar zenith angle of 15°, 40°, and 60°.Three columns represent the result of archetypal BRDF Nos. 1, 4, and 6.

Figure 10 .
Figure 10.Comparison of the MODIS NDVI and NDVIs retrieved from different archetypal BRDFs.(a-c) refer to the result at solar zenith angle of 15 • , 40 • , and 60 • .Three columns represent the result of archetypal BRDF No. 1, 4, and 6.

Table 1 .
The AFX and six BRDF archetypal parameters in the original (f * ) and normalized (F * ) forms for the red and NIR bands.

Table 2 .
The number of pixels and the percentages of BRDF archetype classes in the MODIS data.

Table 2 .
The number of pixels and the percentages of BRDF archetype classes in the MODIS data.