A statistics-based temporal filter algorithm to map spatiotemporally continuous shortwave albedo from MODIS data

Land-surface albedo plays a critical role in the earth’s radiant energy budget studies. Satellite remote sensing provides an effective approach to acquire regional and global albedo observations. Owing to cloud coverage, seasonal snow and sensor malfunctions, spatiotemporally continuous albedo datasets are often inaccessible. The Global LAnd Surface Satellite (GLASS) project aims at providing a suite of key land surface parameter datasets with high temporal resolution and high accuracy for a global change study. The GLASS preliminary albedo datasets are global daily land-surface albedo generated by an angular bin algorithm (Qu et al., 2013). Like other products, the GLASS preliminary albedo datasets are affected by large areas of missing data; beside, sharp fluctuations exist in the time series of the GLASS preliminary albedo due to data noise and algorithm uncertainties. Based on the Bayesian theory, a statisticsbased temporal filter (STF) algorithm is proposed in this paper to fill data gaps, smooth albedo time series, and generate the GLASS final albedo product. The results of the STF algorithm are smooth and gapless albedo time series, with uncertainty estimations. The performance of the STF method was tested on one tile (H25V05) and three ground stations. Results show that the STF method has greatly improved the integrity and smoothness of the GLASS final albedo product. Seasonal trends in albedo are well depicted by the GLASS final albedo product. Compared with MODerate resolution Imaging Spectroradiometer (MODIS) product, the GLASS final albedo product has a higher temporal resolution and more competence in capturing the surface albedo variations. It is recommended that the quality flag should be always checked before using the GLASS final albedo product.


Introduction
Land-surface albedo refers to the ratio of reflected to incoming solar radiation at the earth's surface over the solar spectral domain (Dickinson, 1983).As one of the fundamental forcing parameters in climate models, surface albedo plays an important role in the earth's radiant energy budget (Barnes andRoy, 2008, 2010;Liang et al., 2010;Manalo-Smith et al., 1998).Spatio-temporal variability in albedo is often associated with environmental change and human activities as well (Bsaibes et al., 2009;Dirmeyer and Shukla, 1994;Jin and Roy, 2005;Moritz et al., 2002).According to the Global Climate Observing System (GCOS), spatiotemporally continuous albedo products with 5 % relative accuracy, 1 km spatial resolution and 1 day temporal resolution are required by climate studies (GCOS, 2011).Several factors (e.g.cloud coverage, seasonal snow cover and sensor malfunction), however, lead to large areas of data gaps in current albedo products.Take the MODerate resolution Imaging Spectro-radiometer (MODIS) standard albedo product (MCD43B3) as an example, 20 % to 40 % of the global land pixels miss valid shortwave black-sky albedo (Table 1).A variety of strategies have been proposed to fill in the gaps of current albedo datasets.One approach is to improve the inputs of the albedo retrieving algorithm.Based on the standard retrieving technique of MODIS nadir BRDFadjusted reflectance products, an adapted method was employed to improve the integrity of reflectance series (Ju et al., 2010).A Kalman filter was used to improve the completeness of BRDF coefficients series (Samain et al., 2008).In the retrieving algorithm of the GlobAbledo albedo product, a regularization method was utilized to generate daily kernel coefficients (Quaife and Lewis, 2010;ESA GlobAbedo ATBD, 2012).Beside the above mentioned strategy, post-processing current albedo products is another solution for filling albedo gaps.An ecosystem curving fitting (ECF) method was performed on the MOD43B3 product to generate spatially complete albedo datasets (Moody et al., 2005).In the context of the ECF method, missing data are interpolated by using pixel-level and regional albedo climatology curves.However, Fang et al.'s (2007Fang et al.'s ( , 2008) ) studies demonstrated that multiyear average is more suitable than regional average for characterizing a pixel's profile.Fang et al. (2007Fang et al. ( , 2008) ) designed a temporal-spatial filter (TSF) to generate gapless albedo and leaf area index products.In the TSF method, the climatology of each pixel is first determined from multiyear observations.Then, both multiyear averages and neighboring observations are averaged with weights to replace missing values.
The objective of this study is to map spatiotemporally continuous daily shortwave albedo for global land surface.The paper is organized as follows: Sect. 2 introduces the GLASS preliminary albedo products.The statistics-based temporal filter algorithm is described in Sect.3. Section 4 analyzes the results and Sect. 5 concludes this study.

The GLASS preliminary albedo products
The Global LAnd Surface Satellite (GLASS) is a project sponsored by the Chinese "State Program for High-Tech Research and Development" program, which aims at providing a suite of key land surface parameter datasets with high temporal resolution and high accuracy for a global change study (Liang et al., 2013).The GLASS albedo product is generated in two steps (Liu et al., 2013).In the first step, albedo is retrieved from remotely sensed reflectance data by the inversion of a radiative transfer model, or its simplification.The resultant albedo retrieved in the first step is called the GLASS preliminary product.In the second step, the GLASS preliminary products are post-processed into intermediate and final products.The GLASS preliminary albedo products are composed of four kinds of global daily 1 km albedo datasets generated from two different algorithms (Angular Bin algorithms, AB1 and AB2).Based on Liang's method, AB algorithms divide the solar/view-geometry hemisphere into many angular bins, and then establish a statistical regression model that links directional reflectance in each angular bin to shortwave albedo (Liang, 2003;Liang et al., 2005;Qu et al., 2013).The AB1 algorithm employs the Aqua or Terra MODIS land surface reflectance products as inputs while the AB2 algorithm uses the MODIS top-of-atmosphere reflectance data as inputs.The GLASS preliminary albedo products are named as GLASS02A2x, where x = 1 is for the AB1+MOD09GA combination, 2 is for AB1+MYD09GA, 3 is for AB2+MOD02 and 4 is for AB2+MYD02.Data contained in the GLASS preliminary albedo products are shortwave black-sky albedo (BSA) at local solar noon, shortwave white-sky albedo (WSA) and quality assessment (QA) flag.The uncertainty of the GLASS albedo is encoded as a four-bit word in the QA flag.Based on the GLASS preliminary products, the GLASS intermediate albedo products GLASS02A0x are obtained by averaging the good-quality GLASS preliminary albedo in a 16 day temporal window on an 8 day cycle.Details of the GLASS preliminary and intermediate albedo products are listed in Table 2.
The combination of AB algorithms and MODIS data has greatly improved the temporal resolution of the GLASS preliminary products.However, there are still limitations in the GLASS preliminary albedo products (Liu et al., 2011).First, the GLASS preliminary albedo products are affected by frequent data gaps which are mainly caused by cloud coverage and seasonal snow.Second, the GLASS preliminary albedo series fluctuate sharply because of data noises and algorithm uncertainties.Finally, the four GLASS preliminary albedo products need to be merged into one albedo product for the convenience of users.Therefore, a post-processing algorithm is needed to generate a spatiotemporally continuous GLASS final albedo product.

Methodology
As illustrated in Fig. 1, the procedure of our STF algorithm comprises three parts.First, the a priori statistics of global surface albedo is calculated.Then, the filter parameters are derived from the a priori statistics.Finally, the GLASS final product is generated by merging and filtering the four GLASS preliminary albedo products.

The formula of the statistics-based temporal filter
Based on the temporal correlation of albedo measurements in neighboring days, it is reasonable to assume that the true albedo α k on the kth day is linearly correlated with the true albedo α k+ k on the (k + k)th day: where both a k and b k are regression model coefficients.
The regression residual e k is assumed to be Gaussian distributed with a zero mean and a ζ k standard deviation (see Sect. 3.2 for the calculation of a k , b k and ζ k ).Then given the GLASS02A2x retrieval α * k+ k and its uncertainty η k+ k (estimated from the GLASS02A2x QA data, see Sect. 2) on the (k + k)th day, the probability density function (PDF) of the true albedo α k on the kth day is in this form: where N denotes the normal distribution.In other words, the predicted albedo on the kth day is a k α * k+ k + b k , and the corresponding prediction uncertainty consists of two components: the regression model uncertainty ζ 2 k and the propagated observational uncertainty a 2 k η 2 k+ k .Specifically, the PDF of α k is N (α * k , η 2 k ) when k equals to zero.Moreover, the a priori PDF of α k is N(µ k |σ 2 k ) characterized by the multiyear mean µ k and standard deviation σ k (see Sect. 3.2 for the calculation of µ k and σ k ).
Assume that the regression residual e k and retrieval uncertainty η k+ k ( k = −K, . . ., K) are independent of each other.Then, given the surrounding GLASS02A2x observations α * k + k and the a priori multiyear mean µ k , the joint PDF of α k can be expressed as Consequently, the maximum likelihood estimate of α k takes the form: where c is the corresponding prediction variance: Equation ( 4) shows that the STF method, similar to the TSF algorithm proposed by Fang et al. (2007Fang et al. ( , 2008)), is essentially a weighted average of statistical and observational values.The STF method has some advantages over the TSF methods.First, the weights in the STF formula depend on both the temporal correlation and the observational errors of the albedos in neighboring days.In contrast, the weights in the TSF method are set to be an empirical function of temporal distance and take no account of observational errors.According to Eq. ( 4), observations with better correlation and less noise contribute more to the filtered result.Second, the STF method can improve both the integrity and the smoothness of the GLASS preliminary albedo series.In the context of the STF method, the observational albedo on the kth day is smoothed if it is valid; or if not, the missing value is filled.Third, the a priori mean and standard deviation are employed to fill in data gaps in case of the absence of all surrounding observations.Finally, with Eq. ( 5), the STF method provides an assessment of the uncertainty of filtered results.

The a priori knowledge database of global surface albedo
The a priori knowledge is of importance in retrieving BRDF/albedo parameters.For example, the retrievals of BDRF/albedo have been greatly improved by adding the a priori knowledge into the kernel-driven inversion model (Li et al., 2001).The back-up algorithm of MODIS BRDF/albedo products employ the a priori BRDF information to address the limited angular sampling issue (Schaaf et al., 2002).In fact, both the regional climatology in the ECF method and the multiyear average in the TSF method are a form of a priori information as well.In our study, the a priori statistics of global surface albedo are calculated from the multiyear MCD43B3 product (from 2000 to 2010).The MCD43 series of albedo product have been generated from MODIS data and distributed by the MODLAND team of NASA since 2000.The highest spatial resolution of the MCD43 albedo product is 500 m, but here the 1kmresolution MCD43B3 product was used to be consistent with the GLASS products.The MCD43 series of albedo product has been extensively validated and recognized as one of the most widely used and best quality global albedo datasets.The MCD43B3 albedo product is the basis of the statistics which derive the STF filter parameters (a k , b k , ζ k , µ k and σ k ) in Eq. ( 4).The a priori statistics derived from the MCD43B3 product include: (1) multiyear albedo mean µ k and standard deviation σ k (k = 1, 9 . . .361), and (2) correlation coefficient ρ k ( k = −32, −24 . . .32) for the albedo in two days with a lag k.In consideration of storage consumption, the a priori knowledge database is derived from the 1 km MCD43B3 data and then aggregated into 5 km resolution.Since part of the MCD43B3 product is contaminated by cloud and takes invalid value, the 5 km aggregated statistics are more stable than the 1 km statistics.When applying the 5 km database to a 1 km GLASS pixel, the database record that is closest in location is used.Some regions, such as tropical areas covered by persistent cloud or polar areas affected by polar night, may miss multiyear statistical means.The following strategy is adopted to complete the a priori database: for polar areas, the average of statistical means on the two days before polar night is used to fill the gaps in the a priori mean; for other regions, albedo means for each IGBP classification at each 10 • latitudinal zone are calculated, and are used to fill the gaps in the a priori mean according to the pixel's IGBP classification.
The temporal resolution of the a priori knowledge database is 8 day.To obtain the daily filter parameters and statistics required by Eq. ( 4), the 8 day a priori statistics are interpolated to a daily resolution and then converted to the filter parameters.To this end, a polynomial interpolation technique is applied to the 8 day statistical mean and standard deviation, and an exponential interpolation function is used for the 8 day correlation coefficients.The curve fitting functions are as follows: where k is the time interval and λ i (i = 1, 2 . . .10) is the model fitting coefficient.Once the 8 day statistical mean, standard deviation and correlation coefficients are interpolated to a daily resolution, the daily statistics are converted to the filter parameters in Eq. ( 4): 4 Results and analysis

Consistency check of the STF inputs
The inputs of the STF algorithm, the four kinds of GLASS02A2x (x = 1, 2, 3, 4) albedo products, may have systematic discrepancies.Therefore, it is necessary to examine the consistency among the GLASS02A2x products before implementing the STF algorithm.Wide field measurements are the best way to examine the consistency between different products.However, such comprehensive datasets are inaccessible since field observations usually cannot represent pixel-level albedo.Therefore, to examine the consistency between different products by using a reference product is an alternative approach.The MODIS standard albedo product (MCD43B3) is employed to examine the consistency between GLASS02A2x products over different land covers (see Table 3 for the classification criteria).
To avoid the influence of spatial and temporal discrepancies between different products, spatiotemporally homogeneous pixels are first selected from the 8 day GLASS02A0x and MCD43B3 products.For a pixel centered in a 3 × 3 spatial window on the kth day, it is considered to be spatiotemporally homogeneous if 1. all the pixels in a 3 × 3 spatial window from day (k − 8) to day (k + 8) have the same land cover; 2. spatial relative differences between the albedo α 0 of center pixel and that of the surrounding pixels α i (i = 1, 2 . . .8) are less than 5 %; 3. temporal relative differences between the average albedo on the kth day and the average albedo on the neighboring days are less than 5 %.
Once the spatiotemporally homogeneous pixels are selected, the statistics, including average, standard deviation, correlation coefficient, bias and root mean square error (RMSE) are computed at every 10 • latitudinal zone.The statistical results provide some useful information about the consistency among albedo products (Table 4).First, the high correlation coefficients indicate a strong correlation between the GLASS02A0x products and the MCD43B3 product.Second, small biases and RMSEs indicate a good agreement between the GLASS02A0x products and the MCD43B3 product as well.Therefore, the systematic discrepancies among various GLASS02A0x products are considered to be negligible.

Analysis of the STF results
The STF technique was performed to post-process the global GLASS02A2x products and the resultant albedo dataset is named as GLASS02A06.Like the GLASS preliminary albedo products, there is a quality assessment flag in the GLASS02A06 product.Figure 2a and b are the comparison between the GLASS02A22 and the GLASS02A06 black-sky albedo of tile H25V05 in the MODIS Sinusoid grid (from day 209-361, 2008, every 8 days).This area is located on the Tibet Plateau and is covered by grassland, bare soil, lake and snow.It can be seen from Fig. 2a that large number of data gaps (gray color) can be found in this area owing to cloud and seasonal snow coverage.After the STF algorithm is implemented, the completeness of the albedo map has been greatly improved (Fig. 2b).The spatial distribution of different land covers is well reflected in the GLASS02A06 albedo maps.Take the albedo map of day 209 as an example, bare soil areas are with red color, grasslands areas are with green color, small lakes are with dark blue color and snow-covered mountain areas are with yellow color.Seasonal variations in albedo The field albedo measurements at three automatic weather stations are selected to evaluate the performance of the STF algorithm.These sites are BJ (31.3687 • N, 91.8987 • E), D105 (33.0643 • N, 91.9426 • E) and MS3478 (31.9262 • N, 91.7147 • E) with homogeneous grassland.Figure 3 depicts the time series of field measurements, the MODIS albedo (MCD43B3), the GLASS preliminary (GLASS02A22) and final (GLASS02A06) albedo at these stations for 2008.The a priori statistical mean and standard deviation series are also plotted in Fig. 3. a ρ is the correlation coefficient between the GLASS02A0x (x = 1, 2, 3 and 4) product and the MCD43B3 product.b b is the bias between the GLASS02A0x (x =1, 2, 3 and 4) product and the MCD43B3 product.c ε is the regression root mean square error between the GLASS02A0x and MCD43B3 products.
For all these stations, large statistical standard deviations are observed during the winter season.This large variation may be associated with the surface heterogeneity caused by seasonal snow.In contrast, statistical standard deviations are much smaller for the growing seasons.Seasonal variations or trends in albedo are well reflected by the time series of statistical mean, field measurements, MODIS, the GLASS preliminary and final results.The albedo time series remain stable for the growing season, increase with snow and decrease with snow melting.
Compared to other series, the fluctuations of the GLASS preliminary albedo series (GLASS02A22) are much larger.This may result from the uncertainties of the AB retrieving algorithm and the noise of input data.Besides, there are several data gaps in the GLASS preliminary albedo series due to the presence of cloud cover.After the STF technique is performed, both the integrity and smoothness of the albedo series (GLASS02A06) are improved greatly.However, this should be interpreted with caution because some true fluctuations in albedo may be treated as noisy data as well.For instance, both the field measurements and the GLASS preliminary product (GLASS02A06) at the MS3478 station show that there is an abrupt increase in surface albedo around day 141.However, this variation cannot be depicted by the GLASS final product since the STF method has smoothed the GLASS albedo time series.
It is interesting to know to what extent the GLASS final albedo is influenced by the a priori statistics.Results from these three stations (Fig. 3) indicate that the deviation of the GLASS final albedo from the statistical mean is very small in the non-winter seasons, but this deviation becomes very large in the winter season.The statistical standard deviation of albedo gives an explanation to this phenomenon: in the non-winter seasons, the variation of land surface albedo is very small and the observation agrees well with the statistical mean; although, in winter there is an intense albedo variation due to seasonal snow, the statistical mean is insufficient to describe the actual state of albedo.The GLASS final albedo product, however, can reflect the albedo change during the snowing/melting process.Therefore, with all the benefit of the a priori statistics to fill gaps and stabilize time series, the major signal of the time series of the GLASS final albedo still comes from remote sensing observation.
The GLASS final albedo time series generally agree well with the MODIS albedo time series.The differences between the GLASS final albedo and the MODIS albedo are very small (less than 0.02) during the growing season.However, the differences are much larger during the winter season.The MODIS albedo is lower than the GLASS final albedo when the surface is probably covered by snow (Fig. 3).This can be associated with the MODIS albedo retrieving strategy for a snow-covered surface: observations of partially snowcovered surface are dropped by the MODIS albedo algorithm to maintain the accuracy of snow-free albedo, which would lead to an underestimation in albedo.The GLASS final albedo time series are much more competent in capturing the drastic variations in land surface.For instance, the field measurements indicate that a snowing and melting process happened at the BJ station from day 28 to day 41.This process is well captured by the GLASS final albedo while it is missed by the MODIS albedo.Some discrepancies between the GLASS final albedo and the field albedo are observed.For the BJ station, the field observations are higher than the GLASS final albedo.However, the field observations for the MS3478 station are lower than the GLASS final albedo.For both the D105 and MS3478 stations, the albedo of a probably snow-covered surface in winter decreases much faster than the GLASS final albedo.Spatial scale difference between 1 km pixel and a weather station footprint may be the cause of this discrepancy.Table 5 gives the comparison results of the field measurements with the GLASS final albedo and the MODIS albedo at these three field stations.Since the temporal resolution of the MODIS albedo product is 8 day, both the daily field measurements and the GLASS final albedo are averaged in the 8 days closest to the center of the MODIS product temporal composition window.For both products, the determinant coefficients (R 2 ) value in winter is larger than the R 2 value in non-winter, but the root mean square error (RMSE) in winter is much larger than that in non-winter seasons.For all seasons, the R 2 between the GLASS albedo and field measurements is larger than that of the MODIS albedo.The RMSE of the GLASS albedo is slightly larger than that of the MODIS albedo in winter, but much smaller in non-winter seasons as well as over the year.This indicates that the GLASS albedo have a better agreement with the field measurements than the MODIS albedo.

Conclusions
The GLASS preliminary albedo products (GLASS02A2x, x = 1, 2, 3 and 4) are with high temporal resolution.Due to factors such as cloud coverage, seasonal snow and other uncertainties, the GLASS02A2x albedo products have large areas of missing data and sharp fluctuations.Based on the Bayesian theory, a statistics-based temporal filtering (STF) algorithm was proposed to fill in data gaps and smooth the time series of the GLASS02A2x products.The essence of the STF method is a weighted average of the albedo in neighboring days and a priori statistics.The STF method proposed in this paper has many common features with the traditional parameter estimation methods, such as Kalman filter or variational methods (Liang, 2004;Samain et al., 2008).For instance, the STF method optimizes the albedo parameters by coupling background and observational information.The estimation of background information and observational error is very critical since it affects the accuracy of results.Compared to these methods, the STF method does not require any complicated optimization procedure, and thus is easy to implement.
The a priori statistics plays an important role in the STF algorithm.When the number of the valid GLASS preliminary albedo in the filtering temporal window is insufficient, the quality of filtered results is mostly determined by the a priori statistics.The current a priori database is calculated from the 11 yr of MCD43B3 product.As new high-quality albedo product has become available, the a priori database can be updated.On the other hand, although the STF algorithm needs a priori statistics, the information of the filtered results essentially comes from the GLASS preliminary albedo products.The influence of the a priori database is negligible when there are enough valid GLASS preliminary albedo products in the temporal composition window.Because the GLASS preliminary albedo products are generated with the AB1 and AB2 algorithms which make use of daily clear-sky Terra/Aqua MODIS observation, they usually provide abundant information for the STF algorithm to generate good quality GLASS final albedo product.However, there are still some occasions in which the STF method over smoothes the time series, especially when the albedo change is accompanied by cloud coverage.The uncertainty of the STF result is indicated in the product quality flag.It is recommended that the quality flag should always be checked when using the final GLASS albedo product.
Results and preliminary analysis show that the STF method has greatly improved the spatial integrity and temporal smoothness of the GLASS albedo product.Seasonal variations and trends in albedo are well depicted by the GLASS final albedo product.Compared with the MODIS product, the GLASS final albedo product is much more competent in capturing surface albedo variations.Although it is an algorithm designed for the GLASS albedo product, the STF algorithm is able to incorporate other albedo products, so long as the newly introduced product is statistically consistent with the current products and the a priori database.It is also possible to apply the STF method to other parameters, such as LAI and soil moisture, which presents a subject for our future study.

Table 2 .
Basic information of GLASS preliminary and intermediate albedo products.

Table 5 .
Comparison of field measurement with GLASS and MCD43B3 albedo.