Theoretical Evaluation of Anisotropic Reﬂectance Correction Approaches for Addressing Multi-Scale Topographic Effects on the Radiation-Transfer Cascade in Mountain Environments

: Research involving anisotropic-reﬂectance correction (ARC) of multispectral imagery to account for topographic effects has been ongoing for approximately 40 years. A large body of research has focused on evaluating empirical ARC methods, resulting in inconsistent results. Consequently, our research objective was to evaluate commonly used ARC methods using ﬁrst-order radiation-transfer modeling to simulate ASTER multispectral imagery over Nanga Parbat, Himalaya. Speciﬁcally, we accounted for orbital dynamics, atmospheric absorption and scattering, direct- and diffuse-skylight irradiance, land cover structure, and surface biophysical variations to evaluate their effectiveness in reducing multi-scale topographic effects. Our results clearly reveal that the empirical methods we evaluated could not reasonably account for multi-scale topographic effects at Nanga Parbat. The magnitude of reflectance and the correlation structure of biophysical properties were not preserved in the topographically-corrected multispectral imagery. The CCOR and SCS+C methods were able to remove topographic effects, given the Lambertian assumption, although atmospheric correction was required, and we did not account for other primary and secondary topographic effects that are thought to significantly influence spectral variation in imagery acquired over mountains. Evaluation of structural-similarity index images revealed spatially variable results that are wavelength dependent. Collectively, our simulation and evaluation procedures strongly suggest that empirical ARC methods have significant limitations for addressing anisotropic reflectance caused by multi-scale topographic effects. Results indicate that atmospheric correction is essential, and most methods failed to adequately produce the appropriate magnitude and spatial variation of surface reflectance in corrected imagery. Results were also wavelength dependent, as topographic effects influence radiation-transfer components differently in different regions of the electromagnetic spectrum. Our results explain inconsistencies described in the literature, and indicate that numerical modeling efforts are required to better account for multi-scale topographic effects in various radiation-transfer components. Our results reveal that the we multi-scale effects, countering results The magnitude of reﬂectance and the correlation structure of properties were not in the Our best results are obtained by using the CCOR and SCS+C methods, although comprehensive atmospheric correction required.


Introduction
Remote-sensing science and technology is essential for studying and understanding the complexities associated with bio-physical/geophysical parameters and landscape dynamics [1,2]. This is especially the case in mountain environments where climate and tectonic forcings govern topography, high-magnitude surface processes, rapid environmental change and geohazards [3][4][5][6]. The inherent coupling of atmospheric, surface and tectonic processes highlights the complexity associated with mountain geodynamics, and the need for quantitative information about surface parameters, process rates and resource availability. Due to logistic, geopolitical, geohazard and physical limitations associated with fieldwork in many mountains, remote sensing represents the most effective approach for producing spatiotemporal quantitative and thematic information [7].
Remote sensing of mountain environments, however, is notoriously difficult, as multi-scale topographic effects govern the anisotropic nature of the radiation-transfer cascade (RTC) including atmospheric scattering, irradiance components, and surface reflectance [1,8]. Research indicates that it is essential to address the anisotropic nature of spectral variation caused by various radiation-transfer parameters that are strongly controlled by topographic effects [1].
Anisotropic-reflectance correction (ARC) research has been ongoing for about 40 years, e.g., [9][10][11][12][13]. Numerous approaches to ARC have been characterized as: (1) spectral-feature extraction methods; (2) statistical-empirical methods; (3) semi-empirical methods; and (4) radiation-transfer modeling [14,15]. The most frequently utilized methods from the first three categories can all be considered empirical in nature, as proxy parameters are used as surrogates for representing radiation-transfer parameters [1]. Furthermore, they are popular and have been compared to each other in various environments, because researchers and practitioners do not usually have access to spatiotemporal data, information and parameters that are required by more rigorous radiation-transfer models.
Commonly used empirical ARC methods used in comparative studies include the Statistical-Empirical Correction, B-Correction, Variable Empirical Coefficient Algorithm, Cosine Correction, C-Correction, Soil-Canopy Correction, SCS+C Correction, and Minnaert Correction. These and other approaches can be utilized based on variations in computing parameters using stratification of topography and land cover conditions in an attempt to optimize results. Given that these parameterization schemes are simplistic and cannot account for topographic anisotropy and the anisotropic nature of irradiance and reflectance, evaluation of these methods is subject to substantial uncertainties. An array of assumptions, exclusion of significant RTC parameters, and inconsistent results regarding best methods and parameterization schemes [16][17][18] strongly suggest that these ARC methods are not being adequately evaluated, given that the nature-of-the-problem is fundamentally a radiometric calibration issue.
Fundamental issues related to the complexity of the problem include the following: (1) Can ARC methods be used to produce reliable quantitative surface parameters? This will require the magnitude of surface reflectance to be valid given the viewing geometry. (2) Can ARC methods be used to process entire scenes of imagery? This will require that methods can address spatial variations in solar geometry, topographic geometry and land cover conditions. (3) Can various methods account for wavelength dependence of topographic effects? (4) Can they account for multi-scale topographic effects? (5) To what degree are methods or approaches dependent upon sensor characteristics? (6) Can they account for landscape complexity, or only work for homogeneous land cover conditions? (7) To what degree are various methods influenced by the Lambertian assumption or by the bidirectional reflectance distribution function (BRDF)? (8) Do these methods result in information compression which may be interpreted as a reduction in topographic effects? These issues have not been adequately addressed by comparative studies, although it is assumed that different methods work better in different mountain environments [16,18].
Consequently, our overall research objective was to utilize radiation-transfer modeling to improve our understanding of ARC issues by evaluating commonly used empirical methods. We generated a first-order approximation of the RTC over the Nanga Parbat Himalaya, Pakistan, and produced spectrally simulated ASTER satellite imagery. The Nanga Parbat massif is ideally suited to investigate issues associated with topographic effects and reflectance variations, as it exhibits extreme relief, the full range of slope and slope-azimuth angles, extreme spatial variability in morphometric properties including topographic shielding, and a wide variety of surface composition and land cover classes [19,20]. Specifically, we accounted for Earth-Sun orbit dynamics, atmospheric absorption and scattering, direct-and diffuse-skylight irradiance, land cover spatial structure, and surface biophysical variations in order to evaluate widely utilized empirical ARC methods for their ability to accurately account for multi-scale topographic effects, and adequately predict the magnitude and spatial variation of surface reflectance. In our research, multi-scale topographic effects were defined to be morphometric (morphology of the topography) constraints and parameters that partially govern radiation-transfer parameters. The analysis of the topography and the generation of topographic parameters and RT parameters were scale dependent with respect to direction and distance (i.e., multi-scale). It was not the purpose of this research to evaluate the influence of sensor-system characteristics on topographic-correction methods, as spatial, spectral and radiometric sensitivity also governs the variance structure in imagery. Limitations of this study include utilizing the Lambertian assumption for surface reflectance, and not accounting for the surface BRDF and the adjacent-terrain irradiance.

Radiation-Transfer Cascade
The ARC problem requires addressing the RT components of the RTC. Specifically, the spatial, temporal, and spectral dimensions of numerous RTC parameters must be accounted for. This includes orbital, solar, atmospheric, topographic, surface composition and sensor-system characteristics.
Solar geometry variations are strongly controlled by orbital dynamics, as orbital parameters such as eccentricity, obliquity, and longitude-of-perihelion govern the Earth-Sun distance, solar declination, and solar-zenith and -azimuth angles [21,22]. The magnitude of the exoatmospheric irradiance at the top of the Earth's atmosphere (E 0 ,[W m −2 µm −1 ]) is a function of the spectral exitance of the sun and the Earth-Sun distance. The atmosphere attenuates the solar beam irradiance via gaseous absorption and molecular and aerosol scattering [22,23]. Atmospheric processes are highly wavelength dependent and controlled by spatiotemporal variations in various atmospheric constituents, e.g., aerosols, water vapor; [22,24].
Atmospheric transmittance accounts for the primary constituents, such that the total downward and upward transmittance (T ↓↑ ) is represented as [25]: where T r represents transmittance due to Rayleigh scatter, T a is transmittance due to aerosol scattering, T g is transmittance due to primary gas absorption, T o is transmittance due to ozone absorption, T w is transmittance due to water vapor absorption, θ s is the solar-zenith angle (θ v should be substituted for θ s for upward transmittance; θ v is the sensor zenith-viewing angle), and λ is the wavelength of light. The total atmospheric transmittance is a function of the total optical depth of the atmosphere that is related to the hypsometry of the topography, such that variations in topographic relief can cause mesoscale variations in the atmospheric optical depth and volumetric scattering of radiation. Many ARC investigations do not account for this topographic effect by applying adequate atmospheric correction procedures, although atmospheric-topographic coupling must be accounted for [1]. Various other topographic effects govern the surface irradiance (E) such that [15]: where E b represents the direct-beam irradiance from the Sun, E d represents the diffuse-skylight irradiance given atmospheric scattering, and E t is the adjacent-terrain irradiance given the irradiance scattered from the surrounding terrain. Each irradiance component is partially regulated by the atmosphere and by multi-scale topographic effects, such that it is very difficult to address topographic effects without accounting for the coupled influence of the atmosphere [26]. These components should be accounted for in mountain environments, as they can significantly influence spectral variability in satellite imagery [8].

Direct Irradiance
The direct-beam irradiance is the dominant irradiance component [8]. It is strongly controlled by local and meso-scale topographic effects. The parameterization scheme for an exact computation of this RTC parameter is: where E 0 is the exoatmospheric irradiance, which is modified by the Earth-Sun distance correction factor, T ↓ is the downward total transmittance, cos i is the cosine of the incidence angle (i), and S is a binary coefficient that accounts for cast shadows being present (S = 0.0) or absent (S = 1.0) on the landscape. The total transmittance can be computed from Equation (1). Numerous other atmospheric constituents can be accounted for using various atmospheric-correction models (e.g., 6S and MOTRAN). Depending upon the climatic and tectonic setting that governs erosion, uplift and relief production, atmospheric transmittance may or may not be a significant RTC parameter for ARC. It is not usually accounted for in representing E b , because cos i is used as a proxy in most empirical ARC methods. It is necessary to address local topographic effects that account for the relationship between the solar and terrain geometry. Specifically, the incidence angle of illumination (i) between the sun and normal to the ground surface is defined as: where θ s is the apparent solar-zenith angle that accounts for altitude variations and atmospheric refraction given the atmospheric temperature and pressure profiles, θ t represents the terrain slope angle, φ s is the solar-azimuth angle, and φ t is the terrain slope-azimuth angle. Calculation of cos i is possible with the use of a digital elevation model (DEM), and negative values must be corrected to 0.0. The solar geometry is governed by orbital parameters that vary over time. The magnitude of variation increases over larger time periods. Most researchers do not compute the solar geometry for every pixel, thereby assuming the small-angle approximation, using one estimate for these parameters obtained from imagery for the entire scene or subscene. Nevertheless, as solar geometry is a function of location and time, the solar geometry varies for each pixel in a scene. It is unclear to what degree the small-angle approximation governs the accuracy of ARC efforts, as researchers usually do not evaluate full scenes of satellite imagery. Finally, the meso-scale relief structure of the topography, coupled with temporally-dependent solar geometry, governs the presence of cast shadows on the landscape. Ray-tracing algorithms can be used to account for this influence on E b [27][28][29][30]. It should be noted that this parameter is independent of the computation of cos i, and many studies do not account for it. It may or may not be an issue in ARC investigations depending upon the time of image acquisition and the magnitude of relief in the direction of the solar azimuth.

Diffuse-Skylight Irradiance
Atmospheric scattering will produce a hemispherical source of irradiance that contains isotropic and anisotropic components [8,31]. The anisotropic nature of E d has been the largest source of error associated with estimating this component, as it is governed by circumsolar brightening due to forward scattering of aerosols, and horizontal brightening due to multiple Rayleigh scattering [31]. An approximation model assumes that: where E r is the Rayleigh-scattering component, E a is the aerosol-scattering component, and E g is the ground/sky-backscattering component caused by multiple interactions between the ground surface and the atmosphere [1,8,24]. One can account for single-scattering or a multiple-scattering Rayleigh atmosphere. The single-scattering albedo is required to account for aerosol scattering which is a function of both wavelength and humidity. Aerosol type and models such as rural, urban, maritime or tropospheric can also be accounted for [24]. Finally, the ground-backscattering component can be modeled by accounting for the zonal-ground reflectance from the direct irradiance, the reflectance from the diffuse irradiance, and the overall sky reflectance [24].
To account for the diffuse irradiance for every pixel, numerous topographic effects must be accounted for and computed. Proy et al. [8] provided the computation solution such that: where L ↓ is the downward radiance from incident directions θ i and φ i , which are the zenith and azimuth angles of the incident energy from the hemisphere, respectively; α s is the solar elevation angle (α s = π/2 − θ s ); and I is the incidence angle of the direction defined by hemisphere and terrain geometry, similar to Equation (4) using incident geometry. Such an exact computation is computationally intensive given that E d is wavelength dependent, and hemispherical-topographic shielding must be accounted for. Multiple approaches to account for topographic shielding are based on the computation of a maximum relief angle θ max for φ i , e.g., [27,32]. The so-called skyview-factor coefficient is usually multiplied by estimates of diffuse irradiance for a horizontal surface to compute E d . Many researchers do not utilize this parameter due to computational intensity, but estimate it based on the slope angle, e.g., [33,34]. It is rare for investigators to account for both cos I and θ max (φ i ), thereby accounting for local and meso-scale topographic effects.

Adjacent-Terrain Irradiance
The aforementioned irradiance components and surface reflectance interact with the surrounding terrain geometry to produce the E t irradiance component. Numerous researchers have indicated that this component must be accounted for in mountains due to snow, vegetation and steep slopes [8,33,35]. It is an extremely complicated parameter to compute, as it is governed by numerous multi-scale topographic effects and surface anisotropic-reflectance conditions. A first-order approximation to E t was formulated by Proy et al. [8], assuming Lambertian reflectance. It is rarely accounted for in empirical ARC studies.
Its exact computation is: where L is the surface-reflected radiance coming from the effective incident direction, ρ brd f is the surface BRDF, θ i is the incident vertical-hemispherical zenith angle to account for terrain radiance above and below a pixel location, T ↓↑ t is the atmospheric transmittance given the optical depth of the atmosphere due to relief and propagation zenith angle through the atmosphere (θ v ) between two locations, I t represents the terrain-incidence angle given the influence of the local terrain geometry in relation to the incident directional geometry, and S t represents terrain blocking of the surface radiance between any two points of the surrounding terrain. The superscript e represents the effective zenith and azimuth angles that account for the influence of the terrain slope and slope-azimuth angles on the incident (θ e i , φ e i ) and viewing geometry (θ e v , φ e v ) of the surface BRDF compared to a horizontal surface. The computation of E t should account for the terrain conditions extending out to approximately 5 km from each pixel location [8]. Another complicating factor is that each pixel exhibits a unique BRDF, as topography and varying land cover structure and biophysical properties govern the BRDF [36][37][38].
Most empirical ARC investigations do not adequately account for E t , as they assume Lambertian reflectance, insignificant atmospheric attenuation, and do not iteratively compute estimates of E t over time (E t is partially governed by E t ). Such assumptions are not valid in many mountain environments, as steep slopes and highly variable biophysical properties caused by high-magnitude surface processes, create anisotropic reflectance conditions and extreme relief [11,35,[39][40][41][42][43][44][45]. Furthermore, topographic effects on the surface BRDF are not fully understood, although we know that the viewing geometry changes for every pixel, given variations in terrain slope and slope-azimuth angles [1].

Bi-Directional Reflectance Distribution Function
The characteristics of surface-reflectance properties can range between diffuse and specular [1,46]. A diffuse or Lambertian surface reflects radiation in all directions equally, and is considered an isotropic reflector. We know, however, that this assumption is not valid in mountains, and the surface BRDF should be accounted for to enable effective ARC [1]. The BRDF is used to describe the anisotropic nature of surface reflectance characteristics. It is a scattering function and can theoretically describe anisotropic reflectance given all input-output angles. It is characterized as: The effective illumination (θ e i , φ e i ) and viewing directions (θ e v , φ e v ) are a function of θ s , φ s , θ t , φ t , θ v , and φ v . Consequently, the topography has a significant effect on the BRDF [1]. The BRDF, however, is difficult to measure accurately, and, in practice, the bidirectional reflectance factor (BRF) is commonly used [1,47]. The application of different BRDF models in ARC studies has been relatively limited, and most studies utilize semi-empirical models because they can be fit to measured data [1]. These models are very flexible and account for a geometric component due to land cover structure and a diffuse volume component. Empirical coefficients can be used to determine the relative strength of these two scattering characteristics. The importance of topographic effects on the BRDF in ARC research cannot be overstated. They partially govern E d , E t , L, and L p (the ground component of additive-path radiance). Clearly, empirical ARC methods are too simplistic to account for the spectral variation of numerous radiation-transfer parameters that are part of the RTC.

ARC
The strong influence that topography has on spectral response as measured by satellite-scanner systems in rugged terrain has long been recognized, e.g., [48]. Various approaches have been introduced in an attempt to reduce or remove these effects. Nevertheless, a topographic-correction method has not emerged that is consistently superior in all cases [16,18]. In addition, commonly applied topographic-correction methods do not provide a fully comprehensive ARC solution. Topographic-correction approaches can be organized into three categories: spectral-feature extraction methods, empirical modeling, and radiation-transfer modeling [49][50][51]. Empirical methods can be divided into those based on a Lambertian assumption of surface reflectance and those based on a non-Lambertian assumption. The following background section discusses empirical topographic-correction approaches, as these represent the most commonly used methods and were evaluated in this study.

Empirical Corrections
In contrast to spectral-feature extraction methods such as band ratios, which involve a linear transformation of spectral values without incorporating ancillary data [52,53], empirical topographic corrections utilize slope and slope-azimuth angles computed from a DEM, and attempt to simulate the direct-beam irradiance using the cosine of the incidence angle. Some empirical methods rely on using statistical analysis to adjust spectral response values based on mean spectral values calculated for shaded and sunlit slopes, e.g., [54], for entire images or for selected cover types. Other methods use regression analysis to derive coefficients in an attempt to represent physical properties such as the diffuse-skylight irradiance, or the BRDF of surface cover types. These methods are all essentially empirical, as their calculations and generation of coefficients depend on the specific spatial and temporal characteristics of individual study areas (e.g., slope and slope azimuth, solar geometry, land cover).
Perhaps the most often used, and reportedly successful method adjusting spectral response based on mean values is the Statistical-Empirical Correction (SEC; see Equation (19b)) [55,56]. The SEC equation has also been referred to as the Teillet regression [14,18,57]. The SEC is a statistically based model that assumes that the mean spectral values on a horizontal surface are equal to those on tilted surfaces [58]. A number of authors have indicated the superior performance of the SEC when comparing topographic-correction algorithms, e.g., [59][60][61][62]. As with other topographic correction methods reported in the literature, different algorithms have been used to represent the SEC, e.g., [13,51,63].
The variable empirical coefficient algorithm (VECA) method (see Equation (22c)) was introduced by Gao and Zhang [14] and has been favorably compared to other topographic-correction methods [12,62]. The VECA was developed to overcome some of the challenges of applying more complex correction methods [14]. Similar to the SEC, the VECA algorithm adjusts the reflectance of pixels on an inclined surface based on a linear regression, incorporating cosi, and using the mean radiance values for pixels on horizontal and tilted surfaces [14,57,62]. The B-correction as proposed by Vincini et al. [64] can be calculated using either a linear isotropic (BLC, see Equation (20c)) or non-linear anisotropic form (BNC, see Equation (21b)). According to Gao and Zhang [14], the VECA and B-corrections proved to be capable of removing topographic effects in a Landsat ETM+ image, although the success of the algorithms were wavelength dependent and they each performed better in different bands [14].
If the Lambertian assumption accurately describes a surface, then topographic correction can be accomplished using a simple Cosine Correction (COSC, see Equation (23)) [39]. Numerous studies have found that the Cosine Correction overcorrected for local topographic effects, particularly on northern slopes, e.g., [65][66][67]. These overcorrections have been attributed to only E b being accounted for, as weakly illuminated areas can actually receive a significant amount of E d [56,68], and the inability of the Cosine Correction to account for BRDF variations of land cover surfaces [65].
Teillet et al. [55] modified the Cosine Correction with an additive term C (CCOR, see Equation (24c)) that the authors stated may emulate the effect of indirect illumination, although, the analogy is not exact. The C coefficient has been used to prevent the overcorrection of images [69], and notwithstanding the ambiguity of the original authors, has often been credited with accounting for E d , e.g., [16,26,61]. As a Lambertian model, the C-Correction does not take into account the BRDF of land cover features. Results from applying the C-Correction in studies comparing topographic correction methods have been mixed, with the correction performing well in some studies, e.g., [61,63], and not as well in others [51,70]. For example, Fan et al. [57] found the C-Correction to be theoretically and empirically inappropriate for forested terrain.
Topographic-correction methods have been applied to correct spectral-reflectance values for single and multiple cover types using sun-terrain-sensor (STS) geometry. To improve the representation of the geotropic nature of trees, Gu and Gillespie [71] introduced a method based on sun-canopy-sensor (SCS) geometry, which incorporates slope into the Cosine Correction (see Equation (25)) [54].
Application of the SCS method has often resulted in overcorrection, e.g., [14,69,72]. To address overcorrections, Soenen et al. [72] introduced the SCS+C method (see Equation (26)) incorporating the additive term C to attempt to account for E d , in a similar fashion to how the C coefficient adjusts the Cosine Correction [69]. Generally, the SCS+C correction has outperformed the SCS method, e.g., [18,69,73].
A commonly applied non-Lambertian topographic correction method that attempts to incorporate a representation of land cover BRDF is the Minnaert Correction (GMC, see Equation (27)). According to Smith et al. [39], the Minnaert constant k [74] is related to surface roughness. Calculation of k varies based on combinations of slope, slope azimuth and land cover, e.g., [9,11]. A global Minnaert k value is calculated per band using reflectance values from multiple cover types across an entire image, e.g., [75]. However, calculating the empirical coefficients C and k over large areas is problematic due to their dependence on land cover variation [57,61,62]. Attempts have been made to calculate k for land cover classes, e.g., [19], and using stratifications based on slope, slope azimuth, and land cover, such as the Modified Minnaert [16] and the Pixel-Based Minnaert (PMC; see Equation (28)) [76]. Minnaert corrections do not incorporate contributions from E d and E t [9].

ARC Evaluation
Effective topographic correction should result in a reduction of spectral variation globally across an image caused by multi-scale topographic effects, and an increase in local spectral variation in poorly illuminated areas due to biophysical property or land cover variations [15]. The performance of topographic-correction methods is dependent upon evaluation strategies [18]. Various approaches have been used to evaluate the effectiveness of topographic-correction methods; however, evaluation results are difficult to compare and there is still not a simple and objective approach for evaluating correction effectiveness [63]. Diagnostic methods generally fall within four categories according to a review of recent articles comparing topographic-correction methods. These categories include visual assessment, statistical methods, classification accuracy, and comparisons to simulated/synthetic images.
The most common and usually initial method for evaluating topographically-corrected images is a visual assessment. Visual assessments have been employed since correction methods were first attempted, e.g., [55] and continue to be applied, e.g., [13]. Usually singular images are evaluated, although image bands have been draped across DEMs to create an orthographic image for interpretation [60]. Even though visual assessment is used extensively as an evaluation method [16], it is fundamentally a subjective approach that is based heavily on the experience and skill of the analyst, e.g., [57] and must be combined with a quantitative assessment [62]. Visual assessment results may differ markedly from quantitative evaluations. For example, corrected images can appear similar and be statistically different [16,60]. In addition, variance compression may occur during topographic correction, which can visually be interpreted as a reduction of topographic effects [19].
Univariate and multi-variate statistical analysis methods have been applied towards evaluating the effectiveness of topographic-correction approaches including intraclass interquartile range reduction [63] and comparisons with in-situ observations of spectral reflectance [70]. Two of the most common methods for comparing imagery before and after correction include assessing a reduction in standard deviation (SD), and comparing the coefficient of variation (CV). Evaluations of the reduction of SD have been made across images [50] or within the same cover types on different slopes and aspects [18,61,73]. Ideally, an effective correction will maintain the mean of the original band or cover type and the SD will be reduced [50,66]. Essentially, these as well as other statistical diagnostic methods are testing for an increase in the homogeneity of reflectance values after correction, e.g., [50,62]. It is not valid, however, to presume spectral homogeneity across an image with different land covers or within the same land cover, such as vegetation with different stand ages, understory, or crown canopies, which can result in different spectral responses, e.g., [18]. Additionally, effective corrections do not always substantially change standard deviations and histogram distributions of imagery data [60]. Therefore, these methods may not account for the spatially-dependent variance structure across an image or within some land cover types.
A range of inferential-statistical methods have been used to evaluate the effectiveness of topographic corrections such as the univariate analysis of variance [69], multivariate analysis of variance [77], and homogeneity of variance [11]. Geo-statistical analyses, such as semi-variogram analysis [19], have also been applied. Perhaps the most common statistical method for evaluating preand post-corrected imagery is based on a linear regression (coefficient of determination, r 2 ) with the spectral values represented as the dependent variable and cos i as the independent variable, e.g., [39]. The spectral values have been represented as digital numbers (DN), e.g., [56,59,60,73] and radiance values, e.g., [14,63,70]. The assumption is that an effective correction will result in the removal of reflectance dependence on cosi. However, in areas where land cover distribution is affected by slope and aspect [61], this assumption is not valid [63]. In mountainous areas with large elevation gradients, spatial patterns of land cover will be dependent upon topographic variables, and a residual correlation would be expected between reflectance values and cosi after correction [18,63]. Additionally, cosi accounts only for local topography and does not incorporate multi-scale topographic effects [15]. Decorrelation of the linear regression is not sufficient to indicate a successful correction [57].
Two of the primary purposes for performing topographic corrections is to improve land cover classification [9,18,78] and biophysical parameter extraction [63]. An improvement in classification accuracy is a common method for determining the effectiveness of topographic correction. Classification methods have been evaluated such as maximum likelihood [26,50,59], artificial neural networks [50,54], and object-based classification [73]. The primary difficulty in using classification accuracy as a diagnostic method is that many factors and potential sources of error are inherent in the classification process [19,50]. Classification analysis is affected by study area characteristics, land cover types, within-class covariance, quality of reference data, sensors, DEMs, atmospheric effects and corrections, geometric corrections, classification algorithms, and topographic effects [26,57,62]. The uncertainties inherent in these factors and their various combinations make it difficult to relate quantitative classification results directly to the quality of topographic corrections, and to compare study results [26,63,79]. Classification provides an indirect evaluation [57], and is not a good indicator of topographic-correction effectiveness [19,50,57,79].
Given the shortcomings and difficulties in comparing results when using traditional evaluation strategies, the development of a more standardized and objective procedure is warranted. The comparison of a corrected image to a synthetic image could provide the basis for this type of approach [18,79]. The synthetic image represents idealized radiance according to specific conditions, and can be generated representing a flat surface or incorporating topography [79]. Quantitative evaluation of the comparison between corrected and synthetic images can be undertaken using a structural-similarity index [79,80], or comparable similarity indices [63,81]. Potential drawbacks to this approach include the complicated process required for producing synthetic images [81].

Study Area
We simulated a first-order approximation of the RTC over the Nanga Parbat, Himalaya in northern Pakistan ( Figure 1). This geographic region is excellent for testing empirical-correction procedures, as it exhibits extreme relief and steep slopes due to relief production caused by high-magnitude erosion and uplift [3,20,82,83]. We specifically simulated the RTC over a 60 km × 60 km area to evaluate the ability of empirical-correction procedures to account for the spatial complexity of an entire scene of ASTER imagery.
Topographic conditions over the region are highly variable due to complex spatiotemporal surface-process dynamics [3] that are governed by climate and topographic forcing. Glaciation, rapid river bedrock incision, and ubiquitous mass movements govern the dynamic nature of the sediment fluxes, and differential denudation coupled with uplift has resulted in significant relief production [83]. Consequently, surface conditions are very spatially variable from material composition and surface/landcover structural perspectives. Such conditions dictate anisotropic-reflectance conditions. Surface-process dynamics and landcover conditions are also regulated by two climatic systems that include the westerlies and the monsoon [20,82]. For more details regarding the study area characteristics, see Bishop and Shroder [20] and Bishop et al. [3].

Data
The ASTER Global Digital Elevation Model Version 2 (ASTER GDEM; [84]) was utilized to represent the topography and account for several multi-scale topographic effects governing the RTC. We spectrally simulated ASTER imagery (bands 1, 2, 3 and 5) for an entire scene which accounted for a 2000 × 2000 grid cell area centered over Nanga Parbat peak. We used a 4000 × 4000 dataset for computation of the cast shadows and topographic-shielding parameters to avoid edge effects. This study also utilized ASTER satellite imagery acquired over the area on 13 September 2004, to provide information about the spatial structure of the land cover characteristics that are described in greater detail below. Mean Exo-atmospheric irradiance values from Bird and Riordan [25] were used in simulating the RTC. Reference reflectance spectra were acquired from the United States Geological Survey Spectral Library Version 7 [85] and the ASTER Spectral Library [86].

Radiation-Transfer Modeling
We simulated ASTER multispectral imagery in order to evaluate frequently used empirical ARC methods. This involves implementing a first-order radiation-transfer model to simulate the RTC from the sun to the at-satellite sensor. We account for multi-scale topographic effects for each radiation-transfer parameter used in our simulation. Although higher-order modeling parameterization schemes can be utilized, they are not required to test and evaluate the most fundamental topographic effects, as the most common ARC methods account for only 1-3 proxy parameters representing the RTC. Furthermore, our approach permits detailed quantitative analysis of the significance of individual parameters and their wavelength dependence.

Exoatmospheric Irradiance and Solar Geometry
The mean exoatmospheric irradiance spectrum of Bird and Riordan [25] was modified using temporally dependent orbital parameters to compute the spatiotemporal dependent exoatmospheric irradiance spectrum (E 0 (λ)). Earth-sun orbital parameters control the magnitude of E 0 ; therefore, orbital forcing must account for eccentricity (e), obliquity ( ) and other parameters such as the longitude of perihelion [21]. We used the Berger [21] amplitudes, rates and phases of trigonometric expansions to predict these orbital parameters for 15 September 2018 to compute the Earth-Sun distance (d), the distance-correction factor ( f r ) and solar geometry at 10: local time, such that: where θ s is the solar-zenith angle.
Orbital parameters are also required to compute the spatial variation in the solar geometry such that: where δ is the solar declination and λ l is the true longitude of the Earth relative to the vernal equinox. The solar-zenith angle (θ s ) and solar-azimuth angle (φ s ) vary for each pixel if spatial and topographic factors are accounted for. The geocentric solar-zenith angle is: where ϕ is latitude, and H is the hour angle of the sun. This parameter, however, must be modified using a parallax correction to compute the apparent solar-zenith angle from the Earth's surface given topography. Parallax correction accounts for the radius of the Earth at a particular latitude, the height relative to the ellipsoid, and the distance from the sun, and is added to the geocentric solar-zenith angle. Finally, we accounted for atmospheric-refraction correction based on the parameterization scheme of U.S. Naval Observatory and Nautical Alamanc Office [87], where we accounted for the atmospheric profile of temperature and pressure of the US standard atmosphere [88]. The solar-azimuth angle for each pixel was computed as: where the sine and cosine components are computed as: and where α s is the solar elevation angle. The solar-azimuth angle was then corrected for grid convergence, which is a function of latitude and the longitude of a pixel location with respect to the central meridian of the projection used (i.e., Transverse Mercator). See Table A1 (Appendix A) for values utilized for modeling solar geometry and Table A2 (Appendix B) for explanation of symbol notation utilized throughout this paper.

Atmospheric Parameters
Atmospheric absorption and scattering processes account for atmospheric attenuation and the additive path-radiance components. These processes are wavelength dependent, and we simulated them from blue to short-wave infra-red wavelengths. We computed the downward and upward total transmittance parameters based on primary atmospheric constituents and processes. Specifically, they include Rayleigh scattering transmittance (T r ) using the parameterization scheme of Gueymard [24] that is wavelength and pressure dependent. Atmospheric pressure was computed based on the scale height of the atmosphere, which accounts for the virtual temperature, and the variation in gravitational acceleration, which is a function of latitude and altitude. Aerosol-scattering transmittance (T a ) was accounted for using the scheme of Bird and Riordan [25] using a rural aerosol model. Atmospheric absorption transmittance of water (T w ), ozone (T o ) and primary gases (T g ) were simulated using the parameterization schemes of Leckner [89]. All parameterization schemes effectively accounted for altitude variations in governing the pressure and relative optical air mass given by Kasten [90]. The atmospheric parameters used in these parameterization schemes to simulate atmospheric transmittance are listed in Table A1. Finally, the total downward (T ↓ ) and upward (T ↑ ) atmospheric transmittances were computed using Equation (1).
The additive path radiance (L p ) due to atmospheric scattering can be computed as [79]: where ρ a is the sky reflectance of the atmosphere. Different sky-reflectance models exist, and we used the parameterization scheme of Gueymard [24]. It accounts for scattering in the sensor field-of-view due to Rayleigh scattering and aerosol absorption and scattering. It also accounts for the single-scattering component of the scattering flux which is dependent upon the aerosol asymmetry factor coefficients found in Table A1. For more details regarding the equations used to compute ρ a , see Gueymard [24]. See Table A1 for all values utilized for modeling atmospheric constituents and properties.

Surface Irradiance
We accounted for the E b and E d in our simulations. We did not account for E t , as an accurate characterization of this component requires knowledge of the BRDF for each location, as the landcover structure and topography influence the BRDF [1]. Researchers have indicated that topographic effects have not been adequately accounted for in BRDF modeling efforts [37]. Although researchers have indicated that E t should be incorporated into ARC efforts in mountain environments [8,91], the uncertainties associated with land cover structural and topographic effects precludes accurate estimation of the magnitude of this component, given the high degree of spatial variability of morphometric properties at Nanga Parbat.
We computed E b using Equation (3), such that altitude, local topographic conditions, and the meso-scale relief structure in the direction of the φ s were accounted for. Specifically, we utilized linear regression as the basis for computing the slope and slope-azimuth angles for the computation of cos i, and implemented a ray-tracing algorithm that accounts for cast shadows (S).
The diffuse-skylight irradiance was computed on the basis of Equation (6), where we assumed isotropic scattering. Specifically, we utilized the parameterization scheme of Bird and Riordan [25] to estimate the diffuse irradiance on a horizontal surface at a particular altitude. Rayleigh and Aerosol scattering terms are utilized, governed by atmospheric transmittance parameters. We then equally partitioned the irradiance value over the hemisphere using one-degree zenith and azimuth angle intervals. The incoming irradiance geometry, in relation to the local topographic conditions were used to compute the local influence of the topography on incoming diffuse irradiance. We then accounted for hemispherical topographic shielding based on the maximum relief angle over a 25 km distance for each hemispherical azimuth direction. Knowing the maximum relief angle enabled us to sum and compute the diffuse irradiance over the appropriate zenith angle range for each azimuth direction. In this way, the full computation of local and mesoscale topographic affects were accounted for. Furthermore, we did this for each wavelength interval, as E d is highly wavelength dependent.

Surface Reflectance
Surface reflectance was simulated based on characterizing the fundamental spatial structure of land cover conditions. To accomplish this, we utilized a Level 1A ASTER scene. Our first objective was to map the basic land cover types that would serve to spatially constrain spectral-mixing modeling efforts to account for surface biophysical properties. These cover types include water, snow, ice, vegetation, rocks and sediment, and desert varnish.
The ASTER 3/4 band ratio was utilized to map snow and ice cover distribution similar to other researchers, e.g., [92,93]. Simple thresholding from 0.18 to 0.35 and >0.35 was used to classify ice and snow, respectively. The optical-band differential index [94] was used for water mapping with a value >−0.025 used to identify surface water. For vegetation, we used the normalized-difference vegetation index (NDVI) for basic mapping [95]. NDVI values greater than 0.1 were used to classify vegetated areas. For desert varnish, we used unstandardized principal-component analysis on all spectral bands in the VNIR, SWIR, and TIR regions of the spectrum. The Principal Component 5 image highlights the spatial distribution of desert varnish very well, located in the arid low-altitude portions of the Indus and Astor river valleys.
Each individual classification map was sequentially merged to produce a composite map, such that the remaining unclassified pixels in the composite map were classified as rock and sediment. A majority filter was iteratively applied to the composite map until reclassification results did not change. This produced a reasonable distribution of relatively homogeneous land cover classes that corresponds well to the land cover structure over Nanga Parbat (Figure 2). Surface reflectance was simulated using a spectral-mixing model based on compositional variations within each land cover class. The multiple end-member linear spectral-mixing analysis (LSMA) model suffices to mix the spectra of multiple kinds of matter [96]: where ρ c s is the composite surface reflectance, f j is the fraction of end member j, ρ j is the reflectance of end-member j, and e j is an error term that we assume is 0.0. All end-member spectra were resampled to a common reference: the wavelengths tabulated by Bird and Riordan [25]. For each land cover class, we accounted for matter composition variation that has been reported in the literature, and randomly assigned the relative percentages or fractions of matter constituents for each pixel. In this way, we ensured that the spatial variation of surface reflectance caused by land cover structure and compositional variations in surface matter are different from the spatial variance structure of topographic effects that govern other radiation-transfer parameters. This is essential to determine the effectiveness of empirical ARC methods to remove only topographic effects and accurately represent spatial variation in surface reflectance, given the possibility of spectral compression or over-correction.
Rock and sediment were both simulated utilizing relative mineral abundances described for Nanga Parbat leucogranites and sediment [97][98][99]. The relative mineral abundances provided by these authors, although not always explicitly stated, were probably determined by petrologic analysis under a microscope. As such, they represent an areal fraction suitable for spectral mixing. Otherwise, the assumption was made that relative abundances by weight or volume approximate areal abundances. Because these studies were primarily interested in rock or river sediment, some extrapolation was necessary to decrease albedo and qualitatively approach the reference satellite imagery. This was primarily accomplished by increasing the fraction of clay, iron oxides, and heavy minerals. Illite was chosen to represent clay minerals, hematite for iron oxides, hornblende for heavy minerals (primarily amphibole minerals), oligoclase for plagioclase feldspars, orthoclase for non-plagioclase feldspars, and quartz for lithics (Table 1). Mineral spectra were acquired from the USGS Spectral Library v7 [85] and resampled via linear interpolation to 0.005 µm. Table 1. End-member minerals used to simulate reflectance values for rocks and sediment and desert varnish in the Nanga Parbat region. End-member fractions were randomly determined from the indicated range in the order shown in the table. The remainder was then divided among the "remainder" end-members. See Table A3 (Appendix C) for the names of USGS Spectral Library files utilized for these simulations.  [101]. Desert varnish composition was interpreted from chemical composition, assuming that relative abundances by weight approximate areal abundance. Spectra were acquired from the USGS Spectral Library v7 [85]. As the spectral library lacked an entry for Birnessite (manganese oxide), the spectra for Desert Varnish found in the database was used to represent the spectra of Birnessite. We did this to introduce compositional variation based on fractional mineral abundance (Table 1).

End-Member Rocks and Sediment Desert Varnish
For vegetation, we did not account for variation of sagebrush, coniferous, deciduous, and tundra plant species, as these are generally not available in most spectral libraries. Instead, we utilized the spectra of the Lodgepole Pine from the USGS Spectral Library v7 [85] as the basis for characterizing vegetation reflectance, regardless of variations in species. We assumed that the magnitude of the NDVI was related to the spatial leaf-area coverage for a pixel, and that the remaining coverage represented sediment with no ground cover vegetation. We randomly assigned the composition of the sediment using the minerals presented in Table 1. In this way, we could account for the spatial distribution of the presence of vegetation, and account for variations in leaf coverage and surface composition. It is important to note that the objective was not to accurately characterize the biophysical properties of vegetation canopies at Nanga Parbat, but to account for spatial biophysical variation in surface reflectance due to leaf-area coverage and different ground matter conditions. Surface water at Nanga Parbat primarily consists of highly turbid streams that originate at the termini of glaciers surrounding the mountain. They drain into the Astor River and the Large Indus River on the north side of Nanga Parbat. Water spectra were acquired from the ASTER Spectral Library v2.0 [86], and we uses the tap water spectra (tapwater.jhu.becknic.spectrum) rather than distilled water because of the lack of data for the latter for wavelengths less than 4 µm. Turbid-water spectra for each pixel is simulated such that the water fraction was randomly set between 0.5 and 0.6 and the remainder was set as a fixed sediment composition (Table 1).
For snow, we modeled the hemispherical albedo (α H ) for non-isotropic scatterers as a function of grain size and wavelength, according to Hapke [102] as follows: where µ 0 = cos θ s , b is the asymmetry factor (second coefficient of the phase function expansion in Legendre polynomials, as shown by Hapke [102]), and ω is the single-scattering albedo. We modeled b and ω as functions of grain size, wavelength, and complex refractive index, using Mie scattering theory [103]. Equation (15a) can also be used for isotropic scatterers if the second term is not considered.
For the purposes of this study, we generated a database of hemispherical albedos using Equation (15a), for both isotropic and anisotropic scatterers. The wavelength ranges from 0.4 to 2.5 µm (with variable steps with values around 0.005 ± 0.002), and grain size ranges from 1 to 1000 µm. Finally, we compare our results with those of Wiscombe [103]. The shape of our curves are extremely similar, with small differences in the values probably due to the use of a slightly different complex refractive index.
To decrease the degrees of freedom for the simulation, the solar-zenith angle was assumed to be 0. For each pixel classified as snow, grain size was randomly selected between 1 and 1000 µm. Ice was simulated using the spectrum for the 1000 µm grain size, and is the only invariant class of land cover for reflectance simulations. All of the reflectance spectra were then sampled to conform to the spectral resolution of the ASTER sensor. We used the spectral response functions of the ASTER spectral bands to produce a weighted-average spectral reflectance value (ρ s ) for each pixel. This is required to produce reference data for validation of ARC methods in reducing topographic effects. Simulated spectra were computed as follows:ρ where n represents the number of wavelengths sampled by the spectral-response function, and w represents the spectral-response coefficient weight that characterize the nature of the spectral bandwidth. The spectrally-weighted reflectance images are presented in Figure 3.

Simulated ASTER Imagery
All of the aforementioned radiation-transfer parameters were simulated from 0.4 to 3.0 µm for each pixel in the scene. We assumed that surface reflectance is Lambertian, and computed the at-satellite radiance (L 0 ) for the full spectrum as: We then applied the ASTER spectral-response functions to each L 0 spectra and used Equation (16) to simulate the ASTER spectral images. The spectral characteristics of the simulated images are presented in Table 2. Although we recognize that surface reflectance characteristics at Nanga Parbat are anisotropic in nature, we did not attempt to simulate the BRDF given a lack of data and our objective to evaluate ARC methods to account for topographic effects in E b , T ↓ , E d , T ↑ , and L p . Table 2. Spectrally Simulated ASTER imagery characteristics. Mean exoatmospheric irradiance at 1 astronomical unit (E 0 ) from Thome et al. [104]. It is important to note that we did not simulate the spatial or radiometric resolution of the ASTER sensor, but forced the spatial resolution to coincide with the resolution of the DEM. We also assumed that the sensor could account for the magnitude of the simulated at-satellite radiance values, even though the ASTER sensor exhibits 8-bit radiometric resolution.

Anisotropic Reflectance Correction
We evaluated a variety of empirical ARC methods. To do so, our evaluation considered an evaluation of whether atmospheric correction is required for suitable topographic correction. Given our RTC simulation, perfect atmospheric correction was applied before topographic correction as follows: The ARC methods we evaluated are as follows, where L n represents normalized imagery after applying the technique: 1. Statistical-Empirical Correction [55,56]. Assuming a linear correlation between at-satellite radiance and local topographic conditions: where β 0 and β 1 represent the y-intercept and slope coefficient of the linear regression, respectively. 2. B-Correction assuming isotropic reflectance conditions [14].
where β 0 and β 1 represent the y-intercept and slope coefficient of the linear regression, respectively, and x represents the residuals of the radiance data to the predicted regression line. 3. B-Correction assuming a non-linear relationship and an anisotropic reflectance model [105].
where β 0 and β 1 represent the y-intercept and slope coefficient of the linear regression, respectively. 4. Variable Empirical Coefficient Algorithm [14].
where C is computed using Equation (24b). 9. Global Minnaert Correction [39]: where e is the exitance angle (e = θ t for nadir viewing), and k is a globally-derived dimensionless Minnaert coefficient that is wavelength dependent and ranges from 0 to 1. It was calculated using least-squares regression on the variables x and y, where x = log(cos i cos e) and y = log(L 0 cos e).
The slope of the regression equation represents k. The correction procedure defaults to the Lambertian assumption when k = 1.0. 10. Pixel-Based Minnaert Correction [26,76]: where k is computed for slope ranges using an interval of five degrees, based on regression analysis using variables x and y.

Validation and Statistical Analysis
The evaluation of ARC methods is notoriously difficult due to subjective approaches that have been typically used. The magnitude and spatial distribution of L n resulting from topographicallycorrected imagery must be correct to ensure accurate biophysical remote sensing and thematic mapping capabilities. Subjective approaches including visualization, global image statistics and other techniques do not make use of reference data and account for magnitude and correlation structure. Therefore, we compared our simulated surface reflectance values,ρ s , to computed surface reflectance values (ρ n ) based on the topographically corrected radiance values (L n ). Given that we have the exact magnitude of the surface irradiance component simulated, under the Lambertian reflectance assumption: Our validation and statistical analysis approach was as follows: 1. First, we evaluated the image magnitude of reflectance differences using the Root Mean Squared Error (RMSE) as: 2. We computed the magnitude of the Pearson-Product Moment correlation coefficient (r), to ensure that the spatial variability in reflectance accounts for reflectance variation due to land cover structure and biophysical variation as: 3. We then utilized a structural-similarity index (SSI) developed and evaluated by Wang et al. [80] to determine the collective influence of reflectance magnitude (l), global variation (c; contrast) and the correlation coefficient (r), on the topographically corrected reflectance values compared to simulated reflectance. It was computed as follows: where C 1 = (K 1 R) 2 , K 1 was set to 0.01, and R is the dynamic range of an image, which we set to 255, given scaled reflectance images. The contrast was computed as: where C 2 = (K 2 R) 2 , K 2 was set to 0.03, and R is the dynamic range of an image. Finally, the SSI was computed as: where the SSI index is a coefficient that ranges from 0 to 1, where a 0.0 value indicates no structural similarity, and a value of 1.0 represents identical values, spatial variance and correlation structure. We computed this for each topographic-correction procedure for each spectral band. In this way, we have an integrated metric which accurately evaluates an entire image. 4. Nevertheless, to better characterize the spatial nature of how each topographic-correction procedure performs under different topographic conditions, we implemented the SSI over an arbitrary window size of 11 × 11 to ensure an adequate sample size for computing the statistic. We generated SSI images to evaluate where the topographic-correction procedures worked well, and where they did not. This provided for better interpretation of the results in terms of what topographic conditions cause problems for various techniques.
Perfect topographic correction of imagery should result in a RMSE value of 0.0, the correlation with the surface-reflectance benchmark should be 1.0, and the SSI value should be 1.0 for the globally and locally computed indices. Evaluating the departures from perfect numbers related to application requirements was not our objective, although our evaluation strategy permitted us to determine the degree to which this radiometric calibration issue is addressed over our geographic region.

RTC Simulation
The direct-beam irradiance image is presented in Figure 4. As expected, the magnitude of E b decreases with longer wavelengths and local topographic effects cause significant spatial variation in irradiance. The northern slopes receive less energy than southern slopes, and cast shadows exist at the highest altitudes on the northern face of Nanga Parbat. Magnitude variations in E b are larger for shorter wavelength regions, and decreases at longer wavelengths. This can be seen when comparing the green and SWIR bands ( Figure 4A,D). The diffuse-skylight irradiance is presented in Figure 5. In general, the magnitude and spatial variability of E d decreases with increasing wavelength. This is related to the significant reduction in atmospheric scattering at longer wavelengths. Meso-scale topographic effects can clearly be seen, as the magnitude of E d is greatest at higher altitudes in glacier accumulation areas, where less topographic shielding occurs. Glacier erosion and river incision generate greater relief that reduces this irradiance component. This can be seen over the Indus River valley in the northwestern portion of the scene.  Figure 6 depicts the atmospheric path-radiance images for ASTER spectral bands. The magnitude of this component is highly wavelength dependent, as greater atmospheric scattering occurs at shorter wavelengths, and lower altitudes increase the optical depth that facilitate greater volumetric atmospheric scattering and higher values of L p . This can be clearly seen in the Indus Valley that exhibits the lowest altitudes in the region ( Figure 6A).
The at-satellite radiance images (L 0 ) are presented in Figure 7. Simulated images for the ASTER spectral bands clearly depict the multi-scale topographic effects that are inherent in satellite imagery acquired over this region. The images are extremely realistic and visually similar to real imagery, depicting extreme topographic and land cover spectral variations. Univariate statistics for all the simulated radiation-transfer parameters are presented in Table 3.

ARC Validation
Validation statistics for characterizing the empirical ARC methods are presented in Table 4. In general, the so-called statistical-empirical methods did not produce reasonable results. For example, the SEC method produced results with relative high RMSE and low r 2 values for all spectral bands. Global SSI values were approximately 0.0, with or without the use of atmospheric correction. Similarly, for the VECA method, global SSI values were 0.0 for all spectral bands, with or without the use of atmospheric correction, although RMSE and r 2 values were exceptionally high, demonstrating that the scaling factor can account for topographic effects given the Lambertian assumption, although the technique does not produce reasonable values from a magnitude perspective. Our results indicate that the B-Correction method performed the best out of this category of approaches. Specifically, the BNC method outperformed the BLC method with lower RMSE values, and higher r 2 and SSI values for all spectral bands. For both methods, improved results were obtained given that imagery were atmospherically corrected before applying topographic correction.
Results for so-called semi-empirical methods assuming Lambertian reflectance are also presented in Table 4. Results were mixed, although several methods generated superior results compared to the other methods evaluated. The COSC and SCS methods performed the most poorly for this category, with relatively high RMSE values, and low r 2 and SSI values for all spectral bands. Atmospheric correction before topographic correction did not significantly affect the results for these methods. Conversely, the CCOR and SCS+C methods performed the best given atmospherically corrected imagery, with relatively low RMSE values and high r 2 and SSI values for all spectral bands. Both techniques remove almost all the multi-scale topographic variation, and the corrected imagery represents the inherent spectral variation related to land cover structure and biophysical variation. It is important to note, however, that atmospheric correction makes a significant difference and accounts for approximately 10-35% of the explained variability due to topographic coupling that governs atmospheric conditions, depending on the spectral band.
The results for so-called semi-empirical methods assuming non-Lambertian reflectance are also presented in Table 4. The GMC method outperformed the PMC method, and atmospheric correction improved topographic correction results for the GMC method. It did not, however, outperform the BNC, CCOR and SCS+C methods. These results are not surprising given that we utilized the Lambertian assumption in our simulations, and that a proxy parameter for E d is not included in the GMC and PMC parameterization schemes.
Statistics characterizing the local SSI index for each method (given atmospheric correction) evaluated are presented in Table 5. The minimum and maximum values of the local SSI index reveal that all methods produced good and bad correction results, suggesting that these methods produce highly variable results over a large geographic region depending upon the spatial variability in morphometric properties. Perhaps the best metric to evaluate is the mean local SSI value that provides an overall assessment of local topographic-correction effectiveness. In this regard, it is clear that the CCOR and the SCS+C methods produce the most superior results, with the SCS+C method technically outperforming the CCOR. Examination of SSI images indicate that both methods effectively correct for multi-scale topographic effects in simulated imagery. It is interesting to note that these methods do not perform as well on the SWIR spectral band, compared to the other bands, as numerous locations exhibit local SSI values <0.5. These locations are associated with relatively high-altitude ridge tops and steep slopes near ridges throughout the scene. This suggests that the methods do not adequately account for the diffuse-skylight irradiance at longer wavelengths, given the use of a proxy parameter that better represents diffuse irradiance at shorter wavelengths given more atmospheric scattering.
Most other methods exhibit relatively low local SSI values throughout the scene. Topographic correction results are highly variable, and in general, topographic correction is superior on sunlit slopes. For example, the SEC method that has been noted for being superior for topographic correction cannot account for multi-scale topographic effects. Figure 8 clearly depicts very low local SSI values on north, northwestern facing slopes (red), whereas sunlit slopes exhibit higher SSI values (green).
Methods that produce higher mean SSI values, such as the BNC and the GMC methods produce results that appear to be wavelength and spatially dependent. For example, the BNC method ( Figure 9) works better for the green and red regions of the spectrum, but the diffuse-irradiance component cannot be accounted for at high altitude on steep slopes. The results for the NIR and SWIR bands are more spatially variable than results in shorter wavelengths and the method does not effectively account for less variation in the diffuse irradiance in these regions of the spectrum. Similarly, Figure 10 depicts high variability in SSI values for all spectral bands from correction using the GMC method, although the spatial variability in accuracy is higher in the NIR and SWIR bands. Figure 8. Local structural-similarity index for SEC corrected imagery given atmospheric correction: (A) SSI for ASTER spectral band 1 (green); (B) SSI for ASTER spectral band 2 (red); (C) SSI for ASTER spectral band 3 (NIR); and (D) SSI for ASTER spectral band 5 (SWIR). All SSI images are overlaid onto the SEC corrected imagery for that spectral band. Figure 9. Local structural-similarity index for BNC corrected imagery given atmospheric correction: (A) SSI for ASTER spectral band 1 (green); (B) SSI for ASTER spectral band 2 (red); (C) SSI for ASTER spectral band 3 (NIR); and (D) SSI for ASTER spectral band 5 (SWIR). All SSI images are overlaid onto the BNC corrected imagery for that spectral band. Figure 10. Local structural-similarity index for GMC corrected imagery given atmospheric correction: (A) SSI for ASTER spectral band 1 (green); (B) SSI for ASTER spectral band 2 (red); (C) SSI for ASTER spectral band 3 (NIR); and (D) SSI for ASTER spectral band 5 (SWIR). All SSI images are overlaid onto the GMC corrected imagery for that spectral band. Table 4. Statistical summary of anisotropic-reflectance correction methods for spectrally simulated ASTER imagery. Methods include the statistical-empirical correction (SEC), B correction (linear form; BLC), B correction (non-linear form; BNC), variable empirical coefficient algorithm (VECA), Cosine Correction (COSC), C correction (CCOR), sun-canopy-sensor correction (SCS), sun-canopy-sensor plus C correction (SCS+C), global Minnaert Correction (GMC), and pixel-based Minnaert Correction (PMC). Surface radiance related parameters include the minimum surface radiance (L min , [W m −2 sr −1 µm −1 ]), maximum surface radiance (L max , [W m −2 sr −1 µm −1 ]), surface radiance (L, [W m −2 sr −1 µm −1 ]), and standard deviation of surface radiance (L σ ). Evaluation metrics include the root-mean-squared-error (RMSE), coefficient of determination (r 2 ), and structural-similarity index (ssi).   Table 5. Univariate statistics for structural-similarity index (SSI) images based on comparison of simulated reflectance and normalized reflectance given atmospheric correction. Statistics parameters include the minimum (Min), maximum (Max), mean (µ), and standard deviation (σ) of local SSI values in the image.

Discussion
Radiation-transfer simulations and our theoretical assessment of empirical ARC methods reveal new insights into the nature of research involving the evaluation of topographic-correction methods. Our results clearly reveal that the methods we evaluated did not reasonably account for multi-scale topographic effects, countering results presented in the literature. The magnitude of reflectance and the correlation structure of biophysical properties were not preserved in the topographically-corrected imagery. Our best results are obtained by using the CCOR and SCS+C methods, although comprehensive atmospheric correction is required.
We are surprised that the CCOR and SCS+C methods accounted for multi-scale topographic effects related to the diffuse-skylight irradiance component. The C parameter is an empirical coefficient that is thought to be able to account for the variation in E d . Although the authors who introduced the CCOR were not definitive in their description of the C parameter [55,72], it has since become commonly accepted in the literature that C accounts for E d , e.g., [16,26,57]. This is not possible, however, because it is a wavelength dependent constant that does not covary with E d . It is used to modify θ s and cos i, such that the ratio-scaling factor is highly correlated to E b + E d , thereby accounting for the primary topographic effects in our simulation. This suggests that θ s accounts for some variation in E d related to diffuse scattering, and that cos i also accounts for variation in topographic shielding given steep slopes, as E b is relatively low with steeper slopes. Nevertheless, we should expect both of these methods to be totally inadequate for removing multi-scale topographic effects in real imagery, given that a multitude of other topographic factors govern E t , the BRDF, and therefore other secondary radiation-transfer parameters, as all of the methods evaluated cannot account for these more complex radiation-transfer parameters due to their limited parameterization schemes.
Our results strongly suggest that these empirical methods cannot be used to address the problem of multi-scale topographic effects, given the simplicity of their parameterization schemes. Other multi-scale topographic effects must be accounted for [8], as well as the anisotropic nature of atmospheric scattering on E d and the BRDF [8,106]. The literature assumes that there is a best correction method for different environments and that investigators should evaluate a variety of methods on a case-by-case basis and use the best one, e.g., [18]. This assumption is most likely false given that high spatial variability in topographic and land cover conditions governs the magnitude and spatial variation in wavelength dependent radiation-transfer parameters, as our simulation demonstrates. Furthermore, this assumption is based on the belief that researchers can accurately evaluate topographic-correction procedures, accounting for magnitudes and the spectral correlation structure of land cover and biophysical variations, which is effectively a radiometric calibration issue. This information over an entire geographic region is generally not available in mountain environments. Therefore, common evaluation procedures are inadequate. Recently, investigators have begun to recognize this fact and have utilized a multi-criteria-based approach where numerous evaluation methods are considered, e.g., [63]. We note, however, that these approaches cannot account for variations in numerous radiation-transfer parameters that are governed by topographic effects.
Fundamentally, the inherent nature-of-the-problem stems from inadequate parameterization schemes and inadequate evaluation approaches. ARC is fundamentally a radiometric calibration issue that requires the removal of spectral variation caused by coupled atmosphere-topography effects. Evaluation procedures must account for various radiation-transfer parameters in order to remove those image information components. The use of previously reported evaluation procedures are highly subjective given unrealistic assumptions. The issues include: 1. The use of the human visualization system to evaluate empirical ARC methods. Decreased spectral variation is thought to signify effective topographic correction. This approach cannot be used to evaluate the degree of information compression given the use of a scaling factor, the correct magnitude of surface reflectance, or the spatial variation structure of surface spectral reflectance. Furthermore, interpretation is subject to biases resulting from a-priori experience, generalization caused by 8-bit display monitors, color bins, and graphic production techniques. 2. The concept of similar surface reflectance conditions for homogeneous landcover is assumed to be the result of effective ARC. This assumption is not realistic in mountain environments because of biophysical property variation related to topography and land cover structure. Environmental dynamics in mountain environments determine the magnitude and extent of various processes which dramatically alter surface irradiance, surface composition and the BRDF. These include processes such as evapotranspiration, gravitational sediment fluxes, erosion, deposition and slope stability that governs plant dispersion and mineralogical composition. In addition, climate forcing and variations in temperature and precipitation govern vegetation density, glaciation and snow cover. ARC methods need to work over large spatial and temporal scales, therefore the operational scale dependencies of processes may vary significantly, and the assumption that homogeneous land cover conditions will exhibit similar reflectance on illuminated and shadowed terrain, which serves as a basis for evaluation, does not account for spatial variation in surface-process dynamics and radiation-transfer theory. 3. The relationship between cos i and L should not be used to evaluate the effectiveness of topographic correction, as it does not account for numerous radiation-transfer parameters (e.g., E d and E t ) that exhibit spectral variation at different scales and magnitudes than cos i. In addition, land cover patterns are dependent upon topographic variables and residual correlations between cos i, and spectral reflectance values would be expected following topographic corrections [18,63].
Furthermore, as we have demonstrated, a global image evaluation procedure cannot be used to account for local effectiveness of a topographic-correction method. 4. It is frequently assumed that spectral homogeneity should increase after topographic correction, e.g., [18]. This assumption is not consistent with the concept of enhancing the moderate-to high-frequency surface reflectance conditions, as lower-frequency topographic information is reduced. This was demonstrated by Bishop and Colby [15], who used semi-variogram analysis to evaluate surface reflectance variability after topographic correction. This is especially the case in more dynamic mountain environments where moisture variations, sediment fluxes, and anthropogenic forcing can cause significant compositional and biophysical variation that can increase spectral variation over various spatial scales. 5. The utilization of classification accuracy as a way to evaluate topographic correction.
This approach is extremely problematic and does not address numerous radiometric calibration issues related to the fundamental radiation-transfer components. The magnitude and spatial variation of reflectance is also not accounted for. It introduces significant uncertainty given various assumptions, and is influenced by overcorrection or information compression, the nature of the algorithm, training, as well as dimensionality of the feature space and specific information components making up that feature space. Therefore, many potential sources of error are intrinsic in the classification process [19,50]. Given high spectral variability in mountain environments, and the availability of high resolution imagery, good classification accuracy can be difficult to obtain given increased spectral variation. Nevertheless, good classification results do not equate to optimized ARC results that need to be representative of surface biophysical property variation to enable the prediction of surface parameters and generation of thematic information. 6. Issues of empiricism associated with using parameter values and thresholds to modify results for a particular geographic region. Often, such attempts are applied over relatively small geographic areas, and their utilization may become questionable over larger regions due to more spatial variability in morphometric properties of the topography and RT parameters. We would argue that, given the paucity of data over large areas by which to formally validate ARC results, investigators cannot effectively select the appropriate magnitude or interval range of parameters in an attempt to optimize results. For example, empirically modifying results based on visual interpretation does not effectively address the removal of radiation-transfer components from the imagery, as the parameterization schemes and modification procedures cannot deterministically account for other topographic effects. It also results in non-repeatable results, as such procedures may not produce consistent results in different environments, at different times, or with different analysts. Consequently, empiricism is also partially responsible for inconsistent results presented in the literature. 7. Research involving the ranking of empirical ARC procedures. We would argue that this should not be a research objective, as traditional statistical and classification evaluation methods do not account for the reduction of radiation-transfer components encapsulated in multispectral imagery. As demonstrated in our simulation, radiometric calibration is required and the results are highly dependent upon the wavelength region, time, and morphometric properties of the topography.
Most of these factors have not been systematically evaluated. Rather, investigators assume that more evaluation procedures are better and produce reliable results when, in fact, the radiometric calibration issues are not being addressed using simplistic empirical ARC procedures. Rankings which are based on different evaluation approaches also lead to inconsistencies in the literature. 8. Mountain topographic and spectral variability. The literature is ubiquitous with statements indicating that investigators are evaluating ARC methods in complex topography. Complexity needs to be semantically defined, as indicated by Bishop and Dobreva [107], who evaluated the term from morphometric property and system perspectives. Clearly, the magnitude and spatial variability of morphometric properties govern the multi-scale topographic effects that partially regulate irradiance and reflectance variations. As demonstrated, even secondary or tertiary parameters can significantly influence results if they are not formally accounted for (i.e., need for atmospheric correction given relief variations). Similarly, the radiation-transfer process dynamics need to be accounted for in order to represent their variation adequately so that it can be removed. Consequently, we could expect that various empirical ARC methods will better function in environments with less topographic and spectral variability. We would argue that the concept needs to be evaluated from a radiation-transfer systems perspective. Therefore, researchers need to account for the morphometric properties of the topography within their study areas. 9. Small angle approximation. Most investigators utilize constant solar geometry parameters for atmospheric correction, the cosine of the incidence angle and ARC methods. Solar geometry parameters vary for each pixel and are a function of time, latitude, longitude, and altitude. The computation of these parameters is more important for processing large regions or scenes compared to small areas. Nevertheless, it is assumed that the variance is insignificant. It is yet unclear as to the impact that this assumption has on ARC, as it is a parameter that partially governs atmospheric transmittance, path radiance, direct irradiance, diffuse-skylight irradiance, adjacent-terrain irradiance and the BRDF. 10. Performance dependence on land cover.
Numerous researchers have indicated that the effectiveness of empirical ARC methods is largely dependent on land-surface types, e.g., [17,81,108]. We would argue that this is not technically correct, as effective ARC depends upon a multitude of factors including: (1) parameterization scheme of technique/model; (2) spatial variability in morphometric properties; (3) landcover type and structure; and (4) sensor characteristics, such as wavelength. Clearly, the morphometric properties of the topography coupled with atmospheric conditions and the land cover surface structure will collectively govern the irradiance fluxes and the BRDF. Sensor-system characteristics will then determine the level of generalization associated with recording the RTC. Empirical ARC methods do not formally account for all of the radiation-transfer components that are caused by multi-scale topographic effects, especially the BRDF. It has yet to be conclusively determined to what degree the BRDF accounts for spectral variation in imagery. In such dynamic environments, areal and intimate mixtures of matter must also be accounted for, and is another factor governing the BRDF.
Our simulation results strongly suggest that inconsistencies in research results can be understood on the basis of environmental complexity, as it governs radiation-transfer parameters. Less complex topography, and more spatially homogeneous land cover conditions, may dictate that a reduced number of radiation-transfer parameters need to be accounted for, and that E b may be the most dominant factor. Under these conditions, various empirical ARC methods may work to some degree. As the environmental complexity increases, however, different approaches may need to be used. Nevertheless, ARC investigations do not generally account for a multitude of multi-scale topographic effects and parameters related to atmospheric conditions, cast shadows, E d , E t , or the BRDF. Although empirical methods that assume non-Lambertian reflectance have been evaluated, the parameters that govern the BRDF are not fully accounted for. Furthermore, parameterization scheme errors have propagated in the literature, and researchers may have utilized different or incorrect equations for some ARC methods. To complicate matters, researchers have demonstrated that results are temporally variable, e.g., [54].
The aforementioned issues demonstrate the multi-faceted complexity of ARC. Numerous radiation-transfer components are collectively governed by atmospheric conditions, multi-scale topographic effects and landcover dynamics. Characterizing and removing these variance components from multispectral imagery is notoriously difficult and will require radiation-transfer modeling. We should not expect the use of a few proxy parameters to address the nature-of-the-problem given the anisotropic nature of irradiance and reflectance. Furthermore, more research is warranted on developing and evaluating new parameterization schemes that better account for variations in the components of the RTC. More research investigating topographic effects on the BRDF is sorely needed. Only in this way will we be able to better understand which radiation-transfer components must be accounted for in ARC methods and models.

Conclusions
Remote sensing of mountain environments is challenging, as multi-scale topographic effects govern the anisotropic nature of the radiation-transfer cascade. Although anisotropic-reflectance correction investigations have been ongoing for approximately 40 years, inconsistent research results caused by inadequate evaluation approaches have resulted in little progress on effective topographic correction of multispectral imagery in mountain environments. Consequently, we diagnostically evaluated commonly used empirical ARC methods using spectrally simulated ASTER imagery over the Nanga Parbat, Himalaya. Specifically, we used first-order radiation-transfer modeling to account for a variety of spectral, spatial and temporal factors including orbital dynamics, atmospheric absorption and scattering, direct-and diffuse-skylight irradiance, land cover structure, and surface biophysical variations.
As expected, we found that the dominance and spatial variability of various radiation-transfer components were wavelength dependent. Our results clearly reveal that most of the empirical ARC methods we evaluated could not effectively or reasonably account for the multi-scale topographic effects at Nanga Parbat. The magnitude of reflectance and the correlation structure of biophysical properties were not preserved in the topographically-corrected multispectral imagery. Using a Lambertian assumption, better results were obtained using the CCOR and SCS+C methods, although comprehensive atmospheric correction was required, and we did not account for other primary and secondary topographic effects that are thought to significantly influence spectral variation in imagery acquired over mountains (adjacent-terrain irradiance and the BRDF). Evaluation of structural-similarity images revealed spatially variable results that are wavelength dependent. Collectively, our simulation and evaluation procedures strongly suggest that empirical ARC methods cannot be used to address the problem of multi-scale topographic effects in mountain environments, given the simplicity of the parameterization schemes and their inability to adequately account for anisotropic irradiance and reflectance.
Characterizing and removing radiation-transfer parameter variance components from multispectral imagery is notoriously difficult and will require radiation-transfer modeling efforts. More research is warranted on developing and evaluating new parameterization schemes that better account for variations in the components of the RTC. More research investigating topographic effects on the BRDF is sorely needed. The next phase of our research will be to incorporate more primary, secondary and tertiary parameters of the RTC to evaluate ARC methods.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: