Reflectance Anisotropy from MODIS for Albedo Retrieval from a Single Directional Reflectance

Abstract

Surface reflectance anisotropy and insufficient multi-angular observations are the main challenges in albedo estimation from satellite observations. Numerous studies have been developed for albedo retrieval from a single directional reflectance by associating the anisotropy information extracted from coarse-resolution bidirectional-reflectance distribution function (BRDF) data. The contribution of land-cover type (LCT) and the Normalized Difference Vegetation Index (NDVI) in distinguishing reflectance anisotropy in these methods remains controversial. This study first proposed an approach to extracting a priori BRDF (F) from the MODIS BRDF/albedo product by considering the distribution characteristics of the model parameters. LCT- and NDVI-based F were also extracted from the corresponding subset. Then, the F-based albedo was derived from simulated or satellite directional reflectance and the anisotropic information of F. Finally, the directional reflectance and F-based albedo were compared with the MODIS albedo or ground measurement, in order to show the ability of F to compensate for the effect of reflectance anisotropy in the albedo retrieval process. The method was fully validated by the global and time-series MODIS BRDF data. The results showed that reflectance anisotropy has an aggregated distribution pattern, and F can represent the reflectance anisotropy of most pixels within a tile. The improvement of LCT and NDVI only occurs when the tile contains a large area of vegetated and barren ground. With the exception of the hotspot and large viewing-zenith-angle area in the forward hemisphere, the F-based shortwave albedo has high consistency with the MODIS albedo product. A comparison with the ground measurements and MODIS albedo showed that the F-based albedo from a single directional reflectance generally achieves an absolute accuracy requirement, with a root-mean-square error (RMSE) of 0.027 and 0.036.

1. Introduction

Land-surface albedo, which is a critical geophysical variable controlling Earth’s energy budget, must be determined with acceptable accuracy [1]. Satellite remote sensing provides a practical way to perform surface albedo mapping and monitoring [2]. Over the past two decades, the global surface albedo at various spatial and temporal resolutions has been derived from the data from various satellite sensors [2,3,4,5,6,7]. Satellite-derived surface albedos have been widely used as inputs to various climate models [8,9], and the absolute accuracy of the albedo required for model studies has been estimated at ±0.02–0.05 [3,10]. However, existing satellite products suffer from several problems, including the lack of a long-term, fine-spatial-resolution albedo. This is mainly due to a lack of information on surface reflectance anisotropy, the scanning limit of fine-resolution sensors, and contamination by atmospheric aerosols and clouds [2,11,12].

For most natural surfaces, the reflection of incoming solar radiation is anisotropic [13]. The effect of surface anisotropy on remotely sensed observations has been investigated by many researchers over the past two decades [2,14,15]. It is widely accepted that large uncertainties may be introduced to an albedo if the effect of reflectance anisotropy is not appropriately considered in the retrieval process [16]. Usually, the bidirectional-reflectance distribution function (BRDF) is used to describe surface reflectance anisotropy [13]. Various empirical and semi-empirical BRDF models have been developed to convert directional reflectance to albedo. Most BRDF models can accurately perform the BRDF/albedo retrieval inversion given a sufficient number of observations (i.e., >7) that are well distributed in the viewing angle [12]. However, problems will occur when angular samplings are sparse or are not well-distributed [17,18].

The capacity of sensors to acquire multi-angular measurements is most probably limited by sensor scanning configurations and the platform’s orbital characteristics [11]. Since coarse-resolution satellite sensors—such as the Moderate Resolution Imaging Spectroradiometer (MODIS) [2] and the Polarization and Directionality of the Earth’s Reflectances (POLDER) [19]—have wide swaths (approx. 2400 km for both sensors), large maximum satellite zenith angles, and short periods of orbital repeat cycle, they are able to observe a location from a variety of angles in a relatively short period of time [2]. However, most middle-spatial sensors (such as the Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), Operational Land Imager (OLI) [20,21,22], and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [23]) and high-spatial sensors (such as the Satellite Pour l’Observation de la Terre (SPOT) [24] and GaoFen (GF) [25,26]) provide a narrow swath-width (e.g., Landsat: 185 km and GF-2: 45 km). As such, orbits from consecutive days do not overlap, making them unable to obtain multi-angular observations. This makes it difficult to derive BRDF directly from such middle- and high-resolution sensors.

Various BRDF/albedo products have been derived from existing coarse-resolution satellite data [2,19]. The coarse-resolution BRDF/albedo product has been used to normalize directional reflectance into a common viewing and illumination geometry [27,28,29], and to improve albedo retrieval from fine-resolution reflectance [20,21,22,30]. The influence of spatial resolution on the angular variation patterns of optical reflectance retrieved from various pixel-scale measurements has been typically analyzed [31,32]. There are two main inversion strategies for albedo retrieval from directional reflectance. One is based on directional reflectance at the Top of Atmosphere (TOA) by utilizing linear functions [21,33,34], and the other strategy is based on surface directional reflectance and BRDF knowledge estimated from existing BRDF products [20,22]. In the present study, we focus on the latter and mainly explore how to extract proper BRDF knowledge from a historical MODIS BRDF product for albedo retrieval from a single directional reflectance.

A priori information is widely used in albedo retrievals from poor angular sampling [35,36,37,38]. Early studies based on field measurements revealed that reflectance anisotropy was associated with different homogeneous LCTs [14,39]. Strugnell and Lucht pointed out that the BRDF shape may have a strong relationship with the LCT, and a back-up algorithm was developed for MODIS BRDF products when the observations could not meet the requirement of the full inversion [40]. Varying LCT-based or NDVI-based BRDF knowledge, which was extracted from MODIS BRDF products, has been applied to Landsat or Sentinel surface reflectance to generate albedo, and a comparison with field measurements showed that the retrieved fine-resolution albedo generally satisfied the precision requirement [20,22,41]. However, a different perspective was presented based on POLDER data, stating that the variation in directional signatures is more significant in intra-class than inter-class types [42]. MODIS BRDF and the Anisotropic Flat Index (AFX) theory suggest that reflectance anisotropy does not necessarily have a dependency on LCT and NDVI; therefore, global MODIS BRDF archetypal maps have been generated [43]. It is important to understand the cause of these potentially conflicting conclusions for a deeper understanding of reflectance anisotropy and its effect on albedo retrieval.

Recently, we noticed that the distribution pattern of the model parameters is clustered. A study based on AFX theory and simulated MODIS reflectance revealed that the intermediate archetypal BRDF can improve the accuracy of albedo retrieved from insufficient observations [44]. Combining the two results, we are curious as to whether a priori BRDF extracted from MODIS BRDF can represent the dominant reflectance anisotropy of the underlying surface to convert single directional reflectance to albedo. In the present study, we first explore the distribution pattern of reflectance anisotropy with four MODIS tiles based on the scatter diagram of normalized BRDF model parameters. A priori BRDF, which is expected to represent the reflectance anisotropy of the most pixels, is extracted from a high-quality MODIS BRDF product or from subsets separated by LCT or NDVI in a tile. Then, a priori BRDF is applied to a simulated single directional reflectance to retrieve albedo using the magnitude-inversion method. The accuracy of albedo is evaluated using the high-quality MODIS albedo product. The effect of illumination and observation geometry and the ability of LCT and NDVI to distinguish reflectance anisotropy for albedo retrieval are also discussed. Finally, a fine-spatial-resolution albedo is also generated using a priori BRDF and Landsat reflectance, and the accuracy of the Landsat albedo is validated using ground measurements.

2. Materials and Methods

2.1. Data Used and Study Area

2.1.1. MODIS Data

The MODIS BRDF/albedo product has a temporal and spatial resolution of daily and 500 m, respectively. It can satisfy the accuracy requirements of model studies [12,45,46], and has been widely used in various reflectance anisotropy studies [20,22,43]. To guarantee the precision of the result in prior BRDF extraction and directional reflectance simulation, only the snow-free high-quality MODIS BRDF/albedo product [12] was engaged in this study.

This paper focuses on the analysis of four MODIS tiles (h12v04, h30v12, h20v05, and h18v07), which have different surface structures, to exhibit the distribution characteristics of reflectance anisotropy. The International Geosphere-Biosphere Programme (IGBP) land-cover map in 2010 and the NDVI map on DOY 170 in the four tiles are shown in Figure 1a–h, respectively. On DOY 170, the path of the sun in the sky is farthest north in the Northern Hemisphere. The solar zenith angle at local solar noon is about 22°, 56°, 11°, and 8° around the center of the four tiles. The land cover data were obtained from MCD12Q1, and the NDVI was recalculated from the simulated nadir reflectance based on the high-quality MODIS BRDF product and RTLSR BRDF model at local solar noon. H12v04 is in the northeast of the United States, and h30v12 is in the southeast of South Australia. Both tiles contain various kinds of forests, shrublands, savannas, and grasslands. The large value of the NDVI in h12v04 suggests that it has denser vegetation than h30v12. H20v05 is in the eastern Mediterranean Sea. The northern region of the tile is mainly covered by forests, savannas, and grasslands, while the southern region is covered by barren land. H18v07 is in North Central Africa; similarly, the northern region is covered by barren land, while the southern region is mainly covered by grasslands. Due to the low-quality of the reflectance in the south region, the main land cover of the tile is barren land at this time. Beyond that, the time-series data of MODIS BRDF from 2008 to 2010 over the four tiles, and the global MODIS BRDF data on DOY 170, 2010, were used in the analysis of the temporal and spatial performance of a priori BRDF for albedo retrieval from directional reflectance. The MODIS BRDF data in h12v04, h10v05, h08v05, h09v04, and h11v04 during 2000–2010, which are temporally and spatially consistent with seven surface-radiation budget network (SURFRAD) stations, were also used to extract a priori BRDF for the retrieval of albedo from Landsat reflectance.

To further verify the potential of using a priori BRDF for albedo retrieval from a single directional reflectance, multi-angular observations from MODIS (MOD09GA) during a 16-day composite period (2010.162-177) in h20v04 were also used in this study. Only the first layer of the product was of high quality during the 16-day period and was extracted. Figure 2 shows an example of the MODIS angular samplings during the period. The solar zenith angle ranges from 26° to 32°, the largest viewing zenith angle is around 60°, and the relative azimuth angle is around 45° (240°) in the backward (forward) hemisphere.

2.1.2. Landsat Satellite Data and Field-Measured Albedo

The United States Geological Survey (USGS) provides calibrated and atmospherically corrected Landsat surface reflectance data via Earth Explorer (https://earthexplorer.usgs.gov/, accessed on 1 September 2020.). In this study, all available TM and ETM+ land-surface reflectance data from 2001 to 2010 corresponding to the seven ground sites were downloaded (

Table 1. Study areas represented by MODIS tile number and Landsat path/row, as well as the site of SURFARD network.

MODIS Tile No.	Lat./Lon.	Dominant LC	Landsat		SURFRAD Station
			Path/Row	No.
h12v04	40–50N/93–65W	Forest, Mixed Forests, Savannas, and Croplands	16/32	62	PSU
h30v12	30–40S/138–169E	Forest, Shrublands, Savannas, and Grasslands	-	-	-
h20v05	30–40N/26–35E	Grasslands, Croplands, and Barren	-	-	-
h18v07	10–20N/0–10E	Grasslands, Shrublands, and Barren	-	-	-
h10v05	30–40N/104–80W	Forest, Savannas, Croplands, and Mosaics	23/36	235	GWN
h08v05	30–40S/104–130E	Barren or Sparsely Vegetated	40/35	277	DRA
h09v04	40–50S/105–140E	Grasslands and Evergreen Needleleaf Forests	34/32	244	TBL
h11v04	40–50S/78–109E	Croplands and Deciduous Broadleaf Forests	23/32	195	BON
			35/26	296	FPK
			29/30	114	SXF

). The Landsat 5 TM and 7 ETM+ nominal radiometric calibration accuracy is reported as 5% [47,48]. To ensure the precision of the retrieved albedo, any data gaps in ETM+ caused by the failure of the scan-line corrector (SLC) were omitted in the albedo retrieval process.

SURFRAD has been providing accurate field observation of broadband albedo since 1995 [22]. The dominant LCTs for these sites are grassland and agricultural. The measurements of these sites are spatially representative of the conditions of the surrounding area, with an uncertainty range from ±2% to ±5% [49]. The ground-measured albedo at the seven SURFRAD sites is used to evaluate the accuracy of the retrieved albedo. The 10 m high pyranometer has a field of view of 81°, and covers a radius of ~63 m on the ground, which is equal to a 3 by 3 Landsat 30 m pixels [21,22]. The ground-measured broadband albedo within 30 min before and after the satellite overpass time was averaged to represent the ground albedo [21,22]. The shortwave (0.25–2.5 um) broadband reflectance, which can be taken as the albedo under Lambertian assumption, is calculated from the spectral reflectance via narrow- to broadband conversion coefficients [50]. It is also directly compared with ground measurement to show the improvement in F for albedo retrieval.

2.2. A Priori BRDF (F) Extracted from MODIS BRDF

The RTLSR BRDF model has the following form [2,5,51]:

$$R\left({\theta }_s,{\theta }_v,\Phi ,\lambda \right)=f_{iso}\left(\lambda \right)+f_{vol}\left(\lambda \right)K_{vol}\left({\theta }_s,{\theta }_v,\Phi ,\lambda \right)+f_{geo}\left(\lambda \right)K_{geo}\left({\theta }_s,{\theta }_v,\Phi ,\lambda \right)$$

(1)

where R is the surface bidirectional reflectance in the wavelength λ; f_iso, f_vol, and f_geo are the spectrally dependent BRDF model parameters; K_vol and K_geo are kernel functions of the volumetric [52,53] and geometric optical [15,51] scattering, respectively; θ_s and θ_v denote the illuminating and viewing zenith angles; and Φ is the relative azimuth angle.

The characteristics of reflectance anisotropy are complex, and the BRDF model parameters provide a feasible method for the quantitative evaluation of reflectance anisotropy. The difference in spectral reflectance amplitude causes the BRDF shapes to be mixed together and makes them difficult to distinguish one from the other. To remove the effect of spectral reflectance amplitude and make the surface reflectance anisotropy comparable, the BRDF parameters are simply normalized by a scale factor α/f_iso [43], where α is set to 0.5 to limit the range of normalized BRDF shapes to 0–1.0 when drawing them. The surface reflectance can be rewritten as:

$$R\left({\theta }_s,{\theta }_v,\Phi ,\lambda \right)=({F^{\prime }}_{iso}\left(\lambda \right)+{F^{\prime }}_{vol}\left(\lambda \right)K_{vol}\left({\theta }_s,{\theta }_v,\Phi ,\lambda \right)+{F^{\prime }}_{geo}\left(\lambda \right)K_{geo}\left({\theta }_s,{\theta }_v,\Phi ,\lambda \right)/\left(\frac{\alpha }{f_{iso}}\right)$$

(2)

where F′ is the normalized BRDF model parameters with F′_iso = α, F′_vol = αf_vol/f_iso, and F′_geo = αf_geo/f_iso.

Figure 3 shows the density distribution of the original (f) and normalized (F′) volumetric and geometric optical scattering parameters of a tile. The scatter diagram of the two model parameters is divided into S_max × L_max grids with a size of k. S_max and L_max refer to the sampling number of columns and rows. Values of 0.005, 260, and 60 were set for k, S_max and L_max. The distribution pattern of both the original and normalized model parameters is aggregative. Figure 3a also shows the distribution of the normalized model parameters (F’) in four grids (A–D). The pixels that are located in the same grid in the original mode are separated into different grids in the normalized mode. The shapes of the original BRDFs in a grid are almost parallel with each other due to the differences in f_iso, while the shapes of the normalized BRDFs in a grid are nearly the same since the F′_iso is a constant. BRDF normalization provides a more feasible way to define the reflectance anisotropy in a grid. Furthermore, in the albedo retrieval process (Section 2.3), a multiplicative factor is applied to generate a proper BRDF from a priori BRDF (Equations (6)–(9)), and the BRDF normalization process does not affect the albedo retrieval. Therefore, a priori BRDF is extracted from the normalized BRDFs rather than the original BRDFs.

The integration processes of BRDF can weaken the effect of reflectance anisotropy on the albedo, and an intermediate archetypal BRDF can be used to improve albedo retrieval from insufficient observations [44]. The model parameters within one MODIS tile shows an aggregated distribution (Figure 3). We are curious as to whether a priori BRDF extracted from MODIS BRDF can represent the dominant reflectance anisotropy of the underlying surface for converting single directional reflectance to albedo. A priori BRDF (F) is extracted from the MODIS BRDF product within a MODIS tile based on the probability-weighted method. For column i and row j in the grid, the number of pixels n(i, j) in the grid can be counted. The central value F’(i, j) is taken as the reflectance anisotropy of the grid (Equation (3)). When the grid contains fewer than 10 pixels (gray color in Figure 3), it will not be considered in the process of extracting F. This strategy can filter out some noise pixels that have exception parameters. The probability of each unit can be calculated using the ratio of n(i, j) to the total number of pixels (n) that participated in the operation. The F value can be calculated through the sum of the products of F’(i, j) and n(i, j)/n (Equation (4)). The F value, which is expected to represent the dominant reflectance anisotropy, is taken as the a priori BRDF in this study.

$$F^{\prime }\left(i,j\right)=\left(k\ast i, k\ast j\right)-\frac{k}{2}$$

(3)

$$F={\displaystyle \sum }_{i=1}^{s_{max}}{\displaystyle \sum }_{j=1}^{L_{max}}{F}^{\prime }\left(i,j\right)\times n\left(i,j\right)/n$$

(4)

To learn the ability of LCT and NDVI in order to distinguish reflectance anisotropy, the high-quality MODIS BRDF product is further segmented by LCT or NDVI, and F is extracted from each specific BRDF subset. F, F_L and F_N are used to represent the a priori BRDF extracted from all MODIS BRDF products, and the LCT- and NDVI-based high-quality ones.

2.3. Albedo Inversion Method and Evaluation

2.3.1. F-Based Albedo from a Single Directional Reflectance

The method used to adjust F to fit single directional reflectance follows the MODIS magnitude-inversion algorithm [12]. For a given single directional observation (ρ), the simulated directional reflectance ρ_s, which has the same viewing geometry as ρ, can be simulated according to RTLSR BRDF model and F (Equation (1)). Suppose that the ratio of directional reflectance and reflectance anisotropy and the albedo should be similar [18,40]; the strength of the BRDF effect could be determined by the ratio (a) between ρ_s and ρ at given viewing and solar geometries (Equation (5)).

$$a=\rho /{\rho }_s$$

(5)

Using the RTLSR model, the Black Sky Albedo (BSA) at any solar zenith angle and White Sky Albedo (WSA) of the pixel under consideration can be calculated (Equation (6)). The BSA and WSA are, then, combined as a function of the fraction of diffuse skylight (S) to calculate the blue-sky albedo (α) (Equation (7)) [54]. The albedo refers to blue-sky albedo in the subsequent section. In this study, an aerosol optical depth (AOD) of 0.2 is used to calculate S when learning the anisotropic reflectance characteristics based on the simulated MODIS directional reflectance. When calculating the Landsat albedo at the SURFRAD stations, S is measured using different pyranometers:

$$BSA\left(\theta \right)=a{{\displaystyle \sum }}_{i=1}^3{F}_i{h}_k\left(\theta \right);WSA=a{\sum }_{i=1}^3{F}_i{H}_k$$

(6)

$$\mathsf{\alpha }\left(\theta \right)=\mathrm{S}\left(\theta ,{\tau }_{550\mathrm{nm}}\right)BSA\left(\theta \right)+\left(1-\mathrm{S}\left(\theta ,{\tau }_{550\mathrm{nm}}\right)\right)\mathrm{WSA}$$

(7)

where h_k and H_k are the directional-hemispherical and bi-hemispherical integrals of the BRDF model kernels, do not depend on the observations, and can be precomputed and stored in a look-up table; and τ_550nm is the AOD at 550 nm.

2.3.2. The Feasibility Analysis Based on Simulated MODIS Directional Reflectance

The feasibility of using F for albedo retrieval from a single directional reflectance is explored. The RTLSR model in the forward mode (Equation (1)) is used to simulate shortwave single directional reflectance (ρ) for each pixel based on the high-quality MODIS BRDF product. The simulated MODIS blue-sky albedo (A_F) can be calculated from ρ and F using Equations (6)–(8). The relative difference (Rd) and absolute difference (Ad) between A_F and the high-quality MODIS albedo product (A_m) can be calculated for each pixel using Equations (8) and (9). The percentage of pixels that have an absolute difference range from −0.02 to 0.02 (P_0.02) is used to quantitatively evaluate the performance of F (Equation (10)). The larger the P_0.02, the more representative it is of F for albedo retrieval from a single directional reflectance.

$$\mathrm{Rd}=\frac{A_F-A_m}{A_m}\times 100\%$$

(8)

$$Ad=A_F-A_m$$

(9)

$$P_{0.02}=\frac{n(\left|Ad\right|<0.02)}{n}$$

(10)

To study the impact of viewing geometry, first, the image of MODIS shortwave-band reflectance in nadir—at a forward and backward angle of 30° with a solar zenith angle at local solar noon, and F-based albedos on DOY 170, 2010 in h12v04—are used to provide a visual evaluation. Second—based on the four MODIS tiles on DOY 170, 2010—nadir reflectance with varying solar zenith angles, and the reflectance in any direction in the entire viewing hemisphere with a solar zenith angle at local solar noon are simulated; then, the albedo is retrieved from a directional reflectance by considering F, F_L, and F_N, respectively. Finally, the performance of F is fully validated by the spatial and temporal MODIS BRDF product. From a spatial perspective, the global albedo on DOY 170, 2010 is retrieved from F and nadir reflectance. From a temporal perspective, the time-series albedo over the four tiles from 2008 to 2010 is retrieved from F, and any directional reflectance in the principal plane (PP) with a solar zenith angle of the local solar noon. The retrieved albedo is compared with the MODIS albedo product, and the P_0.02 in each direction is used to evaluate the effect of viewing geometry and the improvements in LCT and NDVI for albedo retrieval.

2.3.3. Albedo Retrieval from F and Real Satellite Observation

We explore the feasibility of using the MODIS-based F to retrieve albedo from MODIS and Landsat observation. The Fs, which are spatially and temporally consistent with the observation, are extracted from the high-quality MODIS BRDF product based on the probability-weighted method (Section 2.2). To avoid the uncertainty caused by the BRDF model and the atmosphere [45,55,56], the observation with a solar zenith angle greater than 65°, or with a ratio of diffuse skylight measured at ground sites greater than 0.3 are omitted. The F-based Landsat and MODIS albedos are retrieved from 1423 Landsat scenes and MODIS observations during 162–177, 2010 in h20v04.

3. Results

3.1. A Priori BRDF (F) from MODIS

The MODIS BRDF product over H12V04, h20v05, h18v07, and h30v12 on DOY 170, 2010 is used as an example to exhibit the distribution characteristics of reflectance anisotropy. In the case of the four tiles, the number of high-quality pixels is more than one million in the spectral and shortwave broadband, and the number of data can serve as sufficient material for the analysis of the distribution of reflectance anisotropy. Figure 4 shows the scatter diagrams of F′_vol and F′_geo in the four MODIS tiles. To keep a consistent style, these figures use an X-axis range from 0 to 1 and a Y-axis range from 0 to 0.3.

Although each band contains varying reflectance anisotropy, the distribution pattern of reflectance anisotropy presents an aggregated distribution. The range and gathering center of reflectance anisotropy differ between tiles, which may be caused by ground-surface structures. When the surface is dominated by dense vegetation (e.g., h12v04) and the surface structure is supposed to be complex, the distribution of reflectance anisotropy is more diffuse. When the surface is dominated by desert and the surface structure is simple (e.g., h18v07), the distribution range of reflectance anisotropy is reduced remarkably. The most dense area in the scatter diagrams is close to zero, especially for the geometric optical scattering parameter. When the surface contains a large proportion of both vegetated and barren land (e.g., tile h20v05), the scatter diagrams may have multiple gathering centers. When the surface is covered by vegetation with various densities (e.g., h30v12), the percentage of pixels that have a medium density increases significantly.

Figure 5 shows the scatter diagrams of F′_vol and F′_geo for different subsets. Because of the small number of LCTs and relatively stable range of NDVI in h18v07, the tile is not considered in this part. The pentacle represents the location of the F_L or F_N. Although the number of available pixels for the subset is reduced, these scatter diagrams of reflectance anisotropy in different subsets show an aggregated distribution pattern. In h12v04, the distribution pattern of reflectance anisotropy may be influenced by the number of pixels; however, the gathering center remains consistent for different subsets. In h20v05, the vegetated LCTs have a similar gathering center and the ENF has a relatively strong F′_vol. The non-vegetated barren area, which usually has NDVI values less than 0.2, is more Lambertian, as the value of F′_vol and F′_geo is small. When the NDVI ranges from 0.2 to 0.7, the gathering center remains consistent. In h30v12, the gathering center of EBF and woody savannas is located in a F′_vol-dominated area, while open shrubland has a weak F′_vol. The other three LCTs have a similar distribution pattern. With an increase in NDVI, the gathering center gradually moves towards the F′_vol dominated area. The ability of LCT and the NDVI to distinguish reflectance anisotropy is limited, and a visible difference only appears when there is a surface change from non-vegetated (Barren or NDVI < 0.2) to vegetated (NDVI > 0.2).

Figure 6 shows the values of F, F_L, and F_N that were extracted from the high-quality MODIS BRDF dataset or from the subset divided into LCT or NDVI in the three tiles. The values of F_L and F_N are extracted when the number of pixels in the subset is greater than 20,000. The variation regularity of F_L and F_N is not obvious. The woody, shrub, and herbaceous plants, which have obvious differences in structure, could even have a similar F_L values. In the three tiles, the forest (ENF, MF, EBF) usually has a larger F_vol, while the barren and open shrubland usually have the smallest F_vol. The F_vol values in the near-infrared and shortwave broadband have a decreasing trend with the increase in NDVI, while the other bands have an increasing trend in h12v04. In h20v05 and h30v12, F_vol has an increasing trend when the NDVI is less than 0.3, but the momentum appears to wane when the NDVI is greater than 0.3. This may be because the surface structure of the two tiles is not as complex as in h12v04. The F_geo is always small in the three tiles, and its change with NDVI and LCT is not significant.

3.2. Albedo Retrieval from F and Nadir Reflectance

3.2.1. Impact of the BRDF on Albedo Retrieval from Nadir

The effect of reflectance anisotropy on albedo retrieval is first investigated based on F and nadir reflectance. Figure 7c shows the BRDF shape of F extracted from h12v04 on DOY 170, 2010. As a comparison, geometric-optical-scattering-dominated (F_g) and volumetric-scattering-dominated (F_v) BRDF shapes are also shown together. The F_g and F_v values are determined according to the distribution of the model parameter (Figure 4a). F_g has a relatively large value of F_geo, while F_v has a relatively large value of F_vol, and the other parameter is set to 0. The 3D BRDF shapes are plotted using a visualization tool provided by Dong [57]. The F_g shows a dome shape, while the F_v shows a bowl shape. The shape of F is in-between.

Nadir reflectance is simulated according to the RTLSR BRDF model and high-quality MODIS BRDF product in h12v04 on DOY 170, 2010. The F-based albedo (A_F) is retrieved from the nadir reflectance and F. The Ad and Rd between the A_p and MODIS albedo product (A_m) are calculated for each pixel. The distribution of Rd in the scatter diagrams and the accumulative distribution of Ad for the shortwave broadband are shown in Figure 8. The results show that when the pixels have similar reflectance anisotropy to F, the relative difference is close to 0. The relative difference gradually increases with increasing distance between F and the real reflectance anisotropy. In conjunction with the distribution of reflectance anisotropy (Figure 4), F is located near the high-density region of reflectance anisotropy, and the retrieved albedo of most pixels has high precision and reliability. F_g and F_v, which are away from the high-density region, may introduce large errors in the albedo for most of the pixels.

The percentage of pixels that have an Rd between −10% and 10% is 13.4%, 7.5% and 61.3% for F_g, F_v, and F, respectively. For the albedo based on F, the percentage of pixels that have an Ad within −0.02 and 0.02 (P_0.02) is 94% (Figure 8d), while the P_0.02 is 63% and 28% for F_g and F_v, respectively. The P_0.02 is 81%, 78%, and 69% for h30v12, h18v07, and h20v05, respectively, in shortwave broadband. The small value of P_0.02 may be caused by the small solar zenith angle (h18v07) and multiple gathering centers of reflectance anisotropy (h20v05). The details of the result will be discussed in the subsequent sections.

3.2.2. Impact of the Solar Zenith Angle on Albedo Retrieval from Nadir

In order to explore the effect of the solar zenith angle on albedo retrieval, F, F_L, and F_N in the four tiles are fitted to the corresponding simulated MODIS nadir reflectance, with different solar zenith angles, to generate the albedo. As a comparison, the nadir reflectance is also compared with the MODIS albedo. Figure 9 shows the change in the P_0.02 with the solar zenith angle. The nadir reflectance greatly differs from the albedo, and the retrieved albedo is improved by employing the constraint of F, F_L, or F_N, especially when the solar zenith angle is small or large. This phenomenon can be explained by the weak fitting ability of the RTLSR BRDF model around the hotspot area [58] and large solar zenith angle condition [45]. When the solar zenith angle is around 40°, both the nadir reflectance and the F-, F_L-, or F_N-based albedos are close to the MODIS albedo, with a P_0.02 greater than 90%. In h20v05, the F_L or F_N can increase the value of P_0.02 by 10–20 percent for the small or large solar zenith angle conditions. This could be caused by the large proportion of both vegetated and barren land in h20v05. In h30v12, when the solar zenith angle is less than 40°, the performance of F, F_L, and F_N are similar. The difference appears when the solar zenith is greater than 45°, with an increase of 10 percent in P_0.02 for F_L and F_N. In the other two tiles, the F_L and F_N have similar performance when compared with F.

3.2.3. Spatial Validation with Nadir Reflectance

The spatial variations in reflectance anisotropy are analyzed based on global MODIS shortwave BRDF data on DOY 170, 2010. Figure 10a,b show the model parameters of F. For this period, the northern hemisphere and north-central part of South America usually has a relatively strong F_vol, and the maximum value appears in h20v03, h20v04, and h21v03; this is covered by dense crop, mixed forests, and grass. The relatively strong F_geo appears around the Americas and the South Africa, and the maximum value appears on the west coast of the Americas, which are mainly covered by shrublands and savannas. For the desert regions in north Africa, the values of both F_vol and F_geo are small.

Like the results shown in Figure 9, Figure 10c,d show that the nadir reflectance in the tiles that have a small solar zenith angle greatly differs from the MODIS albedo caused by the hotspot effect. When the solar zenith angle is around 40°, the nadir reflectance has the best consistency with the MODIS albedo. For most of the tiles, F has an appositive impact on albedo retrieval and increases the value of P_0.02; notably, the number of tiles that have a P_0.02 greater than 80% changes from 56 to 116 compared with the Lambertian assumption result. It also shows that the improvement is more significant for the tiles that have a small or large solar zenith. A small P_0.02 appears in h21v06, h22v06, h22v07, h23v07, and h18v05, which may be caused by a large percentage of both vegetation and barren surface, and a small solar zenith angle. The small P_0.02 in h15v01 may be caused by snow and ice since the NDVI of most pixels is negative and the mean f_iso is greater than 0.5.

Figure 10e shows that for most of the tiles, the difference between F_L- and F-based P_0.02 is between −1% and 3%, and only six tiles have a difference greater than 10%; this indicates that the improvement of F_L is limited in most cases. The improvement appears to be more significant when the solar zenith angle is small, and the largest improvement (20%–25%) appears in h18v05 and h22v06, which contain a large proportion of both vegetation and barren land. There is no obvious distinction between the effect of F_N and that of F_L for most of the tiles, which have a difference range from −3% to 3% (Figure 10f).

3.2.4. Comparison of Ground Albedo with F-Based Landsat Albedo

The F extracted from a MODIS tile and the Landsat reflectance are combined to generate 30 m-resolution albedo. The F is spatially and temporally consistent with the Land-sat reflectance around the SURFRAD sites during 2001–2010. A comparison with the ground measurements shows that F performs well in albedo retrieval from the TM and ETM+ data (Figure 11). The validation of Landsat shortwave albedo against ground measurements showed a bias of 0.04 and an RMSE of 0.027 for the seven SURFRAD sites. Specifically, Landsat albedo shows a slight underestimation at BON, GWN, and PSU, with a bias of −0.016, −0.017, and −0.023. At the Desert Rock (DRA) and TBL sites, surface albedo has been consistently overestimated compared with ground measurements, with bias ranging from 0.03 to 0.034. The results confirm the findings reported at the SURFRAD sites in previous studies [20,21,22,30].

The Landsat shortwave surface reflectance (Lambertian albedo) is also compared with ground measurements. In general, the surface Lambertian assumption may introduce great uncertainty into the albedo, with RMSEs ranging from 0.022 to 0.042. By applying the concurrent F to the reflectance, the RMSE will be obviously decreased, except for that of the TBL station. For all the data, the RMSE can be reduced from 0.032 to 0.027. This suggests that the F extracted from the MODIS BRDF product can improve the accuracy of the 30 m-resolution albedo retrieved from a single directional reflectance. When the solar zenith angle is small or large, unlike the result based on simulated reflectance (Figure 9), the difference between the reflectance and measured albedo does show an obvious increasing trend. Meanwhile, at the Desert Rock (DRA) and TBL sites, both the Landsat reflectance and albedo have consistently underestimated surface albedo. Both problems may be caused by the uncertainty in the reflectance, since the proposed method to estimate the albedo from single directional reflectance is sensitive to noise, and the uncertainty in the reflectance will be completely transferred to the albedo.

3.3. F-Based Albedo from Any Single Directional Reflectance

3.3.1. The Potential of Using F for Albedo Retrieval from Any Single Directional Reflectance

The reflectance in any direction for each pixel is simulated based on the MODIS BRDF product and RTLSR BRDF model. The F, F_L, and F_N and the directional reflectance are then used to retrieve the albedo. The P_0.02 in each direction is shown in Figure 12 for the four tiles. For h12v04 and h18v07, except for the directional reflectance around the hotspot area and the large viewing-zenith-angle area in the forward directions, P_0.02 is always greater than 90%. For h20v05 and h30v12, the directions that have low P_0.02 increase obviously; these are mainly concentrated around the hotspot and the large viewing-zenith-angle area and in the forward hemisphere. This could be explained by the multi-gathering center of reflectance anisotropy in h20v05, and the relatively large solar zenith angle in h30v12. The performance of F_L and F_N for albedo retrieval from directional reflectance is similar. There is little or no improvement in h12v04 and h18v07 since F has good performance. For h20v05, visible improvement appears around the hotspot and large viewing-zenith area. For h30v12, the improvement mainly occurs in the forward hemisphere, but the accuracy is relatively low.

3.3.2. Temporal Validation with Reflectance in PP

The three-year time series (2008 to 2010) of F, F_L, and F_N in the shortwave band are extracted from four MODIS tiles. The directional reflectance in PP for each pixel is simulated based on the MODIS BRDF product and the RTLSR BRDF model. The surface albedos are retrieved based on F, F_L, and F_N and the simulated single directional reflectance in PP. The P_0.02 in each direction is also used to evaluate the accuracy of the albedo. The model parameters of the time series of F (right y-axis), the P_0.02 in PP, and the solar zenith angle for the four tiles are shown in Figure 13. Some of the F that has low accuracy is plotted with a red circle, because the number of high-quality pixels in the tile is insufficient (i.e., <10⁵). The results show that the inter-annual variability and trend of reflectance anisotropy are similar, which helps to reveal the regularity of the surface anisotropy variation and proves that the probability-weighted method for extracting prior BRDF captures the dominant surface reflectance anisotropy. With the underlying surface change from forest (h12v04) to savanna (h20v05) to semi-desert (h18v07), the degree of variation in F_vol and F_geo gradually decreases. There are some significant changes in F_vol, especially in the forest areas, while the change in F_geo is not obvious.

Like the results shown in Figure 12, for h12v04 and h18v07, except for the hotspot and large viewing-zenith-angle area, more than 90% of pixels have a similar albedo to MODIS. A notable feature is that during the winter time, in h20v04, F-based albedo has a low accuracy; this is possibly caused by the effect of snow. For h20v05, the accuracy of the F-based albedo decreases obviously; this is because both the vegetation and barren ground account for a relatively large proportion within the tile. For the four tiles, when the solar zenith angle is greater than 50°, the performance of F starts to show a decline in the forward direction. Apart from these, despite only a priori BRDF being used for one MODIS tile, the F-based albedo in most directions has a high consistency with the MODIS albedo over the three years. Figure 13e,f shows the P_0.02 in PP based on the F_L and F_N. The directions that have high accuracy are increased, especially in h20v05. The accuracy around the forward direction is obviously improved when the solar zenith angle is greater than 50°. There is no significant difference between the F_L- and F_N-based albedo retrieved from directional reflectance.

3.3.3. F-Based Albedo from Simulated and Real MODIS Directional Reflectance

Figure 14a–c show three simulated shortwave reflectance maps based on the MODIS BRDF product of tile h12v04 and the RTLSR BRDF model with a solar zenith angle of 22°. Compared with the MODIS shortwave albedo product (Figure 12d), the reflectance in nadir and at a backward angle of 30° is larger, while the reflectance at a forward angle of 30° is similar, with a mean difference of 0.014, 0.027 and −0.007, respectively. The albedo of all pixels, except for the water and fill pixels, is simulated based on F and directional reflectance. All the F-based albedos have a high consistency with the MODIS albedo. In the estimation of surface albedo, large errors may occur without considering the anisotropy reflectance effect. However, F, which represents the dominant reflectance anisotropy within a MODIS tile can improve the accuracy of albedo retrieval from a single directional reflectance.

Figure 15 shows the comparison of MODIS WSA with directional reflectance or with F-based WSA. The results indicate that the directional reflectance may overestimate or underestimate the WSA, while the F-based albedo retrieval from the directional reflectance can reduce the effect of reflectance anisotropy with RMES change from 0.051 to 0.036.

4. Discussion

Fine-resolution albedo products are very helpful in various applications such as agricultural and forest management, and ecological environment monitoring and assessments. Additionally, it is crucial for validating coarse-resolution satellite albedos, particularly over heterogeneous land surfaces. With the rapid development of remote-sensing technology, many fine-resolution satellites have been launched or will be launched soon. Therefore, the need to pursue research on new methodologies for the derivation of fine-resolution land-surface albedos is increasing. In this study, we investigated the distribution the characteristics of reflectance anisotropy based on the MODIS BRDF product, and proposed a probability-weighted method to extract a priori BRDF from a MODIS tile. A priori BRDF was used to improve albedo retrieval from a single directional surface reflectance.

The F used in this study was extracted from the normalized model parameters, which is one of the main differences between this method and others. The normalization process removes the effect of spectral reflectance amplitude and makes the BRDF comparable. The corresponding F can better represent the dominant reflectance anisotropy. Another difference is that pixels that have a highly aggregative state, as opposed to “pure” MODIS pixels, were used to extract F. The reason for this is that the pure pixel in the MODIS scale is only relative, and the pure pixels used in the BRDF assignation procedure are only spectrally similar and may not necessarily be similar in surface anisotropy. Furthermore, sometimes, the chance of finding “pure” MODIS pixels for a specific class is lower. This study considered the distribution characteristics of reflectance anisotropy and neglected some discrete pixels in the process of extracting F, which could filter out some exception parameters and help to obtain a more representative F. It should be emphasized that a priori BRDFs can only be used for the retrieval of albedo since albedo is integral to BRDF over the viewing and illuminating hemispheres, since the integration process of BRDF weakens the effect of reflectance anisotropy on albedo. However, this a priori knowledge is not so accurate in retrieving other vegetation structures, such as the clumping index [59,60], particularly when a specific directional reflectance is needed.

Among studies on a priori BRDF extraction from historical MODIS BRDF for albedo retrieval, surface land-cover type and NDVI are usually taken as the basis for classifying reflectance anisotropy. The anisotropy of the reflectance is closely related to the structure of the target, in theory; however, both land-cover type and the NDVI mainly reveal the ground spectral information [61]. Based on AFX and archetypal BRDFs, we have proven that the land-cover type and NDVI are not necessarily appropriate for identifying reflectance anisotropy [43,62]. Although the albedo retrieved from LCT- or NDVI-based a priori BRDF shows a varying degree of consistency with ground measurements [40,42], this does not mean the LCT or NDVI can effectively distinguish reflectance anisotropy. In this study, a notable improvement in LCT or NDVI-based F only appears when both vegetated and barren ground occupy a large proportion in a tile. In most situations, similar results also can be obtained based only on F extracted from all the data. All the a priori BRDFs obtained via different methods captured the dominant reflectance anisotropy, which represent the majority of reflectance anisotropy of the ground. It is important to note that F is very similar to the intermediate archetype BRDF, which can correct the effect of reflectance anisotropy and obtain a relatively stable albedo [44,63].

Concurrent coarse-resolution BRDF data are an essential item for retrieving fine-resolution albedo from nadir reflectance. However, the BRDF assignation procedure may not be effective when sufficient MODIS BRDF data are not available or the size of the study area is too small. In this case, a priori BRDF would be borrowed from spatially adjacent images on the same day. The proposed algorithm provides another possible approach to borrowing a priori BRDF from the time adjacent images. The F extracted from a MODIS tile, based on the probability-weighted method, showed similar interannual variability. Even if the clear-sky MODIS data on the same day are not available, a priori BRDFs near real-time or in a different year may be able to provide an accurate estimation of reflectance anisotropy for albedo retrieval. Another potential advantage is that the present method reduces computational complexity, while providing accurate reflectance anisotropy for albedo retrieval. One specific a priori BRDF can be applied to the entire MODIS tile area for most viewing directions and ensure that most of the pixels have high accuracy. For the observations of most of the polar orbital satellites located around the cross-principal plane (CPP), F usually has a good performance. A priori BRDF extracted from a MODIS tile may be unreasonable, especially for the tiles that have multiple gathering centers. This study only provides a feasible method for a priori BRDF extraction and carries out some simple validation using simulated and real satellite data. More efforts are needed to further refine the method and find the best balance between reflectance anisotropy properties and spatial scale.

Comparing the Landsat surface reflectance and the F-based albedo with ground measurements indicates that the Lambertian assumption could introduce significant errors into the albedo, while after the correction of a priori BRDF, the uncertainty in albedo is obviously reduced. Surface reflectance anisotropy is one of the many factors that affect albedo retrieval from directional reflectance. Since the proposed method needs a single directional reflectance as the most important input, an accurate atmospheric correction is desperately required to generate reliable surface reflectance from the TOA reflectance [64]. However, in general, accurate aerosol loadings around the satellite overpass time are not available for most fine-resolution satellite data; therefore, the atmospheric correction may sometimes be unreliable for obtaining an accurate albedo. Beyond that, the F-based albedo of some pixels shows a relatively large difference compared with the MODIS albedo product, and this problem requires further study for an improved understanding of the characteristics of reflectance anisotropy.

Scale issues in the retrieval of fine-resolution albedo should not be neglected since the MODIS BRDF and Landsat reflectance have different spatial resolutions. Studies have shown that albedo can be significantly impacted by the use of BRDFs of different resolutions [31,34]. However, because of the lack of BRDF data at a fine resolution, in most studies, the reflectance anisotropy information, used for fine-resolution reflectance anisotropy correction or albedo retrieval, is extracted from a coarse-resolution BRDF product [20,21,22,28,30]. The results showed that these coarse-resolution a priori BRDFs can yield reasonable results when applied to fine-resolution reflectance. For the purpose of this study, we expected to find a specific a priori BRDF which can represent the reflectance anisotropy of most pixels for converting directional reflectance to albedo. Following the strategy used in earlier studies, this study does not consider the scale effect, and the extracted BRDF can still yield accurate results when used as a priori knowledge to retrieve the albedo.

5. Conclusions

In this study, the distribution characteristics of reflectance anisotropy over one MODIS tile were discussed, and a priori BRDF was extracted using the probability-weighted method. The effect of reflectance anisotropy on albedo retrieval was learned by adjusting F to nadir reflectance. The possibility of using F to retrieve albedo from various viewing zenith angles, the ability of LCT or NDVI to distinguish reflectance anisotropy, and the improvement of F_L and F_N for albedo retrieval from a single directional reflectance were also discussed. The possibility of using F to retrieve albedo from a single directional reflectance was validated using global (on DOY 170, 2010) and time-series (over the four tiles) MODIS BRDF data. The F extracted from MODIS was also applied to the Landsat surface reflectance to retrieve the fine-resolution albedo, and the Landsat albedo was validated using ground measurements.

Based on the four tiles, the surface reflectance anisotropy, represented by normalized BRDF model parameters in each spectral band and different tiles, showed a pattern of aggregated distribution. Except for the tiles that had a large percentage of both vegetated and non-vegetated surface and had multiple gathering centers, most of the tiles only had one gathering center. In the shortwave band, the reflectance anisotropy characteristics of different land-cover types or with a certain range of NDVI also presented aggregated distribution patterns. Significant differences appeared only when the NDVI changed from greater than 0.2 to less than 0.2, or when the ground changed from vegetated to barren. For the tiles that had a single anisotropy gathering center, F could capture the dominant reflectance anisotropy of most pixels.

The albedo retrieved from F and nadir reflectance had a high consistency with the MODIS albedo. When the solar zenith angle was around 40°, the nadir reflectance and F-based albedo were close to the MODIS albedo. With an increase or decrease in the solar zenith angle, the difference between nadir reflectance and MODIS albedo was gradually increased. However, the F-based albedo always kept a relatively high accuracy. An improvement in F_L and F_N occurred when the surface contained a large proportion of both vegetated and barren ground. The global MODIS BRDF data showed that the accuracy of the albedo retrieved from nadir reflectance can be improved efficiently by introducing F. The improvement in F_L and F_N was not obvious when compared with F, and only seven tiles that had a P_0.02 increased more than 10%. A priori BRDF can also be used to generate a reliable fine-resolution albedo from the directional reflectance. The fine-resolution albedo retrieved from the Landsat surface reflectance was highly consistent with ground measurements, with an RMSE of 0.027; however, the influence of spatial resolution on the angular variation patterns of optical reflectance between Landsat and MODIS cannot be considered at this stage. This approach could improve the accuracy of albedo compared with the Lambertian assumption.

For the entire viewing hemisphere, except the hotspot and large viewing-zenith-angle area, the F-based albedo from directional reflectance had high consistency with the MODIS albedo, with more than 90% of pixels having an absolute difference between −0.02 and 0.02. However, when the tile contained a large proportion of both vegetated and barren ground, the directions that had low accuracy would increase. When the solar zenith angle was large, the accuracy around the forward direction would also decrease. The accuracy of F_L- and F_N-based albedo around these directions could be improved. The time-series BRDF product over the four MODIS tiles also proved these conclusions. Beyond that, F showed obvious seasonal features, and the interannual variability of F showed remarkable regularity. F-based albedo retrieval from simulated or real directional reflectance also had high consistency with the high-quality MODIS albedo product.