PLC: A simple and semi-physical topographic correction method for vegetation canopies based on path length correction

Remote Sensing of Environment 215 (2018) 184–198

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

PLC: A simple and semi-physical topographic correction method for vegetation canopies based on path length correction ⁎

T

⁎

Gaofei Yina, , Ainong Lia, , Shengbiao Wub, Weiliang Fanc, Yelu Zengd, Kai Yane, Baodong Xub, Jing Lib, Qinhuo Liub a Research Center for Digital Mountain and Remote Sensing Application, Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu 610041, China b State Key Laboratory of Remote Sensing Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100101, China c School of Environmental and Resources Science, Zhejiang A & F University, Lin'an 311300, Zhejiang, China d Department of Global Ecology, Carnegie Institution for Science, Stanford, CA 94305, USA e School of Land Science and Techniques, China University of Geosciences, Beijing 100083, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Topographic effect Path length correction BRDF correction Radiative transfer

Rugged terrain distorts optical remote sensing signals, and land-cover classification and biophysical parameter retrieval over mountainous regions must account for topographic effects. Therefore, topographic correction is a prerequisite for many remote sensing applications. In this study, we proposed a semi-physically based and simple topographic correction method for vegetation canopies based on path length correction (PLC). The PLC method was derived from the solution to the classic radiative transfer equation, and the influence of terrain on the radiative transfer process within the canopy is explicitly considered, making PLC physically sound. The radiative transfer equation was simplified to make PLC mathematically simple. Near-nadir observations derived from a Landsat 8 Operational Land Imager (OLI) image covering a mountainous region and wide field-of-view observations derived from simulation using a canopy reflectance model were combined to test the PLC correction method. Multi-criteria were used to provide objective evaluation results. The performances were compared to that of five other methods: CC, SCS + C, and SE, which are empirical parameter-based methods, and SCS and DS, which are semi-physical methods without empirical parameter. All the six methods could significantly reduce the topographic effects. However, SCS showed obvious overcorrection for near-nadir observations. The correction results from D-S showed an obvious positive bias. For near-nadir observations, the performance of PLC was comparable to the well-validated parameter-based methods. For wide field-of-view observations, PLC obviously outperformed all other methods. Because of the physical soundness and mathematical simplicity, PLC provides an efficient approach to correct the terrain-induced canopy BRDF distortion and will facilitate the exploitation of multi-angular information for biophysical parameter retrieval over mountainous regions.

1. Introduction Optical remote sensing observation is subject to terrain-induced radiometric distortion (Proy et al., 1989; Schaaf et al., 1994). Generally, rugged terrain has two major effects on the observed radiance (Dymond et al., 2001; Li et al., 2012). First, topography modifies the illumination conditions of the surface through redistribution of the incoming irradiance (Kobayashi and Sanga-Ngoie, 2008; Schulmann et al., 2015; Wen et al., 2009). Illumination is the upper boundary condition of the radiative transfer process within the inclined canopy (Mousivand et al., 2015; Yin et al., 2017). Second, topography incurs variations in canopy structure (Fan et al., 2014; Luisa et al., 2008), and the free path length

⁎

Corresponding authors. E-mail addresses: [email protected] (G. Yin), [email protected] (A. Li).

https://doi.org/10.1016/j.rse.2018.06.009 Received 11 April 2018; Received in revised form 25 May 2018; Accepted 7 June 2018 0034-4257/ © 2018 Published by Elsevier Inc.

of the photon during the radiative transfer process within the canopy is distorted (Yin et al., 2017). Therefore, inclined and horizontal surfaces present different Bidirectional Reflectance Distribution Functions (BRDFs), exemplified by the asymmetry of the BRDF over an inclined surface (Schaaf et al., 1994). Topographic effects hinder the retrieval of biophysical variables (Gonsamo and Chen, 2014; Pasolli et al., 2015), land-cover classification (Hantson and Chuvieco, 2011) and applications that aim to detect land surface changes (Tan et al., 2013; VicenteSerrano et al., 2008). Therefore, topographic correction has become a prerequisite for remote sensing applications over mountainous areas. Corresponding to the two main topographic effects, topographic correction methods can be divided into illumination correction and BRDF

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

et al., 2008). Path length has been applied to characterize 3-D canopy structure for leaf area index (LAI) retrieval (Hu et al., 2014, 2018; Yan et al., 2016) and is a critical variable influencing the radiative transfer process within canopy (Hapke, 1981; Ross, 1981). The slope significantly distorts the angular distribution of the path length. Surface topography stretches the path length in the up-slope direction and compresses the path length in the down-slope direction (Duursma et al., 2003; Luisa et al., 2008). Path length distortion is the critical factor causing BRDF distortion. Yin et al. (2017) demonstrated the potential of path length correction (PLC) in reflectance simulation for sloping terrain. The objective of this paper is to develop a mathematically simple and physically sound topographic correction method based on PLC. The PLC method is by nature a BRDF correction, which can convert BRDF for an inclined surface to its counterpart over a horizontal surface. Near-nadir observations from a Landsat 8 Operational Land Imager (OLI) image covering a mountainous region and wide field-of-view observations derived from simulation using a canopy reflectance model were utilized to test the PLC correction method.

correction. Redistribution of the incoming irradiance by rugged terrain has been extensively examined. The most straightforward treatment for this topographic effect is to convert the direct irradiance to the inclined surface based on the cosine law, forming the basis for the cosine correction method (Teillet et al., 1982). The cosine correction method ignores the contribution of diffuse irradiance from the sky and the surrounding terrain and is often criticized for its overcorrection under low illumination conditions (Hantson and Chuvieco, 2011; Richter et al., 2009). To account for the diffuse irradiance, the CC correction (Teillet et al., 1982) method was proposed in somewhat empirical manners. More physically sound illumination correction methods accounting for radiative transfer in the atmosphere and land surface were also developed (Li et al., 2015; Li et al., 2012; Wen et al., 2015). These methods represent the influence of rugged terrain on sky obstruction and surrounding irradiance in detail yet mostly ignore the BRDF distortion caused by surface orientations. Illumination simulation and correction over rugged terrain are generally well-studied. However, the terrain-induced variation in canopy BRDF is still far from being well understood (Combal et al., 2000; Fan et al., 2014; Yin et al., 2017). One popular treatment is based on simple normalization of the sun-target-sensor (STS) geometry from the horizontal datum to the local terrain by coordinate rotation. This STS paradigm underlies many topographic correction methods. For example, Li et al. (2012) proposed a physics-based correction method that accounts for the illumination and BRDF variations, in which the BRDF variation is solved by straightforward coordinate rotation. This treatment implicitly assumes that vegetation grows perpendicular to the local slope, but trees over sloping terrain may still be erect because of the geotropic control (Gu and Gillespie, 1998; Soenen et al., 2005). Therefore, this underlying assumption of STS methods may cause uncertainty in sloping BRDF characterization (Fan et al., 2014; Yin et al., 2017). Sun-canopy-sensor (SCS) correction is a pioneering method in BRDF correction (Gu and Gillespie, 1998). It assumes that the sunlit canopy is the main factor contributing to the pixel-level reflectance and implements the topographic correction through normalizing the sunlit canopy area (Gu and Gillespie, 1998). One disadvantage of SCS correction is that it does not consider the viewing geometry effects; this may limit its usability for wide field-of-view sensors such as GF-1 WFV cameras with a view zenith angle from 0° to 24° (Feng et al., 2016). Dymond and Shepherd (1999) proposed a BRDF correction method, referred to as the D-S method in this paper, in which the sun and viewing geometries are both explicitly involved. However, the canopy structure is not accounted for in the D-S method, causing unsatisfactory topographical correction in some cases (Richter et al., 2009). The vegetation structure at the scale of the individual plant should be independent from the terrain because of the geotropic nature of vegetation (Wen et al., 2018). However, the topography significantly influences the topological relationship between different plants. This terrain-induced variation in canopy structure is the main determinant of the BRDF distortion over the inclined surface (Combal et al., 2000; Schaaf et al., 1994). This factor led to the development of canopy reflectance models for sloping terrain (Fan et al., 2014; Yin et al., 2017). Based on these sloping reflectance models, fully physical BRDF correction methods using model inversion have been developed. These fully physical BRDF correction methods have solid physical basis and outperform empirical methods (Couturier et al., 2013; Soenen et al., 2008; Yin et al., 2017). However, the ill-posed nature of the inverse problem counteracts their advantage in physics (Quan et al., 2015). A good method should be as physically sound as possible to ensure its transferability and as mathematically simple as possible to facilitate its operational application (Yin et al., 2015). A more advantageous semiphysical BRDF correction method is required to fulfill the above requirements. Path length is defined as the distance between the bottom and the top of the canopy along a direction relative to canopy height (Luisa

2. Method development 2.1. Theory background The proposed correction method is based on the classic radiative transfer theory that assumes the canopy is composed of homogeneously distributed leaves bounded by two parallel sloping planes: a sloping bottom soil surface and a sloping canopy top surface. The radiance (I) distribution function can be given by the radiative transfer equation (Ross, 1981):

dI (z , Ω2 ) 1 + μG (Ω2 ) I (z , Ω2) = ds π

∫4π

I (z , Ω1)Γ(Ω1 → Ω2 )dΩ1

(1)

where ds is the path length parallel to direction Ω2; z is the relative optical height inside the canopy layer along the direction of gravity, which, by convention, ranges from −1 at the bottom to 0 at the top of the canopy; μ is the leaf area volume density (m2/m3); G(Ω2) is the total leaf area projected on a plane perpendicular to the direction Ω2 by leaves of all orientations; Γ(Ω1 → Ω2) is the leaf scattering phase function quantifying the fraction of the intercepted energy from Ω1 that is scattered into a unit solid angle around Ω2. G(Ω2) and Γ(Ω1 → Ω2) can be parameterized by the canopy structure and leaf optical properties. For their detailed expressions, (Knyazikhin et al., 1992; Myneni et al., 1991) can be referenced. 2.2. PLC formulation Derivation of the PLC method consists of two steps: First, the radiative transfer equation (Eq. (1)) was simplified based on some assumptions, and the canopy reflectance was expressed as an explicit function of the path length. Based on the simplified radiative transfer function, a BRDF correction formula was deduced. Second, the path lengths over horizontal and sloping terrains were formulized. We obtained the corrected reflectance after substituting these path lengths in the BRDF correction formula derived in the first step. The main symbols used in formulation of the PLC are summarized in Table 1. 2.2.1. Radiative transfer equation simplification We assume the radiative transfer process within a sloping canopy meets the following assumptions: Assumption I. The canopy is illuminated by collimated light, i.e., the sky diffuse and surrounding reflected radiation is negligible. Assumption II. The radiance collected by the sensor is only from the single-scattering of leaves, i.e., the contributions from the soil reflectance and the multiple-scattering of leaves are negligible. 185

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

directions over flat terrain, respectively, and St(Ω1) and St(Ω2) are their counterparts over sloping terrain. The underlying physics of Eq. (7) is that the BRDF distortion over mountainous areas is mainly caused by the terrain-induced path length variation. Therefore, the reflectance from sloping terrain can be converted to its horizontal surface counterpart through path length correction (PLC).

Table 1 Main symbols used in the derivation of PLC. Symbol

Explanation

Ω1 Ω2 θ1 φ1 θ2 φ2 S St α β i

Solar direction Viewing direction Solar zenith angle Solar azimuth angle Viewing zenith angle Viewing azimuth angle Path length over flat terrain Path length over sloping terrain Slope angle Aspect angle of the slope Incident zenith angle between the sun and the surface normal, which can be derived from Eq. (10) Observation zenith angle between the sensor and the surface normal, which can be derived from Eq. (10)

o

2.2.2. Path length formula A key component in implementing PLC correction is calculating the path length over horizontal and sloping terrains. For a canopy over flat terrain, the path length is characterized by azimuthal symmetry and can be simply computed as a function of the solar/viewing zenith angle (θ) (Yin et al., 2017),

For sloping terrain, the path length is dependent on the solar/ viewing direction (with zenith and azimuth angles of θ and φ, respectively) and on the topographic features (with a slope and aspect of α and β, respectively). Two typical path length expressions in the literature can be potentially used; one is (Gonsamo et al., 2011; Walter and Torquebiau, 2000):

The radiance collected by the sensor can be formulated as (Hapke, 1981; Liangrocapart and Petrou, 2002; Yin et al., 2017):

I (0, Ω2) =

0

∫−1

F (z , Ω2) exp[−L (z ) G (Ω2 ) S (Ω2)]/cos(θ2 ) dz

(2)

where L(z) is the accumulated LAI from the top of the canopy to a height z, S(Ω2) is the path length along the viewing direction Ω2, θ2 is the viewing zenith angle, and F(z, Ω2) is the source function (Liangrocapart and Petrou, 2002), and

F (z , Ω2) =

J Γ(Ω1 → Ω2) exp(−L (z ) G (Ω1) S (Ω1)) π

St (θ, φ, α, β ) = 1/cos θ′

cos θ′ = cos(θ) cos(α ) + sin(θ) sin(α ) cos(φ − β ) (3)

St (θ, φ, α, β ) =

1 cos θ (1 − tan α cos(φ − β ) tan θ)

(11)

The key difference between Eqs. (9) and (11) lies in the definition of the reference direction. The path length calculated from Eq. (9) was normalized according to the normal of the sloping surface, and the minimum value is 1.0, when the solar/viewing direction is perpendicular to the slope. Eq. (11) calculates the path length over a horizontal surface and the unity value appears in the vertical direction. Because the radiative transfer equation in our study (Eq. (1)) is defined over a horizontal reference plane, we used Eq. (11) to ensure physical consistency between the path length formulation and radiative transfer simplification.

(4)

In the case of dense vegetation, 1 - exp[-LG(Ω1)S(Ω1) - LG(Ω2)S(Ω2)] in the numerator is approximated to 1.0. Therefore, Eq. (4) can be further reduced to:

Γ(Ω1 → Ω2) 1 cos(θ1) cos(θ2 ) G (Ω1) S (Ω1) + G (Ω2 ) S (Ω2 )

(10)

The other is (Duursma et al., 2003; Luisa et al., 2008):

ρ=

ρ=

(9)

where θ′ is the zenith angle relative to the local surface normal, calculated as,

where J is the incoming irradiance from the sun. According to the definition of the Bidirectional Reflectance Factor (BRF) (Schaepman-Strub et al., 2006), the canopy reflectance can be calculated as:

I (0, Ω2 ) π J cos(θ1 ) Γ(Ω1 → Ω2) {1 − exp[−LG (Ω1 ) S (Ω1) − LG (Ω2 ) S (Ω2 )]} = cos(θ1) cos(θ2 ) G (Ω1) S (Ω1) + G (Ω2 ) S (Ω2 )

(8)

S (θ) = 1/cos θ

(5)

3. Method evaluation

We further assumed that, The proposed PLC correction method was evaluated in two main ways: First, the PLC method was implemented with a Landsat 8 OLI image to examine its performance for near-nadir viewing observations. Second, a canopy reflectance model dedicated for sloping terrain was employed to simulate the reflectance over sloping and horizontal surfaces. The horizontal reflectance served as the reference value of the corresponding sloping reflectance with identical canopy structure and observation geometry. The second evaluation method focused on assessing the potential of PLC in topographic correction for wide field-ofview observations.

Assumption III. The leaf inclination distribution function is spherical, thus, G ≡ 0.5. The canopy reflectance can be further formulized as:

ρ=

2Γ(Ω1 → Ω2) 1 cos(θ1) cos(θ2 ) S (Ω1) + S (Ω2 )

(6)

Because of the geotropic nature of plant growth, the first term of the right hand equation in Eq. (6), i.e., 2Γ(Ω1 → Ω2) / cos(θ1)cos(θ2), has the same form for both horizontal and sloping surfaces (Wang et al., 2002; Yin et al., 2017). Therefore, in Eq. (6), only the term of 1 / (S (Ω1) + S(Ω2)) is related to the topography. The canopy reflectance defined over the horizontal reference plane (ρt, the image observed reflectance) and its counterpart defined over the sloping reference plane (ρPLC, the targeted reflectance) can both be formulated as Eq. (6). We obtained the following formula for topographic correction:

ρPLC = ρt

S (Ω1) + S (Ω2 ) St (Ω1) + St (Ω2 )

3.1. Topographic correction for Landsat 8 OLI image 3.1.1. Study area A mountainous region at the junction of Sichuan, Shanxi and Gansu Provinces in China was selected as the study area. The area is centered at approximately 32° 96′ N and 105° 56′ E (Fig. 1(a)) with altitudes from 500 m to 2500 m (Fig. 1(b)). The following factors make our study area a satisfactory candidate for evaluation of reflectance correction: First, the shadows caused by sun obstruction by the surrounding terrain are limited, so we can focus on the topographic effects caused by

(7)

where S(Ω1) and S(Ω2) are the path lengths along solar and viewing 186

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 1. Study area. (a) Landsat8-OLI image from July 15, 2015 (the R, G and B color space correspond to bands 5, 4 and 3, respectively). (b) Elevation map of the study area. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(1) Visual inspection: The visual inspection checked if the reflectance difference between the slopes facing toward and away from the sun was reduced. (2) Correlation analysis: Correlation analysis between the cosine of the local incident angle (cos(i), see Table 1) and the reflectance is one of the most widely used quantitative evaluation methods. The efficiency of the correction methods can be quantified by the R2 and the slope of their linear regression. The ideal correction method should result in R2 and slope values of zero. (Yin et al., 2017). (3) Aspect-dependence of reflectance: The topography causes the reflectance on the sunny aspect to appear higher than that on the shady aspect. Theoretically, a good correction method should eliminate this aspect-dependence of reflectance (Balthazar et al., 2012). (4) Radiometric stability: Theoretically, the maximum/minimum reflectance in the original image before correction should appear in the sunny/shady slope and will be reduced/increased after topographic correction. Consequently, the reflectance range will be reduced by a successful correction method. Moreover, the median reflectance is relatively stable and should not change significantly after correction (Sola et al., 2016). (5) Physical soundness: Assuming the surrounding-reflected irradiance is negligible, the slope will exert no influence on the reflectance when the solar and viewing directions are perpendicular to the aspect. In this special case, the canopy structure over flat terrain is the same as that over sloping terrain (Yin et al., 2017). Therefore, the reflectance should not change after topographical correction over such illumination-observation geometries. Landsat 8 OLI is a near-nadir viewing sensor, and its viewing direction can be assumed to be perpendicular to all the aspects. Therefore, we analyzed the reflectance over aspects of approximately 23.7° and 203.7°, which are perpendicular to the solar direction. A physically sound correction method should not change the reflectance over these aspects.

surface orientation per-se rather than those caused by the surrounding terrain, which are more difficult to correct. Second, the land-cover types are relatively limited and are dominated by forests. Therefore, the variation in reflectance incurred by topography may exceed that caused by cover heterogeneity. 3.1.2. Data A cloud-free Landsat 8 Operational Land Imager (OLI) image (path 129 and row 37) acquired on July 15, 2017 was used to provide the canopy reflectance over the study area. The solar zenith and azimuth angles were 23.5° and 113.7°, respectively. The image was downloaded through the Earth Resources Observation and Science Center Science Processing Architecture (ESPA) on-demand interface (USGS, 2017b), accessible at https://espa.cr.usgs.gov/. The ESPA provides canopy reflectance through an atmospheric correction algorithm based on the Second Simulation of the Satellite Signal in the Solar Spectrum Vectorial (6SV) model (Vermote et al., 2016). The 6SV model has been integrated in the Landsat 8 Surface Reflectance Code (LaSRC) used to generate the high-level Landsat data products (USGS, 2017a). In addition to the canopy reflectance, the LaSRC also delivers pixel-wise illumination-observation products (including solar zenith angle, solar azimuth angle, viewing zenith angle and viewing azimuth angle), which are important inputs in our correction method to calculate the path length (see Eqs. (8) and (11)). In this study, the red and near infrared (NIR) bands, which are the most frequently used bands in biophysical parameter retrieval, were used to assess the efficiency of the topographic correction methods. The topographic factors used in the topographic correction (i.e., slope and aspect) were derived from the Advanced Spaceborne Thermal Emission and Reflectance Radiometer (ASTER) Global Digital Elevation Model (GDEM) version 2 product (Tachikawa et al., 2011). ASTER GDEM is a product of the Ministry of Economy, Trade, and Industry (METI) in Japan and the National Aeronautics and Space Administration (NASA) in the United States. The spatial resolution is 1 arc-second (~30 m), approximately the same as the Landsat8-OLI image. The overall accuracy of the GDEM version 2 product was reported to be approximately 17 m at the 95% confidence level (METI and NASA, 2011).

3.2. Evaluation using radiative transfer model simulation Canopy reflectance is anisotropic, so the successful mitigation of topographic effects may not lead to similar reflectance for canopies with identical structures but different observation geometries. This scenario makes it more difficult to evaluate the performance of topographic correction methods for wide field-of-view observations than for near-nadir viewing observations.

3.1.3. Evaluation strategies Many strategies were used to assess the performance of the topographic correction methods (Sola et al., 2014, 2016). To obtain an objective evaluation result, five different methods were selected: 187

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

A canopy reflectance model dedicated for sloping terrain is an alternative to evaluate PLC for wide field-of-view observations. Sloping reflectance and horizontal reflectance with identical structure and observation geometry can be simultaneously simulated, and the latter reflectance can serve as a reference value for the sloping reflectance after correction. Therefore, direct comparison between the corrected and reference reflectance is possible using canopy reflectance model simulation.

Table 3 Expressions of topographic correction methods for comparison. A and B are the slope and intercept of the fitted line between cosi and ρt, C is the ratio between A and B, and ρt is the mean of the original reflectance. For the meanings of other parameters, see Table 1. Correction method

Expression

Reference

CC

cos α + C cos i + C cos α cos θ1 ρSCS = ρt cos i cos α cos θ1 + C ρSCS + C = ρt cos i + C

Teillet et al. (1982)

ρSE = ρt − (A cos i + B ) + ρt

Teillet et al. (1982)

SCS

3.2.1. Simulation of sloping and horizontal reflectance The canopy reflectance model developed by Yin et al. (2017) was selected to simulate concomitant sloping and horizontal reflectance. The effects of sloping terrain on single-order and diffuse scatterings are accounted for by path length correction and modification of the fraction of incoming diffuse irradiance, respectively. The performance of this model for sloping surface was comprehensively validated using Monte Carlo and image simulations (Yin et al., 2017). The PLC validation dataset was obtained through model simulation with the following input parameters: The hot spot size parameter and fraction of diffuse irradiance to the total irradiance were set as 0.05 and 0, respectively; the solar zenith and azimuth angles were specified as the same as those of the Landsat 8 OLI image (i.e., 23.5° and 113.7°, respectively); three types of leaf inclination distribution functions were used, spherical, planophile, and erectophile; the observation zenith and azimuth angles were specified in the range of [0°, 40°] and [0°, 330°], respectively, with steps of 5° and 30°, respectively, resulting in 193 observation geometries; the LAI ranged from 1 to 5 with a step of 1; and the slope and aspect were specified in the range of [5°, 40°] and [0°, 330°], respectively, with steps of 5° and 30°, respectively, resulting in 192 sloping surfaces. Detailed descriptions of the above parameters can be found in the literature (Yin et al., 2017). The Red and NIR bands were selected; the optical properties of the leaf and soil components are listed in Table 2, as recommended by the RAMI exercise (Pinty et al., 2001), and 2895 structure-observation cases of horizontal reflectance (3 leaf inclination distribution functions × 193 observation geometries × 5 LAI values) were simulated with the slope and aspect set to 0°, and each horizontal reflectance corresponds to 192 sloping reflectance cases.

SCS + C SE D-S

NIR band

Leaf reflectance Leaf transmittance Soil reflectance

0.0546 0.0149 0.1270

0.4957 0.4409 0.1590

Soenen et al. (2005)

Dymond and Shepherd (1999)

3.3. Comparison of PLC and other topographic correction methods Five commonly used topographic correction methods were selected for comparison. Table 3 shows their expressions. CC, SCS + C and SE were selected because they perform excellently in most cases (Hantson and Chuvieco, 2011; Sola et al., 2016). However, they all rely on empirical parameters. Therefore, SCS and D-S were also compared. Similar to our PLC method, SCS and D-S were both derived according to the physical mechanism of remote sensing imaging and no empirical parameter was involved. 4. Results 4.1. Topographic correction for near-nadir observations from Landsat 8 OLI image 4.1.1. Visual inspection Fig. 2 shows a visual comparison of the uncorrected image and corrected images using different methods. The images in Fig. 2 were all rendered by applying a linear method with a 2% clip on both ends of the displayed data. There were obvious terrain effects in the uncorrected image, as indicated by the bright regions and dark shadows facing toward and away from the sun, respectively. The reflectance variations caused by terrain were significantly reduced in the corrected images. However, the SCS method generated some artifacts. Many black clusters appeared in the SCS-corrected image, and the image became blurred. Two different types of hues exist in the corrected images: SCS and D-S generated darker results than the original uncorrected image, whereas the CC, SCS + C, SE and PLC corrected images retained a similar hue as the original image. A detailed visual comparison of the enlarged area (white box in Fig. 2) is shown in Fig. 3. The SCS and D-S corrected images represent somewhat overcorrection. For example, the shady slopes in the northwestern part of this sub-area became too bright after SCS correction. This excessive brightness of the same shady slopes can also be found in the D-S corrected image with less magnitude. In addition, the ridges became very dark after CC, SE and D-S correction, causing the corresponding corrected images to appear extremely fragmented. Generally,

Table 2 Optical properties of the leaf and soil components. Red band

ρD − S =

cos θ1 + cos θ2 ρt cos i + cos o

Gu and Gillespie (1998)

In addition, we analyzed the influence of the observation angle on PLC performance as follows: We simulated the sloping reflectance before and after correction and the corresponding horizontal reflectance in the solar principal plane (with observation azimuth angles of 113.7° and 293.7°). During the simulation, the LAI was set at a moderate value of 3, and a spherical leaf inclination distribution function was used. Other input values were specified according to the PLC validation dataset described in Section 3.2.1. The reflectance variation for a specified observation geometry is exclusively caused by surface orientation. The reflectance variations in different observation zeniths were analyzed to reveal the influence of observation angle.

3.2.2. Evaluation strategies The radiative transfer simulation enables a direct comparison between topographically corrected and reference reflectance. In addition to visual inspection, two criteria were adopted to quantitatively assess the correction results. The first criterion is the relative root mean square error (R-RMSE) between the topographically corrected and reference reflectance. The second criterion is the coefficient of variation (MCV) of the corrected reflectance, averaged over the 2895 structure-observation cases. The R-RMSE was used to quantify whether the correction method can reconstruct horizontal reflectance; a low R-RMSE indicates a high absolute accuracy. The MCV was designed to assess if the correction method can reduce the topographic effect, irrespective of its absolute accuracy. The reflectance variation in each structure-observation case is totally caused by surface orientation, so a satisfactory correction method will have a low MCV. However, a correction method with low MCV does not necessarily have a high absolute accuracy because it may systematically overestimate or underestimate the corresponding horizontal reflectance.

Variables

ρCC = ρt

188

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 2. False color composite images before and after topographic correction. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

amplified (Li et al., 2015; Li et al., 2012). After correction using the SCS method, the R2 decreased, indicating that the topographic effect was weakened. However, the slopes became negative (see Table 4), i.e., SCS overcorrected the original image. The finding of overcorrection problem of the SCS method is consistent with other studies (Kane et al., 2008; Soenen et al., 2005). The D-S method obtained satisfactory results but exhibited a slight overcorrection in the NIR bands with a slope of −0.007 (see Table 4). In addition, the scatterplots for the D-S method show a shift to the high reflectance values compared to other correction methods, i.e., the D-S method may cause a systematic bias in the reflectance. Consistent with existing studies (Hantson and Chuvieco, 2011; Sola

the SCS + C and PLC methods generated the superior corrected results without artifacts, e.g., overcorrection in the SCS and D-S results and fragmentation in the CC, SE and D-S results. 4.1.2. Correlation analysis The correlation between the red band reflectance and the cosine of the local solar incidence angle (cos(i)) is shown in Fig. 4, and the detailed regression results for the red and NIR bands are shown in Table 4. Due to the topographic effect, significant positive correlations were observed for the uncorrected image, especially for the near infrared band (R2 = 0.384, slope = 0.249, see Table 4), because the sky diffuse radiation in this band is very small and the topographic effect is

Fig. 3. Enlargement of the inset shown in Fig. 2. 189

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 4. Density scatterplots between the red band reflectance and the cosine of the local solar incidence angle (cos(i)) before and after topographic correction.

empirical methods in the red band, but PLC's performance (4.1%) surpassed that of SE (4.5%) in the NIR band.

Table 4 The regression results between the reflectance (y) and the cosine of the local solar incidence angle (x). Lower values of the slope and R2 of the fitted line indicate a better corrected result. RF represents regression function. Correction method

Original image CC SCS SCS + C SE D-S PLC

Red band

4.1.4. Radiometric stability To analyze the radiometric stability, we generated box plots of the red band reflectance from the uncorrected and the corrected images (Fig. 6). Theoretically, the reflectance ranges after correction should be contained in their counterparts before correction (Sola et al., 2016). The SCS + C, SE and PLC corrected reflectance met this theoretical distribution. The upper bound of the CC results was slightly higher than that of the original case. The upper bound of the SCS results was also higher than that of the original case, with a greater magnitude due to its overcorrection problem. The D-S method results showed an obvious shift to high values, revealing the systematically positive bias of the D-S corrected reflectance. We further list the median in the red and NIR bands in Table 6. Generally, CC, SCS and D-S overestimated the original median. SCS + C, SE and PLC obtained radiometrically stable results and could reproduce the original median reflectance very well.

NIR band

RF

R2

RF

R2

y = 0.019x + 0.011

0.089

y = 0.249x + 0.154

0.384

y = 0.027 y = −0.010x + 0.034 y = 0.005x + 0.022 y = 0.005x + 0.021 y = 0.030 y = 0.007x + 0.020

0.0002 0.024 0.006 0.007 0.0002 0.015

y = 0.008x + 0.372 y = −0.146x + 0.476 y = 0.057x + 0.307 y = 0.058x + 0.300 y = −0.007x + 0.413 y = 0.086x + 0.278

0.0005 0.173 0.036 0.033 0.0003 0.078

et al., 2016), CC, SCS + C and SE reduced the correlation between the reflectance and cos(i) significantly. This result can be expected because they all utilize the regression between the reflectance and cos(i) (see Table 3). The PLC correction obtained a similar level of de-correlation. This result is significant considering that there is no empirical parameter involved in its expression (Eq. (7)).

4.1.5. Physical soundness The density scatterplots between the topographically corrected and original reflectance in the red band over aspects perpendicular to the solar direction are represented in Fig. 7. The statistics including the R2, relative RMSE and relative bias in red and NIR bands are shown in Table 7. The points in the CC derived scatterplot were clustered around the 1:1 line but the fitted line showed a slight shift to high values. The scatterplot for the SCS method fitted the 1:1 line very well without bias. However, SCS represented an obvious “bivalve” structure, indicating that overestimation and underestimation might both appear. The scatterplot for SCS + C inherited the “bivalve” structure from SCS. The SCS + C corrected reflectance exhibited more consistency with the original ones in this aspect plane, due to the introduction of an empirical parameter. The D-S method showed a more dispersed distribution. In addition, most of the points were above the 1:1 line, i.e., the reflectance was systematically overestimated by the D-S method. The SE reproduced the original reflectance satisfactorily with very slightly underestimation. Superior to other methods, the PLC corrected reflectance remained nearly the same as the original reflectance for aspects of approximately 23.7° and 203.7°. Therefore, we can conclude that PLC is physically sound.

4.1.3. Aspect-dependence of reflectance Fig. 5 reveals the reflectance variations across different aspect angles. In the original image, the reflectance over the sun-facing aspects (approximately 113.7°, the red point in each plot) was remarkably larger than those over the opposite directions. Compared with the original image, SCS generated an opposite trend due to overcorrection. After the topographic correction of other methods, the reflectance variations over different aspects were significantly reduced. The relationship between the reflectance and aspect expressed in polar coordinates exhibits circles centered at ordinate origins. However, a closer analysis of Fig. 5 revealed that the area of the D-S derived circle was obviously larger than those of other correction methods. In theory, an ideal topographic correction should not significantly alter the area of the circle because the radius in the sunny and shady aspects will be decreased and increased, respectively. Therefore, the reflectance after D-S correction was overestimated. To quantify the aspect distribution of the reflectance, Table 5 lists the coefficient of variation (CV) (%) of the reflectance across different aspects. CC performed the best in both the red (2.6%) and NIR (3.2%) bands. The PLC resulted in a slightly higher CV (3.5%) than the 190

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 5. Mean reflectance across aspect angles. The radius represents the magnitude of the reflectance with zero in the center, and the polar angle represents the aspect. The red point in each plot represents the solar azimuth. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 5 The coefficient of variation (%) of reflectance across different aspects.

Table 6 The median of reflectance over the study area.

Correction method

Red band

NIR band

Correction method

Red band

NIR band

Original image CC SCS SCS + C SE D-S PLC

10.2 2.6 7.5 2.7 3.4 2.8 3.5

11.5 3.2 6.4 3.3 4.5 4.0 4.1

Original image CC SCS SCS + C SE D-S PLC

0.0259 0.0277 0.0264 0.0256 0.0250 0.0301 0.0251

0.3526 0.3819 0.3602 0.3541 0.3475 0.4105 0.3474

4.2. Evaluation for wide field-of-view observations through radiative transfer model simulation 4.2.1. Comparison with simulated horizontal reflectance The comparison between the topographically corrected red band reflectance and the corresponding horizontal reflectance simulated by the canopy reflectance model is exhibited in Fig. 8, and the statistics quantifying the performance of each correction method are shown in Table 8. Before correction, each horizontal reflectance corresponded to many sloping reflectance caused by different surface orientations, resulting in a significantly dispersive scatterplot (Fig. 8(a)). Empirical correction methods (CC, SCS + C and SE) generated similar results: The MCV quantifying the reflectance variations caused by topography decreased from 11.5% before correction to approximately 6.0% in the red band, i.e., approximately half of the topographic effect was removed. Similar performances were obtained in the NIR band for the three methods. In addition, the empirical methods can reconstruct the absolute horizontal reflectance, and the R-RMSE also decreased by nearly half. Of the three selected empirical methods, SCS + C performed the best. This may be because SCS + C combines the physical mechanism of SCS and the empirical parameter of C (see Table 3). SCS also mitigated the topographic effect and reconstructed the absolute horizontal reflectance with a relatively lower magnitude than empirical methods. Regarding MCV in Table 8, the performance of D-S in reducing the topographic effect was comparable with that of SCS + C. However, there was an obvious positive bias in the D-S corrected reflectance,

Fig. 6. Box plots of the reflectance from the uncorrected (UNCORR) and corrected images in the red band. The boxes stretch from the 25th percentile to the 75th percentile. The horizontal solid line within each box represents the median. The bars correspond to ± 1.5σ. Outliers are not displayed. The horizontal dashed lines represent the median and median ± 1.5σ of the original reflectance.

191

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 7. Density scatterplots between the topographically corrected and original reflectance over aspects perpendicular to the solar direction.

rapidly than that of PLC with increasing observation zenith. When the observation zenith was larger than 10°, the CV of PLC was much smaller than that of SCS + C. Therefore, PLC is superior for slant observations.

causing the highest R-RMSE. Among the six methods, PLC performed the best in reducing the topographical effect and reconstructing the absolute horizontal reflectance. Two factors may explain its good performance for wide fieldof-view observations: First, the influence of topography on the radiative transfer process within a sloping canopy is accounted for through path length correction. Second, the observation geometry is explicitly presented (see Eq. (7)).

5. Discussion 5.1. Comparison between PLC and existing correction methods In addition to PLC, five other methods were selected for comparison. CC, SCS + C and SE were selected because they rely on empirical parameters and perform excellently in most cases (Hantson and Chuvieco, 2011; Sola et al., 2016). SCS and D-S were chosen as representatives for semi-physical methods without empirical parameter, similar to the proposed PLC method. The topographic correction for Landsat 8 OLI image found that SCS overcorrected the topographic effect for near-nadir observations. D-S successfully reduced the topographic effect but yielded a systematic positive bias. Therefore, D-S is more suitable for land-cover classification which focus more on relative accuracy, rather than biophysical parameter retrieval which requires high absolute reflectance accuracy. Empirical parameter-based methods such as CC, SCS + C and SE in this study performed the best. Without any empirical parameter, PLC obtained comparable results to empirical parameter-based methods and was more physically sound. For wide field-of-view observations, PLC outperformed all the other methods (see Fig. 8 and Table 8). An additional analysis using a canopy reflectance model simulation showed

4.2.2. Sensitivity to the observation angle The two best performing methods in topographic correction for wide field-of-view observations, i.e., SCS + C and PLC, were selected for the sensitivity analysis to the observation angle. Fig. 9 shows the changes in the red band reflectance with observation zenith angle in the solar principal plane. Because of different surface orientations, each observation angle has 192 sloping and corrected reflectance. For comparison, the horizontal reflectance at each observation zenith is also shown in Fig. 9. The reflectance variations caused by topography were reduced after both SCS + C and PLC correction. Compared to the sloping reflectance, the corrected reflectance aggregated more closely around the horizontal reflectance. PLC performed better than SCS + C, especially for slant observations. We also analyzed the CV of the reflectance for each observation zenith (see Fig. 10). For near-nadir observations (with viewing zenith angle ranging from −10° to 10°), SCS + C obtained better results than PLC. However, the CV of SCS + C increased more

Table 7 The regression results between the corrected (y) and original reflectance (x) over aspects perpendicular to the solar direction. RF represents regression function. Correction method

CC SCS SCS + C SE D-S PLC

Red band

NIR band 2

RF

R

y = 1.031x + 0.001 y = 0.996x y = 0.995x y = 0.985x y = 1.072x + 0.002 y=x

0.957 0.983 0.984 0.981 0.915 1.000

R-RMSE

R-BIAS

RF

R2

R-RMSE

R-BIAS

0.106 0.047 0.044 0.062 0.187 0.006

0.073 0.000 0.000 −0.030 0.150 0.000

y = 1.010x + 0.0215 y = 0.947x + 0.019 y = 0.950x + 0.018 y = 0.963x + 0.005 y = 1.128x + 0.008 y = 0.997x + 0.001

0.773 0.880 0.897 0.866 0.639 0.998

0.103 0.047 0.043 0.058 0.188 0.006

0.071 0.000 0.000 −0.022 0.150 0.000

192

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 8. (a) Scatterplot between canopy reflectance model simulated horizontal and sloping reflectance. (b)–(c) Scatterplots between horizontal and corrected reflectance. For brevity, only the red band is shown.

implements the topographic correction through normalizing the sunlit canopy area (Gu and Gillespie, 1998). The volume scattering of the leaves within the canopy is neglected, which is a key drawback of the geometric optical model (Chen et al., 2000). Although the D-S method was derived from the fully validated Hapke model (Dymond and Shepherd, 1999; Hapke, 1981), it implemented the derivation over the local coordinate system, and resulted in a relative rough assumption, i.e., reflectance is inversely proportional to the sum of the cosine of the incidence angle and the cosine of the exitance angle (Dymond and Shepherd, 1999). Unlike SCS and D-S, the influence of terrain on the radiative transfer process within the canopy is fully accounted for in PLC. The solid physical basis of PLC avoids the appearance of unreasonable correction results, e.g., the overcorrection in SCS and the radiometric instability in D-S. In contrast to PLC, which is derived from radiative transfer equation simplification, some researchers recommended using canopy reflectance models directly to implement topographic correction. For example, based on the Li-Strahler Geometric Optical Mutual Shadowing model (Li and Strahler, 1992; Schaaf et al., 1994) and Discrete Anisotropic Radiative Transfer model (Gastellu-Etchegorry et al., 2004), Soenen et al. (2008) and Couturier et al. (2013) developed the MFMTOPO and MFM-3-D correction methods, respectively. In MFM-TOPO and MFM-3-D, the canopy reflectance models were used in an explicit manner: They firstly retrieved the biophysical parameters using the canopy reflectance models (in inverse mode) using remote sensing image and DEM data. Then, they simulated the corrected reflectance

Table 8 The R-RMSE (%) and MCV (%) of the corrected reflectance. See section 3.2.2 in the main text for the definition of MCV. Correction method

UNCORR CC SCS SCS + C SE D-S PLC

Red band

NIR band

R-RMSE

MCV

R-RMSE

MCV

12.0 6.9 7.6 5.7 6.7 8.1 3.3

11.5 5.6 6.6 5.0 6.4 4.9 2.5

8.9 4.9 7.2 4.1 4.7 7.8 2.9

8.5 4.0 7.3 3.6 4.2 4.8 2.2

that although empirical parameter-based methods obtained slightly better results than PLC for near-nadir observation, PLC was extremely superior to parameter-based methods for slant observations. Therefore, we recommend using PLC to implement topographic correction for multi-angular observations collected from wide field-of-view sensors. One of the obvious advantages of PLC is its physical soundness, and no empirical parameter is involved. In contrary to the SCS and D-S methods, which are also without empirical parameter, the PLC method was derived from the classic radiative transfer theory. The SCS method is based on a geometric optical approach. It assumes that the topographic effect is from the regulation of the sunlit canopy area and

193

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 9. SCS + C and PLC corrected reflectance in the solar principal plane in the red band. For comparison, the sloping and horizontal reflectance simulated by the canopy reflectance model are also shown. The view zenith angles ranged from −40° to 40° with an interval of 5°, and the negative view zenith angles represent backward directions relative to the solar incidence.

sloping surface. To account for the difference in the reference frames of the path length calculation and radiative transfer equation simplification, a geometrical adjustment factor of cos(θ1)/cos(i) should be multiplied by the corrected reflectance. This adjustment is similar to the normalization of the LAI to the horizontal surface in the case of in-situ LAI measurements using Eq. (9) (Gonsamo et al., 2011; Walter and Torquebiau, 2000). The final correction results after this adjustment are shown in Fig. 11(b) and (d). Visual inspection shows that Eq. (9) can obtain excellent correction results. The performances of Eqs. (9) and (11) under different conditions are worthy of further comparison. However, we argue that Eq. (11) may be preferable because of the underlying consistent definition of the path length and corrected reflectance. 5.3. Assumptions in the formulation of the PLC method In the derivation of the PLC method, several assumptions were involved which may incur uncertainties for the topographic correction: Firstly, only direct radiation from the sun is considered, and the diffuse component from the sky and the surrounding terrain is neglected. Neglecting these factors may cause considerable uncertainty in forward radiance modeling, especially for pixels with steep slopes (Mousivand et al., 2015). This uncertainty will propagate to the topographical correction results. Here, we preliminarily analyzed the influence of diffuse irradiance using the following two methods: During the topographic correction for the Landsat 8 OLI image, we took the sky view factor as proxy of the diffuse component. The sky view factor is defined as the relative proportion of the solid angle of the sky that is viewed from a pixel and is not blocked by the surrounding topography (ranges from 0 to 1) (Zaksek et al., 2011). A small sky view factor indicates a large proportion of diffuse component in the incoming irradiation above the canopy, and vice versa. The R2 between the PLCcorrected reflectance and the cosine of the local solar incidence angle was calculated, and the change in R2 with the sky view factor was analyzed (Fig. 12). Generally, R2 decreased with an increasing sky view factor, demonstrating that neglecting the diffuse component decreases the correction efficiency for near-nadir observations by a factor of 2.0 to 3.0 in a highly obstructed case. We also analyzed the influence of neglecting diffuse irradiance through canopy reflectance model simulation. We generated another 20 PLC validation datasets with diffuse irradiance fractions from 0.05 to 1.0 with a step of 0.05, with the same input parameters as in section 3.2.1. The change in MCV with diffuse irradiance fraction is presented in Fig. 13. As in the Landsat 8 OLI, the ability of PLC to mitigate the topographic effect decreases as the diffuse irradiance fraction increases. In the extreme case of a diffuse irradiance

Fig. 10. Coefficient of variation (CV) for the PLC and SCS + C corrected red band reflectance across view zenith angles in the solar principal plane.

using the models with the retrieved parameters but with a flat terrain (slope = 0° and aspect = 0°). The model inversion procedure makes the MFM-TOPO and MFM-3-D methods time-consuming. In addition, the ill-posed nature of the inverse problem may undermine their robustness (Quan et al., 2015). In contrary to MFM-TOPO and MFM-3-D, the PLC method uses the canopy reflectance model implicitly. We established PLC by simplifying the classic radiative transfer theory. Therefore, PLC inherits the physical soundness of a physical model with a concise form.

5.2. Effects of path length formulation To ensure physical consistency between the path length formulation and the radiative transfer equation simplification, we used Eq. (11) rather than Eq. (9) to calculate the path length. Eq. (9) has also been extensively utilized, especially for in-situ LAI measurements over sloping terrain (Gonsamo et al., 2011; Walter and Torquebiau, 2000). Therefore, it may be worthy to test its potential in topographic correction. Fig. 11(a) and (c) represents the corresponding corrected results for the Landsat 8 OLI image and the enlarged area, respectively. Contrary to expectation, the topographic effect was not effectively removed. Eq. (9) calculates the path length relative to the normal of the sloping surface (Gonsamo et al., 2011; Walter and Torquebiau, 2000). Therefore, the radiative transfer equation should be transferred to the same 194

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Fig. 11. PLC-corrected images with path length formulation from Eq. (9) before (a) and after geometric adjustment (b). The enlarged areas in (a) and (b) are shown in (c) and (d), respectively.

Fig. 12. R2 between the PLC corrected reflectance and the cosine of the local solar incidence angle for different sky view factors.

Fig. 13. MCV of the PLC corrected reflectance as a function of the diffuse irradiance fraction.

fraction of 1.0, the MCV increased by a factor of approximately 3.0 for wide field-of-view observations. The analysis of topographic correction for Landsat 8 OLI and model simulated reflectance both demonstrated that neglecting diffuse irradiance affects the performance of PLC. However, it does not significantly reduce its suitability because our method is focused on the correction of BRDF distortion caused by

sloping terrain. For a highly diffuse irradiance-illuminated pixel, the sensor-observed radiance may be changed more by illumination variation by surrounding pixels than BRDF variation from surface orientation. Therefore, for this case, an effective illumination correction method can be combined with our PLC BRDF correction method to further improve the topographic correction. 195

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Ni-Meister et al., 2010; Pinty et al., 2004). The effective LAI can be estimated by multiplying the true LAI by a clumping index (Ryu et al., 2010). Our study further confirms the validity of the 1-D model framework in BRDF simulation and correction.

The second assumption involved in PLC is that the background reflectance and the multiple-scattering of leaves are negligible. For dense vegetation, e.g., when the LAI is greater than 3.0, the reflectance signal is saturated, and neglecting the background reflectance would not cause an obvious influence. When the vegetation is sparse, the influence of the background is determined by its reflectance anisotropy. A Lambertian background, which is a common assumption in canopy reflectance modeling, would not affect the efficiency of PLC. However, the anisotropic background, e.g., snow, may significantly influence the canopy BRDF (Manninen and Stenberg, 2009). The sloping canopy reflectance model used for PLC assessment was based on an isotropic background (Yin et al., 2015), so the influence of background anisotropy on PLC performance was not specifically analyzed in this study. The multiple-scattering of leaves may be safely neglected in topographic correction because the topography-induced distortion for multiple-scattering is much smaller than that for single-scattering. This is also why the influence of surface orientation on multiple-scattering is often neglected in sloping reflectance modeling (Combal et al., 2000; Schaaf et al., 1994; Yin et al., 2017). Thirdly, the leaf inclination distribution function is assumed to be spherical, in which case the leaf area projected function G is always 0.5 for all incidence and observation directions (Ross, 1981). The spherical leaf inclination distribution is often implied in many applications. However, Pisek et al. (2013) found that many tree species are characterized by planophile or plagiophile leaf angle distributions rather than spherical. To examine the influence of leaf area distribution, we compared the MCV for six typical leaf angle distributions: spherical, planophile, erectophile, plagiophile, extremophile and uniform. The results show that leaf angle distribution has a negligible influence on PLC performance (Fig. 14). In addition, the hot spot effect was not accounted for in the derivation of the PLC method. However, it would not significantly influence the correction results because the hot spot effect can easily be accounted for by multiplying Eq. (4) by a hot spot function (Dymond et al., 2001; Hapke, 1981; Liangrocapart and Petrou, 2002). If the hot spot function is not sensitive to terrain (Yin et al., 2017), then it can be canceled when calculating the ratio between ρRT and ρt during the derivation of Eq. (7). Finally, the PLC method was derived from a 1-D model framework and assumes a random distribution of leaves within canopy. The 1-D model framework fails to correctly represent many vegetation canopies that are characterized by various levels of clumping (Widlowski et al., 2005). However, a 1-D model can mimic the reflectance of 3-D canopy if the effective LAI is used rather than the true LAI (Duthoit et al., 2008;

5.4. Prospects for future studies Illumination-observation geometry is explicitly involved in the PLC method, so it can be used for wide field-of-view sensors. It is more difficult to evaluate the performance of topographic correction methods for wide field-of-view observations than for near-nadir viewing observations because canopy reflectance is by nature anisotropic, and successful mitigation of the topographic effects may not lead to similar reflectance for canopies with identical structure but different observation geometries. Therefore, a canopy reflectance model dedicated for sloping terrain (Yin et al., 2017) was used to simulate wide field-ofview observations. The advantage of this evaluation method is that reference values for the corrected reflectance can be extracted from model simulations with identical canopy structures and illuminationobservation geometries but with flat terrain (both the slope and aspect were set to 0). The performance of PLC for wide field-of-view observations was primarily assessed via this method. The applicability of PLC for real wide field-of-view images, e.g., the newly launched GF series (Feng et al., 2016), will be analyzed in a future paper. Existing studies on biophysical parameter retrieval for mountainous areas mostly used canopy reflectance models developed for horizontal surfaces (Gonsamo and Chen, 2014; Pasolli et al., 2015). In these studies, the topographic effect was accounted for by simple normalization of the illumination-observation geometry from a horizontal datum to the local terrain by coordinate rotation. Such treatment implicitly assumes that vegetation grows perpendicularly to local slopes and violates the geotropic nature of vegetation, inducing large uncertainties in the retrieved parameters (Soenen et al., 2010). PLC will facilitate biophysical parameter retrieval over mountainous areas because it can transform the BRDF values over sloping terrain to their horizontal counterparts with identical biophysical parameters. Thus, the canopy reflectance models developed for horizontal surfaces can be safely used without violating the geotropic nature. In addition, multi-angle data, an important data source for biophysical parameter retrieval (Chen et al., 2003), can be readily exploited using PLC because illumination-observation geometry is explicitly involved. It is worth noting that the geotropic nature of vegetation growth, which is an important vegetation property considered in our PLC method, may be questionable. On slopes, the tree crowns are often asymmetric. For example, crowns may expand down-slope to receive more solar radiation (Walter and Torquebiau, 2000). Therefore, the actual vegetation orientation over sloping surfaces may be in between the STS and SCS cases. This explains why the STS paradigm can also obtain satisfactory topographic correction results in some scenarios. The respective suitable scenarios for STS and STS should be clarified in future work. Topography affects the BRDF of a canopy (Combal et al., 2000; Fan et al., 2014; Schaaf et al., 1994; Yin et al., 2017) and also significantly changes the illumination condition (Kobayashi and Sanga-Ngoie, 2008; Shepherd and Dymond, 2003; Tan et al., 2013). PLC implements topographic correction in terms of BRDF correction. Several mountain radiative transfer simulation frameworks have been developed (Li et al., 2012; Mousivand et al., 2015; Sandmeier and Itten, 1997; Shepherd and Dymond, 2003). Due to the simplicity of the PLC formula, it can be easily integrated into these frameworks to provide an integrated radiometric correction for topographic and atmospheric effects. Illumination correction should be included in future research. 6. Conclusions

Fig. 14. Influence of leaf angle distribution on PLC.

From a simplification of the classic radiative transfer equation, this 196

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

paper proposed a semi-physical and simple topographic correction method for vegetation canopies based on path length correction (PLC). A Landsat 8 OLI image (representing near-nadir observations) and sloping reflectance simulated using a canopy reflectance model (representing wide field-of-view observations) were combinedly used to test the PLC correction method, and a multi-criteria analysis was used in the evaluation. The results were compared with empirical parameterbased methods, i.e., CC, SCS + C and SE, and similar semi-physical methods without empirical parameter, i.e., SCS and D-S. Generally, all methods significantly reduced the topographic effect. However, SCS showed obvious overcorrection for near-nadir observations. The D-S corrected reflectance showed a systematic positive bias. Empirical parameter-based methods (CC, SCS + C and SE) performed excellently for near-nadir observations. PLC also generated satisfactory correction results for near-nadir observations without overcorrection and radiometric instability. For wide field-of-view observations, PLC outperformed all other methods. PLC provides an efficient approach to correct the terrain-induced canopy BRDF distortion and will facilitate biophysical parameter retrieval over mountainous areas, excluding selfshadow and cast-shadow areas.

based on subpixel Sun-canopy-sensor geometry. Remote Sens. Environ. 64, 166–175. Hantson, S., Chuvieco, E., 2011. Evaluation of different topographic correction methods for Landsat imagery. Int. J. Appl. Earth Obs. Geoinf. 13, 691–700. Hapke, B., 1981. Bidirectional reflectance spectroscopy. 1. Theory. J. Geophys. Res. 86, 3039–3054. Hu, R., Yan, G., Mu, X., Luo, J., 2014. Indirect measurement of leaf area index on the basis of path length distribution. Remote Sens. Environ. 155, 239–247. Hu, R., Yan, G., Nerry, F., Liu, Y., Jiang, Y., Wang, S., Chen, Y., Mu, X., Zhang, W., Xie, D., 2018. Using airborne laser scanner and path length distribution model to quantify clumping effect and estimate leaf area index. IEEE Trans. Geosci. Remote Sens. 56, 3196–3209. Kane, V.R., Gillespie, A.R., McGaughey, R., Lutz, J.A., Ceder, K., Franklin, J.F., 2008. Interpretation and topographic compensation of conifer canopy self-shadowing. Remote Sens. Environ. 112, 3820–3832. Knyazikhin, Y.V., Marshak, A.L., Myneni, R.B., 1992. Interaction of photons in a canopy of finite-dimensional leaves. Remote Sens. Environ. 39, 61–74. Kobayashi, S., Sanga-Ngoie, K., 2008. The integrated radiometric correction of optical remote sensing imageries. Int. J. Remote Sens. 29, 5957–5985. Li, X.W., Strahler, A.H., 1992. Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy - effect of crown shape and mutual shadowing. IEEE Trans. Geosci. Remote Sens. 30, 276–292. Li, F.Q., Jupp, D.L.B., Thankappan, M., Lymburner, L., Mueller, N., Lewis, A., Held, A., 2012. A physics-based atmospheric and BRDF correction for Landsat data over mountainous terrain. Remote Sens. Environ. 124, 756–770. Li, A., Wang, Q., Bian, J., Lei, G., 2015. An improved physics-based model for topographic correction of Landsat TM images. Remote Sens. 7, 6296–6319. Liangrocapart, S., Petrou, M., 2002. A two-layer model of the bidirectional reflectance of homogeneous vegetation canopies. Remote Sens. Environ. 80, 17–35. Luisa, E.M., Frederic, B., Marie, W., 2008. Slope correction for LAI estimation from gap fraction measurements. Agric. For. Meteorol. 148, 1553–1562. Manninen, T., Stenberg, P., 2009. Simulation of the effect of snow covered forest floor on the total forest albedo. Agric. For. Meteorol. 149, 303–319. METI, NASA, 2011. ASTER GDEM 2 README. Mousivand, A., Verhoef, W., Menenti, M., Gorte, B., 2015. Modeling top of atmosphere radiance over heterogeneous non-Lambertian rugged terrain. Remote Sens. 7, 8019–8044. Myneni, R.B., Marshak, A.L., Knyazikhin, Y.V., 1991. Transport-theory for a leaf canopy of finite-dimensional scattering centers. J. Quant. Spectrosc. Radiat. Transf. 46, 259–280. Ni-Meister, W., Yang, W.Z., Kiang, N.Y., 2010. A clumped-foliage canopy radiative transfer model for a global dynamic terrestrial ecosystem model. I: theory. Agric. For. Meteorol. 150, 881–894. Pasolli, L., Asam, S., Castelli, M., Bruzzone, L., Wohlfahrt, G., Zebisch, M., Notarnicola, C., 2015. Retrieval of leaf area index in mountain grasslands in the Alps from MODIS satellite imagery. Remote Sens. Environ. 165, 159–174. Pinty, B., Gobron, N., Widlowski, J.L., Gerstl, S.A.W., Verstraete, M.M., Antunes, M., Bacour, C., Gascon, F., Gastellu, J.P., Goel, N., Jacquemoud, S., North, P., Qin, W.H., Thompson, R., 2001. Radiation transfer model intercomparison (RAMI) exercise. J. Geophys. Res.-Atmos. 106, 11937–11956. Pinty, B., Gobron, N., Widlowski, J.L., Lavergne, T., Verstraete, M.M., 2004. Synergy between 1-D and 3-D radiation transfer models to retrieve vegetation canopy properties from remote sensing data. J. Geophys. Res.-Atmos. 109. Pisek, J., Sonnentag, O., Richardson, A.D., Mottus, M., 2013. Is the spherical leaf inclination angle distribution a valid assumption for temperate and boreal broadleaf tree species? Agric. For. Meteorol. 169, 186–194. Proy, C., Tanre, D., Deschamps, P.Y., 1989. Evaluation of topographic effects in remotely sensed data. Remote Sens. Environ. 30, 21–32. Quan, X.W., He, B.B., Li, X., 2015. A Bayesian network-based method to alleviate the illposed inverse problem: a case study on leaf area index and canopy water content retrieval. IEEE Trans. Geosci. Remote Sens. 53, 6507–6517. Richter, R., Kellenberger, T., Kaufmann, H., 2009. Comparison of topographic correction methods. Remote Sens. 1, 184–196. Ross, J., 1981. The Radiation Regime and Architecture of Plant Stands. Junk, The Hague. Ryu, Y., Nilson, T., Kobayashi, H., Sonnentag, O., Law, B.E., Baldocchi, D.D., 2010. On the correct estimation of effective leaf area index: does it reveal information on clumping effects? Agric. For. Meteorol. 150, 463–472. Sandmeier, S., Itten, K.I., 1997. A physically-based model to correct atmospheric and illumination effects in optical satellite data of rugged terrain. IEEE Trans. Geosci. Remote Sens. 35, 708–717. Schaaf, C.B., Li, X.W., Strahler, A.H., 1994. Topographic effects on bidirectional and hemispherical reflectances calculated with a geometric-optical canopy model. IEEE Trans. Geosci. Remote Sens. 32, 1186–1193. Schaepman-Strub, G., Schaepman, M.E., Painter, T.H., Dangel, S., Martonchik, J.V., 2006. Reflectance quantities in optical remote sensing-definitions and case studies. Remote Sens. Environ. 103, 27–42. Schulmann, T., Katurji, M., Zawar-Reza, P., 2015. Seeing through shadow: modelling surface irradiance for topographic correction of Landsat ETM plus data. ISPRS J. Photogramm. Remote Sens. 99, 14–24. Shepherd, J.D., Dymond, J.R., 2003. Correcting satellite imagery for the variance of reflectance and illumination with topography. Int. J. Remote Sens. 24, 3503–3514. Soenen, S.A., Peddle, D.R., Coburn, C.A., 2005. SCS + C: a modified sun-canopy-sensor topographic correction in forested terrain. IEEE Trans. Geosci. Remote Sens. 43, 2148–2159. Soenen, S.A., Peddle, D.R., Coburn, C.A., Hall, R.J., Hall, F.G., 2008. Improved topographic correction of forest image data using a 3-D canopy reflectance model in multiple forward mode. Int. J. Remote Sens. 29, 1007–1027.

Acknowledgements This work was funded by the National Natural Science Foundation of China (41631180, 41601403, 41571373, 41401418), the National Key Research and Development Program of China (2016YFA0600103), the GF6 Project (30-Y20A03-90030-17/18), the CAS “Light of West China” Program and the Youth Talent Team Program of the Institute of Mountain Hazards and Environment, CAS (SDSQB-2015-02). We appreciate Professor Liangyun Liu for his helpful comments on our method. References Balthazar, V., Vanacker, V., Lambin, E.F., 2012. Evaluation and parameterization of ATCOR3 topographic correction method for forest cover mapping in mountain areas. Int. J. Appl. Earth Obs. Geoinf. 18, 436–450. Chen, J., Li, X., Nilson, T., Strahler, A., 2000. Recent advances in geometrical optical modelling and its applications. Remote Sens. Rev. 18, 227–262. Chen, J.M., Liu, J., Leblanc, S.G., Lacaze, R., Roujean, J.L., 2003. Multi-angular optical remote sensing for assessing vegetation structure and carbon absorption. Remote Sens. Environ. 84, 516–525. Combal, B., Isaka, H., Trotter, C., 2000. Extending a turbid medium BRDF model to allow sloping terrain with a vertical plant stand. IEEE Trans. Geosci. Remote Sens. 38, 798–810. Couturier, S., Gastellu-Etchegorry, J.P., Martin, E., Patino, P., 2013. Building a forwardmode three-dimensional reflectance model for topographic normalization of highresolution (1–5 m) imagery: validation phase in a forested environment. IEEE Trans. Geosci. Remote Sens. 51, 3910–3921. Duthoit, S., Demarez, V., Gastellu-Etchegorry, J.-P., Martin, E., Roujean, J.-L., 2008. Assessing the effects of the clumping phenomenon on BRDF of a maize crop based on 3D numerical scenes using DART model. Agric. For. Meteorol. 148, 1341–1352. Duursma, R.A., Marshall, J.D., Robinson, A.P., 2003. Leaf area index inferred from solar beam transmission in mixed conifer forests on complex terrain. Agric. For. Meteorol. 118, 221–236. Dymond, J.R., Shepherd, J.D., 1999. Correction of the topographic effect in remote sensing. IEEE Trans. Geosci. Remote Sens. 37, 2618–2620. Dymond, J.R., Shepherd, J.D., Qi, J., 2001. A simple physical model of vegetation reflectance for standardising optical satellite imagery. Remote Sens. Environ. 75, 350–359. Fan, W.L., Chen, J.M., Ju, W.M., Zhu, G.L., 2014. GOST: a geometric-optical model for sloping terrains. IEEE Trans. Geosci. Remote Sens. 52, 5469–5482. Feng, L., Li, J., Gong, W.S., Zhao, X., Chen, X.L., Pang, X.P., 2016. Radiometric crosscalibration of Gaofen-1 WFV cameras using Landsat-8 OLI images: a solution for large view angle associated problems. Remote Sens. Environ. 174, 56–68. Gastellu-Etchegorry, J.P., Martin, E., Gascon, F., 2004. DART: a 3D model for simulating satellite images and studying surface radiation budget. Int. J. Remote Sens. 25, 73–96. Gonsamo, A., Chen, J.M., 2014. Improved LAI algorithm implementation to MODIS data by incorporating background, topography, and foliage clumping information. IEEE Trans. Geosci. Remote Sens. 52, 1076–1088. Gonsamo, A., Walter, J.M.N., Pellikka, P., 2011. CIMES: a package of programs for determining canopy geometry and solar radiation regimes through hemispherical photographs. Comput. Electron. Agric. 79, 207–215. Gu, D., Gillespie, A., 1998. Topographic normalization of landsat TM images of forest

197

Remote Sensing of Environment 215 (2018) 184–198

G. Yin et al.

Wang, S., Chen, W.J., Cihlar, J., 2002. New calculation methods of diurnal distribution of solar radiation and its interception by canopy over complex terrain. Ecol. Model. 155, 191–204. Wen, J.G., Liu, Q.H., Liu, Q., Xiao, Q., Li, X.W., 2009. Parametrized BRDF for atmospheric and topographic correction and albedo estimation in Jiangxi rugged terrain, China. Int. J. Remote Sens. 30, 2875–2896. Wen, J.G., Liu, Q., Tang, Y., Dou, B.C., You, D.Q., Xiao, Q., Liu, Q.H., Li, X.W., 2015. Modeling land surface reflectance coupled BRDF for HJ-1/CCD data of rugged terrain in Heihe river basin, China. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 8, 1506–1518. Wen, J., Liu, Q., Xiao, Q., Liu, Q., You, D., Hao, D., Wu, S., Lin, X., 2018. Characterizing land surface anisotropic reflectance over rugged terrain: a review of concepts and recent developments. Remote Sens. 10, 370. Widlowski, J.L., Pinty, B., Lavergne, T., Verstraete, M.M., Gobron, N., 2005. Using 1-D models to interpret the reflectance anisotropy of 3-D canopy targets: issues and caveats. IEEE Trans. Geosci. Remote Sens. 43, 2008–2017. Yan, G.J., Hu, R.H., Wang, Y.T., Ren, H.Z., Song, W.J., Qi, J.B., Chen, L., 2016. Scale effect in indirect measurement of leaf area index. IEEE Trans. Geosci. Remote Sens. 54, 3475–3484. Yin, G.F., Li, J., Liu, Q.H., Fan, W.L., Xu, B.D., Zeng, Y.L., Zhao, J., 2015. Regional leaf area index retrieval based on remote sensing: the role of radiative transfer model selection. Remote Sens. 7, 4604–4625. Yin, G.F., Li, A.N., Zhao, W., Jin, H.A., Bian, J.H., Wu, S.B.A., 2017. Modeling canopy reflectance over sloping terrain based on path length correction. IEEE Trans. Geosci. Remote Sens. 55, 4597–4609. Zaksek, K., Ostir, K., Kokalj, Z., 2011. Sky-view factor as a relief visualization technique. Remote Sens. 3, 398–415.

Soenen, S.A., Peddle, D.R., Hall, R.J., Coburn, C.A., Hall, F.G., 2010. Estimating aboveground forest biomass from canopy reflectance model inversion in mountainous terrain. Remote Sens. Environ. 114, 1325–1337. Sola, I., Gonzalez-Audicana, M., Alvarez-Mozos, J., Torres, J.L., 2014. Synthetic images for evaluating topographic correction algorithms. IEEE Trans. Geosci. Remote Sens. 52, 1799–1810. Sola, I., Gonzalez-Audicana, M., Alvarez-Mozos, J., 2016. Multi-criteria evaluation of topographic correction methods. Remote Sens. Environ. 184, 247–262. Tachikawa, T., Hato, M., Kaku, M., Iwasaki, A., IEEE, 2011. Characteristics of ASTER GDEM version 2. In: 2011 IEEE International Geoscience and Remote Sensing Symposium. IEEE, New York, pp. 3657–3660. Tan, B., Masek, J.G., Wolfe, R., Gao, F., Huang, C.Q., Vermote, E.F., Sexton, J.O., Ederer, G., 2013. Improved forest change detection with terrain illumination corrected Landsat images. Remote Sens. Environ. 136, 469–483. Teillet, P.M., Guindon, B., Goodenough, D.G., 1982. On the slope-aspect correction of multispectral scanner data. Can. J. Remote. Sens. 8, 84–106. USGS, 2017a. Product Guide: Landsat 8 Surface Reflectance Code (LaSRC) Product. USGS, 2017b. User Guide: Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on Demand Interface. Vermote, E., Justice, C., Claverie, M., Franch, B., 2016. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 185, 46–56. Vicente-Serrano, S.M., Perez-Cabello, F., Lasanta, T., 2008. Assessment of radiometric correction techniques in analyzing vegetation variability and change using time series of Landsat images. Remote Sens. Environ. 112, 3916–3934. Walter, J.M.N., Torquebiau, E.F., 2000. The computation of forest leaf area index on slope using fish-eye sensors. C. R. Acad. Sci. III 323, 801–813.

198