Impact of BRDF Spatiotemporal Smoothing on Land Surface Albedo Estimation

: Surface albedo, as a key parameter determining the partition of solar radiation at the Earth’s surface, has been developed into a satellite-based product from various Earth observation systems to serve numerous global or regional applications. Studies point out that apparent uncertainty can be introduced into albedo retrieval without consideration of surface anisotropy, which is a challenge to albedo estimation especially from observations with fewer angular samplings. Researchers have begun to introduce smoothed anisotropy prior knowledge into albedo estimation to improve the inversion efﬁciency, or for the scenario of observations with signal or poor angular sampling. Thus, it is necessary to further understand the potential inﬂuence of smoothed anisotropy features adopted in albedo estimation. We investigated the albedo variation induced by BRDF smoothing at both temporal and spatial scales over six typical landscapes in North America using MODIS standard anisotropy products with high quality BRDF inversed from multi-angle observations in 500 m and 5.6 km spatial resolutions. Components of selected typical landscapes were assessed with the conﬁdence of the MCD12 land cover product and 30 m CDL (cropland data layer) classiﬁcation maps followed by an evaluation of spatial heterogeneity in 30 m scale through the semi-variogram model. High quality BRDF of MODIS standard anisotropy products were smoothed in multi-temporal scales of 8 days, 16 days, and 32 days, and in multi-spatial scales from 500 m to 5.6 km. The induced relative and absolute albedo differences were estimated using the RossThick-LiSparseR model and BRDFs smoothed before and after spatiotemporal smoothing. Our results show that albedo estimated using BRDFs smoothed temporally from daily to monthly over each scenario exhibits relative differences of 11.3%, 12.5%, and 27.2% and detectable absolute differences of 0.025, 0.012, and 0.013, respectively, in MODIS near-infrared (0.7–5.0 µ m), short-wave for Grassland. The SD of ENF in three broad bands are NIR 2.20, SW 0.97, and VIS 0.52 which are also higher than NIR 0.35, SW 0.22, and VIS 0.16 for Grassland. This indicates that smoothed BRDFs for ENF have more discreteness and variability than Grassland due to spatial smoothing on the corresponding days. Considering albedo is the hemispherical integration of BRDF, variations in BRDF due to spatial smoothing has potential impact on albedo retrievals.


Introduction
Land surface albedo, defined as the ratio of upwelling to downwelling solar radiation at the land surface over the whole solar spectrum, is an important parameter determining the partition of solar energy between Earth's surface and atmosphere [1,2]. Surface albedo quantifies the solar energy intercepted by components within various ecosystems to trigger more energy cycles, such as the phase transition of water via latent heat and heat transfer through conduction, convection, or radiation within components of biogeophysical and biogeochemical processes [2][3][4][5]. Land surface albedo is a major parameter of land surface models and directly affects the heat and water-vapor exchange between ground and atmosphere [6]. The budget of energy between ground and atmosphere is influenced by surface albedo, thus it affects the atmospheric movement and medium/long-term weather forecast [7]. Albedo changes due to land use/cover change (LUCC) such as deforestation have a profound impact on surface net radiation flux, which should not be ignored in the study of land surface temperature and monsoon [8][9][10]. Therefore, land surface albedo is valuable for studies of land surface models, circulation models, LUCC and the global climate.
The accuracy requirements of albedo keep increasing with the need from various user communities or with increasing understanding of processes of earth systems. An albedo accuracy of 0.02-0.05 is required by the community of global climate change [11,12] but increases to 0.02 for regional investigations [13]. Further, a revised land surface model based on albedo with higher accuracy could modify the predicted net radiation and sensible heat exchange by 5-30% [14]. A 10% error in albedo for agricultural landscape may induce a relative uncertainty in net radiation of 5% [15]. In evapotranspiration investigation, an uncertainty of 10% in albedo can introduce a potential error of 20 W/m 2 in net radiation [16]. Therefore, accurate surface albedo is important to achieve relevant energy estimation for numerous molders or applications in ecosystems, the carbon cycle, and climate communities [14,17].
Satellite remote sensing is the only way to routinely produce regional or global albedo. In the early stage, observations with a single angle were the dominant scheme provided by satellite-based sensors. Later, multi-angle sensors represented by POLDER and MISR were launched with the ability to simultaneously provide multi-angle observations [18][19][20]. In addition, multi-angle observations can be achieved in wide-swath satellite images with a trade-off between angles and time [21][22][23]. Compared with previous signal-angle estimation, modern albedo algorithms rely on multiple-angle observations to first build a BRDF (bi-directional reflectance distribution function), then integrate over incident and view hemispheres to calculate albedo. Studies have pointed out that the relative errors can reach up to 45% in inferring hemispherical reflectance from nadir reflectance without consideration of angle effects [24,25]. Thus, multi-angle observations are essential to inverse the BRDF model and guarantee the estimation of albedo accuracy.
BRDF was defined by Nicodemus in 1977 [26] and has been accepted widely as an intrinsic feature of land cover. Numerous studies have demonstrated that ground, airborne, or spaceborne observations exhibit obvious angle-effect phenomena [27][28][29]. BRDF theoretically describes the variation of directional reflectance in terms of incidentview geometry within both hemispheres [22,[30][31][32]. Multi-angle observations have been collected from ground measurements and several satellite sensors, such as POLDER [33], MISR, and MSG, and MODIS and AVHRR have achieved angle sampling from overlapped swaths within a given period. The standard BRDF products retrieved from multi-angle satellite observations include the MODIS nominal daily anisotropy dataset at spatial scales from 500 m to 5.6 km [34], POLDER at 6 km by 7 km spatial scale, and some other sporadic anisotropic products.
Remote Sens. 2022, 14, 2001 3 of 26 Anisotropy information at fine resolution is needed in the retrieval of surface albedo in the tens of meters spatial resolution. Compared with observations at kilometers, satellites providing observations at fine resolution are required to collect surface details in global or regional areas, such as the Sentinel series of Europe, the Landsat series of the United States, and the HuanJing and GF series of China. However, these observation systems can't provide instantaneous multi-angle observations, nor trade-off angle with time through swath overlapping. Thus, researchers have begun to introduce BRDF priori knowledge into surface albedo retrieval, especially at high spatial resolution due to the lack of multi-angle observations. Shuai et al., 2011 and, proposed the MODIS-Concurrent algorithm, which combined MODIS 500 m BRDF with corresponding Landsat surface reflectance to estimate surface albedo at 30 m resolution [35], followed by an extension to the pre-MODIS era, back to the 1980s, through an anisotropic LUT (look-up table) established from high-quality 500 m MCD43A-BRDF data [36]. Franch used smoothed 5.6 km MODIS BRDF with a further adjustment by NDVI in their VJB method, then decomposed it to estimate 30 m Landsat surface albedo [37,38]; Jiao established a BRDF archetype database by using a clustering algorithm and parameter smoothing of BRDF based on anisotropic flat index (AFX) [39,40]. Zhang Hu estimated surface albedo based on the archetype database classified by AFX [41]. However, there is a lack of understanding of the difference within estimated albedo induced by the anisotropy information smoothed at spatiotemporal scales. This work will focus on this issue to investigate the effects of smoothed BRDF on albedo magnitude through case studies using MODIS operational BRDF products.

Study Sites
Six study sites in the MODIS golden tile H09V05 in the Southwestern United States were selected to cover rich types of land surface as the study region (92.38 • W to 117.49 • W, 30 • N to 40 • N). Details of the six selected sites are shown in Table 1, which includes the IGBP classification type, location, and climate zone.

MODIS Products
We selected MODIS BRDF/Albedo products in three broad bands which are nearinfrared (NIR, 0.7-5.0 µm), short wave (SW, 0.3-5.0 µm), and visible (VIS, 0.3-0.7 µm) at 500 m and 5.6 km resolutions, and land cover products (Table 2), which we downloaded from the NASA official open data pool (https://ladsweb.modaps.eosdis.nasa.gov/search/ (accessed on 30 January 2021)). The operational MODIS BRDF/Albedo algorithm assumes that the land surface does not change significantly during the period of data accumulation (16 days), but with a further enhancement on the observations close to the retrieval day in the nominal daily inversion. The kernel-driven linear model expressed as RossThick-LiSparseR is adopted to determine the kernel parameters f iso f geo f vol respectively for isotropic kernel, K iso ; volumetric scattering kernel, K vol ; and geometric kernel, K geo , to estimate the gridded 500 m and 5600 m BRDF/Albedo at that shown in Equation (1), where λ is the spectral band [22,34]. MODIS land cover products MCD12C1 and MCD12Q1 provide global land cover maps at 5600 m and 500 m, at annual time steps. The MCD12Q1 product is produced using supervised classification of MODIS Nadir BRDF-adjusted reflectance and MCD12C1 is aggregated and reprojected from 500 m MCD12Q1 products [42,43]. MODIS land cover type products had an overall accuracy of 75% by cross-validation accuracy assessment [43] and the majority class confidence can assist to assess majority land cover type in MCD12C1 pixels. IGBP contains 17 types of land cover, and its original data source (AVHRR) and goal of modeling biophysical parameters of land surface is similar to MODIS [42][43][44][45]. Thus, we selected IGBP in MCD12C1 and MCD12Q1 to provide the classification of land cover.

CDL Database
CDL (cropland data layer) annually produced by USDA-NASS (United States Department of Agriculture-National Agricultural Statistics Service) is a geospatial crop-specific open 30 m land cover thematic dataset and serving numerous user communities. The 2019 CDL is generated from Landsat-8 OLI/TIRS, DEIMOS-1, UK2, LISS-3 of ResourceSat-2, SENTINEL-2A/B and downloaded from https://nassgeodata.gmu.edu/CropScape/ (accessed on 16 February 2021). The crop-specific classification map with a spatial resolution of 30 m, covering the continental United States, is produced by satellite imagery. The quality of CDL is high, with the total crop mapping accuracies ranging from 85% to 95% for the major crop categories [46]. The accuracy of CDL in non-agricultural land cover classes ranges from 82-89% depending upon the USGS, National Land Cover Database (NLCD) [47]. We selected CDL to assess the land cover types in MODIS pixels in detail.

Landsat-8 Dataset
To evaluate the spatial heterogeneity of land surface in pixels at coarse spatial scale, the corresponding land surface albedo at higher resolution is needed as the basic data. The Landsat-8 satellite is in a sun-synchronous, near-polar, 705 km circular orbit and can produce scenes within 185 km × 180 km. Compared with MODIS, Landsat-8 OLI (operational land imager) can provide reflectance product of land cover at a spatial resolution of 30 m and can be used to reconstruct the surface albedo at 30 m. Landsat-8 has a repeat cycle of 16 days, which can provide sequential repeated images 22 or 23 times for a fixed location per year [48]. Landsat-8 OLI provides spectral reflectance information of land cover in 9 bands from visible to shortwave infrared. We selected landsat-8 OLI surface reflectance products after radiometric and geometric correction under cloud-free condition in 2019 as the basic data and the datasets for six study sites were downloaded from https://earthexplorer.usgs.gov/ (accessed on 15 August 2021).

Data Preprocessing
A series of preprocessing steps were performed, including the selection and assessment of study areas, re-projection between different products, and screening of high quality BRDF/Albedo. MODIS land cover type products (MCD12C1 and MCD12Q1) were used to determine specific study pixels (500 m). MCD12C1 IGBP layer in 5600 m scale was used to choose study sites with high confidence, followed by reprojection to the 500 m IGBP layer  Figure 1, we can see that a representative site (500 m) is selected from homogeneous underlayer with the pixels of the same land cover type around.

Data Preprocessing
A series of preprocessing steps were performed, including the selection and assessment of study areas, re-projection between different products, and screening of high quality BRDF/Albedo. MODIS land cover type products (MCD12C1 and MCD12Q1) were used to determine specific study pixels (500 m). MCD12C1 IGBP layer in 5600 m scale was used to choose study sites with high confidence, followed by reprojection to the 500 m IGBP layer of MCD12Q1 ( Figure 1). From Figure 1, we can see that a representative site (500 m) is selected from homogeneous underlayer with the pixels of the same land cover type around. To further analyze the internal land cover types in MODIS pixels (5600 m and 500 m), CDL pixels were counted to show the proportions of land cover categories within MODIS pixels (5600 m and 500 m) ( Figure 2). It can be seen from Figure 2a-l that, compared with pixels at 5600 m, the proportions of the main categories at 500 m are significantly higher, which indicates that pixels at 500 m are more homogeneous, which is consistent with conventional cognition (the larger the range, the more complex the land cover). We selected the representative MODIS pixels (500 m) with homogeneous underlayer for the study of temporal smoothing, together with their corresponding MODIS pixels at 5600 m for the study of spatial smoothing.
To focus on the investigation of targets located on the land surface, observations contaminated by snow and cloud were removed from the temporal MODIS data series using the snow and cloud flag. Further, high-quality MODIS BRDF and albedo data labeled by related pixel-based QA (0, 1, 2) were screened out from collected datasets. This could not only ensure the acquisition of sufficient DOYs (days of the year) in 2019 with high-quality, but also eliminate the impact of data with low quality to a great extent.
In addition, Re-projection is required between different datasets. MODIS products with spatial resolution of 5600 m adopt longitude and latitude projection, while MODIS products with spatial resolution of 500 m use Sinusoidal projection (Table 2) [49][50][51]. The CDL uses Albers projection, while landsat-8 is UTM (universal transverse Mercator) projection [46,48]. When conversion was required between different types of datasets, re-projection processing was performed: (1) in order to investigate the classification within MODIS pixels (5600 m), MCD12C1 pixels were re-projected to MCD12Q1 products; (2) in order to investigate the type complexity of MODIS pixels in detail, MCD12C1 and MCD12Q1 pixels were re-projected to CDL; (3) to assess the surface heterogeneity, MCD12C1 pixels were re-projected to landsat-8 OLI images. To further analyze the internal land cover types in MODIS pixels (5600 m and 500 m), CDL pixels were counted to show the proportions of land cover categories within MODIS pixels (5600 m and 500 m) ( Figure 2). It can be seen from Figure 2a-l that, compared with pixels at 5600 m, the proportions of the main categories at 500 m are significantly higher, which indicates that pixels at 500 m are more homogeneous, which is consistent with conventional cognition (the larger the range, the more complex the land cover). We selected the representative MODIS pixels (500 m) with homogeneous underlayer for the study of temporal smoothing, together with their corresponding MODIS pixels at 5600 m for the study of spatial smoothing.

Strategy of Spatiotemporal Smoothing
To obtain BRDF in multi-temporal scales, BRDF parameters of land surface are smoothed in different temporal scales (8 days-16 days-1 month), see Figure 3a. As BRDF To focus on the investigation of targets located on the land surface, observations contaminated by snow and cloud were removed from the temporal MODIS data series using the snow and cloud flag. Further, high-quality MODIS BRDF and albedo data labeled by related pixel-based QA (0, 1, 2) were screened out from collected datasets. This could not only ensure the acquisition of sufficient DOYs (days of the year) in 2019 with high-quality, but also eliminate the impact of data with low quality to a great extent.
In addition, Re-projection is required between different datasets. MODIS products with spatial resolution of 5600 m adopt longitude and latitude projection, while MODIS products with spatial resolution of 500 m use Sinusoidal projection (Table 2) [49][50][51]. The CDL uses Albers projection, while landsat-8 is UTM (universal transverse Mercator) projection [46,48]. When conversion was required between different types of datasets, reprojection processing was performed: (1) in order to investigate the classification within MODIS pixels (5600 m), MCD12C1 pixels were re-projected to MCD12Q1 products; (2) in order to investigate the type complexity of MODIS pixels in detail, MCD12C1 and MCD12Q1 pixels were re-projected to CDL; (3) to assess the surface heterogeneity, MCD12C1 pixels were re-projected to landsat-8 OLI images.

Strategy of Spatiotemporal Smoothing
To obtain BRDF in multi-temporal scales, BRDF parameters of land surface are smoothed in different temporal scales (8 days-16 days-1 month), see Figure 3a. As BRDF can be represented by a kernel-driven linear model (Equation (1)), smoothed BRDF can also be represented by these smoothed model parameters (Equation (2)), where f k is the parameter of the isotropic, volumetric, geometric-optical kernel, scale is 8 days, 16 days, or a month, respectively, and N is the number of scale periods. To maintain the consistency of the temporal range, the number of days in the monthly scale is set to 32 days. It should be noted that the MCD43A1 is named 'daily' though it is a nominal daily period and is the weighted synthesis of surface reflectance under clear sky over a 16-day period.
For BRDF smoothing in multi-spatial scales (Figure 3b), as MODIS provides two BRDF products (MCD43A and MCD43C) at different resolutions (500 m and 5600 m), we directly extracted BRDF parameters from these datasets to represent BRDF before and after spatial smoothing. Here, we obtained BRDF parameters before and after smoothing in multi-spatiotemporal scales and the reflectance under any incident-view geometry could be inversed based on a set of fixed BRDF parameters. In this paper, the research version of the operating AMBLARS algorithm [22,30,52], provided by Boston University, was used to retrieve reflectance under given incident-view geometries based on RossThick-LiSparseR model.
According to the integral mode, MODIS-albedo can be defined as: black sky albedo (BSA), which is the integration of BRDF in the view hemisphere with direct incident light (Equation (3)); white sky albedo (WSA), which is the integration of BRDF in both incident and view hemispheres with isotropic incident light (Equation (4)). In Equations (3) and (4), θ i and ∅ i are zenith and azimuth angles of incident solar direction; θ r and ∅ r are zenith and azimuth angles of view direction, and f r (θ i , ∅ i ; θ r , ∅ r ) is the BRDF function. Both the BSA and WSA are based on a set of fixed kernel coefficients to integrate in hemispherical or bi-hemispherical space. Theoretically, albedo obtained from BRDF parameters at the coarse scale after spatial and temporal smoothing is consistent with the smoothed value of albedo before smoothing (Equation (5), and θ and ∅ are the zenith and azimuth angle in hemispherical or bi-hemispherical space).
the BSA and WSA are based on a set of fixed kernel coefficients to integrate in hemispherical or bi-hemispherical space. Theoretically, albedo obtained from BRDF parameters at the coarse scale after spatial and temporal smoothing is consistent with the smoothed value of albedo before smoothing (Equation (5), and and ∅ are the zenith and azimuth angle in hemispherical or bi-hemispherical space).
As MODIS provides two albedo products (MCD43A3 and MCD43C3) at different resolutions (500 m and 5600 m), we extracted albedo from MCD43A3 and smoothed daily albedo in multi-temporal scales to represent the integration of the corresponding As MODIS provides two albedo products (MCD43A3 and MCD43C3) at different resolutions (500 m and 5600 m), we extracted albedo from MCD43A3 and smoothed daily albedo in multi-temporal scales to represent the integration of the corresponding smoothed BRDF, meanwhile we directly extracted albedo from MCD43C3 to represent the integration of the corresponding smoothed BRDF in multi-spatial scales. The workflow of this paper is shown in Figure 4. smoothed BRDF, meanwhile we directly extracted albedo from MCD43C3 to represent the integration of the corresponding smoothed BRDF in multi-spatial scales. The workflow of this paper is shown in Figure 4.

Differences of BRDF
In the view hemisphere, the principal plane, which is determined by the directions of incident illumination and outgoing reflectance that intersect at the identical target, implies typical anisotropy features, such as variations in reflectance in magnitude and shapes of BRDF ( Figure 5a). In this paper, the difference in BRDF is investigated by comparing the magnitude and shape between BRDF principal planes. The average reflectance with an interval of 1° from −70° (backward) to 70° (forward) was used to describe magnitude of principal plane (Equation (6)).
To describe the change in BRDF shape, we took ±70°, ±45°, LSN (local solar noon), zenith direction as characteristic directions and calculated differences in reflectance compared to smoothed BRDF at characteristic directions ( Figure 5b and Equation (7), where N is the number of principal planes to be compared and j is ±70°, ±45°, LSN, 0°, respectively, and Fixed-Reflectance stands for the smoothed principal planes). Furthermore, the SD (standard deviation) of differences of reflectance at characteristic directions is adopted to measure the variation in the shape of BRDF (Equation (8)). If the distribution of differ-

Differences of BRDF
In the view hemisphere, the principal plane, which is determined by the directions of incident illumination and outgoing reflectance that intersect at the identical target, implies typical anisotropy features, such as variations in reflectance in magnitude and shapes of BRDF (Figure 5a). In this paper, the difference in BRDF is investigated by comparing the magnitude and shape between BRDF principal planes. The average reflectance with an interval of 1 • from −70 • (backward) to 70 • (forward) was used to describe magnitude of principal plane (Equation (6)).
To describe the change in BRDF shape, we took ±70 • , ±45 • , LSN (local solar noon), zenith direction as characteristic directions and calculated differences in reflectance compared to smoothed BRDF at characteristic directions ( Figure 5b and Equation (7), where N is the number of principal planes to be compared and j is ±70 • , ±45 • , LSN, 0 • , respectively, and Fixed-Reflectance stands for the smoothed principal planes). Furthermore, the SD (standard deviation) of differences of reflectance at characteristic directions is adopted to measure the variation in the shape of BRDF (Equation (8)). If the distribution of differences in reflectance in the six characteristic directions is relatively concentrated, it indicates that the differences in reflectance in the six characteristic directions is relatively uniform and the Remote Sens. 2022, 14, 2001 9 of 26 difference between BRDF shapes is small, meanwhile the SD is smaller in this situation. On the contrary, a larger SD means that the difference between the shapes of BRDF is greater.
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 27 ences in reflectance in the six characteristic directions is relatively concentrated, it indicates that the differences in reflectance in the six characteristic directions is relatively uniform and the difference between BRDF shapes is small, meanwhile the SD is smaller in this situation. On the contrary, a larger SD means that the difference between the shapes of BRDF is greater.

Differences of Albedo
Albedo is a hemispherical integration of BRDF, the magnitude and shape of BRDF have potential impacts on albedo. MODIS-albedo products are available in two types: WSA and BSA. Considering that WSA is highly correlated with BSA [53,54], we selected WSA as the albedo before and after BRDF smoothing to serve for the difference comparison.
To explore differences in albedo caused by temporal smoothing of BRDF, as shown in Figure 3a, first, smoothed monthly values of albedo for 12 months were calculated; then, the smoothed albedo values at 16-day and 8-day scales within one month were calculated; Lastly the daily albedo values were compared with the corresponding smoothed albedo values in three broad bands: near-infrared (NIR), short wave (SW), and visible (VIS). We selected the relative and absolute differences in albedo to describe the relationship between them (Equations (9) and (10), where scales are 8 days, 16 days, and 32 days). Through the absolute and relative differences in albedo, we can quantitatively analyze the difference in daily albedo from smoothed albedo in a specific temporal scale.
When discussing the influence of BRDF spatial smoothing on the retrieval of albedo, we extracted albedo at 500 m and 5600 m from MCD43A3 and MCD43C3 products, then calculated the relative and absolute differences between them. Lastly, we counted the

Differences of Albedo
Albedo is a hemispherical integration of BRDF, the magnitude and shape of BRDF have potential impacts on albedo. MODIS-albedo products are available in two types: WSA and BSA. Considering that WSA is highly correlated with BSA [53,54], we selected WSA as the albedo before and after BRDF smoothing to serve for the difference comparison.
To explore differences in albedo caused by temporal smoothing of BRDF, as shown in Figure 3a, first, smoothed monthly values of albedo for 12 months were calculated; then, the smoothed albedo values at 16-day and 8-day scales within one month were calculated; Lastly the daily albedo values were compared with the corresponding smoothed albedo values in three broad bands: near-infrared (NIR), short wave (SW), and visible (VIS). We selected the relative and absolute differences in albedo to describe the relationship between them (Equations (9) and (10), where scales are 8 days, 16 days, and 32 days). Through the absolute and relative differences in albedo, we can quantitatively analyze the difference in daily albedo from smoothed albedo in a specific temporal scale.
When discussing the influence of BRDF spatial smoothing on the retrieval of albedo, we extracted albedo at 500 m and 5600 m from MCD43A3 and MCD43C3 products, then calculated the relative and absolute differences between them. Lastly, we counted the daily differences to explore the magnitudes of the differences (Figure 3b). Albedo difference due to BRDF spatial smoothing can be calculated by Equations (11) and (12), where N is the number of MODIS-pixels at 500 m within one MODIS pixel at 5600 m.
Relative di f f erence spatial = Absolute di f f erence spatial =

Semi-Variogram
The spatial heterogeneity is related to the distribution and complexity of land cover, and the different land cover is related to albedo differences. A semi-variogram model is an effective way to describe heterogeneity of land cover, and its model parameters have a certain corresponding relationship with the landscape pattern. Some studies have used a semi-variogram function to evaluate the homogeneity of land cover and applied it to assess the spatial representativeness of ground tower measurement points [55][56][57][58]. As the landscape pattern may affect the distribution of spatial albedo, we analyzed the correlation between the land surface heterogeneity described by semi-variogram parameters and the differences in albedo due to BRDF spatial smoothing.
In order to evaluate the spatial heterogeneity within MODIS pixels, the surface albedo at 30 m needed to be reconstructed based on landsat-8 images under cloud-free weather in 2019. Among the 22 or 23 images over the whole year, we selected images from leaf-on or leaf-off stages under cloud-free weather to compare the seasonality of the surface. To completely cover MODIS pixels, we made a semi-variogram analysis of land cover within 6 km (slightly greater than 5600 m), centered on the selected study sites. The specific steps were as follows: (1) Pixels corresponding to the study sites were extracted from Landsat-8 OLI surface reflectance products in band 2 to band 7; (2) Narrowband-to-broadband parameters were used to produce albedo in three broad bands (near-infrared, short wave, and visible) from Landsat-8 pixels. As the imaging bands of Landsat OLI are different from the three broadbands in MODIS products, narrowband-to-broadband conversion was required (Equations (13)-(15)), where α is the albedo in the corresponding band [59,60]; (3) Equation (16) was used to calculate half the average-squared-difference between pixels with three broadband albedos. The parameters of the semi-variogram model, such as range and sill, were fitted by sequential half the average-squared-difference (Equation (17)). In Equations (16)- (17), h is the lag distance, N(h) is the number of pixel pairs corresponding to h, Z x is the albedo at location x, C 0 is the nugget, C is the sill, and a is the range.
Range (a) describes the correlation distance of point pairs in study sites. There is a certain correlation between point pairs within range distance, and there is no correlated attribute between point pairs beyond this distance. Range is of great significance to judge spatial representativeness between different scales. Sill (C) describes the intensity of albedo change, and a smaller sill represents a more homogeneous land cover. In this paper, the area around the samples is fixed (6 km, covering 5600 m pixel), and sill (C) is used to analyze the heterogeneity of land cover. α V IS = 0.5621α 2 + 0.1479α 3 + 0.2512α 4 − 0.0015 (13) α N IR = 0.5911α 5 + 0.3155α 6 + 0.0731α 7 + 0.0019 (14) α SW = 0.2453α 2 + 0.0508α 3 + 0.1804α 4 + 0.3081α 5 + 0.1332α 6 + 0.0521α 7 + 0.0011 (15)

Vegetation Index
EVIs shown in Figure 6a were calculated from the zenith reflectance provided by MCD43A4. Then, the seasonal distribution of EVI could be counted according to spring (March to May), summer (June to August), autumn (September to November) and winter (December, January and February) (Figure 6b). EVI (enhanced vegetation index) can reflect the growth and coverage of land cover and has a significant advantage in vegetation saturation, atmospheric impact, and soil background compared with other vegetation indexes [61][62][63]. It can be seen from Figure 6a,b that, generally, the EVIs of land cover in spring and summer are higher than in autumn and winter, which indicates that vegetation has a stronger growth state and higher coverage in spring and summer, which are more suitable for growth. But for EBF, the EVI in spring is slightly lower than that in autumn. This may be attributed to the change in EVI over the whole year being relatively mild for EBF growing in four seasons, and that the EVI in winter is high, resulting in a high EVI in autumn as the transition between summer and winter. For ENF, we only extracted data without snow from June to November, while the EVI in summer is higher than autumn. For DBF, Savannas, Grassland, and Cropland, EVI shows the obvious seasonal characteristic that EVIs in spring and summer are higher. = 0.5911 + 0.3155 + 0.0731 + 0.0019 (14) = 0.2453 + 0.0508 + 0.1804 + 0.3081 + 0.1332 + 0.0521 + 0.0011 (15)

Vegetation Index
EVIs shown in Figure 6a were calculated from the zenith reflectance provided by MCD43A4. Then, the seasonal distribution of EVI could be counted according to spring (March to May), summer (June to August), autumn (September to November) and winter (December, January and February) (Figure 6b). EVI (enhanced vegetation index) can reflect the growth and coverage of land cover and has a significant advantage in vegetation saturation, atmospheric impact, and soil background compared with other vegetation indexes [61][62][63]. It can be seen from Figure 6a,b that, generally, the EVIs of land cover in spring and summer are higher than in autumn and winter, which indicates that vegetation has a stronger growth state and higher coverage in spring and summer, which are more suitable for growth. But for EBF, the EVI in spring is slightly lower than that in autumn. This may be attributed to the change in EVI over the whole year being relatively mild for EBF growing in four seasons, and that the EVI in winter is high, resulting in a high EVI in autumn as the transition between summer and winter. For ENF, we only extracted data without snow from June to November, while the EVI in summer is higher than autumn. For DBF, Savannas, Grassland, and Cropland, EVI shows the obvious seasonal characteristic that EVIs in spring and summer are higher.   Figure 7c,d, which shows that the spatial heterogeneity of Savanna on day 205 of 2019 is lower than that on day 365 of 2019 (DOY 365). The reason can be explained as smaller sill (C) corresponding to lower spatial heterogeneity of land cover, and sill (C) being an effective parameter to describe spatial heterogeneity. × 10 −6 ) due to crop rotation and other factors. Meanwhile, the spatial heterogeneity of samples growing in four seasons is more complex. The sill (C) value of EBF on DOY319 within the senescence season is smaller than on DOY207, which indicates that spatial homogeneity is lower. However, heterogeneity of ENF in late autumn (DOY277) is higher than that in early autumn (DOY247). This shows that variation in land surface heterogeneity with temporal series is very complex. Therefore, we took the DOYs with Landsat8 OLI images to explore the relationship between the albedo difference due to BRDF spatial smoothing and the spatial heterogeneity.  Table 3 shows the parameters of the semi-variogram for each study site in leaf-on and leaf-off seasons, corresponding to the growing and senescence seasons. For DBF, Savanna and Grassland samples with four distinct seasons, sill (C) values in the senescence season are significantly higher than those in the growing season. In the NIR and SW bands, Cropland shows similar characteristic to other samples, but the heterogeneity of the VIS band in winter (DOY352, 842 × 10 −6 ) is lower than that in summer (DOY192, 1128 × 10 −6 ) due to crop rotation and other factors. Meanwhile, the spatial heterogeneity of samples growing in four seasons is more complex. The sill (C) value of EBF on DOY319 within the senescence season is smaller than on DOY207, which indicates that spatial homogeneity is lower. However, heterogeneity of ENF in late autumn (DOY277) is higher than that in early autumn (DOY247). This shows that variation in land surface heterogeneity with temporal series is very complex. Therefore, we took the DOYs with Landsat8 OLI images to explore the relationship between the albedo difference due to BRDF spatial smoothing and the spatial heterogeneity.  Figure 8 shows the principal planes after BRDF smoothing in one month compared with daily BRDF over selected study sites. We took June of 2019, in which the EVI level is high, as an example to show the variation in BRDF smoothing. It can be seen from Figure 8 that, compared with daily BRDF, the magnitudes and shapes of smoothed BRDF have changed after temporal smoothing. In addition, smoothed BRDF is almost in the middle location of daily principal plane curves and the ability to capture features of BRDF is reduced. To quantify the changes in BRDF, the difference in magnitude (Equation (6)) and standard deviation of multi-characteristic directions (Equations (7) and (8)) within principal planes were to be counted ( Figure 9). Figure 9a,b shows that the magnitude difference is largest in the NIR band and lowest in the VIS band for each selected study site, for example, the magnitude differences in the NIR, SW, and VIS bands are 0.63%, 0.27%, and 0.09%, respectively, for ENF in 2019/06. The standard deviation of each study site shows a similar characteristic in that the average SDs in the NIR, SW, and VIS bands are 0.32, 0.20, and 0.14, respectively. This indicates that the dispersion and shape change of BRDF are also gradually decreasing from the NIR and SW to the VIS band in 2019/06, which is consistent with the information shown in Figure 8. Figure 9c,d shows the statistics of average magnitude difference and standard deviation in principal planes in 2019 for each study site. The average magnitude difference and standard deviation are highest in NIR and lowest in VIS for each study site. This shows that the spectral band is an important reason for the variation in BRDF. For land cover types, each study site shows various levels of differences: the average difference in BRDF magnitude in three broad bands is EBF 0.25%, DBF 0.41%, Savanna 0.39%, Grassland 0.35%, Cropland 0.46%, and ENF 0.26%, and the average standard deviation is EBF 0.15, DBF 0.31, Savanna 0.38, Grassland 0.32, Cropland 0.51 and ENF 0.34, respectively. This indicates that variations exist in the magnitude and shape of BRDF after smoothing in multi-temporal scales and the differences are related to land cover types and spectral bands. Figure 10 shows the comparison of BRDF principal planes at 5600 m with the BRDF at 500 m over selected ENF and Grassland. From our study sites, we selected ENF and Grassland to show the variation of BRDF after spatial smoothing as the largest and smallest daily albedo difference occurred in ENF and Grassland respectively and the BRDF variation is representative on the extreme day of albedo difference.    Figure 10 shows the comparison of BRDF principal planes at 5600 m with the BRDF at 500 m over selected ENF and Grassland. From our study sites, we selected ENF and Grassland to show the variation of BRDF after spatial smoothing as the largest and smallest daily albedo difference occurred in ENF and Grassland respectively and the BRDF variation is representative on the extreme day of albedo difference. Figure 10a shows the variations in BRDF shapes are largest in the NIR band and smallest in the VIS band, meanwhile the locations of smoothed BRDF are not in the middle of daily BRDF principal planes strictly of three broad bands for ENF. However, Figure 10b shows variations in BRDF shapes of ENF are relatively stronger than those of Grassland and the locations of smoothed BRDF are almost in the middle of daily BRDF principal planes strictly of three broad bands for Grassland. To quantify the variations in BRDF, the difference in magnitude (Equation (6)) and standard deviation of multi-characteristic directions (Equations (7) and (8)) within BRDF principal planes needed to be counted (Table 4). Table 4 shows that the magnitude differences in ENF are NIR 1.57%, SW 0.85%, and VIS 0.39% on DOY176, DOY180, and DOY180, which are higher than NIR 0.86%, SW 0.37%, and VIS 0.29% on DOY152, DOY153, and DOY260 for Grassland. The SD of ENF in three broad bands are NIR 2.20, SW 0.97, and VIS 0.52 which are also higher than NIR 0.35, SW 0.22, and VIS 0.16 for Grassland. This indicates that smoothed BRDFs for ENF have more discreteness and variability than Grassland due to spatial smoothing on the corresponding days. Considering albedo is the hemispherical integration of BRDF, variations in BRDF due to spatial smoothing has potential impact on albedo retrievals.  Figure 10a shows the variations in BRDF shapes are largest in the NIR band and smallest in the VIS band, meanwhile the locations of smoothed BRDF are not in the middle of daily BRDF principal planes strictly of three broad bands for ENF. However, Figure 10b shows variations in BRDF shapes of ENF are relatively stronger than those of Grassland and the locations of smoothed BRDF are almost in the middle of daily BRDF principal planes strictly of three broad bands for Grassland. To quantify the variations in BRDF, the difference in magnitude (Equation (6)) and standard deviation of multi-characteristic directions (Equations (7) and (8)) within BRDF principal planes needed to be counted (Table 4).       Table 4 shows that the magnitude differences in ENF are NIR 1.57%, SW 0.85%, and VIS 0.39% on DOY176, DOY180, and DOY180, which are higher than NIR 0.86%, SW 0.37%, and VIS 0.29% on DOY152, DOY153, and DOY260 for Grassland. The SD of ENF in three broad bands are NIR 2.20, SW 0.97, and VIS 0.52 which are also higher than NIR 0.35, SW 0.22, and VIS 0.16 for Grassland. This indicates that smoothed BRDFs for ENF have more discreteness and variability than Grassland due to spatial smoothing on the corresponding days. Considering albedo is the hemispherical integration of BRDF, variations in BRDF due to spatial smoothing has potential impact on albedo retrievals.

Albedo Differences Induced by BRDF Temporal Smoothing
Figure 11a-c shows the relative differences (Equation (9)) and absolute differences (Equation (10)) in albedo due to BRDF temporal smoothing in multi-temporal scales (8 days, 16 days and 32 days) in three broad bands for each selected study site. As the smoothing scale increases from 8 days to 16 days and 32 days, the albedo differences (relative and absolute) compared to smoothed albedo also gradually increase. It indicates that a larger temporal scale produces higher albedo difference, and BRDF at coarser temporal scale has lower temporal representativeness.

Albedo Differences Induced by BRDF Temporal Smoothing
Figure 11a-c shows the relative differences (Equation (9)) and absolute differences (Equation (10)) in albedo due to BRDF temporal smoothing in multi-temporal scales (8 days, 16 days and 32 days) in three broad bands for each selected study site. As the smoothing scale increases from 8 days to 16 days and 32 days, the albedo differences (relative and absolute) compared to smoothed albedo also gradually increase. It indicates that a larger temporal scale produces higher albedo difference, and BRDF at coarser temporal scale has lower temporal representativeness. Figure 11. Statistics of albedo differences due to BRDF temporal smoothing in NIR (a), SW (b), VIS (c) and the average monthly albedo differences (d) for each study site. Figure 11d shows the statistic of average monthly albedo difference in three broad bands for each selected study site. From the perspective of spectral band, the magnitude of relative difference in the VIS band (7.17%) is the highest, followed by NIR (3.36%) and SW (3.12%), while the absolute difference in the NIR band (0.67%) is the highest followed by (0.40%) and VIS (0.27%). This indicates that albedo in VIS band is generally the lowest among the three broad bands while albedo in the NIR and SW bands is relatively larger, which is consistent with the characteristic that vegetation absorbs energy in VIS band and strongly reflects in the NIR band. Thus, albedo differences in vegetation due to BRDF smoothing in multi-temporal scales are highest in the NIR band absolutely, and largest in the VIS band relatively.
From the perspective of land cover type, the relative and absolute differences in EBF in the three broad bands are smaller than DBF, indicating that EBF-BRDF has better temporal representativeness. Thus, EBF which is in low latitude and growing in four seasons Figure 11. Statistics of albedo differences due to BRDF temporal smoothing in NIR (a), SW (b), VIS (c) and the average monthly albedo differences (d) for each study site. Figure 11d shows the statistic of average monthly albedo difference in three broad bands for each selected study site. From the perspective of spectral band, the magnitude of relative difference in the VIS band (7.17%) is the highest, followed by NIR (3.36%) and SW (3.12%), while the absolute difference in the NIR band (0.67%) is the highest followed by (0.40%) and VIS (0.27%). This indicates that albedo in VIS band is generally the lowest among the three broad bands while albedo in the NIR and SW bands is relatively larger, which is consistent with the characteristic that vegetation absorbs energy in VIS band and strongly reflects in the NIR band. Thus, albedo differences in vegetation due to BRDF smoothing in multi-temporal scales are highest in the NIR band absolutely, and largest in the VIS band relatively.
From the perspective of land cover type, the relative and absolute differences in EBF in the three broad bands are smaller than DBF, indicating that EBF-BRDF has better temporal representativeness. Thus, EBF which is in low latitude and growing in four seasons shows better stability than DBF. For DBF, Savanna, and Grassland growing in four seasons, the average relative differences in the three broad bands (NIR, SW, and VIS) are 6.29%, 3.71%, and 2.49%, respectively. Considering Savanna is a collection of forest and grassland, while Grassland brings lower differences than DBF, thus with the increase in deciduous forest, the relative difference could be significantly increased. Both Cropland and Grassland are sites with uniform land cover types and four distinctive seasons. However, the relative and absolute differences in Cropland in three broad bands are higher than those of Grassland. This may be due to the influence of artificial disturbance and crop rotation. ENF, which is at high latitude and is growing in four seasons, only has data in half the year due to the climate. ENF shows the highest relative differences in summer and autumn, meanwhile their relative and absolute differences are significantly higher than those of EBF, indicating that the temporal representativeness of ENF-BRDF is poor.
From the above analysis, BRDF smoothing at coarser temporal scales can produce higher albedo differences, and the albedo differences of vegetation due to temporal smoothing of BRDF are highest in the NIR band absolutely, and largest in the VIS band relatively. Land cover type can significantly affect albedo differences. EBF, which is in low latitude and growing in four seasons shows better temporal representativeness compared to other forest samples (such as DBF and ENF), meanwhile the proportion of forest and artificial disturbance are also influencing factors for albedo differences. Figure 12a-d shows the monthly relative and absolute albedo differences which were counted according to spring (March-May), summer (June-August), autumn (September-November), and winter (December, January, and February). While ENF only has data from June to November with snowless coverage, only data in summer and winter were counted. It can be seen from Figure 12a-c that the seasonal characteristics and trends of relative difference and absolute difference are consistent for each study site. We combined the data in Figure 12a-c to obtain Figure 12d. The data for each season within Figure 12d comes from the average values of the six study sites in Figure 12a-c. Figure 12d shows that albedo differences due to temporal smoothing of BRDF in spring and summer are significantly higher than those in autumn and winter in the three broadbands (NIR, SW, and VIS). As EVI curves (Figure 6a) show that land covers in spring and summer have stronger growing states and higher vegetation coverage for study sites with obvious growing cycles. This indicates that the growing state of vegetation could significantly affect albedo differences of land cover. Spring and summer, which are suitable for vegetation growth, are more likely to lead to albedo differences, and BRDF is less representative in spring and summer. The relative albedo differences of NIR 4.57% and SW 3.96% in spring are higher than those in summer with NIR 3.73% and SW 3.63%. The absolute differences of NIR 1.08% and SW 0.59% in spring are also higher than those in summer with NIR 0.76% and SW 0.46%. This means that the temporal representativeness of BRDF is lowest in spring in the NIR and SW bands. However, both the relative albedo difference and absolute albedo difference in summer, with values of 10.67% and 0.37%, are higher than those in spring with 7.13% and 0.35% for the VIS band. This indicates that the temporal representativeness of BRDF is the lowest in summer in the VIS band, which may be due to the highest vegetation index in summer and strong absorption by photosynthesis. Autumn and winter are insensitive to albedo differences due to temporal smoothing, BRDF in autumn and winter are more representative. Thus, albedo differences show obvious seasonal characteristics due to BRDF temporal smoothing. Table 5 shows the extreme monthly albedo differences of albedo due to BRDF temporal smoothing at all study sites. The smoothed BRDF from daily to monthly phase in three broad bands (NIR, SW, and VIS) could produce relative differences of 11.3%, 12.5%, and 27.2%, and absolute differences of 0.025, 0.012, and 0.013. Thus, albedo differences due to temporal smoothing of BRDF are significant and should not be ignored.   Table 5 shows the extreme monthly albedo differences of albedo due to BRDF temporal smoothing at all study sites. The smoothed BRDF from daily to monthly phase in three broad bands (NIR, SW, and VIS) could produce relative differences of 11.3%, 12.5%, and 27.2%, and absolute differences of 0.025, 0.012, and 0.013. Thus, albedo differences due to temporal smoothing of BRDF are significant and should not be ignored.  Figure 13 shows the daily relative and absolute albedo difference due to spatial smoothing of BRDF from 500 m to 5600 m for each selected site (Equations (11) and (12)). Figure 14 shows statistics of the relative and absolute differences at annual level in the three broad bands (NIR, SW, and VIS). Combining the data from the two figures, the average absolute difference in NIR (1.33%) is higher than SW (0.81%) and VIS (0.51%), while the average relative difference in VIS (12.53%) is higher than NIR (6.70%) and SW (6.37%). This is consistent with the characteristic of green vegetation with higher albedo in the NIR band and lower albedo in the VIS band as photosynthesis needs to absorb visible light,   Figure 13 shows the daily relative and absolute albedo difference due to spatial smoothing of BRDF from 500 m to 5600 m for each selected site (Equations (11) and (12)). Figure 14 shows statistics of the relative and absolute differences at annual level in the three broad bands (NIR, SW, and VIS). Combining the data from the two figures, the average absolute difference in NIR (1.33%) is higher than SW (0.81%) and VIS (0.51%), while the average relative difference in VIS (12.53%) is higher than NIR (6.70%) and SW (6.37%). This is consistent with the characteristic of green vegetation with higher albedo in the NIR band and lower albedo in the VIS band as photosynthesis needs to absorb visible light, which also indicates that albedo differences due to spatial smoothing are affected by the spectral band.

Figure 14.
Statistics of albedo differences due to BRDF spatial smoothing at annual level for each site in three broad bands. Table 6 shows the maximal daily relative and absolute albedo differences for each selected sites due to spatial smoothing of BRDF. The highest relative differences of 36.5%, 37.1%, and 94.7% and absolute differences of 0.037, 0.024, and 0.018 occurred in NIR, SW, and VIS after spatial smoothing of BRDF. Figure 14. Statistics of albedo differences due to BRDF spatial smoothing at annual level for each site in three broad bands.
Regarding land cover types, both the annual relative and absolute albedo differences of Grassland are smallest in all study sites. Thus, Grassland-BRDF had the best spatial representativeness. Compared with Grassland, albedo differences for the three forest sites (EBF, DBF, and ENF) are complex and higher. With increasing latitude, the relative albedo differences become greater. This shows that spatial representativeness of EBF-BRDF is higher than DBF and ENF.
From DBF, to Savanna, to Grassland, with the decrease in forest proportion, the relative and absolute differences gradually decrease. This indicates that forest effectively increases the albedo difference. However, with similar roughness and uniform land cover types, the albedo differences for Cropland are significantly higher than Grassland due to the disturbance of human activities. Thus, albedo differences due to spatial smoothing of BRDF are influenced by characteristics of vegetation and land cover types. Table 6 shows the maximal daily relative and absolute albedo differences for each selected sites due to spatial smoothing of BRDF. The highest relative differences of 36.5%, 37.1%, and 94.7% and absolute differences of 0.037, 0.024, and 0.018 occurred in NIR, SW, and VIS after spatial smoothing of BRDF. Figure 15 shows the statistical analysis of sill (C) values combined with albedo differences and RMSE on the DOYs with coverage by Landsat-8 images. There is a high correlation between sill values and absolute differences or RMSE, while the overall correlation coefficients in the three broad bands is 0.8876 and 0.8674, respectively. Thus, the heterogeneity of land cover is related to albedo difference. The higher the surface heterogeneity, the higher the Sill value, the lower the spatial representativeness of BRDF at coarser scales, and the larger the albedo difference.  Figure 15 shows the statistical analysis of sill (C) values combined with albedo differences and RMSE on the DOYs with coverage by Landsat-8 images. There is a high correlation between sill values and absolute differences or RMSE, while the overall correlation coefficients in the three broad bands is 0.8876 and 0.8674, respectively. Thus, the heterogeneity of land cover is related to albedo difference. The higher the surface heterogeneity, the higher the Sill value, the lower the spatial representativeness of BRDF at coarser scales, and the larger the albedo difference.

Discussion
It is a feasible approach to introduce BRDF priori knowledge into the estimation of surface albedo at fine resolutions, but there is a potential for apparent uncertainty existing in the derived albedo values, especially when BRDF priori knowledge at coarse scales is

Discussion
It is a feasible approach to introduce BRDF priori knowledge into the estimation of surface albedo at fine resolutions, but there is a potential for apparent uncertainty existing in the derived albedo values, especially when BRDF priori knowledge at coarse scales is used in fine resolution albedo calculations. This work focused on several typical case studies to investigate its variation. Our results show that BRDFs of different land cover types can produce captured differences into the retrieval of albedo after spatiotemporal smoothing. Moreover, the differences have the seasonal characteristic that differences in spring and summer are higher in temporal smoothing and heterogeneity so that the overall correlation coefficient between sill values and absolute differences is also high. In spatial smoothing, the direct causes are changes in land surface condition in the spatial and temporal scale, such as the spatial continuity of surface biophysical characteristics, various rhythmic cycles of land surface biomes with the change of seasons, and as well as residual uncertainties left in the nominal corrected input observations [40,[64][65][66]. Though observations within a 16-day period are enhanced on the most recent day of the retrieval date, it should be noted that the time-angle trade-off scheme preserves land surface signatures from other days in the nominal daily BRDF retrieval.
Regarding analysis of albedo differences due to BRDF temporal smoothing, the temporal smoothing of BRDF shows obvious scale effect and seasonality. With the increase in smoothing scale, the albedo differences gradually increase (Figure 11a-c). This conclusion is applicable to the three broad bands (NIR, SW, and VIS) for each study site. Generally, the absolute differences caused by temporal smoothing of BRDF in the NIR band are the largest, while the relative differences in the VIS band are the largest (Figure 11d). This may be related to the biophysical characteristic of vegetation that it needs to absorb light in the VIS band, meanwhile the light in NIR band could be strongly scattered by leaf structure. In addition, the impact of BRDF temporal smoothing shows certain seasonal characteristics, with elevated variation in spring and summer due to the vigorous growth of vegetation (Figure 12d), and limited variation in autumn and winter due to lesser temporal changes over land surface.
The impact of spatial BRDF smoothing on albedo implies certain spectral features accompanied by connections with spatial heterogeneity. The absolute differences caused by spatial BRDF smoothing show the highest values in the NIR band and the highest relative differences in the VIS band ( Figure 14). This is consistent with the biophysical characteristics of the vegetation analyzed above, that photosynthesis needs to absorb visible light and produce higher reflectance in the NIR band. Compared with flat grassland, forest effectively increased the differences. In addition, we reconstructed the land surface albedo at 30 m based on Landsat8-OLI images to analyze the spatial heterogeneity of land covers. The results show that the spatial heterogeneity of land cover is related to the albedo differences due to spatial smoothing of BRDF ( Figure 15).

Conclusions
Several studies have introduced BRDF priori knowledge into albedo retrievals smoothed at different scales. This paper investigated the effects of smoothed BRDF on albedo differences through case studies over six North American regions using operational MODIS-BRDF/Albedo products. Our results show that: (1) as the BRDFs smoothed temporally from daily to monthly, and spatially from 500 m to 5600 m, variations in the magnitude and shapes during BRDF smoothing can be captured, with potential relationships with spectral band, land cover types, and characteristic of vegetation. (2) Temporally smoothed BRDF in NIR, SW, and VIS could lead to apparent relative differences of estimated albedo to smoothed values of 11.3%, 12.5%, and 27.2% and detectable absolute differences of 0.025, 0.012, and 0.013. Further, albedo differences show an obvious seasonal characteristic that differences in spring and summer are significantly higher than those in autumn and winter in the three broadbands. (3) Spatially smoothed BRDF from 500 m to 5.6 km in NIR, SW, and VIS bands could lead to albedo achieving apparent relative differences to smoothed values of 36.5%, 37.1%, and 94.7% and detectable absolute differences of 0.037, 0.024, and 0.018, while albedo differences due to BRDF spatial smoothing are related to heterogeneity of land cover and the overall correlation coefficient between sill values and absolute differences in the three broad bands was 0.8876. The introduction of BRDF priori knowledge of land cover is an important way to produce albedo at high resolution. This work demonstrated that smoothed BRDF after temporal and spatial smoothing can introduce variations in the magnitude and shape of BRDF, and the uncertainty propagating into albedo retrieval is obvious. Thus, it is necessary to avoid the smoothing process in quantitative remote sensing communities, especially when immediate anisotropy retrievals are available at the required spatiotemporal scale.