A Multi-Scale Validation Strategy for Albedo Products over Rugged Terrain and Preliminary Application in Heihe River Basin , China

The issue for the validation of land surface remote sensing albedo products over rugged terrain is the scale effects between the reference albedo measurements and coarse scale albedo products, which is caused by the complex topography. This paper illustrates a multi-scale validation strategy specified for coarse scale albedo validation over rugged terrain. A Mountain-Radiation-Transfer-based (MRT-based) albedo upscaling model was proposed in the process of multi-scale validation strategy for aggregating fine scale albedo to coarse scale. The simulated data of both the reference coarse scale albedo and fine scale albedo were used to assess the performance and uncertainties of the MRT-based albedo upscaling model. The results showed that the MRT-based model could reflect the albedo scale effects over rugged terrain and provided a robust solution for albedo upscaling from fine scale to coarse scale with different mean slopes and different solar zenith angles. The upscaled coarse scale albedos had the great agreements with the simulated coarse scale albedo with a Root-Mean-Square-Error (RMSE) of 0.0029 and 0.0017 for black sky albedo (BSA) and white sky albedo (WSA), respectively. Then the MRT-based model was preliminarily applied for the assessment of daily MODerate Resolution Imaging Spectroradiometer (MODIS) Albedo Collection V006 products (MCD43A3 C6) over rugged terrain. Results showed that the MRT-based model was effective and suitable for conducting the validation of MODIS albedo products over rugged terrain. In this research area, it was shown that the MCD43A3 C6 products with full inversion algorithm, were generally in agreement with the aggregated coarse scale reference albedos over rugged terrain in the Heihe River Basin, with the BSA RMSE of 0.0305 and WSA RMSE of 0.0321, respectively, which were slightly higher than those over flat terrain.


Introduction
Land surface shortwave albedo is defined as the fraction of incident solar irradiance reflected by Earth's surface over the shortwave band (0.3-3 µm) in the whole solar spectrum [1].It is a key climate-regulating parameter that determines the amount of solar radiation absorbed by the land surface at regional and global scales [2,3].Remote sensing satellites provide a practical method to estimate land surface albedos because of their large spatial scale coverage and a high revisit frequency [4].However, the retrieved albedos suffer from large uncertainties due to the inherent complexity of the physical processes and their parameterization of retrieval algorithms [5].Thus, it is critical to evaluate the performance of the retrieved albedos prior to their wide application.
Many scientists have focused on assessing the accuracy of albedo products in recent decades over flat and homogeneous land surfaces.Taking the MODIS albedo products validation as an example, the validation results showed that the MODIS albedo products displayed high accuracy with an uncertainty (Root-Mean-Square-Error, RMSE) below 0.03 at the snow-free land covers and 0.07 at the snow-covered land surface, respectively, when the validation activities occurred in the homogeneous land surface or in sites with high spatial representativeness [6][7][8][9][10][11][12][13][14][15][16][17][18][19].Generally, the albedo product can be assessed by direct comparison with in situ albedos at the sites where the land surface is sufficiently homogeneous [6,9,10,20].However, in the heterogeneous land surface, the in situ albedo cannot be directly compared with albedo products because of the scale mismatching between the in situ albedo and the albedo products, unless that in situ albedo can be considered with high spatial representativeness over the sampled area [10,16,20,21].The scale mismatching will result in about 15% disagreement between the MODIS albedo and in situ albedo [6,9,22].As the sample sites with limited spatial representativeness, multi-points albedo observing is generally adapted to capture the spatial distribution characteristics of the albedo over the sampled area.The simplest and most efficient method is to average these albedos within the area and as the reference truth to compare with the albedo products [23].Alternatively, the multi-scale validation strategy provides a solution to deal with the scale mismatching over a heterogeneous land surface by introducing fine scale albedo products (e.g., the Enhanced Thematic Mapper plus (ETM+) or China HJCCD (HJ) albedo) as an upscaling bridge between in situ albedo and coarse scale satellite albedo products (e.g., MODIS and Global LAnd Surface Satellite (GLASS) albedo) [8,17,24,25].In situ albedos are used to calibrate the fine scale albedo.Then, the calibrated fine scale albedos are aggregated to the coarse scale and for albedo validation.Previous studies have indicated that the upscaling of albedo from fine scale to coarse scale is highly linear over flat terrains [8,16,17,[24][25][26].Therefore, in the case of flat terrain, the linear weighted average model was considered as a good performance model for upscaling the fine scale albedo to a coarse scale.
As a special heterogeneous land surface, the topography has vast effects on the land surface albedo [24,26].Topographic slope, aspect, shadow, and solar location influence albedo values and their spatial distribution when compared with that over flat terrain [13,[27][28][29][30].The coarse scale albedo decreased with the increase of the slope facing away from the sun and increased when facing toward the sun [30,31], and generally showed larger values over the slope facing toward the sun than that facing away from the sun, especially, in the shadowing case [27,30,32].The complex topography leads to the intensive scale effects on albedo products among different spatial resolutions over the rugged terrain [17,33].However, neglecting the scale effect caused by the complex topography in albedo products results in unreliable validation results [33][34][35].Peng et al., (2014) assessed the MODIS products by using the multi-scale validation strategy with the HJ albedo as the bridge to aggregate the in situ albedo linearly to the MODIS pixel scale.The uncertainty distribution analysis showed that the largest scaling uncertainty was at the pixels over rugged terrain and its uncertainty of MODIS was as high as 0.07, when neglecting the scale effects at the upscaling progress over rugged terrain [17].Therefore, neither the direct comparison nor the linear upscaling model in the multi-scale validation strategy were suitable for the coarse scale albedo products validation over rugged terrain [33].The albedo spatial scale issue caused by topography should be emphasized in the multi-scale validation strategy over rugged terrain.
The objective of this paper was to develop an upscaling method in the procedure of the multi-scale validation strategy for albedo products validation over rugged terrain.Simulated data with different mean slopes and solar zenith angle over nine Digital Elevation Models (DEM) were used to validate the albedo upscaling method.Based on the proposed upscaling method, the aggregated HJ albedo, which had been validated by the in situ albedo over the Heihe River Basin, was used as the reference truth for the MODIS albedo products preliminary validation.The paper is organized as follows: Section 2 describes the multi-scale validation strategy, including the upscaling method and the fine scale albedo products retrieval algorithm; Section 3 describes the experimental area and validation dataset; Section 4 shows the performance of the albedo upscaling model and the preliminary validation results for MODIS albedo products.The discussions are summarized in Section 5. Finally, a brief conclusion is drawn in Section 6.

Multi-Scale Validation Procedure over Rugged Terrain
The general multi-scale validation strategy includes three key procedures over rugged terrain.The first one is the retrieval and accurate evaluation of fine scale albedo.The second is the calibration of fine scale albedo products.Finally, the third process is to scale up the fine scale albedo products to the coarse scale [17,24].The multi-scale validation strategy over rugged terrain has similar procedures than those over flat terrain, which is shown in Figure 1.The DEM was used here to calculate the topographic factors (e.g., slope, aspect), and to couple with fine scale reflectance to retrieve fine scale albedo.The fine scale albedo was assessed by direct comparison with in situ albedos and was calibrated for reducing its uncertainties.An albedo upscaling model, which was based on the mountain radiation transfer theory (MRT-based albedo upscaling model), was proposed to aggregate the fine scale albedo to a coarse spatial scale.Consequently, the aggregated albedos could be directly compared with the coarse scale albedo products.To implement a successful multi-scale validation strategy over rugged terrain, two key issues should be solved including the albedo upscaling method and the fine scale albedo-generated algorithm on sloping surfaces.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 23 reference truth for the MODIS albedo products preliminary validation.The paper is organized as follows: Section 2 describes the multi-scale validation strategy, including the upscaling method and the fine scale albedo products retrieval algorithm; Section 3 describes the experimental area and validation dataset; Section 4 shows the performance of the albedo upscaling model and the preliminary validation results for MODIS albedo products.The discussions are summarized in Section 5. Finally, a brief conclusion is drawn in Section 6.

Multi-Scale Validation Procedure over Rugged Terrain
The general multi-scale validation strategy includes three key procedures over rugged terrain.The first one is the retrieval and accurate evaluation of fine scale albedo.The second is the calibration of fine scale albedo products.Finally, the third process is to scale up the fine scale albedo products to the coarse scale [17,24].The multi-scale validation strategy over rugged terrain has similar procedures than those over flat terrain, which is shown in Figure 1.The DEM was used here to calculate the topographic factors (e.g., slope, aspect), and to couple with fine scale reflectance to retrieve fine scale albedo.The fine scale albedo was assessed by direct comparison with in situ albedos and was calibrated for reducing its uncertainties.An albedo upscaling model, which was based on the mountain radiation transfer theory (MRT-based albedo upscaling model), was proposed to aggregate the fine scale albedo to a coarse spatial scale.Consequently, the aggregated albedos could be directly compared with the coarse scale albedo products.To implement a successful multiscale validation strategy over rugged terrain, two key issues should be solved including the albedo upscaling method and the fine scale albedo-generated algorithm on sloping surfaces.

MRT-Based Albedo Upscaling Model over Rugged Terrain
Coarse scale albedo can be defined as the ratio of the reflected solar radiant flux and incident solar radiant flux.Assuming that the coarse scale pixel is horizontal overall and the topographic effects between the coarse scale pixels can be ignored, the incident solar radiant flux of a coarse scale pixel can be expressed as the sum of direct and diffuse radiant flux at the micro-slope.Therefore, the coarse scale reflected radiant flux was calculated by summing up the micro-slope reflected radiant flux.If the micro-slope albedo is known, the coarse scale albedo can be expressed as:

MRT-Based Albedo Upscaling Model over Rugged Terrain
Coarse scale albedo can be defined as the ratio of the reflected solar radiant flux and incident solar radiant flux.Assuming that the coarse scale pixel is horizontal overall and the topographic effects between the coarse scale pixels can be ignored, the incident solar radiant flux of a coarse scale pixel can be expressed as the sum of direct and diffuse radiant flux at the micro-slope.Therefore, the coarse scale reflected radiant flux was calculated by summing up the micro-slope reflected radiant flux.If the micro-slope albedo is known, the coarse scale albedo can be expressed as: where A c is the coarse scale albedo; φ ↑ c , φ ↓ c are the coarse scale pixel's reflected and incident solar radiant flux, respective; N is the amount of micro-slope within the corresponding area of coarse scale pixel; ∂ k and dφ T↓ k are the albedo and incident radiant flux of the kth micro-slope, respectively.φ ↓ sc , φ ↓ dc are the coarse scale direct solar radiant flux and diffuse radiant flux, respectively.They can be expressed as: where E s , E d are the incident direct solar irradiance and diffuse irradiance on the horizontal plane; Area C is the area of the coarse scale pixel, which can be expressed as the sum of the projected area (dA th ) of the micro-slope on horizontal plane.
For each micro-slope, dφ T↓ k can be expressed as the sum of direct radiant flux, diffuse radiant flux, and terrain radiant flux reflected by adjacent terrain (Figure 2) [32,36].The incident direct and diffuse solar radiant flux of a coarse scale pixel can be expressed as the sum of the micro-slope's solar direct and diffuse radiant flux on the projected horizontal plane, respectively.Hence, Equation (1) can be displayed as: , are the coarse scale direct solar radiant flux and diffuse radiant flux, respectively.They can be expressed as: (2) where E E are the incident direct solar irradiance and diffuse irradiance on the horizontal plane; C Area is the area of the coarse scale pixel, which can be expressed as the sum of the projected area ( th dA ) of the micro-slope on horizontal plane.
For each micro-slope, T k d φ ↓ can be expressed as the sum of direct radiant flux, diffuse radiant flux, and terrain radiant flux reflected by adjacent terrain (Figure 2) [32,36].The incident direct and diffuse solar radiant flux of a coarse scale pixel can be expressed as the sum of the micro-slope's solar direct and diffuse radiant flux on the projected horizontal plane, respectively.Hence, Equation ( 1) can be displayed as: The superscript 'T' means that the pixel is on a sloping surface; the subscript 'k' means that it is the th k fine scale micro-slope; , , are the reflected direct solar radiant flux, diffuse sky flux, and terrain radiant flux of the th k micro-slope, respectively.They can be written as: The superscript 'T' means that the pixel is on a sloping surface; the subscript 'k' means that it is the kth fine scale micro-slope; dφ T↑ sk , dφ T↑ dk , dφ T↑ ak are the reflected direct solar radiant flux, diffuse sky flux, and terrain radiant flux of the kth micro-slope, respectively.They can be written as: where cos i sk , cos θ sk are the cosine of the relative incident angle on sloping surface and the incident zenith angle on horizontal surface; Θ k is the binary coefficient in the direction of the sun which is used to show whether the pixel is shadowed by the terrain; E h sk , E h dk are the incident solar direct and diffuse irradiance on horizontal surface; dA tk is the area of the kth micro-slope surface; V d is the sky-view factor, which is the sky portion seen from a specific surface [37]; and V t is the portion of adjacent terrain seen from a surface (the methods for calculating these parameters are listed in Appendix A). α k is the slope of kth micro-slope surface; a T bk , α T wk are the micro-slope's directional-hemisphere reflectance (also as black-sky albedo, BSA) and bi-hemispherical reflectance (white-sky albedo, WSA) on a sloping surface, respectively.Inserting Equations ( 6)-( 9) into (5), the coarse scale albedo can be rewritten as: We define a parameter S as the sky diffuse ratio factor, which is the proportion of diffuse solar irradiance to the total radiation.It can be calculated by the ratio of the diffuse radiance on a horizontal surface and the downward global solar radiance and exposes values between 0 and 1.
Inserting Equation ( 11) into (10), the coarse spatial scale albedo can be simplified as: Similar to the descriptions of the blue sky albedo of MODIS [38,39], the aggregated coarse scale albedo can also be approximated through a linear combination of BSA and WSA, weighted by the sky diffuse ratio factor.Thus, when S equals zero, it means that it has no diffuse skylight and the value of A C can be considered as the coarse scale BSA (BSA C ).However, when S is 1, the A C is the coarse scale WSA (WSA C ): It is obvious that the BSA and WSA have intense topographic effects as they depend on topographic factors and the solar incident angle [33].Thus, they have scale effects in the validation of albedo over rugged terrain.Furthermore, Equations ( 13) and ( 14) provide the albedo upscaling method and function as the reference truth in albedo validation over rugged terrain.

Fine Scale Albedo Retrieval Algorithm over Rugged Terrain
Currently, the Bidirectional Reflectance Distribution Function-based (BRDF-based) albedo retrieval algorithms cannot be applied directly to fine scale albedo retrieval because of the lacking of enough fine scale multi-angle observations [40,41].A feasible method is the direct retrieval algorithm such as the Angular Bin (AB) algorithm, which estimates the surface broadband albedo based on a single-date/angular observation [40][41][42][43].In the AB algorithm, the incident and observing hemisphere were divided into several angular bins.The POLarization and Directionality of Earth Reflectance-BRDF (POLDER-BRDF) data were used as prior knowledge in this algorithm for extracting anisotropy reflectance information to build a Look-Up Table (LUT) of the land surface albedo at each angular bin [41].Then, the multi-variant linear regression models were established at each angular bin, which linked the narrowband surface directional reflectance with broadband albedo, specifically, the shortwave WSA and BSA corresponding to the solar angle at local noon [40,41].The multi-variant linear regression relationship can be expressed as: where α AB represents the land surface fine scale albedo including BSA and WSA; C i (θ s , ϕ s , θ v , ϕ v ) is the regression coefficient; and ρ i (θ s , ϕ s , θ v , ϕ v ) is the surface directional reflectance at band i, which are the functions of solar/view angles (θ s , ϕ s , θ v , ϕ v ) are the solar zenith angle, solar azimuth angle, view zenith angle, and view azimuth angle).The AB algorithm has the advantages of simple computation, fewer input parameters, and consideration of surface bidirectional and spectral characteristics [41].Preliminary validation of the AB algorithm showed good applicability for albedo retrieval with an absolute error of 0.009 for vegetation, 0.012 for soil, and 0.030 for snow/ice [41][42][43].However, the current AB algorithm had to be improved before being applied over rugged terrain as the regression models and the LUT were built on horizontal land surface.Therefore, the coefficient and the LUT were not suitable for the albedo retrieval over rugged terrain.One accessible approach to improve the AB algorithm over rugged terrain was to rebuild the regression models and LUT on the sloping surface.Under the assumptions that the BRDF shape was the same over the land cover, the BRDF shape for the slope was the BRDF for the rotated angles.Sloping surface incident/observation hemispheres were re-divided and rotated to the sloping surface at the sloping coordinate system.Thus, it can be re-written on a sloping surface as: where α ABT is the fine scale sloping surface albedo, i s , φ s , i v , φ v are the relative solar zenith angle, solar azimuth angle, view zenith angle and azimuth angle on the sloping coordinate system, respectively.
is the land surface reflectance on the sloping surface.C i (i s , φ s , i v , φ v ) is the regression coefficient on a sloping coordinate system.The sloping surface reflectance can be retrieved by the Coupled-BRDF mountain radiation transfer model, which was developed by Wen in 2015 [44].

Simulated Coarse Scale and Fine Scale Albedos
The simulated coarse scale albedos were used as reference data for comparison with the aggregated coarse scale albedos, which were upscaled from a simulated fine scale albedo by using the MRT-based upscaling model.The BSA can be simulated by integrating directional reflectance over the exitance hemisphere, and the WSA can be obtained by integrating the directional reflectance over all viewing and irradiance reflectance directions.The coarse scale directional reflectance can be simulated according to the radiosity theory [33,34,45], which coupled with the DEM and micro-slope reflectance.
In the coarse scale reflectance simulation, the micro-slope or the fine scale directional reflectance can be directly simulated by the PROSAIL model [46], which couples the PROSPECT leaf optical properties model [47] with the SAIL canopy reflectance model [48] and has been widely validated and applied to reflectance modeling studies [49].For directional reflectance simulation over rugged terrain, the leaf inclination distribution function (LIDF) was assumed as a spherical type, and the incident/observation geometrics were corrected to the sloping coordinate system to assess the photon path length alteration.Table 1 illustrates the parameters of the inputted PROSAIL model for land surface reflectance simulation over the sloping surface.Nine DEMs (Table 2) were generated with different Gaussian height distributions to provide the various slope and aspect of micro-slope (shown in Figure 3).Specifically, the nine DEMs had the same reference template, which was simulated by a Gaussian random distribution model with one unit mean elevation and 0.25 unit standard error within the area of 100 units × 100 units.By linking with the real distance of 30 m, the simulated micro-slope could be considered with the 30 m resolution [33].Additionally, the vertical elevation was exaggerated in 1, 10, and 20, respectively.Through appropriate exaggeration in the elevation and different Gaussian smoothing filtering parameters, the micro-slopes were generated with different slope and aspect.However, only the central part 510 m × 510 m (17 × 17 grid cells) of the above DEM was considered as a coarse scale pixel coverage to avoid the ambiguous calculation errors at the edge of the DEM.The mean slopes listed in Table 2 were defined as the average slope of the grids within the 17 × 17 grid cells.Since the BSA depends on the solar geometry (including SZA and SAA), the SZA varied from 0 • to 60 • with a 10 • interval, and SAA varied from 0 • to 360 • with a 30 • interval.Therefore, 819 BSAs and nine WSAs from nine DEM files were simulated.

In Situ Albedo Measurements
The upper stream of the Heihe River Basin (HRB) was selected as the study area for coarse scale albedo validation over rugged terrain.HRB is a typical inland river in the arid region of northwest China (97.1°E-102.0°E,37.7°N-42.7°N).This region was selected as the core experimental watershed of Hi-WATER, because of the abundant accumulations of long temporal scale records since 1995 [50,51].There are mountains in this region where the altitude ranges from 1025 m to 5076 m.The dominated land covers are cropland, Gobi Desert, grassland, forest, and high-latitude meadows.
A prototype watershed observation system has been established since 2009 [52].Fifteen Automatic Weather Stations (AWS) were mounted in this observation system (Figure 4), which recorded the essential parameters every ten minutes including upwelling solar shortwave radiance (USR), and downwelling solar shortwave radiance (DSR).In situ albedo can be calculated as the ratio of USR and DSR at the local solar noon (11:30-12:30).Figure 4 illustrates the study area and the location of the selected AWSs.
The 15 AWSs in HRB are temporary sites where the measurements were only recorded at the periods of the Hi-WATER experiments.Two CNR pyranometers were mounted back to back between 2 m and 2.5 m above the top of the canopy at every AWS.The measurement footprints were circular and 80% of the signal came from a region with a diameter between 34.8 m and 43.5 m, which could be easily calculated as per Sailor's work [53].The measurements recorded at 15 AWSs were considered with the high representativeness of mean albedo within the 30 m HJ pixel.Table 3 summarizes the details of the albedo measured sites used for validation.

In Situ Albedo Measurements
The upper stream of the Heihe River Basin (HRB) was selected as the study area for coarse scale albedo validation over rugged terrain.HRB is a typical inland river in the arid region of northwest China (97.1 • E-102.0 • E, 37.7 • N-42.7 • N).This region was selected as the core experimental watershed of Hi-WATER, because of the abundant accumulations of long temporal scale records since 1995 [50,51].There are mountains in this region where the altitude ranges from 1025 m to 5076 m.The dominated land covers are cropland, Gobi Desert, grassland, forest, and high-latitude meadows.
A prototype watershed observation system has been established since 2009 [52].Fifteen Automatic Weather Stations (AWS) were mounted in this observation system (Figure 4), which recorded the essential parameters every ten minutes including upwelling solar shortwave radiance (USR), and downwelling solar shortwave radiance (DSR).In situ albedo can be calculated as the ratio of USR and DSR at the local solar noon (11:30-12:30).Figure 4 illustrates the study area and the location of the selected AWSs.
The 15 AWSs in HRB are temporary sites where the measurements were only recorded at the periods of the Hi-WATER experiments.Two CNR pyranometers were mounted back to back between 2 m and 2.5 m above the top of the canopy at every AWS.The measurement footprints were circular and 80% of the signal came from a region with a diameter between 34.8 m and 43.5 m, which could be easily calculated as per Sailor's work [53].The measurements recorded at 15 AWSs were considered with the high representativeness of mean albedo within the 30 m HJ pixel.Table 3 summarizes the details of the albedo measured sites used for validation.Ten AWSs were mounted in the grassland, three AWSs were located in the cropland, one site in the desert, and one site in the forest.The slope of the AWSs and the mean slope of coarse scale pixels (500 m in this paper) were obtained from the DEM were listed in Table 3.According to the mean slope, if the mean slope was less than 5°, the coarse scale pixel was thought of as being of a gentle slope.If the mean slope was greater than 5° and less than 10°, the coarse scale pixel was a relatively rugged terrain.When the mean slope was greater than 10°, the coarse scale pixel was considered as steep terrain.The DEM were collected from the Shuttle Radar Topography Mission (SRTM) DEM data [54,55], and had a 30 m spatial resolution and UTM projection [56].Slope, aspect, shaded factor, and other terrain parameters were derived from the DEM referenced to the algorithm in Dozier's work [37].

Satellite Imagery
The 500 m Collection V006 MODIS albedo product (MCD43A3 C6) was selected as the coarse scale albedo products for validation.The semi-empirical, kernel-driven BRDF model was the primary algorithm to retrieve surface BRDF and albedo at 16-day time periods [57].A back-up algorithm was  Ten AWSs were mounted in the grassland, three AWSs were located in the cropland, one site in the desert, and one site in the forest.The slope of the AWSs and the mean slope of coarse scale pixels (500 m in this paper) were obtained from the DEM were listed in Table 3.According to the mean slope, if the mean slope was less than 5 • , the coarse scale pixel was thought of as being of a gentle slope.If the mean slope was greater than 5 • and less than 10 • , the coarse scale pixel was a relatively rugged terrain.When the mean slope was greater than 10 • , the coarse scale pixel was considered as steep terrain.The DEM were collected from the Shuttle Radar Topography Mission (SRTM) DEM data [54,55], and had a 30 m spatial resolution and UTM projection [56].Slope, aspect, shaded factor, and other terrain parameters were derived from the DEM referenced to the algorithm in Dozier's work [37].

Satellite Imagery
The 500 m Collection V006 MODIS albedo product (MCD43A3 C6) was selected as the coarse scale albedo products for validation.The semi-empirical, kernel-driven BRDF model was the primary algorithm to retrieve surface BRDF and albedo at 16-day time periods [57].A back-up algorithm was employed for deriving the albedo in a situation of insufficient angular sampling because of cloud cover or orbital constraints [28].The MCD43A3 products provided shortwave broadband black-sky albedo (BSA) and white-sky albedo (WSA) at a 500 m resolution [38].
The fine scale remote sensing data were collected from the China HJ-1/CCD sensor which provides an opportunity to record earth land surface reflectance with high spatial resolution and a broad coverage in China [44].The HJ-1/CCD can record 4-band images with a 30 m spatial scale and a revisiting circle of less than two days [17].The HJ images were geo-rectified for matching with the DEM file.

Assessment of the MRT-Based Upscaling Model
A high-quality albedo upscaling model is of great importance for the application of multi-scale validation strategy.The accuracy of the MRT-based albedo upscaling model was assessed through a comparison of the aggregated fine scale albedo with the reference coarse scale albedo.Moreover, to show the performance of the MRT-based upscaling method, the line weighted average model was also used to scale up the fine scale albedo to a coarse spatial scale, which is the most commonly-used model in the multi-scale validation method over flat terrain [8,9,24,58].The bias, root-mean-square error (RMSE), mean absolute percent error (MAPE), and coefficient of determination (R 2 ) were used to show the accuracy of the aggregated coarse scale albedos, which are expressed as follows: where p i , p i are the albedo observation and average albedo observation; o i , o i are the reference albedo measurement and average of reference albedo measurement; and n represents the amount of simulated albedo measurements.

Accuracy of Using MRT-based Albedo Upscaling Model
Figure 5 shows the scatter plots between the reference coarse scale albedo and aggregated coarse scale albedo as well as the histogram of bias distributions.The aggregated coarse scale BSAs showed significant agreements with the reference simulated coarse scale albedo with an RMSE of 0.0029, Bias of −0.0001, MAPE of 0.94%, and R 2 of 0.9962 (Figure 5A).The histogram showed that the bias reflected a normal distribution with the maximum bias less than 0.015 (Figure 5B). Figure 5C shows that the aggregated coarse scale WSAs had a similar accuracy with the BSAs.The bias was 0.0007, the MAPE was 0.79%, the RMSE was less than 0.0017, and the R 2 is as high as 0.9984.The bias in Figure 5D showed that the maximum bias among the nine simulated WSAs was less than 0.003 following the change of the mean slope.Overall, these results showed that the aggregated coarse scale albedos had fewer discrepancies with the reference coarse scale albedos.The BSA depends on the incident solar geometry and is intensely affected by complex terrain, therefore, the accuracy of aggregated coarse scale albedo varied with different SZAs and terrain slopes.Figure 6 shows the comparison between the aggregated coarse scale albedo and simulated coarse scale albedo at the solar zenith angles from 0 • to 60 • with an interval of 20 • .It was obvious that the aggregated coarse scale BSAs had the significant agreement with the reference coarse scale BSAs at the four selected SZAs.With the increase of SZA, the RMSE values showed a slight increase from 0.0024 at SZA of 0 • to 0.0025 at SZA of 40 • .The Bias, MAPE, and R 2 had a little change when the solar zenith angle varied from 0 • to 40 • .Even though the SZA increased to 60 • , the aggregated coarse scale BSAs still had significant agreement with the reference simulated coarse scale albedo with a bias of 0.0021, MAPE of 1.51%, RMSE of 0.0049, and R 2 of 0.9971.The bias distribution histogram also demonstrated that the maximum bias of the two types of BSAs was less than 0.015 at the four selected SZAs.
Figure 7 shows the comparison between the simulated coarse scale BSAs and reference coarse scale BSAs over four selected DEM types with different mean slopes.The mean slope of the DEM in Figure 7A was lower than 10 • , which was considered as the flat terrain.The results showed that the aggregated coarse scale BSAs had great agreement with the reference simulated coarse scale albedo with an RMSE of 0.0001.When the slope is varied from greater than 10 • and less than 20 • , where the terrain can be considered as a the gentle slope, the RMSE was increased to 0.0019 (Figure 7B).When the mean slope was greater than 20 • and less than 30 • (where the land surface can be considered some rugged), the RMSE increased to 0.0034 (Figure 7C).However, when the slope was larger than 30 • , the RMSE still reminded as 0.0041 (Figure 7D).Though, the regression coefficient and R 2 were much closer to 1 at these four situations, the MAPE had an obvious increase from 0.05% to 2.04%.From the histograms of the mean slopes less than 30 • (Figure 7A-C), we could see that more than 70% of the biases were distributed at the intervals less than 0.01.The maximum bias was smaller than 0.02 at the four mean slopes of the DEM.These four figures showed that the aggregated coarse scale BSAs had the great agreement with the simulated coarse scale BSAs at the four selected mean slopes.Figure 7 shows the comparison between the simulated coarse scale BSAs and reference coarse scale BSAs over four selected DEM types with different mean slopes.The mean slope of the DEM in Figure 7A was lower than 10°, which was considered as the flat terrain.The results showed that the aggregated coarse scale BSAs had great agreement with the reference simulated coarse scale albedo with an RMSE of 0.0001.When the slope is varied from greater than 10° and less than 20°, where the terrain can be considered as a the gentle slope, the RMSE was increased to 0.0019 (Figure 7B).When the mean slope was greater than 20° and less than 30° (where the land surface can be considered some rugged), the RMSE increased to 0.0034 (Figure 7C).However, when the slope was larger than 30°, the RMSE still reminded as 0.0041 (Figure 7D).Though, the regression coefficient and R 2 were much closer to 1 at these four situations, the MAPE had an obvious increase from 0.05% to 2.04%.From the histograms of the mean slopes less than 30° (Figure 7A-C), we could see that more than 70% of the biases were distributed at the intervals less than 0.01.The maximum bias was smaller than 0.02 at the four mean slopes of the DEM.These four figures showed that the aggregated coarse scale BSAs had the great agreement with the simulated coarse scale BSAs at the four selected mean slopes.Figure 7 shows the comparison between the simulated coarse scale BSAs and reference coarse scale BSAs over four selected DEM types with different mean slopes.The mean slope of the DEM in Figure 7A was lower than 10°, which was considered as the flat terrain.The results showed that the aggregated coarse scale BSAs had great agreement with the reference simulated coarse scale albedo with an RMSE of 0.0001.When the slope is varied from greater than 10° and less than 20°, where the terrain can be considered as a the gentle slope, the RMSE was increased to 0.0019 (Figure 7B).When the mean slope was greater than 20° and less than 30° (where the land surface can be considered some rugged), the RMSE increased to 0.0034 (Figure 7C).However, when the slope was larger than 30°, the RMSE still reminded as 0.0041 (Figure 7D).Though, the regression coefficient and R 2 were much closer to 1 at these four situations, the MAPE had an obvious increase from 0.05% to 2.04%.From the histograms of the mean slopes less than 30° (Figure 7A-C), we could see that more than 70% of the biases were distributed at the intervals less than 0.01.The maximum bias was smaller than 0.02 at the four mean slopes of the DEM.These four figures showed that the aggregated coarse scale BSAs had the great agreement with the simulated coarse scale BSAs at the four selected mean slopes.

Accuracy Analysis by using the Linear Weighted Average Upscaling Model
Figure 8 shows the comparison between the reference coarse scale albedos and the aggregated coarse scale albedos, which averaged the fine scale albedo within a coarse scale pixel.The linear averaged albedos showed a big discrepancy with the simulated reference coarse scale albedo with an RMSE of 0.0744, and a bias of −0.0481 overall, for the different slopes and incident solar geometries.The R 2 was reduced as low as 0.0223.Though 80% of the bias was distributed symmetrically in the interval of −0.10 to 0.15, the maximum bias increased to nearly 0.25.Additionally, the comparison with WSA showed similar results as the BSA with a bias of −0.0383, and RMSE of 0.0532, respectively.Figure 8D shows the changes of biases between the aggregated coarse scale WSAs and the reference simulated coarse scale WSAs following the increase of mean slope.Furthermore, it shows the obviously increase of the bias following the change of the mean slope with the maximum absolute bias greater than 0.12 at the place that the mean slope is increased to 47.18.The results showed that the linear weighted average model was not very suitable to deal with the scaling mismatching for coarse scale albedo validation over rugged terrain.
The statistics metrics, which are summarized in Table 4, record the comparison results between the reference coarse scale albedo and the aggregated coarse scale albedo which was aggregated by the MRT-based upscaling model and linear weighted average model, respectively.Overall, the MRTbased upscaled coarse scale BSAs had a great agreement with the simulated reference coarse scale BSAs.Following the increase of the solar zenith angle, the uncertainty of the MRT-based upscaled coarse scale BSAs had a slight increase as the RMSE increased from 0.0022 to 0.0049, and the MAPE increased from 0.71% to 1.50%.Meanwhile, the RMSE and MAPE also showed an obviously increase trend following the increase of mean slope of the DEM where the RMSE varied from 0.0001 to 0.0041, and the MAPE varied from 0.05% to 2.04%.However, the aggregated coarse scale BSAs upscaled by the linear average model had large uncertainty when compared with the reference simulated coarse scale BSAs.Following the increasing of SZA and mean slope of DEM, the bias, MAPE, and RMSE had a similar increasing trend as that of the aggregated coarse scale BSA upscaled by the MRT-based model.When the mean slope was smaller than 10°, the coarse scale BSAs aggregated by the linear average model had great consistency with the reference simulated coarse scale BSAs with an RMSE of 0.0043, bias of 0.0002, MAPE of 1.26%, and R 2 of 0.9028.When the slope is increased to larger than 30°, the absolute of the bias, RMSE, MAPE, and R 2 also increased as larger as 0.0860, 0.0809, 37.37%, and 0.1288, respectively.

Accuracy Analysis by using the Linear Weighted Average Upscaling Model
Figure 8 shows the comparison between the reference coarse scale albedos and the aggregated coarse scale albedos, which averaged the fine scale albedo within a coarse scale pixel.The linear averaged albedos showed a big discrepancy with the simulated reference coarse scale albedo with an RMSE of 0.0744, and a bias of −0.0481 overall, for the different slopes and incident solar geometries.The R 2 was reduced as low as 0.0223.Though 80% of the bias was distributed symmetrically in the interval of −0.10 to 0.15, the maximum bias increased to nearly 0.25.Additionally, the comparison with WSA showed similar results as the BSA with a bias of −0.0383, and RMSE of 0.0532, respectively.Figure 8D shows the changes of biases between the aggregated coarse scale WSAs and the reference simulated coarse scale WSAs following the increase of mean slope.Furthermore, it shows the obviously increase of the bias following the change of the mean slope with the maximum absolute bias greater than 0.12 at the place that the mean slope is increased to 47.18.The results showed that the linear weighted average model was not very suitable to deal with the scaling mismatching for coarse scale albedo validation over rugged terrain.
The statistics metrics, which are summarized in Table 4, record the comparison results between the reference coarse scale albedo and the aggregated coarse scale albedo which was aggregated by the MRT-based upscaling model and linear weighted average model, respectively.Overall, the MRT-based upscaled coarse scale BSAs had a great agreement with the simulated reference coarse scale BSAs.Following the increase of the solar zenith angle, the uncertainty of the MRT-based upscaled coarse scale BSAs had a slight increase as the RMSE increased from 0.0022 to 0.0049, and the MAPE increased from 0.71% to 1.50%.Meanwhile, the RMSE and MAPE also showed an obviously increase trend following the increase of mean slope of the DEM where the RMSE varied from 0.0001 to 0.0041, and the MAPE varied from 0.05% to 2.04%.However, the aggregated coarse scale BSAs upscaled by the linear average model had large uncertainty when compared with the reference simulated coarse scale BSAs.Following the increasing of SZA and mean slope of DEM, the bias, MAPE, and RMSE had a similar increasing trend as that of the aggregated coarse scale BSA upscaled by the MRT-based model.When the mean slope was smaller than 10 • , the coarse scale BSAs aggregated by the linear average model had great consistency with the reference simulated coarse scale BSAs with an RMSE of 0.0043, bias of 0.0002, MAPE of 1.26%, and R 2 of 0.9028.When the slope is increased to larger than 30 • , the absolute of the bias, RMSE, MAPE, and R 2 also increased as larger as 0.0860, 0.0809, 37.37%, and 0.1288, respectively.The discrepancies between the MRT-based upscaled albedo and the reference albedo raised slightly following the increasing of the incident solar zenith angle and the mean slope, especially when the SZA was larger to 60° and the mean slope was larger than 30°.However, even though the  The discrepancies between the MRT-based upscaled albedo and the reference albedo raised slightly following the increasing of the incident solar zenith angle and the mean slope, especially when the SZA was larger to 60 • and the mean slope was larger than 30 • .However, even though the SZA is increased to 60 • and the mean slope increased to larger than 30 • , the MRT-based albedo upscaling model showed the great performance over rugged terrain with a maximum bias and RMSE smaller than 0.0032 and 0.0049, respectively.

Accuracy Assessment of HJ Albedo
The accuracy and uncertainty of fine scale HJ blue-sky albedos have an essential influence on coarse scale albedo validation in the multi-scale validation strategy [24].Therefore, the HJ blue-sky albedos were assessed firstly with in situ albedos.The HJ blue sky albedo can be approximated through a linear combination of BSA and WSA, weighted by the fraction of actual direct to diffuse skylight.Figure 9A shows the comparison between the in situ albedo and HJ blue-sky albedo.It was obvious that the two albedos had a less discrepancy with a bias of 0.0028, RMSE of 0.0272, MAPE of 10.0%, and R 2 of 0.7470.Figure 9B shows the error histogram of bias.It was easily found that the biases were normally distributed significantly.Though 80% of the bias was distributed within the interval of −0.05 to 0.05 and the maximum bias was less than 0.07.The validation results indicated that within the 95% confidence interval (shown by the p-value, which is less than 0.005), the regression model could be used to describe the relationship between the fine scale HJ albedos and in situ albedos.Additionally, it can be used as the calibration model to calibrate the HJ albedos.
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 23 SZA is increased to 60° and the mean slope increased to larger than 30°, the MRT-based albedo upscaling model showed the great performance over rugged terrain with a maximum bias and RMSE smaller than 0.0032 and 0.0049, respectively.

Accuracy Assessment of HJ Albedo
The accuracy and uncertainty of fine scale HJ blue-sky albedos have an essential influence on coarse scale albedo validation in the multi-scale validation strategy [24].Therefore, the HJ blue-sky albedos were assessed firstly with in situ albedos.The HJ blue sky albedo can be approximated through a linear combination of BSA and WSA, weighted by the fraction of actual direct to diffuse skylight.Figure 9A shows the comparison between the in situ albedo and HJ blue-sky albedo.It was obvious that the two albedos had a less discrepancy with a bias of 0.0028, RMSE of 0.0272, MAPE of 10.0%, and R 2 of 0.7470.Figure 9B shows the error histogram of bias.It was easily found that the biases were normally distributed significantly.Though 80% of the bias was distributed within the interval of −0.05 to 0.05 and the maximum bias was less than 0.07.The validation results indicated that within the 95% confidence interval (shown by the p-value, which is less than 0.005), the regression model could be used to describe the relationship between the fine scale HJ albedos and in situ albedos.Additionally, it can be used as the calibration model to calibrate the HJ albedos.

Comparison of the Aggregated Coarse Scale Albedo with MODIS Albedo Products
The MRT-based albedo upscaling model was used for aggregating the calibrated HJ albedo to a coarse spatial scale.Figure 10 displays the validation results by comparing the coarse scale MCD43A3 C6 albedos with the coarse scale aggregated HJ albedos.The results showed that MCD43A3 C6 BSAs with full inversion had agreement with the aggregated HJ BSAs, with an RMSE of 0.0305, bias of 0.008, MAPE of 13.50%, and R 2 of 0.6517.While, the MCD43A3 C6 BSAs with magnitude inversion showed slight differences with the aggregated HJ BSAs with an RMSE of 0.0511 and bias of 0.0097, MAPE of 19.17%, and R 2 of 0.3817.The MCD43A3 WSAs showed the similar accuracy and uncertainty as the BSAs, with an RMSE of 0.0321 for full inversion and an RMSE of 0.0531 for magnitude inversion.The results showed that the uncertainty of MCD43A3 C6 over rugged terrain was slightly larger than 0.02, which has been shown in many studies [15,16,18,20,22,[59][60][61].The two bias histograms showed that the biases were distributed symmetrically at the interval zero and the maximum bias were smaller than 0.1.

Comparison of the Aggregated Coarse Scale Albedo with MODIS Albedo Products
The MRT-based albedo upscaling model was used for aggregating the calibrated HJ albedo to a coarse spatial scale.Figure 10 displays the validation results by comparing the coarse scale MCD43A3 C6 albedos with the coarse scale aggregated HJ albedos.The results showed that MCD43A3 C6 BSAs with full inversion had agreement with the aggregated HJ BSAs, with an RMSE of 0.0305, bias of 0.008, MAPE of 13.50%, and R 2 of 0.6517.While, the MCD43A3 C6 BSAs with magnitude inversion showed slight differences with the aggregated HJ BSAs with an RMSE of 0.0511 and bias of 0.0097, MAPE of 19.17%, and R 2 of 0.3817.The MCD43A3 WSAs showed the similar accuracy and uncertainty as the BSAs, with an RMSE of 0.0321 for full inversion and an RMSE of 0.0531 for magnitude inversion.The results showed that the uncertainty of MCD43A3 C6 over rugged terrain was slightly larger than 0.02, which has been shown in many studies [15,16,18,20,22,[59][60][61].The two bias histograms showed that the biases were distributed symmetrically at the interval zero and the maximum bias were smaller than 0.1.Figure 11 illustrates the MCD43A3 C6 products validation results at different mean slopes around the AWSs.This showed that the MCD43A3 C6 products had good agreements with the aggregated HJ albedos at the surface where the slope was lower than 5°.The bias, MAPE, RMSE and R 2 for the MCD43A3 full inversion BSAs were 0.0108, 10.88%, 0.0244, and 0.6136, respectively (Figure 11A).When the slope increased to greater than 5° and less than 10°, the RMSE increased to 0.0315 for BSA and 0.0308 for the WSAs of full inversion (Figure 11C,D).However, when the slope was larger than 10°, there were large discrepancies between the MCD43A3 C6 full inversion BSAs and the aggregated HJ BSAs with an RMSE of 0.0365 and MAPE of 19.71% (Figure 11E).The WSAs of full inversion showed a similar increase trend, with a bias of 0.0126, MAPE of 11.39% and RMSE of 0.0293 (Figure 11B), when the mean slope was lower than 5°.When the mean slope was raised to larger than 10°, the bias, MAPE, and RMSE increased to −0.0133, 15.45% and 0.0313, respectively (Figure 11F).
Generally, the MCD43A3 C6 albedos of magnitude inversion had lower quality than the albedos of full inversion.Following the increase of mean slope, the RMSE and MAPE displayed an outstanding rise for the two inversions MCD43A3 C6 BSAs, as well as the magnitude inversion WSAs, especially where the mean slope was greater than 10°.Over rugged terrain, the land surface albedos were seriously influenced by the complex interactions among the topography and the incident radiation from the sun location, the shadow, the atmosphere, and surface scattering irradiance from the adjacent terrain.The topographic effects led to more uncertainty in albedo remote sensing retrieval.Neglecting these effects, the albedo algorithm may lead to lower quality albedos.Figure 11 illustrates the MCD43A3 C6 products validation results at different mean slopes around the AWSs.This showed that the MCD43A3 C6 products had good agreements with the aggregated HJ albedos at the surface where the slope was lower than 5 • .The bias, MAPE, RMSE and R 2 for the MCD43A3 full inversion BSAs were 0.0108, 10.88%, 0.0244, and 0.6136, respectively (Figure 11A).When the slope increased to greater than 5 • and less than 10 • , the RMSE increased to 0.0315 for BSA and 0.0308 for the WSAs of full inversion (Figure 11C,D).However, when the slope was larger than 10 • , there were large discrepancies between the MCD43A3 C6 full inversion BSAs and the aggregated HJ BSAs with an RMSE of 0.0365 and MAPE of 19.71% (Figure 11E).The WSAs of full inversion showed a similar increase trend, with a bias of 0.0126, MAPE of 11.39% and RMSE of 0.0293 (Figure 11B), when the mean slope was lower than 5 • .When the mean slope was raised to larger than 10 • , the bias, MAPE, and RMSE increased to −0.0133, 15.45% and 0.0313, respectively (Figure 11F).
Generally, the MCD43A3 C6 albedos of magnitude inversion had lower quality than the albedos of full inversion.Following the increase of mean slope, the RMSE and MAPE displayed an outstanding rise for the two inversions MCD43A3 C6 BSAs, as well as the magnitude inversion WSAs, especially where the mean slope was greater than 10 • .Over rugged terrain, the land surface albedos were seriously influenced by the complex interactions among the topography and the incident radiation from the sun location, the shadow, the atmosphere, and surface scattering irradiance from the adjacent terrain.The topographic effects led to more uncertainty in albedo remote sensing retrieval.Neglecting these effects, the albedo algorithm may lead to lower quality albedos.Figure 11 illustrates the MCD43A3 C6 products validation results at different mean slopes around the AWSs.This showed that the MCD43A3 C6 products had good agreements with the aggregated HJ albedos at the surface where the slope was lower than 5°.The bias, MAPE, RMSE and R 2 for the MCD43A3 full inversion BSAs were 0.0108, 10.88%, 0.0244, and 0.6136, respectively (Figure 11A).When the slope increased to greater than 5° and less than 10°, the RMSE increased to 0.0315 for BSA and 0.0308 for the WSAs of full inversion (Figure 11C,D).However, when the slope was larger than 10°, there were large discrepancies between the MCD43A3 C6 full inversion BSAs and the aggregated HJ BSAs with an RMSE of 0.0365 and MAPE of 19.71% (Figure 11E).The WSAs of full inversion showed a similar increase trend, with a bias of 0.0126, MAPE of 11.39% and RMSE of 0.0293 (Figure 11B), when the mean slope was lower than 5°.When the mean slope was raised to larger than 10°, the bias, MAPE, and RMSE increased to −0.0133, 15.45% and 0.0313, respectively (Figure 11F).
Generally, the MCD43A3 C6 albedos of magnitude inversion had lower quality than the albedos of full inversion.Following the increase of mean slope, the RMSE and MAPE displayed an outstanding rise for the two inversions MCD43A3 C6 BSAs, as well as the magnitude inversion WSAs, especially where the mean slope was greater than 10°.Over rugged terrain, the land surface albedos were seriously influenced by the complex interactions among the topography and the incident radiation from the sun location, the shadow, the atmosphere, and surface scattering irradiance from the adjacent terrain.The topographic effects led to more uncertainty in albedo remote sensing retrieval.Neglecting these effects, the albedo algorithm may lead to lower quality albedos.

Discussion
The topography affects the land surface albedo.Compared with the validation over flat terrain, the albedo validation over rugged terrain should consider the scale effects.This paper proposed an MRT-based albedo upscaling model for dealing with the scale effects in the multi-scale validation strategy over rugged terrain.The simulated coarse scale albedos were used as the reference truth to assess the MRT-based albedo upscaling model.
The assessments showed that the coarse scale albedo upscaled by the MRT-based model had the great agreement with the reference coarse scale albedo in different solar zenith angle or terrain with different mean slopes (Table 4).The uncertainty of the aggregated coarse scale albedos had the slight variation following the change of SZA and mean slope.In particular, the coarse scale albedos aggregated by the MRT-based model still showed the great consistency with the reference coarse scale albedo when the mean slope was larger than 30° or the solar zenith angle was up to 60° (Figures 6 and 7).Adversely, the coarse scale albedo aggregated by the linear average model showed a larger discrepancy with the reference coarse scale albedo (Figure 8).Following the increase of solar zenith angle and the mean slope of the terrain, the discrepancy had an obviously increasing trend (Table 4).The MAPE and RMSE were raised to 37.37% and 0.0809, following the increase of the mean slope, respectively.Comparisons between the validation results of the two upscaled coarse scale albedos showed that the MAPEs increased obviously following the increase of mean slope (Table 4).When the mean slope was larger than 30°, the difference of MAPE increased larger than 35%.However, when the mean slope was lower than 10°, the MAPE did not show an obvious difference.Similarly, the difference between the two upscaled models also increased following the increase of the solar

Discussion
The topography affects the land surface albedo.Compared with the validation over flat terrain, the albedo validation over rugged terrain should consider the scale effects.This paper proposed an MRT-based albedo upscaling model for dealing with the scale effects in the multi-scale validation strategy over rugged terrain.The simulated coarse scale albedos were used as the reference truth to assess the MRT-based albedo upscaling model.
The assessments showed that the coarse scale albedo upscaled by the MRT-based model had the great agreement with the reference coarse scale albedo in different solar zenith angle or terrain with different mean slopes (Table 4).The uncertainty of the aggregated coarse scale albedos had the slight variation following the change of SZA and mean slope.In particular, the coarse scale albedos aggregated by the MRT-based model still showed the great consistency with the reference coarse scale albedo when the mean slope was larger than 30 • or the solar zenith angle was up to 60 • (Figures 6 and 7).Adversely, the coarse scale albedo aggregated by the linear average model showed a larger discrepancy with the reference coarse scale albedo (Figure 8).Following the increase of solar zenith angle and the mean slope of the terrain, the discrepancy had an obviously increasing trend (Table 4).The MAPE and RMSE were raised to 37.37% and 0.0809, following the increase of the mean slope, respectively.Comparisons between the validation results of the two upscaled coarse scale albedos showed that the MAPEs increased obviously following the increase of mean slope (Table 4).When the mean slope was larger than 30 • , the difference of MAPE increased larger than 35%.However, when the mean slope was lower than 10 • , the MAPE did not show an obvious difference.Similarly, the difference between the two upscaled models also increased following the increase of the solar zenith angle.When the solar zenith angle was 60 • , the difference between the two MAPEs was larger than 20%.Notably, the two MAPEs also showed a large difference when the solar zenith angle was 0 • .
The simulated results showed that the SZA and mean slope had a significant influence on the scale effects of albedo over rugged terrain.When the land surface was flat or the mean slope was lower than 10 • , the scale effects could be thought as very gentle, and either the MRT-based albedo upscaling model or the linear average upscaling model could be used to upscale fine scale albedo to a coarse scale.When the land surface was steep or the mean slope larger than 10 • , the difference caused by scale effects was larger than 15% of MAPE.The solar zenith angle coupled with complex terrain can cause the mutual shadowing and re-distribution of download solar radiation.These characteristics resulted in the land surface being more heterogeneous and larger albedo scale effects.
The MRT-based upscaling model was applied to MCD43A3 C6 albedo products validation in HRB.Fine scale albedos were retrieved from the HJCCD sensor by the improved AB algorithm on a sloping surface and used as a bridge to deal with the scale mismatching between the in situ albedos and the MCD43A3 C6 albedos.The HJ albedos were evaluated and calibrated by in situ albedos for reducing the systemic errors and the accidental errors.Furthermore, they were upscaled to the coarse scale and were functioned as the reference truth over albedo validation.The validation results showed the MCD43A3 C6 products of full inversion had high agreements at the sites where the mean slope was lower than 5 • .Following the increase of mean slope, the MCD43A3 C6 BSAs of full or magnitude inversion showed the outstanding increase of RMSE, with an RMSE varied from 0.0244 to 0.0365 for the full inversion BSA and 0.040 to 0.0601 for the magnitude inversion, respectively.The MCD43A3 of the full WSAs products also showed a similar result with the RMSE varying from 0.0293 to 0.0313 following the mean slope increasing from 0 • to greater than 10 • .
However, the multi-scale validation strategy used for albedo validation over rugged terrain also still faces challenges.The geometric misalignment between DEM and fine scale albedo had an intense influence on the application of the MRT-based albedo upscaling method.In this paper, the HJ reflectance data were geometrically rectified by the geo-referenced Landsat images and rechecked manually.Second, the MRT-based albedo upscaling model was built under the assumption that the coarse scale pixel was overall horizontal.To use the MRT-based method where the coarse scale pixel was not overall horizontal, the simple way was to change the angular effects in model development.Therefore, the MRT-based upscaling model can also be suitable for using in coarse scale albedo validation.The quality of fine scale albedo and the multiple land covers were also affected the validation results.Thus, in the future, we plan to improve the method's applicability and assess the accuracy of MCD43A3 C6 products comprehensively under multiple land covers and long temporal periods over rugged terrain at the global scale.

Conclusions
This paper proposed a multi-scale validation strategy for coarse scale albedo products validation over rugged terrain, where the fine scale albedos on sloping land surfaces were used as a bridge to upscale in situ albedo to coarse scale.In this paper, the MRT-based albedo upscaling method was developed and preliminarily applied for aggregating the fine scale albedos to coarse scale.The applicability and performance of the MRT-based upscaling model were assessed by comparing the upscaled coarse scale albedo with the simulated reference coarse scale albedo.The multi-scale validation strategy with an MRT-based albedo upscaling model was applied to assess the accuracy of MCD43A3 products over rugged terrain in the Heihe River Basin, China.The HJ satellite albedo was selected as the fine scale albedo and was upscaled to the coarse scale by the MRT-based albedo upscaling model.
The validation results showed that the coarse scale albedo validation are suffered the scale effects over rugged terrain.The scale effects increased following the increase of mean slope and the solar zenith angle over rugged terrain.Furthermore, it cannot be overlooked when the mean slope was larger than 10 • .The MRT-based albedo upscaling model had the great performance in dealing with the issue of scale mismatching for albedo validation over rugged terrain.
The linear average model showed a similar performance with the albedo validation over flat terrain.However, over rugged terrain, the coarse scale albedos upscaled by the linear average model showed a bigger discrepancy than the simulated coarse scale albedos, which indicated that the linear average model had worse applicability for albedo validation over rugged terrain, especially, when the slope was larger than 10 • .
The MRT-based albedo upscaling model was used to assess the MCD43A3 products.It showed great performance in aggregating the fine scale HJ albedo to coarse scale.The validation results showed that in the Heihe River Basin, the MCD43A3 BSA products with full inversion algorithm had an increased uncertainty following the increase of the mean slope, with an RMSE varying from 0.0244 to 0.0365, and a MAPE that varied from 10.88% to 19.71.Meanwhile, the MCD43A3 WSA products with the full inversion algorithm had an RMSE varying from 0.0293 to 0.0313 and MAPE from 11.39% to 15.49%.

Figure 1 .
Figure 1.The coarse albedo validation procedure over rugged terrain.

)
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 23 where c A is the coarse scale albedo; c φ ↑ , c φ ↓ are the coarse scale pixel's reflected and incident solar radiant flux, respective; N is the amount of micro-slope within the corresponding area of coarse scale pixel; k ∂ and T k d φ ↓ are the albedo and incident radiant flux of the th k micro-slope, respectively.

Figure 2 .
Figure 2. The contributions of incident and outgoing radiant flux over a coarse scale pixel; direct, diffuse and terrain irradiances.

Figure 2 .
Figure 2. The contributions of incident and outgoing radiant flux over a coarse scale pixel; direct, diffuse and terrain irradiances.

Figure 4 .
Figure 4. Study area: (A) the overview of the study area; (B) the location of the Automatic Weather Stations (AWSs) in the DEM imageries; and (C) one example of the AWSs.

Figure 4 .
Figure 4. Study area: (A) the overview of the study area; (B) the location of the Automatic Weather Stations (AWSs) in the DEM imageries; and (C) one example of the AWSs.

Figure 5 .
Figure 5. Scatter plots and the histogram evaluate the aggregated coarse scale albedo (using a Mountain-Radiation-Transfer (MRT)-based model) against the reference simulated coarse scale albedo for (A) the aggregated coarse scale BSA vs. the simulated coarse scale BSA; (B) the bias histogram of the aggregated coarse scale BSA minus the reference coarse scale BSA; (C) the aggregated coarse scale WSA vs. the simulated coarse scale WSA; and (D) the bias between the nine aggregated coarse scale WSAs and the reference simulated WSAs following the increase of the mean slope.

Figure 7 .
Figure 7. Validation of the MRT-based upscaling model at different slopes (A) the slope is smaller than 10 degree; (B) the slope is greater than 10 degree and less than 20 degree; (C) the slope is between 20 and 30 degrees; and (D) the slope is larger than 30 degrees.

Figure 7 .
Figure 7. Validation of the MRT-based upscaling model at different slopes (A) the slope is smaller than 10 degree; (B) the slope is greater than 10 degree and less than 20 degree; (C) the slope is between 20 and 30 degrees; and (D) the slope is larger than 30 degrees.

Figure 8 .
Figure 8. Scatter plots and the histogram evaluate the fine scale albedo against simulated coarse scale albedo for: (A) the aggregated fine scale BSA vs. the simulated coarse scale BSA; (B) the histogram of the aggregated fine scale BSA minus the coarse scale BSA; (C) the aggregated fine scale WSA vs. the simulated coarse scale WSA; and (D) the bias distribution of the aggregated fine scale WSA and the coarse scale WSA following the increase of the mean slopes.

Figure 8 .
Figure 8. Scatter plots and the histogram evaluate the fine scale albedo against simulated coarse scale albedo for: (A) the aggregated fine scale BSA vs. the simulated coarse scale BSA; (B) the histogram of the aggregated fine scale BSA minus the coarse scale BSA; (C) the aggregated fine scale WSA vs. the simulated coarse scale WSA; and (D) the bias distribution of the aggregated fine scale WSA and the coarse scale WSA following the increase of the mean slopes.

Figure 9 .
Figure 9. Scatter plots and bias histogram between fine scale albedos and the in situ albedos: (A) in situ albedo vs. the fine scale albedo; (B) Bias histogram of fine scale blue-sky albedo minus the in situ albedo The colors refer to the density of points (from highest (red) to lowest (blue).

Figure 9 .
Figure 9. Scatter plots and bias histogram between fine scale albedos and the in situ albedos: (A) in situ albedo vs. the fine scale albedo; (B) Bias histogram of fine scale blue-sky albedo minus the in situ albedo The colors refer to the density of points (from highest (red) to lowest (blue).

Figure 10 .
Figure 10.Coarse scale MCD43A3 C6 albedos validation by comparison with the aggregated coarse scale albedos: (A) the aggregated coarse scale BSAs vs. the MCD43A3 C6 BSAs; (B) the aggregated coarse scale WSAs vs. the MCD43A3 C6 WSAs.The colors refer to the density of points (from highest (red) to lowest (blue).

Figure 10 .
Figure 10.Coarse scale MCD43A3 C6 albedos validation by comparison with the aggregated coarse scale albedos: (A) the aggregated coarse scale BSAs vs. the MCD43A3 C6 BSAs; (B) the aggregated coarse scale WSAs vs. the MCD43A3 C6 WSAs.The colors refer to the density of points (from highest (red) to lowest (blue).

23 Figure 10 .
Figure 10.Coarse scale MCD43A3 C6 albedos validation by comparison with the aggregated coarse scale albedos: (A) the aggregated coarse scale BSAs vs. the MCD43A3 C6 BSAs; (B) the aggregated coarse scale WSAs vs. the MCD43A3 C6 WSAs.The colors refer to the density of points (from highest (red) to lowest (blue).

Figure 11 .
Figure 11.Coarse scale MCD43A3 C6 products validation by comparing with the aggregated fine scale albedo at different mean slopes.The colors refer to the density of points (from highest (red) to lowest (blue).

Figure 11 .
Figure 11.Coarse scale MCD43A3 C6 products validation by comparing with the aggregated fine scale albedo at different mean slopes.The colors refer to the density of points (from highest (red) to lowest (blue).

Table 1 .
The parameters and the values of the PROSAIL model.

Table 2 .
Mean slopes of the DEMs.

Table 3 .
The description of the selected AWSs.

Table 3 .
The description of the selected AWSs.

Table 4 .
Statistical metrics from the validation between coarse scale BSA and aggregated fine scale BSA aggregated by different upscaling methods (α* means the slope of the DEM).

Table 4 .
Statistical metrics from the validation between coarse scale BSA and aggregated fine scale BSA aggregated by different upscaling methods (α* means the slope of the DEM).