A Method for Retrieving Daily Land Surface Albedo from Space at 30m Resolution

Land surface albedo data with high spatio-temporal resolution are increasingly important for scientific studies addressing spatially and/or temporally small-scale phenomena, such as urban heat islands and urban land surface energy balance. Our previous study derived albedo data with 2–4-day and 30-m temporal and spatial resolution that have better spatio-temporal resolution than existing albedo data, but do not completely satisfy the requirements for monitoring high-frequency land surface changes at the small scale. Downscaling technology provides a chance to further improve the albedo data spatio-temporal resolution and accuracy. This paper introduces a method that combines downscaling technology for land surface reflectance with an empirical method of deriving land surface albedo. Firstly, downscaling daily MODIS land surface reflectance data (MOD09GA) from 500 m to 30 m on the basis of HJ-1A/B BRDF data with 2–4-day and 30-m temporal and spatial resolution is performed: this is the key step in the improved method. Subsequently, the daily 30-m land surface albedo data are derived by an empirical method combining prior knowledge of the MODIS BRDF product and the downscaled daily 30-m reflectance. OPEN ACCESS Remote Sens. 2015, 7 10952 Validation of albedo data obtained using the proposed method shows that the new method has both improved spatio-temporal resolution and good accuracy (a total absolute accuracy of 0.022 and a total root mean squared error at six sites of 0.028).


Introduction
Due to the importance of land surface albedo, in recent years a great deal of research has been carried out on obtaining land surface albedo at moderate to low scales of spatio-temporal resolution by using several types of remote sensing data obtained from different sensors [1][2][3][4][5][6].However, the low spatio-temporal resolution of these albedo data restricts the high-frequency monitoring of land surface changes caused by natural and social events operating at local scales.Thus, retrieval of land surface albedo with better spatio-temporal resolution has been a subject of increasing research attention, and a number of related studies have been conducted.Liang et al. made use of radiative transfer simulation to construct the empirical relation between reflectance and albedo based on a land cover spectral library [7,8].Recently, an improved empirical algorithm introduced the Polarization and Directionality of Earth's Reflectances Bidirectional Reflectance Distribution Function (BRDF) as prior knowledge to consider land surface reflective anisotropy [9]; using this method, topographic effects have been considered [10].Additionally, Shuai et al. (2011) fused Landsat 5 30-m reflectance and MODerate Resolution Imaging Spectroradiometer (MODIS) BRDF at 500-m resolution to retrieve 30-m LANDSAT albedo for a snow-free area [11].Franch et al. (2014) proposed a method to derive the Landsat surface albedo using known BRDF parameters from the MODIS Climate Modeling Grid product using the VJB method [12].Another study retrieving daily MODIS land surface albedo based on the improved VJB method was performed by Franch et al. (2014) [12,13].In addition, some research has focused on fusing albedo products to improve the spatio-temporal resolution [14].These studies improved only the spatial or the temporal resolution of albedo.To acquire albedo data with better spatio-temporal resolution, our previous work combined the reflectance observed by HJ-1A/B at a spatio-temporal resolution of 30 m and 2-4 days with MODIS BRDF at a spatio-temporal resolution of 8 days and 500 m to obtain albedo data at a spatio-temporal resolution of 30 m and 2-4 days [15].Compared to other previous studies, that work retrieved albedo data with an improved spatio-temporal resolution, however, the albedo spatio-temporal resolution was still limited by the spatio-temporal resolution of the HJ-1A/B data used.All of these studies retrieved an image of albedo based on an image of reflectance using the assumption that the BRDF of a single land cover classification is stable.
In general, the spatio-temporal resolution of the albedo data retrieved by the above methods is restricted by the spatio-temporal resolution of the input reflectance data.Improvements to the spatio-temporal resolution of albedo are limited unless reflectance data with better spatio-temporal resolution can be obtained.Downscaling technology provides an opportunity for acquiring reflectance data with better spatio-temporal resolution.
Recently, downscaling of remote sensing data has been the focus of many research studies, particularly for coarse-resolution satellite data [16].These include studies of precipitation [17][18][19][20], soil moisture [20,21], and land surface temperature [22,23].Sorooshian et al. (2011) demonstrated the use of downscaling as a key technology in the application of satellite data [24].For downscaling land surface reflectance, Gao et al. (2006) proposed a spatial and temporal adaptive reflectance fusion model (STARFM) which fuses MODIS reflectance and Landsat reflectance to predict daily Landsat reflectance by considering spectral weight, distance weight, and temporal weight [25].Enhanced STARFM was used to improve the algorithm for predicting daily Landsat reflectance in heterogeneous regions based on STARFM.Song and Huang (2013) combined MODIS data and Landsat data by using high-pass modulation [26].
The method of downscaling reflectance (e.g., Gao and Song methods) and the empirical method of retrieving albedo (e.g., Shuai and Franch methods) have not yet been combined to obtain an albedo product with higher spatial and higher temporal resolution than the currently available albedo products [11,12,25,26].To address this issue, this study describes a new method to retrieve daily land surface albedo from space at 30-m resolution by fusing HJ-1A/B and MODIS.The proposed method first downscales Terra MODIS reflectance (MOD09GA) to HJ-1A/B spatial resolution based on HJ-1A/B BRDF, which includes corrections for illuminating and viewing geometry.Thus, albedo data with higher spatio-temporal resolution are retrieved by downscaled reflectance, which is similar to Shuai's method.Compared to other methods of downscaling reflectance, the proposed method possesses two improvements.Firstly, the proposed method initially converts reflectance from narrowband to shortwave broadband (NTOB) for the MODIS and HJ-1A/B data, which theoretically eliminates the error from the band mismatch between MODIS and HJ-1A/B.Secondly, an empirical relation is constructed between MODIS reflectance and HJ-1A/B reflectance with the same illuminating and viewing geometry as MODIS based on HJ-1A/B BRDF, which effectively eliminates the error resulting from illuminating and viewing geometry mismatches between MODIS and HJ-1A/B.The proposed method also has several differences from the Shuai and Franch methods [11][12] for retrieving albedo: firstly, pre-processing of the NTOB reflectance of MODIS and HJ-1A/B data effectively corrects the error resulting from band mismatch between sensors; secondly, the proposed method combines the methods of downscaling reflectance and retrieving albedo based on prior knowledge and, thus, can retrieve albedo data with better spatio-temporal resolution than the Shuai or Franch methods.
Using the empirical method from Shuai's work, we use downscaled land surface reflectance and MODIS BRDF data to derive the daily 30-m land surface albedo.Sections 2 and 3 introduce the theoretical basis of and describe the procedure for deriving daily 30-m albedo based on downscaled reflectance data and prior MODIS BRDF.Section 4 provides an overview of the study area and experimental data.Section 5 validates these results and describes the evaluation of accuracy.Section 6 lists the conclusions and outlines the disadvantages of the proposed method.

Land Surface Albedo and BRDF
Land surface albedo is the ratio of the upwelling irradiance to the incident irradiance over upper hemispherical direction.It can be expressed as 2 2 2 /2 0 0 0 0 2 /2 00 ( , ) ( , , , ) cos sin cos sin ( , ) cos sin where R(θi, θv, Фv, Фi) is the BRDF; F(θi, Фi) is the downwelling solar irradiance; θ is the zenith, Ф is the azimuth and the subscripts i and v indicate "illuminating" and "viewing", respectively.From Equation (1), we infer that the land surface BRDF is essential for obtaining land surface albedo.In this paper, the kernel-driven model of the BRDF model of MCD43A1 product's official algorithm is adopted [27][28][29]: ( , , , , ) ( , , , ) ( , , , ) where fvol is the coefficient of the scattering kernel kvol, fgeo is the coefficient of the geometrical optical kernel kgeo and fiso is the coefficient of isotropy.

Spatio-temporal Downscaling of Land Surface Reflectance
The principle of downscaling of remote sensing parameters is fixing the relation between the parameter that requires downscaling and its most relative parameter that has a higher spatio-temporal resolution.Land surface reflectance is a fundamental remote sensing parameter and has the following characteristics.Firstly, land surfaces with the same land cover and structure have the same BRDF; based on Tobler's First Law, the nearest land surface with the same land cover should have the same or familiar structure.The theories of radiative transfer and geometric optics show that the reflectivity of the land surface is mainly dependent on the structure of the land surface.Secondly, the BRDF and the magnitude of land surface reflectance depend on the land cover; thus, the land surface reflectance changes seasonally and regularly.A mathematical description of land surface reflectance is: where Ra(θi, θv, Фv, Фi) is the radiance reflected by the land surface in the direction of sensor observation (θv, Фv).In a region with the same area as a pixel of moderate spatial resolution (such as a MODIS pixel), the F(θi, Фi) values of all locations are approximately equal; therefore, the land surface reflectance of a moderate-spatial-resolution pixel includes the average status of all land surface reflectance rhr of the high-spatial-resolution pixels (such as HJ-1A/B) within the moderate-spatial-resolution pixel.In the meantime, considering the reflective anisotropy of the land surface, the land surface reflectance from different satellites must be transformed into that with the same illuminating and viewing geometry (θi, θv, Фv, Фi) as Equation (4).A mathematical expression for this is: ( , , , ) ( , , , ) where rmr is the reflectance of the moderate-spatial-resolution pixel, rhr is the reflectance of high-spatial-resolution pixels in the moderate-resolution pixel, n is the number of rhr within the moderate-resolution pixel and (θi, θvmr, Фvmr, Фi) is the illuminating and viewing geometry of rmr.
Equation ( 4) holds when rmr and rhr are observed by different sensors at the same illuminating and viewing geometry.In reality, the reflectance data observed by different sensors have different illuminating and viewing geometries, and the reflectance changes with variations in the illuminating and viewing geometry.In Equation ( 4), the BRDF corresponding to the reflectance rhr of high-spatial-resolution pixels is required to convert the illuminating and viewing geometry of the reflectance in those pixels.A mathematical expression for this is: where BRDFhr is the BRDF of the land surface reflectance rhr in the high-spatial-resolution pixels and BRDFhr and rhr have the same spatio-temporal resolution.The value of rhr for any illuminating and viewing geometry (θi, θr, Фr, Фi) can be calculated from BRDFhr.
In fact, rmr is not equal to the average of rsr, because of the observed time difference of rmr and rhr between different sensors.If the land surface is stable and its BRDF shifts only in magnitude within the observed time difference, rmr obeys a linear relation coh with the average of rhr.Thus, Equation ( 4) can be adjusted to: Usually rmr has a higher temporal resolution than BRDFhr/rhr, and more than one rmr value is observed in the revisit period of rhr.Thus, from Equation ( 6), we infer that more than one coh value can be acquired with the spatio-temporal resolution of rmr in the revisit period of rhr.In addition, rhr multiplies coh to derive more than one land surface reflectance rh with the same spatial resolution as rhr and the same temporal resolution as rmr in the revisit period of rhr:

Method
Based on Section 2.1, downscaling land surface reflectance by fusing MODIS and HJ-1A/B data ensures that daily 30-m land surface reflectance for use as input data can be acquired after atmospheric correction and the conversion of narrowband to broadband (NTOB) shortwave surface reflectance [7,19,30,31].Then, incorporating the classification of MODIS pixels as pure or mixed, two empirical methods were adopted to combine the downscaled reflectance and MODIS BRDF and, thus, retrieve daily 30-m albedo data [15].Within the MODIS pure pixel, the albedo is first derived based on the downscaled reflectance and MODIS BRDF; the derived albedo and the corresponding downscaled reflectance are introduced into a lookup table (LUT) to construct the ratio of the downscaled reflectance and albedo to be derived in the MODIS mixed pixel.From HJ-1A/B 30-m land cover data, some classifications occupy more or less than 60% of a MODIS 500-m pixel, these MODIS pixels are regarded as pure or mixed pixels, respectively.In our previous study [15], a specific discussion of the thresholds was carried out to demonstrate that 60% is a reasonable threshold value.A flowchart for data retrieval is provided in Figure 1.

NTOB Reflectance
For HJ-1A/B observations at the top of the atmosphere (TOA), an atmospheric correction must be performed based on Second Simulation of a Satellite Signal in the Solar Spectrum (6S).Due to the lack of information on corresponding atmospheric parameters (e.g., aerosol optical depth and water vapor), cloud masking is initially carried out using the four band thresholds of HJ-1A/B based on the prior knowledge, to ensure that the reflectance at TOA with no cloud was selected for further processing.Subsequently, 6S runs were carried out based on the HJ-1A/B parameters (e.g., spectral response function and illuminating and viewing geometry) and the empirical atmospheric parameters for the atmospheric correction for HJ-1A/B band reflectance.
After atmospheric correction, the spectral difference of data from different sensors must be considered.Spectral difference hinders fusing of data of different sensors, which could introduce errors in the course of multi-source application.Thus, the reflectance of different sensors must be normalized before running the algorithm to retrieve albedo.Table 1 shows the spectral differences between MODIS and HJ-1A/B.Considering the spectral range of a sensor, it is possible to obtain the shortwave albedo from the band albedo: where wl is the band weight and l is the number of bands.Tables 2 and 3 list the weights used to convert band to shortwave albedo for HJ-1A/B data.In MODIS official products (MCD43B1), the shortwave BRDF/albedo can be calculated by weighting the band BRDF/albedo using the weights for converting band albedo to shortwave albedo.Thus, in this study, the weights converting HJ-1 band albedo to shortwave albedo are used for converting HJ-1 band reflectance to shortwave reflectance (Table 2); the weights for NTOB MODIS are used for NTOB MODIS reflectance (Table 3).

Downscaling Reflectance by Fusing MODIS and HJ-1A/B Data
The spatio-temporal resolution of the derived albedo depends on the resolution of the remote sensing data that are used for obtaining the albedo.The intended use of the proposed method is to derive the daily 30-m land surface albedo, which requires information about the daily 30-m land surface reflectance.However, current remote sensing products do not possess 30-m spatial resolution and 1-day temporal resolution.Thus, obtaining daily 30-m land surface reflectance is essential for the proposed method.The proposed method solves this problem by fusing high-scale spatial data (HJ-1A/B reflectance) and high-scale temporal data (MODIS reflectance) to downscale land surface reflectance.
We assume that land use type is unchanged over an 8-day period.The unchanged land use means that the BRDF shape of the land surface is unchanged, and the BRDF for a single location shows only a shift in magnitude over an 8-day period.Therefore, there must be one image of the HJ-1A/B reflectance Rhj/BRDF BRDFhj that is the closest time and location to each image of MODIS reflectance Rmodis over an 8-day period.For each MODIS reflectance image, the corresponding HJ-1A/B BRDF at 30 m is used to calculate the HJ-1A/B reflectance with the same illuminating and viewing geometry as the corresponding MODIS based on Equation (5).Subsequently, within a MODIS pixel, a coefficient is calculated that is the ratio of the MODIS reflectance to the average value of the corresponding HJ-1A/B reflectance derived from HJ-1A/B BRDF: this coefficient has the same spatio-temporal resolution as the MODIS reflectance (based on Equation ( 6)).The daily 30-m land surface reflectance Rh is calculated from the coefficient multiplied by the corresponding reflectance (Equation ( 7)).The equation for Rh is: where Rmodis is the reflectance observed by MODIS at the illuminating and viewing geometry (θi, θvmodis, Фvmodis, Фi, Tmodis); BRDFhj is the known HJ-1A/B BRDF in the 30-m subpixel of Rmodis; thj is the time when Rhj is obtained; Tmodis is the time when Rmodis is observed and n denotes the number of HJ-1A/B pixels within the MODIS pixel.

Retrieving Daily 30-m Land Surface Albedo for Pure MODIS Pixels
Within a MODIS pure pixel, the MODIS BRDF parameters and the downscaled daily 30-m land surface reflectance have the same land surface reflective distribution shape, except for some shifts in magnitude caused by the optical differences between the products from different sensors over the 8-day measurement period (the temporal resolution of the MODIS BRDF product).Within one MODIS pixel, the MODIS BRDF provides the BRDF shape of every 30-m pixel during these eight days, but cannot provide information on the magnitude of the shift.Thus, for a HJ-1A/B pixel within a pure MODIS pixel, the fitted linear relation between the MODIS BRDF and the daily land surface BRDF in the pixel is given by Equation ( 10): ( , , , , ) / ( , , , , ) where Rh is the downscaled daily 30-m land surface reflectance; BRDFm is the calculated reflectance that is the value of the MODIS BRDF at (θi, θvh, Фvh, Фi) and is the same as Rh; Tm indicates the 8-day period of the MODIS BRDF product and th is the frequency of HJ-1A/B observations during Tm.Thus, the daily 30-m BRDF can be calculated using the following equation: , , , ) After the daily BRDF for each 30-m pixel has been retrieved, BRDFh is integrated to derive the daily 30-m albedo on the basis of Equation (1).

Retrieving the Daily 30-m Land Surface Albedo for Mixed MODIS Pixels
The LUT is constructed on the basis of retrieved daily 30-m albedo data and the corresponding downscaled daily 30-m reflectance data within the pure MODIS pixel, which is used for deriving the daily 30-m albedo within mixed MODIS pixels.For every 30-m pixel in a mixed MODIS pixel, a window (3 × 3 km) is delimited around the 30-m pixel, and all retrieved values of the albedo and the corresponding reflectance of the 30-m pixels with the same land cover as the target 30-m pixel in this image subset are extracted to construct the LUT.In the LUT, every pair of values for retrieved albedo and corresponding land surface reflectance can be assigned a ratio, con, from Liang et al. [8].
For eliminating error from a singular value, averaging all con in LUT can provide a uniform ratio comixed: As the daily 30-m reflectance rhj and the corresponding comixed are known, the albedo Albedo_HJ for a 30-m pixel within a mixed MODIS pixel is:

Study Area and Study Data
The Heihe river basin is an inland river basin in China that contains various types of ecosystems and land cover.There have been many scientific experiments conducted in this area, providing large numbers of ground measurements; therefore, this region was selected as the study area.

Ground Measurements
Ground measurements from the HiWATER project obtained from seven sites during the 2012 growing season were selected to validate the retrieved results (Figure 2 and Table 4) [33,34].In addition to the large amount of available data from previous studies, selection of this study area and the corresponding ground measurements facilitates comparison between the results of the proposed method and the results of our previous study in 2014, allowing a determination of whether the proposed method is able to improve the temporal resolution of albedo data with guaranteed precision.

Study Data
Land surface reflectance includes two types of remote sensing data: HJ-1A/B band reflectance and MODIS reflectance (MOD09GA).
The HJ-1 satellites were launched on 6 September 2008.The HJ-1 series includes two satellites (HJ-1A and HJ-1B) that bear a wide-coverage multispectral charge-coupled device (CCD) camera (HJ-A and HJ-1B), a hyperspectral imager (HSI, HJ-A only) and an infrared camera (HJ-1B only) for dynamic modeling of eco-environmental changes over a wide range of wavelengths [35].The HJ-1A/B CCD data are used in this work.Each of the HJ-1A/B satellites has two CCDs containing four shortwave bands (Table 1) with 30-m spatial resolution.The synchronous observations of both satellites (HJ-1A and HJ-1B) result in the CCD data having a 2-4-day temporal resolution.The China Centre for Resources Satellite Data and Application published HJ-1A/B CCD data that include only the geolocated DN, the radiative calibration parameter, and the illuminating and viewing geometry.Based on the above information, batch processing was performed to convert HJ-1A/B DN to reflectance using the IDL language.After converting radiance to reflectance, atmospheric correction was necessary to acquire band land surface reflectance from the band reflectance at the TOA by 6S.Because HJ-1 has no corresponding aerosol data, we carried out atmospheric correction based on empirical atmospheric parameters; however, we selected images with no cloud to ensure accuracy.
The land surface reflectance of Terra MODIS (MOD09GA) has 1-day temporal resolution and 500-m spatial resolution, and is used for downscaling land surface reflectance in this study.
On the basis of the International Geosphere-Biosphere Program [36], the HJ-1A/B land cover data includes land cover type subdivisions of crop and neglected land, which are not present in the Heihe river basin [37].The HJ-1A/B land cover data was validated in the Dahuofang area of Liaoning Province and the Heiquan area of Gansu Province, showing very high mapping accuracy of 95.76% and 83.78% and high Kappa coefficients of 0.9423 and 0.8165, respectively.HJ-1A/B land cover data are used for distinguishing pure and mixed MODIS pixels.
The MODIS BRDF (MCD43A1), which constitutes prior knowledge, includes the coefficients of the shortwave kernel-driven BRDF model with 8-day temporal resolution and 500-m spatial resolution.
The HJ-1A/B BRDF data with a temporal resolution of 2-4 days used in this study were acquired between 1st July and 30th September 2012.The HJ-1A/B BRDF parameters can be obtained based on the MODIS BRDF and HJ-1A/B reflectance of our previous work, which are used for downscaling reflectance [4].

Results
Using the proposed method, we derived a daily 30-m albedo parameter for the Heihe river basin.The analytical results are presented and discussed below.

Comparison of Daily Reflectance at 30 m with That of MODIS and HJ-1A/B
Using the derived daily 30-m reflectance, a comparison of the derived daily 30-m reflectance with reflectance from MODIS and HJ-1A/B is performed.A comparison between the downscaled daily 30-m reflectance and MODIS reflectance on 30th September 2012 in the study area is shown in Figure 3.
From Figure 3, the data of the downscaled daily 30-m reflectance and the MODIS reflectance have similar spatial distributions.Both can correctly represent the spatial distribution of land cover in the study area, but the downscaled daily 30-m reflectance appears to show more detail of the land cover than the corresponding MODIS reflectance.A quantitative comparison of the red transect line illustrated in Figure 4 demonstrates that the daily 30-m land surface reflectance data display a larger dynamic range (the downscaled reflectance has a good correlation with MODIS, and a greater standard deviation than MODIS in Figure 4).This means that the downscaled daily 30-m reflectance has the potential for precise detection of small-scale emergency changes in the land surface.In addition to the higher spatial resolution, the derived daily 30-m reflectance scale also has a higher temporal resolution.The time series of the derived daily 30-m reflectance and the HJ-1A/B land surface reflectance from 1 July to 30 September 2012 at six sites are extracted, and a comparison between the extracted time series is performed to demonstrate the higher temporal resolution of the derived daily 30-m reflectance (Figure 5).Due to the higher temporal resolution, the derived daily reflectance contains more data than the HJ-1A/B reflectance for the same period, which means a higher possibility that the derived daily reflectance can detect natural emergencies and human activities.An absolute error between the two datasets is present because the derived daily reflectance includes information from the MODIS data.Statistical analysis on the data of Figure 5 shows that R 2 of all sites is 0.73 and the root mean squared error (RMSE) of all sites is 0.029.Considering the three-month short phenological period and the stable land use, these statistical results demonstrate that the numerical differences between the "Daily" and "HJ-1A/B" reflectance values fall in an acceptable range.From Figure 5, the daily reflectance has a higher temporal resolution than the HJ-1A/B reflectance.Figures 3 and 4 show that the downscaled daily 30-m land surface reflectance has better spatio-temporal resolution than the MODIS and HJ-1A/B reflectance, which can monitor small-scale changes of land surface albedo caused by natural events and human activities.

Comparison of the Retrieved Albedo Values with Those of MODIS and HJ-1A/B
Figure 6 shows the derived albedo and the corresponding MODIS albedo for three periods between July and September 2012 and corresponding 2D scatter plots comparing albedo values in regions of interest.A transect line is randomly selected to define regions of interest that contain different land cover types, from which data are extracted for further analysis (Figure 7).  Figure 6 compares the derived albedo using the corresponding MODIS albedo as a benchmark for the Heihe river basin.The consistency of spatial distribution and the accuracy are both good.In addition, compared with the MODIS albedo product, the derived albedo describes the distribution of land cover types more precisely (Figure 6).The derived albedo has the potential to be used for precise detection of small-scale changes in land surface albedo caused by natural events and human activities.In addition, the spatial correlation between the daily albedo and the MODIS albedo appears to be not very significant, which is due to the different spatial resolutions of the two datasets.The larger standard deviation of daily albedo and the smaller standard deviation of MODIS albedo can be observed in the 2D scatter plot of Figure 6, which shows that the slopes of the linear fitting equations are less than 1. Figure 7 demonstrates, quantitatively, that both albedo datasets are in good agreement.Land surface albedo is lowest between June and August because of vegetation growth, and albedo increases with vegetation withering and crop harvesting from September onwards.The derived albedo has a larger dynamic range because of its higher resolution.
To demonstrate further the good temporal characteristics of the derived daily albedo, the number of time series of the derived daily 30-m albedo and of HJ-1A/B albedo retrieved using Shuai's method from 1 July to 30 September 2012 at six sites are extracted (Table 5).Deleting singular values caused by cloud or other reasons, the temporal series of the derived albedo has more virtual values than the HJ-1A/B albedo retrieved using Shuai's method.This demonstrates that the proposed method is able to derive albedo data with a higher temporal resolution than other currently used methods (Shuai's work and Franch's work).

Validation of the Retrieved Albedo
Based on the ground measurements, validation of the retrieved albedo from July to September 2012 was carried out for the study area.The time series and scatter plots of the retrieved albedo and the ground observations at different sites are illustrated in Figure 8.Comparison of the time series shows a reasonable absolute difference in magnitude, but also a seemingly unrealistic correlation between the retrieved albedo and that observed on the ground.Two possible reasons for why the correlation is not good are the shorter period of data validation and the use of an atmospheric correction based on prior atmospheric parameters.The atmospheric correction based on the previous knowledge of HJ-1A/B reflectance is carried out only on sunny days, and the retrieved albedo values show a fluctuation around the ground measurements, which means that atmospheric correction can be used to ensure reasonable accuracy.The other, more important, reason is that the period of data validation is short, only three months, during which the land surface shows no obvious variations in magnitude.Other studies have conducted validation on the basis of several years' worth of ground measurements.Although this study conducted validation using results for three months, it also includes more comparative samples (because of the temporal advantage of the retrieved albedo data).Comparison of the results from the proposed method and MODIS albedo (Table 6) shows that the derived daily 30-m land surface albedo has reasonable retrieval accuracy.The proposed method combines downscaling of reflectance and empirical albedo retrieval, and is able to obtain albedo data with better temporal resolution and accuracy than other currently used methods (Table 7).

Test of Applicability and Robustness of the Proposed Method
From the above analysis, the proposed method can retrieve daily albedo data with 30-m resolution that have acceptable precision.However, there are some issues that may influence the accuracy of the proposed method.For example, there are few pure MODIS pixels in the study area, which may limit the applicability of the proposed method.In addition, the proposed method may have lower precision for mixed MODIS pixels than for pure MODIS pixels.To address this, the applicability and robustness of the proposed method is tested.
Pure/mixed MODIS pixel mapping based on two different thresholds (50% and 60%) is illustrated in Figure 9.There should be more pixels using the 60% threshold than using the 50% threshold.More than 76% of the entire study area is pure MODIS pixels for a threshold of 60%.In the center of the study area both densely populated and agricultural zones are present, which are more complex land use types than the natural watershed.Thus, the proposed method can extract sufficient "pure" MODIS pixels for albedo retrieval; and is suitable for albedo retrieval in other natural watersheds.Further validation of the proposed method for mixed MODIS pixels is performed.Firstly, classification of pure and mixed MODIS pixels was performed based on the two thresholds (50% and 60%).These two thresholds can meet the demand for retrieval precision, according to our earlier work.Secondly, the proposed method was carried out for both thresholds (50% and 60%), and two daily albedo maps were constructed (Figure 10).During this analysis, some pixels were pure pixels when the threshold was 50% and mixed pixels when the threshold was 60%. Figure 10 shows the difference between the results for pure/mixed MODIS pixels.
The difference between Figure 10a and 10b is difficult to detect visually.Thus, a comparison of Figure 10a and 10b is provided in Figure 10C.Figure 11 shows that there is a very good correlation and a small RMSE for the results for pure and mixed pixels, demonstrating that the accuracy of the proposed method is comparable for both pure and mixed MODIS pixels (Figure 11).

Conclusions
The proposed algorithm makes use of downscaling technology to acquire daily 30-m land surface reflectance with good spatio-temporal resolution that cannot be acquired by any single-sensor observation.The downscaled daily 30-m land surface reflectance and MODIS BRDF data are combined empirically to first derive the daily 30-m land surface albedo.The proposed method of albedo retrieval first introduces the downscaling technology and the idea of data fusion, and improves the retrieval accuracy and the spatio-temporal resolution of albedo data.The accuracy of the derived daily 30-m albedo is comparable to or better than those of existing methods and has better spatio-temporal resolution than the currently available data.The newly developed method is of practical use and is, thus, likely to be adopted.In general, the proposed method combines downscaling of reflectance and empirical retrieval albedo, which provides albedo data with a reasonable precision and a higher spatio-temporal resolution than other methods.

Figure 1 .
Figure 1.Flowchart describing retrieval of daily 30-m land surface albedo data.

Figure 4 .
Figure 4. Land surface reflectance obtained from the two datasets along the line shown in Figure 3. "Daily" means the downscaled daily 30-m land surface reflectance; "MODIS" means the MOD09GA product ("MODIS" is resampled to 30-m spatial resolution to yield the same resolution as "Daily").
(a-left) (a-right) 2D scatter plot of a-left and a-right (b-left) (b-right) 2D scatter plot of b-left and b-right (c-left) (c-right) 2D scatter plot of c-left and c-right

Figure 6 .
Figure 6.Derived and MODIS land surface albedo in the Heihe river oasis (from 1 July to 30 September 2012).(a-c, left) Derived daily albedo, (a-c, right) MODIS albedo.The MODIS albedo is resampled to have the same 30-m spatial resolution as the daily albedo.

Figure 7 .
Figure 7.The spatial series of albedo derived from MODIS and the proposed method along the transect regions in Figure 6: (a) data from Figure 6a,d; (b) data from Figure 6b,e; (c) data from Figure 6c,f.

Figure 10 .
Figure 10.Retrieved albedo data.(a) 50% threshold, (b) 60% threshold.(c) Spatial distribution of critical pixels, which are pixels that is pure when the threshold is 50% and mixed when the threshold is 60%.Other pixels are noncritical pixels.

Figure 11 .
Figure 11.2-D scatter plot of daily albedo values obtained using the proposed method within critical MODIS pixels.

Table 1 .
Spectral distributions of MODIS and HJ-1A/B CCD.

Table 4 .
Descriptions of the study sites and instrumental setup.

Table 5 .
The ratio of the number of time series of the derived daily 30-m albedo to HJ-1A/B albedo.

Table 6 .
Comparison of the albedo values obtained using different methods

Table 7 .
Comparison of the spatio-temporal features of albedo data obtained using different methods.