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

Abstract

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. 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).

1. 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. Section 2 and Section 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.

2. Theoretical Basis

2.1. 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

$$a_{shortwave}=\frac{{\displaystyle \underset{0}{\overset{2\mathrm{\pi }}{\int }}}{\displaystyle \underset{0}{\overset{\mathrm{\pi }2}{\int }}}{\displaystyle \underset{0}{\overset{2\mathrm{\pi }}{\int }}}{\displaystyle \underset{0}{\overset{\mathrm{\pi }/2}{\int }}}F({\mathrm{\theta }}_i,{\mathrm{\varphi }}_i)R({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\varphi }}_v,{\mathrm{\varphi }}_i)\text{cos}{\mathrm{\theta }}_v\text{sin}{\mathrm{\theta }}_v\text{cos}{\mathrm{\theta }}_i\text{sin}{\mathrm{\theta }}_{i}d{\mathrm{\theta }}_{v}d{\mathrm{\varphi }}_{v}d{\mathrm{\theta }}_{i}d{\mathrm{\varphi }}_i}{{\displaystyle \underset{0}{\overset{2\mathrm{\pi }}{\int }}}{\displaystyle \underset{0}{\overset{\mathrm{\pi }/2}{\int }}}F({\mathrm{\theta }}_i,{\mathrm{\varphi }}_i)\text{cos}{\mathrm{\theta }}_i\text{sin}{\mathrm{\theta }}_{i}d{\mathrm{\theta }}_{i}d{\mathrm{\varphi }}_i}$$

(1)

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]:

$$R({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i,t)=f_{iso}+f_{vol}{k}_{vol}({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i)+f_{geo}{k}_{geo}({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i)$$

(2)

where f_vol is the coefficient of the scattering kernel k_vol, f_geo is the coefficient of the geometrical optical kernel k_geo and f_iso is the coefficient of isotropy.

2.2. 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:

$$r=\frac{\pi Ra({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i)}{F({\mathrm{\theta }}_i,{\mathrm{\phi }}_i)}$$

(3)

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 r_hr 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:

$$r_{mr}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmr},{\mathrm{\varphi }}_{vmr},{\mathrm{\varphi }}_i)=\frac{{\displaystyle \sum _n}{r}_{hr}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmr},{\mathrm{\varphi }}_{vmr},{\mathrm{\varphi }}_i)}{n}$$

(4)

where r_mr is the reflectance of the moderate-spatial-resolution pixel, r_hr is the reflectance of high-spatial-resolution pixels in the moderate-resolution pixel, n is the number of r_hr within the moderate-resolution pixel and (θ_i, θ_vmr, Ф_vmr, Ф_i) is the illuminating and viewing geometry of r_mr.

Equation (4) holds when r_mr and r_hr 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 r_hr 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:

$$r_{hr}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmr},{\mathrm{\varphi }}_{vmr},{\mathrm{\varphi }}_i)=BRD{F}_{hr}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmr},{\mathrm{\varphi }}_{vmr},{\mathrm{\varphi }}_i)$$

(5)

where BRDF_hr is the BRDF of the land surface reflectance r_hr in the high-spatial-resolution pixels and BRDF_hr and r_hr have the same spatio-temporal resolution. The value of r_hr for any illuminating and viewing geometry (θ_i, θ_r, Ф_r, Ф_i) can be calculated from BRDF_hr.

In fact, r_mr is not equal to the average of r_sr, because of the observed time difference of r_mr and r_hr between different sensors. If the land surface is stable and its BRDF shifts only in magnitude within the observed time difference, r_mr obeys a linear relation co_h with the average of r_hr. Thus, Equation (4) can be adjusted to:

$$c{o}_h=r_{mr}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmr},{\mathrm{\varphi }}_{vmr},{\mathrm{\varphi }}_i)/\frac{{\displaystyle \sum _n}BRD{F}_{hr}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmr},{\mathrm{\varphi }}_{vmr},{\mathrm{\varphi }}_i)}{n}$$

(6)

Usually r_mr has a higher temporal resolution than BRDF_hr/r_hr, and more than one r_mr value is observed in the revisit period of r_hr. Thus, from Equation (6), we infer that more than one co_h value can be acquired with the spatio-temporal resolution of r_mr in the revisit period of r_hr. In addition, r_hr multiplies co_h to derive more than one land surface reflectance r_h with the same spatial resolution as r_hr and the same temporal resolution as r_mr in the revisit period of r_hr:

$$r_h=r_{hr}c{o}_h$$

(7)

3. 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.

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

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

3.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.

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

Band	1 (nm)	2 (nm)	3 (nm)	4 (nm)	5 (nm)	6 (nm)	7 (nm)
MODIS	620–670	841–876	459–479	545–565	1230–1250	1628–1652	2105–2155
HJ-1A/B	430–520	520–600	630–690	760–900	-	-	-

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

Band	1 (nm)	2 (nm)	3 (nm)	4 (nm)	5 (nm)	6 (nm)	7 (nm)
MODIS	620–670	841–876	459–479	545–565	1230–1250	1628–1652	2105–2155
HJ-1A/B	430–520	520–600	630–690	760–900	-	-	-

Considering the spectral range of a sensor, it is possible to obtain the shortwave albedo from the band albedo:

$$a_{shortwave}={\displaystyle \sum _l{w}_{l}a(\text{λ})}$$

(8)

where w_l is the band weight and l is the number of bands. Table 2 and Table 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).

Table 2. Weights used for NTOB HJ-1 reflectance (referencing to Liu (2012)) [32].

Band	1	2	3	4
HJ1ACCD1	−0.5878	0.6543	0.4183	0.3891
HJ1ACCD2	−0.5138	0.5175	0.4838	0.3865
HJ1BCCD1	−0.6066	0.6175	0.4703	0.3925
HJ1BCCD2	−0.5149	0.5403	0.4592	0.3898

Table 2. Weights used for NTOB HJ-1 reflectance (referencing to Liu (2012)) [32].

Band	1	2	3	4
HJ1ACCD1	−0.5878	0.6543	0.4183	0.3891
HJ1ACCD2	−0.5138	0.5175	0.4838	0.3865
HJ1BCCD1	−0.6066	0.6175	0.4703	0.3925
HJ1BCCD2	−0.5149	0.5403	0.4592	0.3898

Table 3. Weights used for NTOB MODIS reflectance (referencing to Liang et al. (1999)) [7].

Band	1	2	3	4	5	6	7	Const
MODIS	0.3973	0.2382	0.3489	−0.2655	0.1604	−0.0138	0.0682	0.0036

Table 3. Weights used for NTOB MODIS reflectance (referencing to Liang et al. (1999)) [7].

Band	1	2	3	4	5	6	7	Const
MODIS	0.3973	0.2382	0.3489	−0.2655	0.1604	−0.0138	0.0682	0.0036

3.2. 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 R_hj/BRDF BRDF_hj that is the closest time and location to each image of MODIS reflectance R_modis 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 R_h is calculated from the coefficient multiplied by the corresponding reflectance (Equation (7)). The equation for R_h is:

$$R_h=R_{hj}\text{x(}{R}_{modis}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmodis},{\mathrm{\varphi }}_{vmodis},{\mathrm{\varphi }}_i,T_{modis})/\frac{{\displaystyle \sum _{n}BRD{F}_{hj}({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vmodis},{\mathrm{\varphi }}_{vmodis},{\mathrm{\varphi }}_i,t_{hj})}}{n}\text{)}$$

(9)

where R_modis is the reflectance observed by MODIS at the illuminating and viewing geometry (θ_i, θ_vmodis, Ф_vmodis, Ф_i, T_modis); BRDF_hj is the known HJ-1A/B BRDF in the 30-m subpixel of R_modis; t_hj is the time when R_hj is obtained; T_modis is the time when R_modis is observed and n denotes the number of HJ-1A/B pixels within the MODIS pixel.

3.3. 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):

$$c=R_h({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vh},{\mathrm{\phi }}_{vh},{\mathrm{\phi }}_i,t_h)/BRD{F}_m({\mathrm{\theta }}_i,{\mathrm{\theta }}_{vh},{\mathrm{\phi }}_{vh},{\mathrm{\phi }}_i,T_m)$$

(10)

where R_h is the downscaled daily 30-m land surface reflectance; BRDF_m is the calculated reflectance that is the value of the MODIS BRDF at (θ_i, θ_vh, Ф_vh, Ф_i) and is the same as R_h; T_m indicates the 8-day period of the MODIS BRDF product and t_h is the frequency of HJ-1A/B observations during T_m. Thus, the daily 30-m BRDF can be calculated using the following equation:

$$\begin{array}{c}BRD{F}_h({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i,t_h)=c{R}_m({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i,T_m)=c{f}_{iso}(T_m)\\ +c{f}_{vol}(T_m)k_{vol}({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i,T_m)+c{f}_{geo}(T_m)k_{geo}({\mathrm{\theta }}_i,{\mathrm{\theta }}_v,{\mathrm{\phi }}_v,{\mathrm{\phi }}_i,T_m)\end{array}$$

(11)

After the daily BRDF for each 30-m pixel has been retrieved, BRDF_h is integrated to derive the daily 30-m albedo on the basis of Equation (1).

3.4. 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, co_n, from Liang et al. [8].

For eliminating error from a singular value, averaging all co_n in LUT can provide a uniform ratio co_mixed:

$$c{o}_{mixed}=\frac{1}{n}{\displaystyle \sum _{n}c{o}_n}$$

(12)

As the daily 30-m reflectance r_hj and the corresponding co_mixed are known, the albedo Albedo_HJ for a 30-m pixel within a mixed MODIS pixel is:

$$Albedo\_HJ=\frac{r_{hj}}{c{o}_{mixed}}$$

(13)

4. 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.

4.1. 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.

Figure 2. Representative images of the study area (Upper left: HJ-1B true-color image on 25 September 2012; lower right: HJ-1A/B land cover; lower left: HJ-1B false-color image on 25 September 2012).

Figure 2. Representative images of the study area (Upper left: HJ-1B true-color image on 25 September 2012; lower right: HJ-1A/B land cover; lower left: HJ-1B false-color image on 25 September 2012).

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

Site	No. 1	No. 7	No. 15	No. 17	Gobi	Shenshawo
Location (degree)	100.3582E, 38.8932N	100.3652E, 38.8767N	100.3722E, 38.8555N	100.3697E, 38.8451N	100.3042E, 38.9149N	100.4933E, 38.7891N
LC	Vegetable	Corn	Corn	Orchard	Gobi	Desert
Radius of Coverage	10 m	10 m	20 m	10 m	10 m	10 m

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

Site	No. 1	No. 7	No. 15	No. 17	Gobi	Shenshawo
Location (degree)	100.3582E, 38.8932N	100.3652E, 38.8767N	100.3722E, 38.8555N	100.3697E, 38.8451N	100.3042E, 38.9149N	100.4933E, 38.7891N
LC	Vegetable	Corn	Corn	Orchard	Gobi	Desert
Radius of Coverage	10 m	10 m	20 m	10 m	10 m	10 m

4.2. 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].

5. 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.

5.1. 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.

Figure 3. Downscaled reflectance (a) and MODIS reflectance (b) of the Heihe river basin on 30 September 2012.

Figure 3. Downscaled reflectance (a) and MODIS reflectance (b) of the Heihe river basin on 30 September 2012.

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”).

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”).

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² 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.

Figure 5. Time series of the downscaled daily 30-m reflectance and HJ-1A/B reflectance at six sites (a–f); “Daily” denotes the downscaled daily 30-m reflectance; “HJ-1A/B” denotes the HJ-1A/B reflectance.

Figure 5. Time series of the downscaled daily 30-m reflectance and HJ-1A/B reflectance at six sites (a–f); “Daily” denotes the downscaled daily 30-m reflectance; “HJ-1A/B” denotes the HJ-1A/B reflectance.

From Figure 5, the daily reflectance has a higher temporal resolution than the HJ-1A/B reflectance. Figure 3 and Figure 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.

5.2. 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. 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 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 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).

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 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.

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

Site	No. 1	No. 7	No. 15	No. 17	Gobi	Shenshawo	All
Proposed/HJ-1A/B	2	1.96	1.88	2.45	2.08	2	>2.0

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

Site	No. 1	No. 7	No. 15	No. 17	Gobi	Shenshawo	All
Proposed/HJ-1A/B	2	1.96	1.88	2.45	2.08	2	>2.0

5.3. 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).

Figure 8. Time series (a–f, left) and scatter plots (a–f, right) of derived daily 30-m land surface albedo and measurements made on the ground at six sites. ”Daily” denotes the derived daily 30-m land surface albedo; ‘Observed’ denotes ground measurements.

Figure 8. Time series (a–f, left) and scatter plots (a–f, right) of derived daily 30-m land surface albedo and measurements made on the ground at six sites. ”Daily” denotes the derived daily 30-m land surface albedo; ‘Observed’ denotes ground measurements.

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).

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

Product	Temporal Resolution (day)	Spatial Resolution(m)	Total RMSE	Total Absolute Error
MODIS	8	500	0.029	0.023
Proposed method	1	30	0.028	0.022

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

Product	Temporal Resolution (day)	Spatial Resolution(m)	Total RMSE	Total Absolute Error
MODIS	8	500	0.029	0.023
Proposed method	1	30	0.028	0.022

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

Method	Temporal Resolution (day)	Spatial Resolution (m)
Shuai or Franch	2–4	30
Proposed method	1	30

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

Method	Temporal Resolution (day)	Spatial Resolution (m)
Shuai or Franch	2–4	30
Proposed method	1	30

5.4. 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.

Figure 9. Classification map of “pure/mixed” MODIS pixels; (a) 50% threshold, (b) 60% threshold.

Figure 9. Classification map of “pure/mixed” MODIS pixels; (a) 50% threshold, (b) 60% threshold.

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 Figure 10b is difficult to detect visually. Thus, a comparison of Figure 10a and Figure 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).

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 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. 2-D scatter plot of daily albedo values obtained using the proposed method within critical MODIS pixels.

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

6. 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.

Although the proposed method effectively improves the spatio-temporal resolution and accuracy of albedo, there are still some deficiencies. As well as inaccurate atmospheric correction caused by a lack of, and inaccuracy in, atmospheric parameters (e.g., aerosol optical depth and water vapour), uncertainties due to the temporal-scale mismatch of the MODIS BRDF (8-day scale) and HJ data (2–4-day scale) may introduce error. In the course of downscaling the land surface reflectance, HJ-1A/B BRDF and reflectance data may not capture the variation of land cover for every day because of the 2–4-day temporal resolution, but the daily MODIS reflectance can achieve this. When natural and social emergencies occur, albedo may change rapidly. In this situation, the HJ-1A/B and MODIS information may be inconsistent, possibly causing errors.

In conclusion, the proposed method first retrieves the daily albedo at 30 m, but some work needs to be done to further improve the proposed method in the future; one is how to acquire the atmospheric parameters corresponding to the HJ-1A/B data for atmospheric correction; another one is how to eliminate the error from the temporal mismatch among the multisource data (e.g., HJ-1A/B data, MODIS reflectance, MODIS BRDF).