Remote Sensing Derivation of Land Surface Albedo at High Resolution by Combining Hj-1a/b Reflectance Observations with Modis Brdf Products

Land surface albedo is an essential parameter for monitoring global/regional climate and land surface energy balance. Although many studies have been conducted on global or regional land surface albedo using various remote sensing data over the past few decades, land surface albedo product with a high spatio–temporal resolution is currently very scarce. This paper proposes a method for deriving land surface albedo with a high spatio–temporal resolution (space: 30 m and time: 2–4 days). The proposed method works by combining the land surface reflectance data at 30 m spatial resolution obtained from the charge-coupled devices in the Huanjing-1A and-1B (HJ-1A/B) satellites with the Moderate Resolution Imaging Spectroradiometer (MODIS) land surface bidirectional reflectance distribution function (BRDF) parameters product (MCD43A1), which is at a spatial resolution of 500 m. First, the land surface BRDF parameters for HJ-1A/B land surface reflectance with a spatial–temporal resolutions of 30 m and 2–4 day are calculated on the basis of the prior knowledge from the MODIS BRDF product; then, the calculated high resolution BRDF parameters are integrated over the illuminating/viewing hemisphere to produce the white-and black-sky albedos at 30 m resolution. These results form the basis for the final land surface albedo derivation by accounting for the proportion of direct and diffuse solar radiation arriving at the ground. The albedo retrieved by this novel method is compared with MODIS land surface albedo products, as well as with ground OPEN ACCESS Remote Sens. 2014, 6 8967 measurements. The results show that the derived land surface albedo during the growing season of 2012 generally achieved a mean absolute accuracy of ±0.044, and a root mean square error of 0.039, confirming the effectiveness of the newly proposed method.


Introduction
Land surface albedo is defined as that fraction of the incident solar radiation in the 0.4-2.5 μm spectrum domain which is reflected by the land surface [1]; it is an important variable that is used for global/regional climatic modeling and calculation of the net shortwave radiation [2,3].Land surface albedo data with high spatio-temporal resolution allows the monitoring of sudden natural processes or human activities at local scales (e.g., changes due to snowfall, snowmelt, precipitation, crop damage by disaster, and crop harvest).This data can also be used to evaluate the surface radiation budget to monitor the drastic variations in the natural environment, and has been used to monitor the urban heat island effect to guide the fight against climate change in urban areas.
In the past few years, satellite observations have led to the production of several land surface albedo products, including the Moderate Resolution Imaging Spectroradiometer (MODIS) albedo product (with a spatial resolution of 500-1000 m and a temporal resolution of 16 days) [4][5][6], the land surface albedo product from the Polarization and Directionality of the Earth's Reflectance (POLDER) satellite (with a 1/12° spatial resolution and a 10 day temporal resolution) [7,8], the Medium Resolution Imaging Spectrometer (MERIS) albedo product (with a temporal resolution of 16-30 days and an spatial resolution of approximately 0.05° climate modeling grid (CMG)) [9], and the Meteosat and Meteosat Second Generation (MSG) product (with a characteristic time scale of 1 day by a iterative scheme or 10 days by composition of 1 day products) [10][11][12].The spatial resolution of these typical land surface albedo products range from 500 m to 5000 m, and their temporal resolution range from one day to one month.The relatively coarse resolution of these products significantly limits their application to many sectors.
Recently, special attention has been paid for developing algorithms to retrieve land surface albedo data with a higher spatio-temporal resolution.For instance, Liang et al. (1999Liang et al. ( , 2005) ) proposed an albedo retrieval method on the basis of the radiative transfer simulation, usually called the direct-estimation algorithm [13][14][15].It utilized the MODTRAN package [16] to simulate the top of atmosphere (TOA) reflectance under the Lambertian assumption for the land surface; a relation between the MODIS TOA reflectance and the broadband surface albedo was then obtained by training a neural network.This method could produce a daily albedo product using just one scene of MODIS data.Later, this method was further improved by Liang (2005) to estimate the daily land surface albedo of the Greenland ice sheet from MODIS data.Unfortunately, the studies mentioned above did not consider the anisotropy of the land surface.To this end, an improved direct-estimates method that considers the anisotropy using the POLDER Bidirectional Reflectance Distribution Function (BRDF) has been reported using the MODIS dataset [17], which can generate daily land surface albedo with a spatial resolution of 1000 m.Additionally, a method to estimate the snow-free albedo at a spatial resolution of 30 m has also been proposed on the basis of a fusion of the Landsat 5 reflectance data and the classification-based MODIS 500 m anisotropic BRDF information [18].This product includes the snow-free land surface albedo with a temporal resolution of 16 day, and was validated over the Maryland, DC, and eastern Virginia areas of the US.Particularly, generating land surface albedos with high spatio-temporal resolutions over a large area by making full use of currently existing satellite observations is necessary.HJ-1A/B launched by China in 2009, have four charge-coupled devices (CCDs) each (with four bands, similar to Landsat 5), and each of these CCDs have a spatial resolution of 30 m and temporal resolution of 2-4 days [19].Thus, producing a land surface albedo product with high spatio-temporal resolution on the basis of the HJ-1A/B data is possible if an appropriate retrieval algorithm could be devised.To this end, this paper proposes an alternative algorithm based on Shuai's method [18], according to which MODIS pixels (500 m) are classified as either "pure" or "mixed."When the MODIS pixels are pure, an empirical relation between the MODIS BRDF product and the HJ-1A/B data can be used to calculate the BDRF parameters of the HJ-1A/B pixels (30 m) within the pure pixel.On the other hand, if the MODIS pixel is mixed, then an empirical relation between the HJ-1A/B CCD reflectance data and the collocated derived albedo is constructed.In this fashion, the albedos of the cases of both pure and mixed MODIS pixels can be calculated.
The remainder of this paper is organized as follows.Section II provides a generic mathematical description and specific implementation of the improved algorithm for retrieving albedo by combining the HJ-1A/B CCD reflectance data with the MODIS BRDF.Section III presents a description of the study area and the ground measurement data.Section IV gives the results of this experiment, and presents a validation of these results.Section V draws the conclusion and outlines the objectives for further future research.

Theoretical Basis
The land surface albedo is a measure for the reflectivity of the earth's surface.According to its definition (see introduction), a mathematical expression of the shortwave land surface albedo can be given as where R(θ i , θ v , Ф v , Ф i ) is the BRDF, which is defined as the ratio between the reflected radiance in a particular direction and the irradiance incident from the illuminating direction in a solid angle.F(θ i ,Ф i ) is the downward solar irradiance, θ i , and θ v are the solar zenith angle and view zenith angle, respectively, and Ф i and Ф v are the solar azimuth and view azimuth angle, respectively.Considering the fact that the surface downward solar irradiance comprises the direct solar irradiance (F dir ) and the diffuse solar irradiance (F dif ), the shortwave land surface albedo can be written as Equation (2) in the course of a specific calculation [20]: where a dh is the directional hemispherical albedo (black-sky albedo) obtained by integrating the BRDF over all viewing angles at any desired solar angle.a dh is defined as the albedo in the absence of a diffuse solar irradiance, and is a function of the solar zenith angle: , ( A further integration of a dh over all illumination angles results in the bi-hemispherical albedo (white-sky albedo), a bh , under isotropic illumination, which is defined as the albedo in the absence of a direct solar irradiance: To characterize the BRDF in Equations ( 3) and ( 4), a kernel-driven BRDF model is adopted, which is a linear composite of the three kernels for the isotropic scattering, the volume scattering, and geometric optics scattering [4,21]: where the scattering term K vol is taken from the Ross-thick kernel [22], the geometrical optical term k geo is from the Li-SparseR kernel [23], f vol is the coefficient of K vol , f geo is the coefficient of K geo , and f iso is the coefficient of isotropy.The land surface actually scatters anisotropically in the shortwave spectral range, which is usually characterized by the BRDF in the literature.Theoretically, the same land surface should have similar anisotropic scattering in the same spectral range.Based on this, if an area is pure homogeneous in terms of its land cover, its BRDF in the shortwave spectral range should be identical even if it is observed by different sensors.However, the fact is that the BRDF derived from different sensors usually show some shift in magnitude because of minor differences in the visiting times of the sensors.If the BRDF derived from one sensor is known, the BRDF observed by another sensor can be obtained by constructing an empirical relation between the land surface reflectances observed by the two sensors at same solar and viewing geometry.The empirical relation can be expressed as where c is an empirical coefficient between the land surface reflectances, R 1 and r 2 observed by two sensors with the same solar and viewing geometry (θ i , θ v2 , Ф v2 , Ф i , t); in the above equation, the r 2 and R 1 both have t implying that they are taken at the same time period.Hence, the BRDF observed by one sensor can be expressed in terms of that observed by another sensor after combining Equations ( 5) and ( 6): Once the BRDF (R 2 ) becomes known, a dh at any desired solar angle can be derived by the integration of BRDF over all viewing angles (Equation ( 3)); further, a bh under isotropic illumination is derived by integration of a dh over all illumination angles (Equation ( 4)).Finally, the land surface albedo can be calculated on the basis of Equation (2).

Derivation of Shortwave Broadband Albedo with High Resolution
The new method for deriving the shortwave broadband albedo with high resolution presented in this paper starts with an atmospheric correction of the HJ-1A/B data [24,25], after which the land surface reflectance for the HJ-1A/B is calculated.After calculating the land surface reflectance, a narrowband to broadband (NTOB) shortwave surface reflectance conversion is performed to ensure the same spectral range for the HJ-1A/B data and the MODIS BRDF product, which will be shown in Subsection 2.2.1.MODIS 500 m pixels are classified as either "pure pixel" or "mixed pixel".When one specific classification covers more than 50% of the MODIS pixel, the corresponding MODIS pixel is deemed as a "pure pixel," otherwise it is classified as a "mixed pixel".If HJ-1A/B pixels (30 m) are located in a pure MODIS pixel (500 m), an empirical coefficient between these two datasets (Equation ( 7)) is used to estimate the BRDF parameters of the HJ-1A/B pixel.If the MODIS pixel is a mixed pixel, an empirical coefficient between the HJ-1A/B CCD reflectance data and the collocated derived albedo is constructed.Finally, the direct-hemispheric albedo, bi-hemispheric albedo, and the land surface albedo at the scale of HJ-1A/B (30 m) are calculated on the basis of the scaled BRDF parameters.Figure 1 gives the scheme of the algorithm.

Converting NTOB Shortwave Land Surface Reflectance of HJ-1A/B CCD
Based on radiative transfer theory, the ground-measured land surface reflectance spectra of different land surface cover types with different solar and viewing geometries are inputted into MODTRAN4.0[16] to simulate both the narrowband land surface reflectance of HJ-1A/B data and the shortwave land surface reflectance.The empirical formula converting shortwave land surface reflectance from narrow-band to broadband can be acquired by a linear fitting of the narrowband and shortwave broadband land surface reflectances (Equation ( 8)): where, ρ(θ i ,θ v ,Φ v ,Φ i ) is the broadband shortwave land surface reflectance, w l is the weighing coefficient of each band, and l is the specific bands used.Table 1 gives the weights for conversion of the HJ-1A/B data, adapted from [26].As mentioned above, the same land surfaces should have identical anisotropic scattering in the same spectral ranges.If an area covered by both MODIS and HJ-1A/B pixels is pure and homogeneous in terms of the land cover type, then the retrieved BRDF parameters based on the data from the two sensors should have a similar shape, with the exceptions of a few shifts in reflectance magnitude due to the optical difference between the sensors or differences in the dates of acquisition.For example, the MODIS BRDF on Day of Year (DOY) 161 in 2012 is derived from the land surface reflectance between DOY 161 and 177 of 2012, which represents an average land surface scattering status during those 16 days, whereas the HJ-1A/B BRDF obtained on DOY 161 of 2012 characterizes the status on that day, but not an average status over 16 days.Within one MODIS pixel, the HJ-1A/B land surface reflectance during those 16 days is assumed to have a fixed empirical relation with the land surface reflectance of the corresponding MODIS pixel that depends on the solar and viewing geometry in each HJ-1A/B pixel.Hence, if an HJ-1A/B pixel (30 m) is located inside a pure MODIS pixel (500 m), an empirical coefficient between these two datasets is used to estimate the BRDF parameters of the HJ-1A/B pixel: where r hj is the land surface reflectance observed by HJ-1A/B, R m is the land surface reflectance calculated by the BRDF of the MODIS pixel at the HJ-1A/B solar and viewing geometries specified by (θ i , θ vhj , Ф vhj2 , Ф i , t), T indicates each 8 day period of the MODIS BRDF product, and t specifies the number of days within each MODIS 8 day period.Thus, the BRDF of the HJ-1A/B pixel can be calculated on the basis of the following relation: Once the BRDF at each HJ-1A/B pixel is retrieved, integrating R hj over all viewing angles results in a directional hemispherical albedo at any desired solar angle (Equation ( 3)).A further integration over all illumination angles results in a bi-hemispherical albedo under isotropic illumination (Equation ( 4)).Then, the land surface albedo for the HJ-1A/B pixel is finally calculated as the linear composition of the directional-hemispherical albedo and the bi-hemispherical albedo weighted by the components of solar direct and sky-diffused irradiance (Equation ( 2)).

Retrieving Land Surface Albedo at Mixed MODIS Pixels
If a MODIS pixel is covered by multiple land cover types (mixed pixel), the method described in Section 2.2.2 cannot be applied because of the heterogeneous property of the MODIS pixel.In this case, the BRDF of such a MODIS pixel is determined by the mixed scattering status of various land cover types.To derive the albedo of an HJ-1A/B pixel within a mixed MODIS pixel, a look-up table (LUT) is constructed on the basis of the results of Section 2.2.2.Specifically, for the MODIS image being studied, all the pure pixels are first collected, and then the HJ-1A/B pixel albedo is calculated on the basis of the MODIS BRDF, according to the method provided in Section 2.2.2 for each pure MODIS pixel.Then, the derived HJ-1A/B albedo, the position of the corresponding MODIS pure pixel, the land cover types, as well as the corresponding HJ-1A/B reflectance (Equation ( 10)) are all arranged into an LUT.One pure MODIS pixel corresponds to one record in the LUT.When an HJ-1A/B pixel is located in a mixed MODIS pixel, an image subset (1.5 km × 1.5 km) around the HJ-1A/B pixel is delimited, and all the possible MODIS pixels possessing the same land cover type as the HJ-1A/B pixel within this area are extracted from the LUT according to their positions.For each pure MODIS pixel, a ratio, co n , can be calculated according to the work of Liang et al. [14] by dividing the HJ-1A/B surface reflectance by the corresponding derived albedo.
Since there maybe more than one MODIS pure pixel within the delimited area, the following equation is used to obtain a uniform ratio: where co mixed is the mean of all co n within the delimited area.According to the above description, if the HJ-1A/B reflectance, r hj , and the corresponding a mixed are provided, the albedo for an HJ-1A/B pixel within a mixed MODIS pixel can be easily derived as follows: where Albedo_HJ is the derived albedo for an HJ-1A/B pixel within a mixed MODIS pixel.

Study Area and Data
The Heihe river basin in the arid region of northwestern China is selected as the study area.The Heihe river basin is a typical inland river basin, and includes various ecosystems and land covers, including forest, alpine grassland, crop land, shrubs, and the Gobi desert.Numerous scientific experiments have been conducted to collect rich ground albedo measurements in the Heihe river basin [27,28]; hence, this region is an ideal study area.Figure 2 shows the location and the land covers of the study area.

Ground Measurements
The ground measurements of the albedo (which are used to validate the retrieval algorithm) are taken from the results of the HiWATER project at seven sites [28] located in the mid reach of Heihe river basin (Figure 3).Table 2 gives the information of land cover and characteristics of measurements of each site.
The sites are generally equipped with a net radiometer mounted on a horizontal platform atop a 6-12 m tall meteorological tower with both a skyward-and a downward-facing pyranometer to measure the shortwave radiation.The downward and upwelling solar radiations are continuously recorded by the skyward facing and downward facing pyranometers, respectively.Shuai et al. (2011) [18] believed that the downward facing pyranometer is most effective when it has an 81° field of view (FOV).Furthermore, our analysis experiment shows that the energy received by pyranometer is over 75% of the energy received by a pyranometer with a 90° FOV when the FOV is 60°.Thus, the ratio of sensor's height to the radius of the effective field of view is 1/1.732.

Satellite Data
Land surface radiometric data collected by the Chinese satellites HJ-1A/B are used for this study.Each of the HJ-1A/B satellites has two CCDs containing four shortwave bands, and has a revisit period of 2-4 days (Table 3).The data used in this study were acquired from 1st June to 30th September in 2012.Combining with the surface condition of Heihe river basin, Land cover classification data was derived from the HJ-1A/B data, by adopting a more delicate classification than the classification of the International Geosphere-Biosphere Program (IGBP) [29,30] (Figure 2), which classified crop more delicately, and neglected non-exist land cover type of Heihe river basin.The land cover data has a temporal resolution of one month.The temporal advantage of the HJ-1/CCD data not only increases the amount information, but also improves the accuracy on the land cover/use classification.In any case, there are still some inaccurate classifications, which may introduce retrieval errors which will be analyzed in Section 4.
The MODIS BRDF (MCD43A1) product supplies the weighting parameters of the Shortwave BRDF model with a temporal resolution of 8 days [31].The shortwave BRDF parameter data used in this study were acquired between 1st June 2012 and 30th September 2012.The MCD43A2 product is the quality data of the MCD43A1 product, and is used to control the quality of the MCD43A1 data in this study.

Results and Discussion
According to the details of the improved method, an experiment is performed in the Heihe river basin, and the land surface albedo is then derived from the data introduced in Section III.The analysis results are presented as follows.

Comparison of Retrieved Albedos with MODIS Product
Figure 4 compares the HJ-1A/B land surface albedo with a concurrent MODIS 500 m product over the study area during the growing season of 2012.The albedos from these two products can be seen to have identical spatial distributions in their overall patterns, except that the HJ-1A/B albedo map contains more detail about the land.For instance, the Heihe river oasis includes a variety of vegetation, each of them having a different land surface albedo.These details are unfortunately smoothed out in the MODIS albedo product, while they can be seen more clearly in the HJ-1A/B product (see Figure 4).Compared to the MODIS albedo product, the land surface albedo derived from the HJ-1A/B land surface reflectance data with 30 m spatial resolution reveals a larger dynamic range (Figure 5).Thus, the 30 m HJ-1A/B retrievals have the potential to precisely detect fine-scale natural and anthropogenic changes in the land surface albedo.
In accordance with the vegetation phenology in the Heihe river basin, the land surface albedo gradually decreases from June to August with the growth of the vegetation (and thus the increase in the amount of chlorophyll).After August, the land surface albedo tends to rise because of the vegetation withering.This seasonal variation of the albedo is well reflected in Figures 4 and 5 for both products.To examine the spatial characteristics of the albedos from the two products, a transect line and a square-shaped sub-area are selected to build regions of interest, from which the data are picked up for further analysis.The transect line is randomly graphed, and the square-shaped sub-area covers the area with most of the ground sites located.The comparison results within the transect line are shown in Figure 6a-d.and indicate that the HJ-1A/B albedo reveals a more dramatic and greater fluctuation than the MODIS product, which implies that the HJ-1A/B data has a better spatial expression capability (Figure 4a-d).To further qualitatively analyze the spatial characteristics of both products (Figure 4), some statistical measurements (Figure 6e) are calculated for these two products within the square-shaped sub-area mentioned above.The left bar for each date in Figure 6e indicates statistics for the HJ-1A/B data, and the right bar indicates those for the MODIS product.The value range of the HJ-1A/B land surface albedo is obviously larger than that of the MODIS product at any date, which further proves that the HJ-1A/B land surface albedo has better spatial expression capability than that of the MODIS product.

Comparison of the Retrieved Albedos with Ground Measurements
To validate the accuracy of the retrieved results, the albedos retrieved by the HJ-1A/B satellites between June and September 2012 at seven sites with different land cover types are organized into a temporal series.Then, a comparison between the retrieved albedo and ground measurements is performed (Figures 7 and 8).Considering that the pyranometers with 12 m altitude can cover an area within a 20 m radius at site No. 15 and those with a 6 m altitude can cover an area within a 10 m radius at the other sites, the retrieved land surface albedo with a 30 m resolution is theoretically considered to be consistent with the ground measurements in terms of spatial scale.From the comparison, we can see that a minor discrepancy between the albedos from the ground measurements and those from the HJ-1A/B data can be detected for at the No. 7, No. 15, Gobi, and Shenshawo sites, while poorer agreement is found at No. 1, No. 17, and Shidi sites.There are two possible reasons for this, which can be explained as follows: (1) Site No. 1 is located in a vegetable field, which is composed of not only vegetables but also a nearby greenhouse (Figure 3).In the classification criteria for the HJ-1A/B data, the greenhouse and vegetables are classified into the same type (crop), so the BRDF parameters of the HJ-1A/B data are actually estimated from a pseudo-MODIS pure pixel.Since the reflectance of greenhouse is generally higher than that of vegetables, the albedo of the pure pixel is greater than that measured at the ground site, who's FOV only includes the vegetables.(2) Atmospheric effects could have been overcorrected while processing the HJ-1A/B data.
Similarly, site No. 17 is located in an area covered by an orchard, around which various land cover types exist, such as short grass, roads, and bare soil.Most signals from the pyranometers are from the vicinity of the orchard with albedo generally less than that of the mixed pixel with vegetable, trees, and bare soil.
Unlike the situations of the two sites analyzed above, the Shidi site is located in the wetland, and corresponds to a pure pixel for the HJ-1A/B data and thus for MODIS.After careful analysis, we find that the BRDF parameters of the MODIS product are significantly underestimated in this area (not shown here).Since the BRDF parameters from the HJ-1A/B data are mainly derived from MODIS product, the lower values for the latter results of the retrieved albedo is underestimated as compared with the ground measurements.
In addition, a comparison between the MODIS albedo product and the in situ measurements over the seven sites is also conducted (Table 4).On average, the retrieved albedo from the HJ-1A/B data is relatively more accurate than the MODIS product for all seven sites (with a RMSE of 0.049 for the HJ-1A/B product, as compared with 0.065 for the MODIS product).Considering the mismatch in spatial coverage between the MODIS footprint and the ground measurements, as well as the similar accuracy of the retrieved albedos for both datasets at the No. 7, No. 15, Gobi, and Shenshawo sites, it can be concluded that the retrieved albedo from the HJ-1A/B data is at least no worse than the MODIS product.This accuracy is also comparable to, and even better than, the results obtained in [17] and [18].Using the method proposed in this paper, land surface albedos with a high spatio-temporal (30 m and 2-4 days) resolution can be definitively derived using the high spatial and temporal resolution HJ observations.

Analysis of Error Sources
In general, any systematic analysis of error sources for ongoing studies definitely provides a basis for further improvement in the future.In this study, the error sources associated with our newly proposed method include the following: (1) Errors resulting from inaccurate atmospheric correction.The performance of atmospheric correction can directly affect the accuracy of the derived surface reflectance.As pointed out in the literature, when appropriate aerosol optical depth (AOD) is not available, the MODIS atmospheric-corrected reflectance can be unreliable [32].In addition, the atmospheric parameter used for atmospheric correction of the HJ-1A/B band reflectance are acquired on the basis of the empirical knowledge referencing ground measurement and MODIS aerosol product, which may induce some uncertainties; (2) Uncertainties due to the variations in the MODIS BRDF within the 8-day scale.In this study, we hypothesized that the land surface state is steady during the 8 days studied (the temporal resolution of the MODIS albedo product).In fact, when several natural processes (e.g., snowfall, snowmelt, precipitation, and vegetation growth) and human activities (e.g., clearing and planting forests, sowing and harvesting of crops, and burning rangeland) occur, the land surface albedo could change rapidly [33]; (3) Misclassification.Image classification error for the HJ-1A/B data can lead to misjudgments in the classification of MODIS pixels as pure or mixed, and thus may introduce errors into the algorithm; and (4) similar to (3), uncertainties arising from unreasonable assignments for the threshold of verifying a pure pixel.
In this study, a threshold of 50% is used to discriminate a pure MODIS pixel.Is this reasonable?Theoretically, a higher threshold could result in a more accurate adjudgement of a pure pixel.While a higher percentage threshold means less pure MODIS pixels, and thus less valid samples that can be used for constructing LUT at mixed MODIS pixels.Hence, a higher value of the threshold doesn't necessarily make for better retrieving accuracy.This depends on the inhomogeneous degree of the study region.For the purposes of demonstrating the reasonability of the 50% threshold we used in this study, a test study based on the different thresholds (50%, 60%, and 70%) in a sub-area (Figure 9) of our study site is further investigated.From Figure 9, different settings for the threshold could result in obvious judgment differences.However, the retrieved white-sky albedos based on different thresholds have a similar spatial distribution (Figure 10), except for some invalid pixels in the northeastern part of albedo map when a 70% threshold is adopted (Figure 10c).Most of these invalid pixels correspond to mixed pixels (Figure 9c), where not enough pure pixels can be used to construct LUT.For quantitative analysis, some statistical measurements are calculated within a marked area in Figure 10, where the discrimination of the pure pixel is obviously affected by the different assignments of threshold.
Based on Figure 11, the statistics show that the retrieved white-sky albedos based on different thresholds have similar possibility density histograms.Analysis reveals that larger thresholds don't necessarily improve the accuracy of the result significantly in this study.Although it cannot affect the conclusion of this study, a higher threshold and more accurate classification are still recommended while considering the general applicability of the proposed method.Investigation of this issue over different study areas is one of our main works for the future.Notably, although some of the errors mentioned above may exist in the novel method in its present version, the validation and comparison results reveal that the accuracy of the new method is still comparable to that of MODIS and other studies.Further, considering its high spatio-temporal resolution, the newly developed method is truly advantageous.

Conclusions
This paper proposed an algorithm for retrieving land surface albedo by combining the HJ-1A/B land surface reflectance data and the MODIS BRDF product.The advantage of this algorithm is that it makes full use of the high spatio-temporal resolution of the HJ-1A/B data and the prior knowledge of MODIS BRDF product.In total, eight HJ-1A/B images from 1st June 2012 to 30th September 2012 were used to test the proposed method.The validation against measurements of seven ground sites with different land cover types shows that on average, the retrieved HJ-1A/B albedos were consistent with those of ground measurements, except for relatively large discrepancies in the No. 1, No. 17, and Shidi sites because of the heterogeneity of the land surface in these sites.In addition, comparisons between the MODIS albedo product and the in situ measurements, as well as between the MODIS product and the HJ-1A/B data, were also conducted.The results reveal that the accuracy of the newly developed method is not worse than that of other methods.Further, considering its high spatio-temporal (30 m and 2-4 days) resolution, the newly developed method is quite practical and thus very attractive.

Figure 1 .
Figure 1.The flowchart of the improved method for retrieving land surface albedo.

Figure 2 .
Figure 2. (Upper right) A false color image (HJ-1B CCD1 band 4, band 3, and band 2 on 28th September 2012) and (Lower panel) a land cover image based on HJ-1A/B data in September 2012 of the Heihe river basin.

Figure 3 .
Figure 3. Experimental sites of the HiWATER project in the Heihe river basin, shown on the false color composition diagram with HJ-1B CCD (4, 3, and 2 band): No. 1 site, No. 7 site, No. 15 site, No. 17 site, Gobi site, Shenshawo site, and Shidi site; and the photos of experimental sites are from: http://westdc.westgis.ac.cn.

Figure 4 .
Figure 4. Land surface albedo of the Heihe river basin from the two datasets during the 2012 growing season: HJ-1A/B (a-d) and MODIS (e-h); the HJ-1A/B results are derived by eight datasets taken between 1st June 2012 and 30th September 2012, where the cloud and singular retrieval pixels are masked out by setting an appropriate zero value; MODIS results are taken from the MCD43A1 product.

Figure 5 .
Figure 5. Histograms of the two datasets shown in Figure 4.

Figure 6 .
Figure 6.Comparison between the land surface albedos from the two datasets along a transect (black lines in Figure 4): (a) data from Figure 4 (a) and (e); (b) data from Figure 4 (b) and (f); (c) data from Figure 4 (c) and (g); (d) data from Figure 4 (d) and (h).Statistics of retrieved albedo from the square-shaped area (covering 7 ground sties) for both datasets in Figure 4 are given in (e) .

Figure 9 .
Figure 9.The statistics of pure/mixed MODIS pixels in the sub-area of the Heihe river basin based on a varying threshold: (a) 50%, (b) 60%, and (c) 70%; combing with (d) land cover type.

Figure 10 .
Figure 10.The retrieved white-sky albedo in a sub-area of the Heihe river basin on the basis of different thresholds: (a) 50%, (b) 60%, and (c) 70%.

Figure 11 .
Figure 11.Histograms of the albedos shown in the marked area of Figure 10.

Table 1 .
The weights for converting band reflectance to shortwave reflectance.

Table 2 .
Descriptions of study site and instruments set-up.

Table 3 .
Introduction of the study site.

Table 4 .
RMSE of the results in Figures7 and 8at the seven sites.