Estimation of Land Surface Albedo from MODIS and VIIRS Data: A Multi-Sensor Strategy Based on the Direct Estimation Algorithm and Statistical-Based Temporal Filter

: Land surface albedo is an important variable for Earth’s radiation and energy budget. Over the past decades, many surface albedo products have been derived from a variety of remote sensing data. However, the estimation accuracy, temporal resolution, and temporal continuity of these datasets still need to be improved. We developed a multi-sensor strategy (MSS) based on the direct-estimation algorithm (DEA) and Statistical-Based Temporal Filter (STF) to improve the quality of land surface albedo datasets. The moderate-resolution imaging spectroradiometer (MODIS) data onboard Terra and Aqua and the Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi-National Polar-orbiting Partnership (NPP) were used as multi-sensor data. The MCD43A3 product and in situ measurements from the Surface Radiation Budget Network (SURFRAD) and FLUXNET sites were employed for validation and comparison. The results showed that the proposed MSS method signiﬁcantly improved the temporal continuity and estimation accuracy during the snow-covered period, which was more consistent with the measurements of SURFRAD (R = 0.9498, root mean square error (RMSE) = 0.0387, and bias = − 0.0017) and FLUXNET (R = 0.9421, RMSE = 0.0330, and bias = 0.0002) sites. Moreover, this is a promising method to generate long-term, spatiotemporal continuous land surface albedo datasets with high temporal resolution.


Introduction
Land surface albedo, defined as a ratio of the radiation reflected from the land surface and the incident solar radiation (0.3-5.0 µm) [1], is a vital parameter for the Earth's radiation and energy budget, which has been used widely in studies of global climate change, hydrology, agriculture, and weather forecasting [2,3]. Over the past few decades, the land surface albedo has changed because of global climate change and human activities [4,5]. In addition, the land surface could change rapidly with many geographic processes, including snowfall and snowmelt, desertification, deforestation, wildfire, and crop harvest [6][7][8]. Therefore, it is important to monitor the dynamics of the land surface albedo with different instruments [3,9]. Although the land surface albedo can be directly measured using ground and tower-based pyranometers, the global or regional land surface albedo cannot be accurately represented because of the limited footprints of measuring instruments and high spatial heterogeneity of land surface albedo [10]. Thus, the ability to monitor the spatiotemporal variations with remote sensing data offers significant advantages [3].
Estimating land surface albedo from multi-sensor with an adequate algorithm is a potential solution to address these issues. Solomon et al. [36] evaluated the performance of estimating bidirectional reflectance distribution function (BRDF) and land surface albedo with a combination of MODIS data onboard Terra and Aqua. The results showed that the data quality, accuracy, and robustness were improved by combining of two sensors' data. Liu et al. [35] proposed a multi-angular and multi-spectral BRDF model for estimating broadband land surface albedo with multi-sensor data. Compared with the mono-sensor method, the higher temporal resolution albedo with the same or even higher-accuracy estimations could be obtained by combining MODIS and Advanced Very High Resolution Radiometer (AVHRR) data. Muller et al. [24] developed an algorithm for generating the Globalbedo dataset with multi-sensor data, such as MERIS, Along Track Scanning Radiometers (ATSR) and Satellite Pour l'Observation de la Terre (SPOT) Vegetation. Wen et al. [37] built a multi-sensor combined BRDF inversion (MCBI) model for retrieving BRDF and albedo on a small timescale. The estimation results showed that the increase of observation number and diversity of angular sampling could provide additional information for obtaining a more robust land surface BRDF retrieval.
In this study, we developed a multi-sensor strategy (MSS) for improving the temporal resolution and temporal continuity of land surface albedo from MODIS and VIIRS data. The paper is organized as follows: Firstly, we present the overall framework of the MSS method, the satellite and in situ measurement data used in this study in Section 2. Secondly, we provide the comparison and validation results of our proposed method in Section 3. Finally, we provide the primary conclusion of this study in Section 4.

Overall Framework
This study provided an algorithm for generating spatiotemporal continuous land surface albedo datasets with high temporal resolution (e.g., daily). For this purpose, we employed the direct-estimation algorithm (DEA) [16,20,21], statistical-based temporal filter (STF) [38], and multi-sensor remote sensing data (e.g., MODIS and VIIRS) in this study. The flowchart (Figure 1) of the MSS method was categorized into the following four steps: (i) We first used the DEA to estimate land surface albedo from the multi-sensor remote sensing data; (ii) then, we used a weight-based fusion method to obtain the consistent land surface albedo dataset from multi-sensor data; (iii) we used the STF to fill the data gaps

Direct Estimation Algorithm
To obtain the high-temporal-resolution land surface albedo dataset, we employed the DEA approach in this study. The DEA approach enabled us to estimate broadband surface albedo directly from a mono-angular observation, which was suitable for monitoring the spatiotemporal dynamics of the land surface albedo. This method was applied to generate broadband surface albedo datasets from MODIS [20,21], VIIRS [16,18], and Landsat TM/ETM/OLI [31] data. In this study, we used an adapted DEA approach to estimate shortwave land surface albedo directly from the atmospheric corrected land surface reflectance. The main steps of the adapted DEA approach were as follows: First, we derived a training dataset of the land surface reflectance and the corresponding broadband surface albedo based on the BRDF database. The BRDF database used in this study was built based on the monthly POLDER-3 BRDF database (https://www.theia-land.fr/en/product/brdf-polder-3-project). It contained 13,879 globally distributed samples (a monthly database shared in 16 International Geosphere-Biosphere Programme (IGBP) classes) for representing the land surface reflectance anisotropy for different locations, seasons, and land cover types. For each sample, it can provide atmospheric corrected and multi-angular land surface reflectance based on the accumulation of one-month observations from different satellite overpasses and viewing angles [40]. The processing methods of data screening, quality control, classification, and band conversion were fully described in our former study [21]. Then, we established the empirical multiple linear relationships between the broadband surface albedo and atmospheric corrected land surface reflectance at different solar/view angular bins based

Direct Estimation Algorithm
To obtain the high-temporal-resolution land surface albedo dataset, we employed the DEA approach in this study. The DEA approach enabled us to estimate broadband surface albedo directly from a mono-angular observation, which was suitable for monitoring the spatiotemporal dynamics of the land surface albedo. This method was applied to generate broadband surface albedo datasets from MODIS [20,21], VIIRS [16,18], and Landsat TM/ETM/OLI [31] data. In this study, we used an adapted DEA approach to estimate shortwave land surface albedo directly from the atmospheric corrected land surface reflectance. The main steps of the adapted DEA approach were as follows: First, we derived a training dataset of the land surface reflectance and the corresponding broadband surface albedo based on the BRDF database. The BRDF database used in this study was built based on the monthly POLDER-3 BRDF database (https://www.theia-land.fr/en/product/brdf-polder-3-project). It contained 13,879 globally distributed samples (a monthly database shared in 16 International Geosphere-Biosphere Programme (IGBP) classes) for representing the land surface reflectance anisotropy for different locations, seasons, and land cover types. For each sample, it can provide atmospheric corrected and multi-angular land surface reflectance based on the accumulation of one-month observations from different satellite overpasses and viewing angles [40]. The processing methods of data screening, quality control, classification, and band conversion were fully described in our former study [21]. Then, we established the empirical multiple linear relationships between the broadband surface albedo and atmospheric corrected land surface reflectance at different solar/view angular bins based on the training dataset. Finally, the regression coefficients for different angular bins were saved in the direct estimation lookup table (LUT) files. When satellite data were available, we were able to estimate the shortwave surface albedo based on the LUT. The multiple linear relationships between the broadband Remote Sens. 2020, 12, 4131 4 of 18 white-sky and black-sky surface albedo and atmospheric corrected land surface albedo can be expressed as follows [20,21]: where α wsa is the broadband white-sky albedo, α bsa (θ bsa (k)) is the broadband black-sky albedo with a solar zenith angle (SZA) of θ bsa (k), which varies from 0 • to 80 • in increments of 4 • , ρ(θ s , θ v , ϕ, Λ i ) is the atmospheric corrected surface reflectance at band of Λ i with SZA of θ s , view zenith angle (VZA) of θ v , and relative azimuth angle (RAA) of ϕ; and m 0 , n 0 (k), m i , and n i (k) are the regression coefficients of multilinear regression functions. We divided the solar and view geometric space into angular bins according to our former study [20]: (1) The central SZA of each angular bin varied from 0 • to 80 • with increments of 2 • ; (2) the central VZA of each angular bin varied from 0 • to 64 • with increments of 2 • ; and (3) the central RAA of each angular bin varied from 0 • to 180 • with increments of 5 • . The details of the DEA approach can be found in our previous studies [20,21]. Note that the main difference between this study and former studies is that the input variable was the atmospheric corrected land surface reflectance instead of the top of atmosphere (TOA) reflectance.

Band Conversion and Fusion of Multi-Sensor Data
Previous studies have demonstrated that estimation accuracy and robustness can be improved by the combined use of multi-sensor data [36,37]. In this study, we employed multi-sensor data to obtain optimal estimations of land surface albedo. To incorporate the DEA and STF approaches, we used the band conversion and weight-based fusion method to handle the multi-sensor data.
Due to the differences in central wavelengths and spectral response functions of the bands, the VIIRS data cannot be used to estimate land surface albedo directly with the MODIS direct-estimation LUTs. To solve this issue, we used the band conversion method to convert the VIIRS data to MODIS data as follows [41]: where ρ a (Λ i ) is the MODIS reflectance at band of Λ i (i = 1~7), ρ b Λ j is the reflectance of other sensors at band of Λ j (j = 1~9), c 0 is the intercept term, and c j is the band conversion coefficient. Among these variables, the band conversion coefficients could be derived by a linear regression method based on the spectral library. Spectral reflectance data can be calculated from the spectral response function and the ground-measured spectral reflectance as follows: where ρ(Λ) is the spectral reflectance at band of Λ, f (λ k ) and ρ(λ k ) are the spectral response function and spectral reflectance data from spectral library at wavelength λ k , respectively. In this study, we used 1232 samples of spectral reflectance data from the U.S. Geological Survey (USGS) spectral library to establish the relationships between MODIS and VIIRS data. The spectral reflectance data contained 286 samples for vegetation, 209 samples for soil, 1276 samples for rock, 370 samples for water and snow/ice, and 259 samples for other materials. Among these samples, we used 70% for the band conversion, and 30% to evaluate the uncertainties of band conversion.
To get optimal estimations from multi-sensor data, the difference in estimation uncertainty for different sensors should be considered. The accuracy of the multi-sensor data may change dynamically because of various external factors (e.g., the accuracy of calibration, quality of data, and data processing algorithms). Thus, to generate a long-term, spatiotemporal continuous land surface albedo dataset, we applied a data fusion method to enhance the continuity and consistency of the satellite data. In this study, we combined the multi-sensor albedo data (MODIS and VIIRS) for a one-year period (from Julian day 1 to 365) by assigning different weights [42]. The fused albedo can be expressed as follows: where α F stands for the fused albedo; α 1 and α 2 denote the estimated albedo from MODIS and VIIRS data, respectively; and w 1 and w 2 are the weights for MODIS and VIIRS data, respectively ( Figure 2), which were derived based on the estimation uncertainties of MODIS and VIIRS data. Then, w 1 can be expressed as the following function [42]: where x is the Julian day of the year and a, b, c, and d are the parameters used to ensure the shape of this function (a = 11, b = 131, c = 235, and d = 355). Then, w 2 can be expressed as follows: Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 19 370 samples for water and snow/ice, and 259 samples for other materials. Among these samples, we used 70% for the band conversion, and 30% to evaluate the uncertainties of band conversion.
To get optimal estimations from multi-sensor data, the difference in estimation uncertainty for different sensors should be considered. The accuracy of the multi-sensor data may change dynamically because of various external factors (e.g., the accuracy of calibration, quality of data, and data processing algorithms). Thus, to generate a long-term, spatiotemporal continuous land surface albedo dataset, we applied a data fusion method to enhance the continuity and consistency of the satellite data. In this study, we combined the multi-sensor albedo data (MODIS and VIIRS) for a one-year period (from Julian day 1 to 365) by assigning different weights [42]. The fused albedo can be expressed as follows: where stands for the fused albedo; and denote the estimated albedo from MODIS and VIIRS data, respectively; and w1 and w2 are the weights for MODIS and VIIRS data, respectively ( Figure 2), which were derived based on the estimation uncertainties of MODIS and VIIRS data. Then, can be expressed as the following function [42]: where x is the Julian day of the year and a, b, c, and d are the parameters used to ensure the shape of this function (a = 11, b = 131, c = 235, and d = 355). Then, can be expressed as follows: Figure 2. Weights for the combination of multi-sensor albedo data.

Statistical-Based Temporal Filter
To obtain a spatiotemporal continuous land surface albedo dataset, we also used the STF approach in this study. Most of the currently used land surface albedo products suffer from issues of spatiotemporal discontinuity and data gaps resulting from cloud obscuration and contamination. Therefore, we carried out the gap filling, or filtering method, to obtain a spatiotemporal continuous dataset.

Statistical-Based Temporal Filter
To obtain a spatiotemporal continuous land surface albedo dataset, we also used the STF approach in this study. Most of the currently used land surface albedo products suffer from issues of spatiotemporal discontinuity and data gaps resulting from cloud obscuration and contamination. Therefore, we carried out the gap filling, or filtering method, to obtain a spatiotemporal continuous dataset.
In the STF approach, we assumed that the surface albedo had a temporal dependency with the albedo of neighboring days. Thus, it was reasonable to assume that the albedo α k on the kth day had a linear relationship with the albedo α k+∆k on the (k+∆k)th day [38] as follows: where a ∆k and b ∆k are the regression coefficients, and e ∆k is the regression residual. We built a priori surface albedo climatology based on an average of the multi-year MODIS land surface albedo product Collection 6 (MCD43A3) [15]. Considering the differences between the snow-free and snow-covered conditions, we treated the priori surface albedo climatology in these two scenarios differently. Then, we derived a priori STF LUT, and reconstructed the albedoα k on the kth day according to the following equation [38]: where ε 2 is the variance of the estimated albedo. We calculated the STF filter parameters including the multi-year albedo mean µ k , standard deviation σ k (k = 1, 2 . . . , 365), and root mean square error (RMSE) ζ ∆k (∆k 0 and ∆k = −5, −4 . . . , 5) from the surface albedo climatology. We set the temporal window in the STF as five days before or after the kth day.

Satellite Data
In this study, we used the daily atmospheric corrected reflectance data of MODIS on board Terra (MOD09GA) and Aqua (MYD09GA), and VIIRS on board the Suomi-National Polar-Orbiting Partnership (Suomi-NPP) (VNP09GA) as multi-sensor data to estimate the shortwave land surface albedo. The MOD/MYD/VNP09GA datasets provided the atmospheric corrected MODIS/VIIRS daily surface reflectance with a spatial resolution of 1 km, and a temporal resolution of 1 day. All these data were downloaded from the Earthdata Search website (https://search.earthdata.nasa.gov). We designed the LUT of DEA based on the MODIS data, the VIIRS data were converted to MODIS bands based on Equation (3). The coefficients and RMSE values of band conversion are listed in Table 1. The RMSEs of the band conversion were acceptable at most of wavelengths, and the RMSEs at shortwave infrared bands were higher than those of visible or near infrared bands because of the differences in central wavelengths and spectral response functions of MODIS and VIIRS data.

Validation Data
We completed the validation with in situ measurements and MCD43A3 data. Because the footprints of satellite observations and the in situ measurements were quite different, it was important to evaluate the spatial representativeness of the in situ measured sites [43]. In this study, we selected three sites from the Surface Radiation Budget (SURFRAD) [44] Network and 11 sites from FLUXNET [45] for validation, which had a high spatial representativeness in the point-to-pixel comparison [46,47].

MODIS Albedo Product
In this study, we used the MCD43A3 product as a reference for comparison with our algorithm. In the comparison procedure, we obtained the blue-sky broadband albedo from a linear combination of the black-sky and white-sky albedo [48,49] as follows: where α(θ s ) is the blue-sky albedo at the local solar noon, D(θ s ) is the diffuse skylight fraction at solar zenith angle θ s , α bs and α ws are the black-sky albedo and white-sky albedo, respectively.

In Situ Measurements
The SURFRAD Network is made up of eight surface radiation and meteorological measurement sites across climatologically diverse regions of the United States [44]. In this study, we used three SURFRAD sites for validation as shown in Table 2, for which spatial representativeness was evaluated with semi-variograms (sill, range, and nugget effects) at different times of the year [14,17,50]. FLUXNET is a global network of micrometeorological tower sites that use eddy covariance methods to measure the exchanges of carbon dioxide, water vapor, and energy between terrestrial ecosystems and the atmosphere [45]. In this study, we used the measurements from 10 sites to validate the land surface albedo derived using the MSS method. The spatial distribution of these sites are shown in Figure 3, and site information is listed in Table 3. The spatial representativeness of these FLUXNET sites was evaluated by Cescatti et al. [46].

Estimations of Land Surface Albedo
In this study, we compared the estimation results of the MSS method with the in situ measurements with high spatial representiveness. The time series of albedo estimated by MSS method, MCD43A3, and in situ measurements at TBL (40. Figure 4. It is evident that the albedo estimated by the MSS method was gap-free and temporally continous, which provided better representation of the temporal dynamics of the land surface albedo at these sites. In this study, we used the albedo value of 0.4 as a threshold [31]. When the albedo value was greater than 0.4, we considered the pixel to be snow-covered; otherwise, we considered the pixel to be snow-free. Under the snow-free conditions, the performance of the albedo estimated by MSS was consistent with the MCD43A3 product in most cases (Figure 4b). During the wintertime, when snowfall or snowmelt events occurred, the temporal dynamics of the albedo were well captured by the albedo derived from the MSS method, including the ephemeral snowfall and snow-covered events (Figure 4a,c).

Validation with SURFRAD Sites
The comparison results between the SURFRAD measurements, albedo estimated by the MSS method, and MCD43 albedo in 2015 are shown in Figure 5. It is evident that the correlation coefficient R of the MSS method (R = 0.9498 for all conditions, and R = 0.7421 for the snow-free condition) was higher than that of MCD43A3 (R = 0.8460 for all conditions, and R = 0.5318 for the snow-free condition), and the RMSE and bias of the MSS method (RMSE = 0.0387, bias = −0.0017 for all conditions; RMSE = 0.0253, bias = 0.0007 for the snow-free condition) were relatively smaller than those for MCD43A3 (RMSE = 0.0656, bias = −0.0040 for all conditions; RMSE = 0.0326, bias = 0.0042 for the snow-free condition). As shown in the figure, the estimation uncertainty of the MCD43A3 product was relatively larger than those derived from the MSS approach for all conditions, whereas the estimation uncertainties of these two methods under snow-free conditions were similar, which suggested that the MSS approach provided more accurate estimations than the MCD43A3 product during the snow-covered period. This accuracy was mainly due to the fact that the MCD43A3 product was developed by fitting and integrating a semi-empirical kernel-driven model with accumulated multi-angular observations within a short period (e.g., 16 days) [14,50]. The estimation results would be acceptable when the land surface changed gradually, however, when the land surface changed rapidly (e.g., snowfall/snowmelt, crop harvest, and forest fire), the temporal variation of albedo could not be well captured by the MCD43A3 product. The MSS approach, however, could retain the abrupt increasing and decreasing signals and provided much more accurate estimations when the land surface albedo changed rapidly [20,21].   spatial footprints between the remote sensing data and in situ measurements. We used in situ measurements as a reference to assess the estimation accuracy of the results derived from remote sensing data in this study. The temporal variation patterns of the land surface albedo at the pixel scale and ground-measured scale were quite different, however, because of the high spatial heterogeneity [51], which suggested that more fieldwork and investigations are needed to improve the estimation and validation of land surface albedo.  The comparison results of the SURFRAD measurements, the albedo estimated by the MSS method, and MCD43A3 albedo at each site for all conditions and the snow-free condition are shown in Figures 6  and 7, respectively. It is evident that the MSS method had a better consistency with the SURFRAD measurements for all three sites during the snow-covered period. In addition, we found that the lower albedo values at FPK site (Figure 4b) cannot be well represented by neither of the MCD43A3 product and MSS approach. Note that the figures (Figures 6 and 7) show a horizontal spread scatter during the snow-free condition. This scatter was mainly caused by the mismatch of spatial footprints between the remote sensing data and in situ measurements. We used in situ measurements as a reference to assess the estimation accuracy of the results derived from remote sensing data in this study. The temporal variation patterns of the land surface albedo at the pixel scale and ground-measured scale were quite different, however, because of the high spatial heterogeneity [51], which suggested that more fieldwork and investigations are needed to improve the estimation and validation of land surface albedo.

Validation with FLUXNET Sites
The comparison results between the FLUXNET measurements, albedo estimated by the MSS method, and MCD43 albedo are shown in Figure 8. Similar to the validation results of SURFRAD, the results showed that the correlation coefficient R of the MSS method (R = 0.9421 for all conditions, and R = 0.8363 for the snow-free condition) was relatively higher than that of MCD43A3 (R = 0.8980 for all conditions, and R = 0.8043 for the snow-free condition), and the RMSE and bias of the MSS method (RMSE = 0.0330, bias = 0.0002 for all conditions, RMSE = 0.0258, bias = 0.0028 for the snow-free condition) were relatively smaller than those of MCD43A3 (RMSE = 0.0429, bias = −0.0033 for all conditions, RMSE = 0.0274, bias = 0.0004 for the snow-free condition), which indicated that the albedo estimated by the MSS method was more consistent with the FLUXNET measurements over different land cover types.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 19 The comparison results between the FLUXNET measurements, albedo estimated by the MSS method, and MCD43 albedo are shown in Figure 8. Similar to the validation results of SURFRAD, the results showed that the correlation coefficient R of the MSS method (R = 0.9421 for all conditions, and R = 0.8363 for the snow-free condition) was relatively higher than that of MCD43A3 (R = 0.8980 for all conditions, and R = 0.8043 for the snow-free condition), and the RMSE and bias of the MSS method (RMSE = 0.0330, bias = 0.0002 for all conditions, RMSE = 0.0258, bias = 0.0028 for the snow-free condition) were relatively smaller than those of MCD43A3 (RMSE = 0.0429, bias = −0.0033 for all conditions, RMSE = 0.0274, bias = 0.0004 for the snow-free condition), which indicated that the albedo estimated by the MSS method was more consistent with the FLUXNET measurements over different land cover types.
The comparison results between the RMSEs of the albedo estimated by the MSS method and MCD43A3 at each FLUXNET site are shown in Figure 9. According to these results, the RMSEs of the albedo estimated by the MSS method were relatively smaller than those of the MCD43A3 product at most of sites. Notably, differences in the RMSEs for FR-Pue and US-SRM sites were not significant because they were not seriously affected by snowfall events.  The comparison results between the RMSEs of the albedo estimated by the MSS method and MCD43A3 at each FLUXNET site are shown in Figure 9. According to these results, the RMSEs of the albedo estimated by the MSS method were relatively smaller than those of the MCD43A3 product at most of sites. Notably, differences in the RMSEs for FR-Pue and US-SRM sites were not significant because they were not seriously affected by snowfall events.

Assessment of the Temporal Continuity
The data gaps of the albedo dataset could be filled efficiently using the multi-sensor and STF approach. A comparison of the number of effective days per year for different sensor combinations and methods is given in Table 4. It is evident that the numbers of effective days per year increased with the numbers of sensors. If we used the MODIS data onboard Terra (MOD) and Aqua (MYD) together (MOD+MYD), the number of valid albedo values increased slightly compared with the results using MOD data only. If we also used the VIIRS data (MOD+MYD+VIIRS), the number of effective days per year increased significantly because the VIIRS data provided more valid observations for CZ-BK1 and US-Los sites. If, however, the STF method was carried out, all of the data gaps were filled, we could derive a spatiotemporally continuous dataset using this method.

Comparison of the Estimation Results Derived by DEA and MSS Approaches
The validation results with SURFRAD and FLUXNET sites for the DEA and MSS approaches are shown in Table 5. The correlation coefficients and RMSEs of the DEA approach with MODIS, VIIRS, and the combination using MODIS and VIIRS data (MODIS + VIIRS) were similar. If we did

Comparison of the Estimation Results Derived by DEA and MSS Approaches
The validation results with SURFRAD and FLUXNET sites for the DEA and MSS approaches are shown in Table 5. The correlation coefficients and RMSEs of the DEA approach with MODIS, VIIRS, and the combination using MODIS and VIIRS data (MODIS + VIIRS) were similar. If we did not consider gap-filled data when no observation was available, the MSS approach provided better estimations with a higher R and a smaller RMSE (R = 0.9025, RMSE = 0.0351), which suggested that the MSS approach provided an improvement in estimation accuracy. If we considered the gap-filled MSS data, the results had a higher R and a larger RMSE (R = 0.9402, RMSE = 0.0385), which suggested that the gap-filled data introduced uncertainties to the estimation results. Because the MSS approach provided temporally continuous and gap-free albedo estimations, its relatively lower estimation accuracy was acceptable when no observation was available.

Discussion
The comparison and validation results showed that the estimation accuracy and temporal continuity of land surface albedo were significantly improved by the proposed MSS method. In this study, we used the MODIS data onboard Terra and Aqua and VIIRS data as an example to test the accuracy of this algorithm. Other remote sensing data, including AVHRR, MISR, MERIS, POLDER, sensors on board Chinese Fengyun(FY)-3 satellites, Landsat TM/ETM/OLI, Sentinel, sensors on board Chinese Huanjing (HJ) and Gaofen(GF) series satellites, and sensors onboard geostationary satellites could also be used for this purpose. In addition, a potential opportunity exists to improve spatial and temporal resolution (e.g., 30 to 250 m, 1 h to 1 day) based on the spatiotemporal data fusion method [52] for agricultural and hydrology applications.
The DEA approaches used in this study and former studies [20,21] had one main difference that is, the input data we used the atmospheric-corrected land surface reflectance instead of TOA reflectance. Because the calibrated TOA radiance and reflectance data (e.g., MOD/VNP02 products) were original swath data with duplicate observations, the data size was much larger than those from atmospheric-corrected reflectance data (MOD/MYD/VNP09GA products) and needed to be processed to generate the standard tile products. Thus, a change in the input data can reduce the huge computations required for the MODIS and VIIRS swath data, which was more suitable and applicable for generating a long-term land surface albedo dataset.
The STF approach used in this study was not equal to the simple filtering method. Better estimation results could be obtained by the STF approach than by simply filtering [38]. When no available observation was available, however, the albedo estimation was much closer to the albedo climatology, which resulted in large errors when extreme climatic events occurred. The only way to address this issue is by adding more valid observations from multi-sensor data, the importance of which has been demonstrated by this study as well as by former studies [35][36][37].
In the MSS approach, the estimation uncertainties can be introduced by the quality of input multi-sensor data, band conversion, DEA approach, and STF method. Because the observations from different sensors (e.g., MODIS and VIIRS) had different observation configurations (i.e., central wavelengths, spectral response functions, and viewing zenith/azimuth angles) and processing methods (i.e., calibration, cloud detection, and atmospheric correction) [17], the estimation uncertainties varied with sensors. Additionally, the band conversion also introduced estimation uncertainties, especially for the shortwave bands. The estimation uncertainties of multi-sensor data and band conversion could be handled by the weight-based fusion method, in which the weights were derived based on the evaluation of estimation uncertainty. These comparison and validation results showed that the DEA approach provided accurate estimation of land surface albedo, especially when the land surface changed rapidly. The validation results demonstrated that the gap-filling method introduced uncertainties in albedo estimation. This was not a significant shortcoming, however, because the gap-filling method could generate temporally continuous and gap-free datasets. In fact, the MSS framework could be applied to other variables (e.g., leaf area index (LAI) and fractional cover of vegetation). If the combined DEA and STF approaches were adapted to other variables, it also could be used to generate various global climate datasets with high temporal resolution and continuity.

Conclusions
In this study, we proposed an MSS method to estimate of the land surface albedo from multi-sensor data (e.g., MODIS and VIIRS). To obtain high temporal resolution and a temporally continuous albedo dataset, we employed the DEA and STF methods in this study. We compared the SURFRAD, FLUXNET, and MCD43A3 data to validate the efficiency and temporal continuity of the MSS method. The main findings of this study are as follows: (1) We obtained more accurate estimations of land surface albedo during snow-covered period using the proposed MSS method. The albedo estimated by the MSS method was consistent with the measurements of SURFRAD (R = 0.9498, RMSE = 0.0387, and bias = −0.0017) and FLUXNET (R = 0.9421, RMSE = 0.0330, and bias = 0.0002) sites. (2) The temporal continuity of the land surface albedo dataset was significantly improved by employing the multi-sensor data and STF. We found that the number of effective days per year increased with the number of valid satellite observations from multi-sensor data, and temporally continuous, gap-free land surface albedo datasets can be obtained using the proposed MSS method. (3) By incorporating the DEA and STF approaches, the MSS method could be used to generate long-term, spatiotemporal continuous land surface albedo datasets with high temporal resolution (e.g., daily). The results of this study demonstrated that this is a promising method for generating global climate datasets with high temporal resolution and spatiotemporal continuity.