Improved Topographic Normalization for Landsat TM Images by Introducing the MODIS Surface BRDF

In rugged terrain, the accuracy of surface reflectance estimations is compromised by atmospheric and topographic effects. We propose a new method to simultaneously eliminate atmospheric and terrain effects in Landsat Thematic Mapper (TM) images based on a 30 m digital elevation model (DEM) and Moderate Resolution Imaging Spectroradiometer (MODIS) atmospheric products. Moreover, we define a normalized factor of a Bidirectional Reflectance Distribution Function (BRDF) to convert the sloping pixel reflectance into a flat pixel reflectance by using the Ross Thick-Li Sparse BRDF model (Ambrals algorithm) and MODIS BRDF/albedo kernel coefficient products. Sole atmospheric correction and topographic normalization were performed for TM images in the upper stream of the Heihe River Basin. The results show that using MODIS atmospheric products can effectively remove atmospheric effects compared with the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model and the Landsat Climate Data Record (CDR). OPEN ACCESS Remote Sens. 2015, 7 6559 Moreover, superior topographic effect removal can be achieved by considering the surface BRDF when compared with the surface Lambertian assumption of topographic normalization.


Introduction
Optical remote sensing has been widely used to quantitatively estimate surface parameters [1][2][3], such as various biogeophysical (e.g., leaf area index (LAI)) and geophysical parameters (e.g., albedo).Landsat imagery has been used extensively for quantitatively evaluating surface parameters [4].However, the satellite-obtained top-of-atmosphere (TOA) radiance is reflected via complex interactions between (1) atmospheric attenuation; (2) incoming solar irradiance at a given surface point, and (3) the surface bidirectional reflectance distribution function (BRDF).In rugged terrain, variations in surface elevation, slope, aspect, the obstruction coefficient, the sky view factor and the topographic configuration factor can significantly impactthese three factors, which produce atmospheric effects, illumination effects, and BRDF effects in optical high-resolution satellite images [5].Therefore, topographic normalization of satellite imagery has become a prerequisite for retrieving surface parameters, which can reduce or remove atmospheric, topographic and BRDF effects in rugged terrain simultaneously.
The quantitative topographic normalization of Landsat imagery has a relatively long history with the aid of a digital elevation model (DEM).The complexity and difficulty lie in the quantitative simultaneous estimation of the incident radiation, atmospheric transmittance, and surface BRDF for each image pixel in mountainous areas.Earlier studies removed atmospheric effects only by estimating atmospheric parameters, such as aerosols and water vapor, using a radiative transfer equation with a Lambertian surface, as reviewed in the literature [6,4].Masek et al. [7] derived the Landsat surface reflectance product by atmospheric correction based on the 6S methodology, which is the basis of the Landsat Climate Data Record (CDR) distributed by USGS.Some researchers recognize that the radiance received by the satellite is also distorted due to the surface BRDF or sun-surface-sensor geometry for each image pixel, and several atmospheric-correction-coupled BRDFs are utilized for Landsat images [8].In recent years, with the development of moderate-resolution satellites, many space-based sensors (such as the Moderate Resolution Imaging Spectroradiometer, MODIS) provide atmospheric products and multi-angle observations, which yield reliable data on land surface BRDFs [9,10].Some researchers use MODIS atmospheric and BRDF model parameters for Landsat correction [11,12].However, in mountainous areas, the topography strongly influences the radiance distribution recorded by optical sensors [13].For the same land surface characteristics, radiance values on a sunny slope are brighter than that those on a shady slope in images [1,5].Since the 1970s, two main categories of topographic correction methods, empirical and physical topographic corrections, have been developed to eliminate or at least reduce topographic influences [14,15].The first category mainly corrects the solar direct radiance to decrease illumination differences caused by terrain, such as ratio models, cosine models, Minnaert models, sun-canopy-sensor (SCS) models, and C models.Some studies compare several empirical topographic correction approaches for Landsat images [16][17][18][19].These empirical methods are extremely simple and, in particular cases, highly effective; however, they are not suitable for retrieving quantitative surface parameters [5].The second category employs a radiative transfer code to obtain deterministic descriptions of the correction of topographic effects and avoids image dependence [20].
Some investigators have developed complex terrain factors to describe solar illumination obtained by slope pixels, such as slope inclination, aspect, sky view factor and shadowing effects for radiation from the sky and terrain [21][22][23][24][25].Moreover, the spatial distribution of atmospheric transmissivity and optical depth vary with the rapidly changing elevation, aerosol optical depth (AOD) and precipitable water (PW); thus, topographic and atmospheric corrections must be performed simultaneously pixel by pixel [15,26,27].In addition, surface anisotropy, which is a function of the view and illumination angles, is evident.The models with Lambertian assumptions only consider the variation in the actual incident angles caused by the terrain and do not eliminate the influence of real sensor viewing angles [28].Some researchers have used an empirical BRDF to reduce the high reflectance values in regions of low illumination [13,27,29,30].In recent years, many quantitative topographic normalization schemes for removing atmospheric, topographic and BRDF effects on Landsat imagery have been proposed using MODIS BRDF parameters [1,2,5].However, Schaaf et al. [31] discovered that sloping and flat surfaces with similar vegetation have different BRDF shapes that are not symmetrical about the principal plane in rugged terrain.Therefore, to retrieve the surface reflectance from optical remote sensing, it is necessary to add a mountain BRDF model to the physical topographic normalization [32,33].
In this study, an improved topographic normalization algorithm for Landsat TM images is introduced by using the slope and aspect in the kernel-driven semi-empirical Ambrals BRDF model with the aid of the MODIS BRDF/albedo product and the MODIS AOD and PW products.The small variation in the BRDF effect due to variations in the sun and view angles (−7.5° to 7.5°) is ignored within a single Landsat scene.
The remainder of this study is organized as follows.In Section 2, an improved algorithm is proposed to correct the TM images.The study area and data are introduced in Section 3. The results and verification are described in Section 4. The conclusions of this study and a discussion are presented in Section 5.

Methods
To retrieve the surface spectral reflectance from Landsat imagery in rugged terrain, topographic normalization must be applied.This normalization simultaneously removes atmospheric, topographic, and surface BRDF effects via an atmospheric transmittance model, the estimation of solar spectral irradiance, and a mountainous BRDF normalized factor, as shown in Figure 1.

Improved Topographic Normalization
The TOA radiance (λ) where (λ,θ ) v T is the atmospheric transmittance from the observing direction, which is a function of the sensor viewing angle, wavelength λ and atmospheric conditions [15].The reflected radiance from the surface is determined by the surface reflectance in the space-based sensor observation direction and surface irradiance received from the sun.The global irradiance (λ) E can be decomposed into four components [34].Therefore, the reflected radiance from the surface can be written as follows: Where (λ denote the direct, isotropic diffuse, anisotropic diffuse and surrounding-reflected spectral irradiance, respectively.ρ (λ) T is defined as the directional-directional reflectance caused by direct and anisotropic diffuse irradiance at the slope surface.
is the hemispheric-directional reflectance caused by isotropic diffuse and surrounding-reflected irradiance [30].In a slope and aspect ( A S , ) pixel, the zenith and azimuth angles in the directions of the solar incident and sensor observations are ( θ ,  ), respectively.However, because of topographic effects, the actual illumination angles between the solar beam and normal in a slope pixel will be modified as ( is the relative azimuth between the sun and a satellite On an inclined surface, the local terrain can produce drastic variations in both the global irradiance and surface BRDF due to the changes in the associated angles, such as the actual zenith angle.Therefore, four actual angles on tiled surfaces must be changed in the new coordinate system based on the DEM, which can be calculated according to the literature [5].Considering the small variation in the Landsat sensor view angle, the observed zenith and azimuth are assumed to be zero for each pixel [8].The calculations for the actual observed zenith and azimuth angles on a given inclined surface can be simplified as: For the same land cover type, the inclined surface directional-directional reflectance ρ (λ)( , , φ )   , which eliminates the topographic effects in mountainous areas.Both reflectances follow the BRDF characteristics, and their relationship can be expressed as where (λ)( , , φ , θ ,θ , )  is defined as a BRDF normalized factor [1,30], which describes the transformation relation of the directional-directional reflectance transform caused by the topography.
In addition, because the hemispheric-directional reflectance is highly complex and the isotropic diffuse components and reflected irradiance are low on clear-sky days, the land surface can be reasonably assumed to obey isotropic characteristics.Therefore, In Equation (1), the path radiance is generated by Rayleigh and aerosol scattering, and ρ (λ) r and ρ (λ) a denote Rayleigh and aerosol reflectances, respectively.The equation can then be expressed as where 0 D is the calibration factor for the earth's orbit, which is a deterministic parameter of sun-earth geometry; and 0 (λ) is the extraterrestrial solar spectral irradiance at the time of the solar zenith angle s θ .A detailed description and discussion of Rayleigh and aerosol reflectances are provided in previous studies [15].Combining Equations ( 2), ( 5) and ( 6) with Equation (1), the corrected reflectance is flat and completely free of topographic effects, rewritten as When (λ) 1   , the surface is assumed to have Lambertian characteristics.Two cases are possible.
The first case indicates that only the atmospheric effect is considered in the reflectance inversion from the images.Equation ( 8) is modified to the following form [35]: The second case indicates that atmospheric and terrain effects influenced by the variation in the illumination angle are removed from the remote sensing images.In this case, the reflectance is estimated by [28]: However, the value of the BRDF normalized factor is often not equal to 1 for the surface of anisotropic reflection characteristics in rugged terrain.Three important input parameters for inverting the sloping reflectivity are described in Section 3.

Comparison of Four Correction Algorithms
To validate the superiority of the improved topographic normalization, four additional algorithms were analyzed, as presented in Table 1.FLAASH is a MODTRAN-based atmospheric correction software package that provides accurate surface reflectances by estimating the amount of atmospheric water vapor and aerosols from the image directly [6].The corresponding Landsat CDR [7] is the surface reflectance product distributed by USGS in the CDR.InEquation (9), the MODIS water vapor and aerosol products are input parameters for sole atmospheric correction, called the "MODIS-based atmospheric correction" (MBAC).Based on Equation (10), the method of topographic normalization by assuming a Lambertian surface is called the "MODIS-based Lambertian algorithm" (MBLA), and the atmospheric parameters are the same as those in MBAC.The topographic normalization proposed in this study is called the "MODIS-based BRDF algorithm" (MBBA), which is based on the MODIS BRDF/albedo kernel coefficient and MODIS atmospheric products.

Atmospheric Transmittance
Atmospheric attenuation is caused by the interaction between solar radiation and the principal atmospheric constituents, which are highly sensitive to aerosol and water vapor concentrations, followed by Rayleigh scattering and the scattering and absorption by the ozone column.More detailed documentation can be found in the literature [36,37].The atmospheric spectral transmittance formulas of the Landsat TM images follow Li's model parameterization schemes [15].Because of the relative stability over space and time, the correction of the atmospheric effects of ozone absorption and Rayleigh scattering is straightforward, as calculated by Li et al. [15].However, in rugged terrain, the spatial distribution of aerosols and water vapor is highly variable over space and time for large regions, and their values are questionable if obtained by a semi-empirical formula or empirically at a single point.With the rapid development of the American NASA Earth Observing System Program (EOS), many remote sensors have provided abundant atmospheric products, such as the MODIS AOD and PW products, which have become reliable data sources for atmospheric parameters.
In this algorithm, the MODIS water vapor product PW (cm) can be directly input as the w parameter.
The angstrom turbidity exponent is provided by the MODIS AOD products.However, the β value is not directly available for MODIS AOD products and is in need of further deduction.The derivation method can be found in the literature [34].The value of β can be expressed by where 0.55 τ is the AOD at a wavelength of 0.55µm, which is derived from the MODIS AOD products.
Moreover, the local surface pressure can be estimated from the surface elevation [34].

Estimation of the Solar Spectral Irradiance
In mountainous areas, the global spectral irradiance can typically be decomposed into three components, i.e., direct, diffuse, and surrounding-reflected irradiances (Figure 2).Considering the direction of the sun, the diffuse irradiance is further divided into isotropic and anisotropic diffuse irradiances.All components are functions of the wavelength  .The global irradiance can be expressed as follows [34]: (3) isotropic diffuse; and (4) surrounding-reflected.
Because of the variation in atmospheric conditions and terrain orientation, the spatial distribution of the irradiance components received on tiled surfaces presents strong spatial heterogeneity.Although the global irradiance and components are wavelength-dependent in Landsat images, formulas can follow the broadband estimation scheme described by Zhang et al. [34], except for the surrounding-reflected irradiance.The detailed calculation of the relevant terrain factors can be referenced in previous studies [26,15].In addition, the TOA solar spectral irradiance of Landsat5 TM is replaced according to the official website [38].
The surrounding-reflected solar radiation can be estimated by where ij F is the topographic configuration factor, i.e., the ratio of the energy reaching the visible pixel to the energy emitted from the object pixel [15] is the surface-reflected radiance measured via a satellite.

Bidirectional Reflectance Distribution Function (BRDF) Model
The BRDF normalized factor (λ)( , , φ , θ ,θ , )  is defined as the ratio of the tilted surface reflectance to the flat surface reflectance [1,30], eliminating the influence of terrain on homogenous land cover.Two types of reflectivity are all directional-directional reflectance, which obey BRDF laws.The BRDF considers surface reflectance properties at a particular wavelength to be a function of the angle of incidence and observation angles, which can correct directional anisotropy caused by shadows, buildings or topography [39].The MODIS BRDF/Albedo model parameter product (MCD43A1) was derived based on a 16-day reflectance dataset from MODIS on Terra and Aqua using the kernel-driven semi-empirical Ambrals BRDF model.This model can be used to calculate the directional reflectance at any desired viewing and solar geometries.The quality is assessed using the MCD43A2 QA (Quality Assessment) layer, which shows that the quality declines sequentially when the quality flags vary from 0 to 4. Mathematically, the Ross Thick-Li Sparse kernel-driven BRDF model is a linear combination of isotropic (defined as a constant value of 1) volume and geometric optical scattering kernels, which are given by Schaaf et al. [31].
where (θ , θ , ) (θ , θ , )   are kernel functions that depend only on the viewing and illumination geometries.The detailed algorithm of the two kernel functions can be found in Schaaf et al. [31].
(λ) iso f , vol f and vol f are the BRDF model empirical parameters provided directly by the MODIS BRDF product to describe the weights of the three kernel functions, i.e., the isotropic, volume and geometric optical scattering parameters, respectively.However, the spatial resolution of the MODIS BRDF parameter C005 is 500 m, which cannot be directly applied to the 30 m TM images.The parameters require downscaling to 30 m according to the land classification maps.According to the land cover types of the Landsat imagery, the BRDF parameters for each class at the MODIS scale are extracted.Based on land-type-specific BRDF correction factors, which hypothesize that the reflectances should be similar at the Landsat and MODIS spatial resolutions, the 30 m scale of the BRDF parameters for each pixel is extracted.Further details of the methodology can be found in the literature [40].However, the average effect of the BRDF parameters using a 500 m MODIS BRDF to correct 30 m Landsat imagery was not addressed here.Generally, the geometric relationship of the sun-surface-sensor is fixed on a flat surface at a given local time; thus, the two kernel functions are constant values.The MODIS BRDF/Albedo product does not explicitly consider slope and aspect effects [12].However, the two kernel functions will change for each pixel in rugged terrain due to the intensive variation in the actual illumination and sensor observation angles with the surface slope and aspect.In areas with significant slopes, the BRDF is no longer symmetric about the principal plane, and it must be corrected by adding the aspect and slope from a DEM [31].Therefore, considering the slope and aspect variations, the slope reflectivity can be calculated by ρ (λ)( , ,

Study Area and Data
The selected study area is located in the upper stream of the Heihe River Basin, a typical inland river basin in Northwest China [41], as shown in Figure 3.The area has hilly terrain and an altitude between 1714 m and 5076 m above sea level.There are five types of land cover in the area: farmland, forest dominated by Picea crassifolia, grass, water bodies and bare surface [42,43].Furthermore, the land cover is largely distributed according to terrain characteristics, particularly the terrain aspect caused by the solar radiation heterogeneity.The forest primarily covers the north-facing slope, farmland and grass cover the flat or gentle south-facing slopes, and bare land is widely distributed on the steep southern slope [34].The dataset for this study consists of a GDEM (30 m) derived from ASTER [44], MODIS PW and AOD atmospheric products (MOD05_L2 and MOD04_L2) with spatial resolutions of 1 km × 1 km and 10 km × 10 km (at nadir), respectively, the MODIS BRDF/Albedo product (MCD43A1) and a Landsat TM image (WRS path/row of 133/033, Level 1T product) acquired on August 11, 2009, at a solar zenith of 31.82° and a solar azimuth of 131.26°.In addition, a corresponding Landsat CDR [7] was applied for validating the atmospheric correction and topographic normalization.These data were obtained from the USGS Earth Resources Observation System (EROS) Data Center and were projected to Universal Transverse Mercator (UTM) coordinates.

Bidirectional Reflectance Distribution Function (BRDF) Normalized Factor
Using high-resolution DEM data, important topographic factors, such as the slope, aspect, obstruction coefficient, and sky view factor, are obtained for each pixel with sloping terrain.The flat surface directional reflectivity (LevelBrdf) and corresponding tilted pixel directional reflectivity (HillBrdf) can be derived using an Ambrals forward model based on the MODIS BRDF/Albedo parameter product.In this manner, the BRDF normalized factor (NormBrdf) is obtained; Figure 4 shows the spatial distributions in band 4.

Validation and Discussion
To achieve topographic normalization, the sloping reflectance derived from the sensor must be converted to a flat reflectance based on atmospheric and topographic correction, as expressed in Equation (8).There are two main methods for validating the merit of the topographic correction: visual inspection and quantitative accuracy analysis.Figure 5 shows the result of the partially enlarged images for a clear visual interpretation of the topographic normalization.The corresponding reflectance composite images of bands 4 (red), 3 (green), and 2 (blue) are also shown in Figure 5b.There were terrain effects in the original imagery of the rugged area, as seen by the dark shadows and bright regions facing the sun.The MBBA and MBLA can eliminate or at least reduce the terrain influence and flatten surfaces compared with the raw images (Figure 5a).Two topographic normalization methods can improve the pixel reflectance of self-shadowing and shadowing areas and reduce the contrast of surface features on shady and sunny slopes.After removing the terrain effects, more detailed features can emerge and improve the amount of information available in images.However, the MBLA overcorrects the reflectance in deep shady locations, such as areas A and B, by assuming isotropic surface reflectance.In contrast, the MBBA can reduce the overcorrection effects by considering the surface BRDF (Figure 5c).Moreover, both the MBLA and MBBA methods generate artifacts over bright areas within the image, possibly because the solar radiation was underestimated due to the AOD and PW overestimation of the MODIS atmospheric products [34].True surface measurements of these spectra are lacking.Moreover, in cold and arid study regions, terrain characteristics determine the spatial distribution of land cover types.It is very difficult to sample suitable areas with uniform vegetation on the sunny and shady sides of hills.Thus, in this test region, unsupervised or supervised classification on the original and corrected images is not feasible for accuracy evaluation.Therefore, other quantitative evaluation methods must be applied.
In general, testing whether the cosine of the solar illumination angle and the surface reflectance are linearly correlated has become a quantitative analysis method to verify the topographic normalization.Therefore, 481 pixels were randomly selected for the quantitative analysis in different terrain conditions and land cover types.Figure 6a shows that the reflectance of the sole atmospheric correction varies with the terrain and the cosine of the illumination angle in bands 2, 3, and 4, which are significantly positively correlated with the cosine of the illumination angle.However, the linear relationships between the pixel reflectance and cosine of the incident angle are eliminated or significantly weakened after considering the terrain effects in the same bands (Figure 6b).
In addition, the spectral curves of the four typical features are extracted to verify the advantage of this method insole atmospheric correction and topographic normalization.First, in flat areas, two types of surfaces are selected to illustrate the effectiveness of the atmospheric correction method, MBAC, MODIS atmospheric products as input parameters.As shown in Figure 7, the MBAC can effectively remove the influence of water vapor and aerosols and restore the true spectra of both water and farmland compared with FLAASH and Landsat CDR.Second, Figure 8a indicates that the spectrum of a shady forest pixel (the northwest-facing slope of 34) was enhanced overall in the MBAC compared with the MBLA and Landsat CDR.The spectrum of a relatively flat grass pixel (the southwest-facing slope of 15) was also restored.Moreover, the reflectance of shady forests was further increased after the terrain correction, and the degree of improvement was larger when considering the surface BRDF, particularly in band 4 (Figure 8a).High relief can lead to topography-related image distortion, whereas the topographic slope and aspect can influence the natural spectral variability.Based on the above comparative analysis, surface features distorted by atmospheric and topographic effects can be realistically reproduced using this algorithm.

Conclusions and Discussion
This study presents a method for topographic normalization with two novel features.First, the BRDF normalized factor is used to convert a sloping pixel reflectance into a flat pixel reflectance.Second, MODIS aerosol and water vapor products are directly used as input parameters to determine the atmospheric attenuation.
However, a careful analysis of Figure 5c shows that the topographic effects are not completely eliminated.In areas of broken terrain, the "topography" still exists on images corrected topographically, and overcorrection and even some noise is still present in low illumination areas.These overcorrection or residual topographic effects are related to several factors.First, the estimation of global solar irradiance is inaccurate in the low-illumination region, which is underestimated or overestimated.The direct solar radiation is relatively simple, but the estimation of diffuse and surrounding-reflected radiation has larger uncertainties.Second, the MODIS water vapor and aerosol products are considered sources of significant error because of the product inversion algorithms [34,45].Third, a DEM spatial resolution of 30 m is not sufficient for TM images for topographic normalization [28].Because the DEMs determine the pixel elevation and crucial terrain factors and ultimately affect the quality of the topographic normalization, the DEM grid size should be smaller than the scale of the original image.However, in addition to data acquisition difficulties, utilizing high-resolution DEMs to obtain the target pixel more accurately than topographic factors at the sub-pixel scale has been difficult; this is important to address in future research.Finally, the information gap caused by temporal and spatial resolutions between the MODIS products and the Landsat imagery is ignored in the current work, possibly leading to discrepancies in the scale of the BRDF parameters, atmospheric conditions, and illumination conditions.Furthermore, for imagery in different seasons, the solar zenith angle should be standardized to reduce the remaining BRDF effect caused by variations in the sun angle [1].

Figure 1 .
Figure 1.Flowchart of the improved algorithm.
can be obtained through radiometric calibration according to the digital number (DN) values, offsets and gains provided by the image header files.This radiation computation consists of two parts, the path radiance (λ)

3 (
spectral transmittance and component spectral transmittance caused by ozone and water vapor absorption, Rayleigh scattering and aerosol scattering.

Figure 3 .
Figure 3. Details of the study area based on a (a) DEM and (b) false-color composite (RGB: band 4, 3 and 2).

Table 1 .
Basic descriptions of the four atmospheric and topographic normalization methods.