Modeling Top of Atmosphere Radiance over Heterogeneous Non-Lambertian Rugged Terrain

Abstract

Topography affects the fraction of direct and diffuse radiation received on a pixel and changes the sun–target–sensor geometry, resulting in variations in the observed radiance. Retrieval of surface–atmosphere properties from top of atmosphere radiance may need to account for topographic effects. This study investigates how such effects can be taken into account for top of atmosphere radiance modeling. In this paper, a system for top of atmosphere radiance modeling over heterogeneous non-Lambertian rugged terrain through radiative transfer modeling is presented. The paper proposes an extension of “the four-stream radiative transfer theory” (Verhoef and Bach 2003, 2007 and 2012) mainly aimed at representing topography-induced contributions to the top of atmosphere radiance modeling. A detailed account for BRDF effects, adjacency effects and topography effects on the radiance modeling is given, in which sky-view factor and non-Lambertian reflected radiance from adjacent slopes are modeled precisely. The paper also provides a new formulation to derive the atmospheric coefficients from MODTRAN with only two model runs, to make it more computationally efficient and also avoiding the use of zero surface albedo as used in the four-stream radiative transfer theory. The modeling begins with four surface reflectance factors calculated by the Soil–Leaf–Canopy radiative transfer model SLC at the top of canopy and propagates them through the effects of the atmosphere, which is explained by six atmospheric coefficients, derived from MODTRAN radiative transfer code. The top of the atmosphere radiance is then convolved with the sensor characteristics to generate sensor-like radiance. Using a composite dataset, it has been shown that neglecting sky view factor and/or terrain reflected radiance can cause uncertainty in the forward TOA radiance modeling up to 5 (mW/m²·sr·nm). It has also been shown that this level of uncertainty can be translated into an over/underestimation of more than 0.5 in LAI (or 0.07 in fCover) in variable retrieval.

1. Introduction

Top-Of-Atmosphere (TOA) radiometric data are traditionally converted to Top-Of-Canopy (TOC) reflectance before being used for vegetation biophysical variable retrieval [1]. This conversion requires corrections for atmospheric, angular and topographic effects. To do retrieval with TOC reflectance it is necessary to separate atmospheric correction and inversion, while the strong mutual interactions between surface and atmosphere makes this approach hard to justify. Although most attention has been given to the retrieval of vegetation biophysical and chemical characteristics from the TOC reflectance [2,3,4,5,6,7,8,9,10], there are studies that highlight the advantages of a TOA radiance approach for variable retrieval [11,12,13].

Retrieval using TOA radiance has the advantage that no correction/conversion (e.g., atmospheric correction, angular effects and topographic (illumination) corrections) is necessary and the actual radiometric measurements can be used for retrieval. This is, in particular, important in the case of a multi-sensor retrieval approach [14], where data acquired by sensors with different characteristics are integrated in the retrieval. Furthermore, the TOA retrieval approach would also allow retrieving surface and atmospheric properties simultaneously if spectral radiance measurements are sensitive to those properties.

Vegetation biophysical variable retrieval using TOA radiance data requires modeling of TOA radiance. Radiative transfer models explain the propagation of radiation through a medium (e.g., atmosphere) and have been frequently and successfully used to simulate reflectance and radiance of soil–vegetation–atmosphere systems [15,16,17]. TOA radiance modeling involves combining radiative transfer models of surface and atmosphere. The radiance coming into the field of view of a sensor is the integrated sum of the surface reflectance, atmospheric effects and target’s surroundings. Therefore, a careful modeling of each contribution is required to accurately calculate TOA radiance. A number of comprehensive simulation tools coupling surface and atmosphere radiative transfer processes have been introduced [18,19,20,21,22,23,24,25,26,27,28]. However, some of these studies are often simplified through various assumptions such as flat terrain, Lambertian assumption for the surface reflectance and some others are more accurate but at very high computational cost that limits their application.

Among them, the four-stream approach proposed by Verhoef and Bach [13,26,27] is a good example of a trade-off between computational efficiency and accuracy. It integrates the Soil–Leaf–Canopy radiative transfer model SLC with the atmospheric radiative transfer model MODTRAN. The approach has been validated against both synthetic and satellite data [13,29] and it has been used for the retrieval of vegetation properties over flat terrain [11,14,30]. Variable retrieval from TOA radiance over rugged terrain may need to account for topographic effects.

Topography has three major effects on the radiance sensed by a sensor. First, the total optical depth (gases and aerosols) decreases with increase in elevation. Second, topography regulates the fraction of direct and diffuse irradiance received on a pixel from the sun and the sky, respectively. Finally, topography changes the orientation of a target and therefore alters the illumination and viewing angles, which affect reflected radiance because of the surface anisotropic reflectance. Topographic effects are determined by the aspect and slope of the surface and the local horizon, and also by the geometry of the surrounding areas [31]. The topographic effects can be as important as the atmospheric effects [32].

The four-stream approach is used for the simulation of TOA radiance. The approach was initially developed for flat terrain and it has been recently extended to account for sloped terrain as well [27]. However, the extended approach only accounts for the slope and aspect of the observed pixel, while the topographic effects caused by the adjacent features (e.g., cast shadow and terrain reflected radiance) are not taken into account. A further concern is related to the extraction of the atmospheric coefficients in the four-stream approach. The atmospheric coefficients are derived from three MODTRAN runs at 0, 0.5 and 1 surface albedos with the assumption of Lambertian flat surface [26,27]. Apart from the computational burden of three calls to the code, the use of zero albedo may lead to some inconsistencies in some of the MODTRAN features [33].

This study aims at addressing three issues: (i) investigating the effect of topography on TOA radiance; (ii) quantifying which slope and aspect requires integrated topography modeling; and (iii) how the topographic effect compares with the variability in TOA radiance due to variability of vegetation properties. In order to accomplish these goals, this paper presents an extension of the four-stream radiative transfer theory over heterogeneous non-Lambertian surfaces, in which topographic, adjacency and atmospheric effects are taken into account. Focus is given to the topography-induced contributions to the top of atmosphere radiance simulation. In this way, the impacts of topography on the irradiance and surface reflectance are precisely modeled using sky view factor and terrain reflected radiance from adjacent slopes. The surface reflectance at every location is modeled using SLC, while the atmosphere is modeled using MODTRAN. The MODTRAN Interrogation Technique (MIT) used for the extraction of atmospheric coefficients has been modified and reformulated to use only two MODTRAN runs instead of the three runs suggested in previous works and also to avoid the use of zero albedo in the calculation.

The remainder of the paper is structured as follows. Section 2 gives a detailed background of TOA radiance modeling using the four-stream theory. In Section 3, the extension of the four-stream theory over rugged terrain is described and topography effects on TOA radiance and surface reflectance are elaborated in details. The modified MODTRAN interrogation technique is also explained in Section 3. Section 4 describes the model implementation and the design of the numerical experiment is given in Section 5. The results of TOA radiance simulations considering different cases of rugged terrains are discussed in Section 6. Finally, Section 7 concludes the paper.

2. Background

The spectral radiance recorded by a sensor is the result of the complex interactions of solar radiation with the Earth’s atmosphere and surface. Following Vermote et al. [34], Verhoef and Bach [13,26], and Schott [35], the main paths through which radiation reaches the sensor over a flat surface are illustrated in Figure 1. Each contribution is briefly explained in Table 1.

Table 1. Different contributions in Figure 1 from the atmosphere and the surface to the sensor plus the corresponding fluxes in Figure 2.

Contribution in Figure 1	Description	Corresponding Flux in Figure 2
Contribution 1	Spectral radiation originating from the sun and, after propagation in the atmosphere, scattered into the field of view of the sensor without reaching the surface.	${\rho }_{so}{E}_s\left(t\right)$
Contribution 2	Direct radiation that goes through the atmosphere without being absorbed or scattered and after reflecting back from the surface it reaches the sensor.	integrated in ${\tau }_{oo}{E}_o(b)$
Contribution 3	Spectral radiation scattered by the atmosphere onto the target and then reflected back towards the sensor.	integrated in ${\tau }_{oo}{E}_o(b)$
Contribution 4 and 5	Spectral radiation that directly or diffusely reach the surrounding areas of the target and then reflected or scattered back into the field of view of the sensor. This effect is so called “adjacency effect” or “blurring effect”.	${\tau }_{do}{E}^{+}(b)$
Contribution 6	Diffuse radiation coming from the adjacent features into the field of view of the sensor. This contribution is incorporated in Contribution 3.	-
Contribution 7	So-called trapping effect and it is a part of the radiation reflected from the surface into the air column above the surface being scattered and ultimately reaches the sensor. This contribution is incorporated in Contribution 3.	-

Table 1. Different contributions in Figure 1 from the atmosphere and the surface to the sensor plus the corresponding fluxes in Figure 2.

Contribution in Figure 1	Description	Corresponding Flux in Figure 2
Contribution 1	Spectral radiation originating from the sun and, after propagation in the atmosphere, scattered into the field of view of the sensor without reaching the surface.	${\rho }_{so}{E}_s\left(t\right)$
Contribution 2	Direct radiation that goes through the atmosphere without being absorbed or scattered and after reflecting back from the surface it reaches the sensor.	integrated in ${\tau }_{oo}{E}_o(b)$
Contribution 3	Spectral radiation scattered by the atmosphere onto the target and then reflected back towards the sensor.	integrated in ${\tau }_{oo}{E}_o(b)$
Contribution 4 and 5	Spectral radiation that directly or diffusely reach the surrounding areas of the target and then reflected or scattered back into the field of view of the sensor. This effect is so called “adjacency effect” or “blurring effect”.	${\tau }_{do}{E}^{+}(b)$
Contribution 6	Diffuse radiation coming from the adjacent features into the field of view of the sensor. This contribution is incorporated in Contribution 3.	-
Contribution 7	So-called trapping effect and it is a part of the radiation reflected from the surface into the air column above the surface being scattered and ultimately reaches the sensor. This contribution is incorporated in Contribution 3.	-

As shown in Figure 1, the top of atmosphere radiance over a non-Lambertian, non-homogeneous Earth’s surface is determined by a combination of the surface reflectance, atmospheric effects, target’s surroundings effects and illumination and viewing geometry. This combination can be effectively modeled using the four-stream radiative transfer theory proposed by Verhoef and Bach [13,26,27]. The four-stream approach is applied in this study. Figure 2 illustrates the relevant fluxes at the top and bottom of the atmosphere considered in the four-stream approach. At the top of the atmosphere, the only downward flux is assumed to be the direct solar irradiance and the only upward flux is the radiance in the direction of the sensor. At the bottom of the atmosphere, four different fluxes are considered: an attenuated direct solar flux ($E_s\left(b\right)$); a surface radiance (times $\text{π}$) in the observation direction ($E_o\left(b\right)$); diffuse downward flux from the sky ($E^{-}\left(b\right)$); and diffuse upward flux ($E^{+}\left(b\right)$) are considered. We follow the notations adopted by Verhoef and Bach [13,26,27] in this study as presented in Table 2.

Figure 1. Various paths from the surface and the atmosphere to the Top-of-Atmosphere (TOA) radiance.

Figure 1. Various paths from the surface and the atmosphere to the Top-of-Atmosphere (TOA) radiance.

Figure 2. Four-stream radiation fluxes at the top (TOA) and the bottom (BOA) of the atmosphere (after Verhoef and Bach [26]).

Figure 2. Four-stream radiation fluxes at the top (TOA) and the bottom (BOA) of the atmosphere (after Verhoef and Bach [26]).

According to the four-stream theory, the TOA radiance in the direction of the sensor can be described as the sum of three components; the atmospherically scattered direct solar radiation (atmospheric radiance); diffuse upwelling flux scattered into the direction of the sensor (adjacency effect) and the reflected radiance from the surface (target radiance), given by

$$E_o\left(t\right)={\rho }_{so}{E}_s\left(t\right)+ {\tau }_{do}{E}^{+}(b)+{\tau }_{oo}{E}_o(b)$$

(1)

where

$$E_s\left(b\right)={\tau }_{ss}{E}_s(t)$$

(2)

$$E^{-}\left(b\right)={\tau }_{sd}{E}_s\left(t\right)+ {\rho }_{dd}{E}^{+}(b)$$

(3)

$$E^{+}\left(b\right)=\overline{r_{sd}}{E}_s\left(b\right)+\overline{r_{dd}}{E}^{-}(b)$$

(4)

$$E_o\left(b\right)=E_s\left(b\right)r_{so}+E^{-}\left(b\right)r_{do}$$

(5)

Equation (1) explains how the TOA radiance is generated from the sum of three components multiplied by atmospheric transmittance and hemispherical albedo. Equation (2) gives the attenuated direct solar flux at the bottom of the atmosphere, i.e., determined by the direct atmospheric transmittance from the sun to the ground. Equation (3) describes the diffuse downward flux (sky irradiance) by diffusely transmitted solar flux and upwelling diffuse flux from the surface that is reflected back by the atmosphere. Equation (4) describes the diffuse upward flux by the attenuated direct solar flux times the average surroundings diffuse reflectance and sky irradiance multiplied by the bi-hemispherical reflectance spatially averaged over surroundings. Finally, Equation (5) gives the non-Lambertian target reflected radiance in the observation direction as the sum of the direct incident flux times the bidirectional reflectance factor $r_{so}$ and of the diffuse incident flux times the directional reflectance $r_{do}$.

Natural surfaces usually exhibit anisotropic reflectance characteristics, which describe the dependency of the reflectance spectra on the illumination and viewing geometry. This is defined by the Bidirectional Reflectance Distribution Function (BRDF). The BRDF of vegetation canopies can be described by canopy radiative transfer models, i.e., SLC in this study. To properly model the radiative interactions between the surface and atmosphere, complete BRDF functions are required, which means three- or four-dimensional functions depending on the undertaken assumption of azimuthal symmetry [26]. The anisotropic reflectance patterns of vegetated surfaces in a non-Lambertian approach have to be described by a minimum of four reflectance factors [26]. In this study, the four directional reflectance factors, $r_{so}$, $r_{do}$, $r_{sd}$, $r_{dd}$, calculated by the SLC describe the anisotropic reflectance and can, therefore, be used to determine radiance reflected in any direction of observation, given the irradiance and the direction of illumination.

By substituting different equations in Equation (1) and solving the diffuse fluxes at the surface, the TOA radiance can be obtained [26].

$$L_{TOA}=\frac{{E_{\text{s}}}^o\cos {\theta }_s}{\text{π}} \left[\stackrel{atmosphere}{\overbrace{{\rho }_{so}}}+\stackrel{adjacency}{\overbrace{\frac{{\tau }_{ss}\overline{r_{sd}}+{\tau }_{sd}\overline{r_{dd}}}{1-\overline{r_{dd}}{\rho }_{dd}}{\tau }_{do}}}+\stackrel{sun}{\overbrace{{\tau }_{ss}{r}_{so}{\tau }_{oo}}} + \stackrel{sky}{\overbrace{\frac{{\tau }_{sd} +{\tau }_{ss}\overline{r_{sd}}{\rho }_{dd}}{1-\overline{r_{dd}}{\rho }_{dd}}{r}_{do}{\tau }_{oo}}}\right]$$

(6)

Four surface reflectance factors plus six atmospheric coefficients are required in the four-stream approach to model TOA radiance (see Table 2 for definitions). The surface reflectance components are described by the SLC radiative transfer model [26] and the atmospheric coefficients are derived from the MODTRAN4 atmospheric radiative transfer model. The multiplicative factor $1-\overline{r_{dd}}{\rho }_{dd}$ in Equation (6) accounts for multiple reflections between the surface and the atmosphere [24].

There is no term in Equation (6) to account for the topography effects. Verhoef and Bach [27] proposed an extension of Equation (6), in which topography effects are considered. The proposed approach, however, does not fully account for topographic impacts, since only the slope and aspect of the observed pixel itself are taken into account, while the obstructing effects of neighboring facets are not.

Table 2. Terminology used in the study (after Verhoef and Bach [13,26,27]).

Factor	Description	SI Unit
rso	surface bidirectional reflectance factor	-
$r_{do}$	surface hemispherical-directional reflectance for diffuse incidence	-
$\overline{r_{sd}}$	spatially averaged directional-hemispherical reflectance over surroundings	-
$\overline{r_{dd}}$	spatially averaged bi-hemispherical reflectance over surroundings	-
${\rho }_{so}$	TOA atmospheric bidirectional reflectance	-
${\rho }_{dd}$	BOA spherical albedo of the atmosphere	-
${\tau }_{ss}$	direct atmospheric transmittance from the sun to the ground	-
${\tau }_{oo}$	direct atmospheric transmittance from the ground to the sensor	-
${\tau }_{sd}$	diffuse atmospheric transmittance from the sun to the ground	-
${\tau }_{do}$	diffuse atmospheric transmittance from the ground to the sensor	-
${E_{\text{s}}}^o$	extraterrestrial solar irradiance on a plane perpendicular to the sunrays	Wm−2
$E_{\text{s}}$	direct solar irradiance on a horizontal plane	Wm−2µm−1
$E_o$	radiance in the viewing direction, times π	Wm−2µm−1
$E^{-}$	downward diffuse irradiance	Wm−2µm−1
$E^{+}$	upward diffuse irradiance	Wm−2µm−1
$\mathrm{\zeta }$	binary coefficient to indicate cast-shadowed pixels (0 or 1)	-
${\text{θ}}_{\text{s}}$	local solar zenith angle	Degree
${\text{θ}}_o$	viewing zenith angle	Degree
$i$	angle between the surface normal and the sun ray (illumination angle)	Degree
$\beta$	surface slope angle	Degree
${\phi }_s$	solar azimuth angle	Degree
${\phi }_t$	slope azimuth angle (aspect)	Degree
$t$	indicates the quantity at the top of atmosphere	-
$b$	indicates the quantity at the bottom of atmosphere	-

Table 2. Terminology used in the study (after Verhoef and Bach [13,26,27]).

Factor	Description	SI Unit
rso	surface bidirectional reflectance factor	-
$r_{do}$	surface hemispherical-directional reflectance for diffuse incidence	-
$\overline{r_{sd}}$	spatially averaged directional-hemispherical reflectance over surroundings	-
$\overline{r_{dd}}$	spatially averaged bi-hemispherical reflectance over surroundings	-
${\rho }_{so}$	TOA atmospheric bidirectional reflectance	-
${\rho }_{dd}$	BOA spherical albedo of the atmosphere	-
${\tau }_{ss}$	direct atmospheric transmittance from the sun to the ground	-
${\tau }_{oo}$	direct atmospheric transmittance from the ground to the sensor	-
${\tau }_{sd}$	diffuse atmospheric transmittance from the sun to the ground	-
${\tau }_{do}$	diffuse atmospheric transmittance from the ground to the sensor	-
${E_{\text{s}}}^o$	extraterrestrial solar irradiance on a plane perpendicular to the sunrays	Wm−2
$E_{\text{s}}$	direct solar irradiance on a horizontal plane	Wm−2µm−1
$E_o$	radiance in the viewing direction, times π	Wm−2µm−1
$E^{-}$	downward diffuse irradiance	Wm−2µm−1
$E^{+}$	upward diffuse irradiance	Wm−2µm−1
$\mathrm{\zeta }$	binary coefficient to indicate cast-shadowed pixels (0 or 1)	-
${\text{θ}}_{\text{s}}$	local solar zenith angle	Degree
${\text{θ}}_o$	viewing zenith angle	Degree
$i$	angle between the surface normal and the sun ray (illumination angle)	Degree
$\beta$	surface slope angle	Degree
${\phi }_s$	solar azimuth angle	Degree
${\phi }_t$	slope azimuth angle (aspect)	Degree
$t$	indicates the quantity at the top of atmosphere	-
$b$	indicates the quantity at the bottom of atmosphere	-

3. Four-Stream Theory over Rugged Terrain

In this section, we present an extension of the four-stream theory with the focus on modeling the topography-induced effects on the TOA radiance. Accordingly, Equation (1) will have to be modified to account for the topography effects. Among different terms in this equation, the reflected radiance from the surface (${\tau }_{oo}{E}_o(b)$) is directly related to the surface topography and to the anisotropy of surface reflectance. In this way, we reformulate $E_o(b)$ to include both the diffuse downward flux and the attenuated direct flux to consider topography effects.

3.1. Irradiance Modeling

The total solar irradiance reaching the ground surface on an inclined surface can be expressed as the sum of direct (from the sun) and diffuse (from the sky and the terrain) irradiance components. Figure 3 depicts the contribution of different components to the incoming irradiance on an inclined surface. The ratio of direct and diffuse fluxes is a function of the surface orientation, elevation and shadows cast by topographic features [36]. A Digital Elevation Model (DEM) is required to model topography effects on irradiance. Multiplying the irradiance components by the relevant non-Lambertian surface reflectance factors will give the total outgoing radiance from an inclined surface headed toward the sensor in the four-stream approach. We express the total surface reflected radiance in the observation direction as the sum of the direct ($E_{o_{dir}}\left(b\right)$), diffuse ($E_{o_{diff}}\left(b\right)$) and terrain ($E_{o_{terr}}\left(b\right)$) radiance components;

$$E_o\left(b\right)=E_{o_{dir}}\left(b\right)+E_{o_{diff}}\left(b\right)+E_{o_{terr}}\left(b\right)$$

(7)

In the following sections, we elaborate on each of these components and discuss how they will interact with different upward and downward fluxes.

Figure 3. Contributions of irradiance on an inclined surface; direct, diffuse and terrain irradiances.

Figure 3. Contributions of irradiance on an inclined surface; direct, diffuse and terrain irradiances.

Direct Irradiance: The direct irradiance is the main contribution to the total solar irradiance under clear sky conditions. This component depends on target’s elevation (which is considered in the atmospheric radiative transfer model), sun–target geometry and the horizon profile that determines shadow distribution. The direct irradiance multiplied by the bidirectional reflectance factor ($r_{so}$) of the surface yields the direct radiance at the ground level reflected towards the sensor.

$$E_{o_{dir}}\left(b\right)=\mathrm{\zeta }{E}_s\left(b\right)\frac{\cos i}{\cos {\theta }_s}{r}_{so}$$

(8)

where $\cos i=\cos {\theta }_s\cos \beta +\sin {\theta }_s\sin \beta \cos ({\phi }_s-{\phi }_t)$.

If $\cos i$ is negative, the pixel is self-shadowed and does not receive direct radiation. A pixel may also be cast-shadowed when the direct path between the pixel and the sun is intercepted by another terrain feature. This binary factor $\mathrm{\zeta }$ (0 or 1) indicates whether the pixel is cast-shadowed by adjacent features [37].

Diffuse Irradiance: The diffuse irradiance is a function of the proportion of sky that is visible from the pixel. A simple assumption is to assume that this quantity is isotropic over the entire sky hemisphere. However, in reality, this quantity exhibits some anisotropic behavior, especially when the slope is greater than 20 degree [38], and a more realistic estimate would be to decompose it into isotropic and anisotropic circumsolar components. The ratio of isotropic and anisotropic components can be derived by means of Hay’s [39] “anisotropy index” ($k$). The index $k$ is defined as the ratio of circumsolar diffuse irradiance to the total sky diffuse irradiance. The direct atmospheric transmittance along the sun–target path is used to define this index [39] such that $k=\frac{B_n}{I_0}\approx {\tau }_{ss}$, where $B_n$ denotes the beam radiation on a surface normal to the solar ray and $I_0$ is the extraterrestrial solar constant. In atmospheric radiative transfer models, this index is simply the direct transmittance from the sun to the ground. The diffuse irradiance multiplied by the hemispherical-directional reflectance factor ($r_{do}$) of the surface gives the diffuse radiance at the ground level reflected towards the sensor.

$$E_{o_{diff}}\left(b\right)=E^{-}\left(b\right)\left(\underset{anisotropic}{\underbrace{k\frac{\cos i}{\cos {\theta }_s}}}+\underset{isotropic}{\underbrace{\left(1-k\right)V_{sky}}}\right)r_{do}$$

(9)

The factor $V_{sky}$, known as sky view factor (Figure 4 left), is defined as the relative proportion of the solid angle of the sky, which is seen from a pixel and is not blocked by the surrounding topography (range between 0 and 1). The most widely used method is the simple trigonometric approximation presented in Liu and Jordan [40] and many others:

$$V_{sky}= \frac{1+\cos \beta }{2}$$

(10)

This expression is not in accordance with the definition of the sky view factor and does not represent it properly [41]. The only variable is the slope angle of the surface itself and the obstructing effects of surrounding topography are not considered. For an infinitely long slant surface or semi-rugged areas this would not lead to significant errors in the calculations, but over rugged terrain the effect of the surrounding topography plays an important role and must be taken into account as well. To account for this effect, the algorithm proposed by Dozier et al. [37,42] is applied in this study that provides a more accurate and analytical approximation of the sky view factor and is computationally efficient. The algorithm calculates the local horizon angles for a number of equally distributed directions. The horizon angle is defined as the angle between the zenith and the horizon surface at a given orientation direction. The local horizon angles have to be calculated over the entire azimuthal range (360°) around the target. However, in practice, a specified number of orientation directions (e.g., 16 directions) are taken to determine the horizon angles for the complete circle [42]. The horizon can result either from “self-shadowing” or from adjacent ridges. Horizon angles are interpolated for all other unsearched directions. Then, the sky view factor is calculated by integrating the horizon angles over the azimuth circle.

$$V_{sky}=\frac{1}{\pi }\underset{0}{\overset{2\pi }{{\displaystyle \int }}}\underset{0}{\overset{h_{\phi }}{{\displaystyle \int }}}\sin \theta \left[\cos \theta \cos \beta +\sin \theta \sin \beta \cos (\phi -{\phi }_t)\right]d\theta d\phi$$

(11)

Figure 4. An example of sky view factor (left) and terrain irradiance from pixel P into pixel M (right).

Figure 4. An example of sky view factor (left) and terrain irradiance from pixel P into pixel M (right).

If N is the number of discrete set of orientation directions over the entire circle and $h_{{\phi }_i}$ represents the horizon angle in direction i, then

$$V_{sky}\approx \frac{1}{N}{\displaystyle \sum _{i=1}^N\cos \beta {\sin }^2{h}_{{\phi }_i}+\sin \beta \cos ({\phi }_i-{\phi }_t)(h_{{\phi }_i}-\sin h_{{\phi }_i}\cos h_{{\phi }_i})}$$

(12)

The sky view factor is an intrinsic characteristic of topography and does not considerably change over time. It is, therefore, adequate to compute this factor once at high azimuthal resolution (i.e., for a large number of orientation directions) before the simulation. In this study, we calculated the sky view factor at each pixel for 64 azimuth directions. This number was selected based on some preliminary analysis that showed increasing directions beyond 64, had a very negligible effect on the accuracy of the sky view factor, at least for the terrain considered in our experiment. Thus, we deemed this number sufficient to give accurate estimates of the sky view factor.

Terrain Irradiance: In rugged terrains, apart from the diffuse radiation from the sky dome, topography can add another diffuse irradiance source to each pixel [38]. In this way, the adjacent slopes visible from the pixel will also contribute to the diffuse irradiance reaching the pixel. For instance, a part of the irradiance received by pixel P in Figure 4 (right) may be reflected towards pixel M. This component is often considered to be relatively insignificant and is neglected in most studies [31,38]. However, if the surface’s cover is highly reflective (e.g., snow and ice) or in mountainous areas with deep valleys, it is important to estimate this component. As an example, Gratton et al. [31], showed that in a glacier basin between 10% and 25% of the total irradiance at any location was terrain irradiance for non-shadowed and large shadowed areas.

An accurate method to calculate irradiance from the adjacent terrain was used in Proy et al. [38]. We hereafter refer to this as the exact method. The method calculates the sum of the irradiance coming from all the pixels that are seen from any observed pixel. The terrain irradiance at a pixel M, coming from pixel P is given by:

$$E_{terr}={\displaystyle \sum _{i=1}^N\frac{\cos T_M\cos T_{P}d{S}_{M}d{S}_P}{R^2}}{L_o}^P\left(b\right)$$

(13)

where $T_M$ and $T_P$ are the angles between normal and the line that connects two pixels (Figure 4 right), $d{S}_P$ and $d{S}_M$ are the area of pixels P and M, $R$ is the distance between the pixels, N is the number of visible pixels from M and ${L_o}^P\left(b\right)$ is the total radiance of pixel P coming towards pixel M that can be calculated using Equation (7). However, there is a challenge in accurately determining the total radiance coming from pixel P, since at the same time ${L_o}^P\left(b\right)$ is both the input for terrain irradiance (Equation (13)) and the output in Equation (7). The solution applied in this study was to ignore the terrain irradiance in the first iteration and calculate ${L_o}^P\left(b\right)$ from the sum of the direct and diffuse terms considering the geometry relationship between pixels M and P, and then iterate this procedure to account for the terrain irradiance as well. The attenuation of the radiance between two pixels is ignored. The effective distance from which the radiance of neighboring pixels arrives at the observed pixel can be considered between 0.5 and 1 km [43], which in this study we have considered a radius of 1 km from the observed pixel.

The exact method is computationally intensive and it is often replaced by an approximation. In this case, the portion of the surrounding pixels seen from the observed pixel (known as terrain view factor) is estimated by a simple trigonometric approximation,

$$V_{terrain}\approx \frac{1+\cos \beta }{2}-V_{sky}$$

(14)

Then, the average radiance scattered off the neighboring areas within the horizon of the observed pixel is multiplied by the terrain view factor to get an approximation of the terrain irradiance at the observed pixel. The major drawback of this trigonometric approximation is that it does not consider the spatial distribution of the neighboring terrain elements with their different radiances towards the target pixel. The integration of radiance contribution of different components leads to:

$$E_o\left(b\right)=\underset{direct}{\underbrace{\mathrm{\zeta }{E}_s\left(b\right)\frac{\cos i}{\cos {\theta }_s}{r}_{so}}}+\underset{diffuse}{\underbrace{E^{-}\left(b\right)\left({\tau }_{ss}\frac{\cos i}{\cos {\theta }_s}+\left(1-{\tau }_{ss}\right)V_{sky}\right)r_{do}}}+\underset{terrain}{\underbrace{E_{terr}{r}_{do}}}$$

(15)

This equation describes the reflected radiance by a non-Lambertian, non-homogeneous rugged terrain at the bottom of the atmosphere. Modeling of the non-Lambertian reflectance is described below. Since terrain-related quantities (e.g., sky view, terrain view, slope, and aspect) are spectrally independent, they can be considered as multiplicative constants. Accordingly, we define the terms below:

$$F_{sun}=\mathrm{\zeta }\frac{\cos i}{\cos {\theta }_s};F_{sky}={\tau }_{ss}\frac{\cos i}{\cos {\theta }_s}+\left(1-{\tau }_{ss}\right)V_{sky}$$

Therefore, Equation (5) on sloped terrain turns into;

$$E_o\left(b\right)=F_{sun}{E}_s\left(b\right)r_{so}+F_{sky}{E}^{-}\left(b\right)r_{do}+E_{terr}{r}_{do}$$

(16)

It has to be noted that the spatial resolution of the DEM plays a crucial role in the estimation of the topography effects. To accurately estimate these effects, a DEM of a resolution equal to or higher than the sensor’s spatial resolution is required [44].

3.2. Reflectance Modeling

Topography influences the orientation of surface and hence alters the illumination and viewing angles of the observed pixels. Apart from the bi-hemispherical reflectance ($r_{dd}$), the other three reflectance factors calculated by the SLC determine additional changes in the reflected radiance depending on the illumination and viewing geometry. To account for surface anisotropic reflectance effect, the SLC calculations have to be done with the illumination and viewing angles considering the topography effect on the orientation. This effect is also considered in the estimation of terrain irradiance from one pixel to another, as described above.

It has to be noted that leaves tend to maintain their LIDF with respect to the gravity field, so in principle the BRDF would stay invariant, but in the case of sloped terrain photon paths inside the vegetation layer change compared to those for a flat terrain case, so some effect on the BRDF is expected. However, this has not been considered in this study.

3.3. Top of Atmosphere Radiance over Rugged Terrain

By substituting Equation (16) in Equation (1) and solving the diffuse fluxes at the surface, the TOA radiance over rugged terrain is obtained by

$$L_{TOA}=\frac{{E_{\text{s}}}^o\cos {\theta }_s}{\text{π}} \left[\stackrel{atmosphere}{\overbrace{{\rho }_{so}}}+\stackrel{adjacency}{\overbrace{\frac{{\tau }_{ss}\overline{r_{sd}}+{\tau }_{sd}\overline{r_{dd}}}{1-\overline{r_{dd}}{\rho }_{dd}}{\tau }_{do}}}+\stackrel{sun}{\overbrace{{\tau }_{ss}{r}_{so}{\tau }_{oo}{F}_{sun}}} + \stackrel{sky}{\overbrace{\frac{{\tau }_{sd}+ {\tau }_{ss}\overline{r_{sd}}{\rho }_{dd}}{1-\overline{r_{dd}}{\rho }_{dd}}{r}_{do}{\tau }_{oo}{F}_{sky}}}+\stackrel{terrain}{\overbrace{E_{terr}{r}_{do}{\tau }_{oo}}}\right]$$

(17)

The definitions of different terms are given in Table 2. The spatial filtering kernel of the adjacency effect is defined by $\overline{r_{dd}}= \frac{1}{n^2}{\displaystyle {\displaystyle \sum }_{i=1}^n}{r}_{dd}$, where $n$ is the number of pixels for the area taken for adjacency effect (≈1 km in this study). However, due to the fact that the scattered photons are spread over a much larger region when the deflection angle increases and due to the strong forward scattering of aerosols, the required filtering kernel should be sharply peaked in the center. Moreover, this filtering kernel would give an approximation of the adjacency effect in the case of sloped terrain since the mean reflectance $\overline{r_{dd}}$ will depend on other factors as well, such as the orientation of adjacent terrain facets. Additionally, the atmospheric reflectance ${\rho }_{so}$ in Equation (17) may change according to the local topography, which in that case, it requires pixel-wise modeling of the atmosphere. In this study, we have assumed a constant atmosphere over the study area and hence ${\rho }_{so}$ is assumed to be constant across the study area.

In order to calculate TOA radiance, the unknown parameters in Equation (17) have to be determined. Amongst the parameters in the equation, four surface reflectance factors are calculated for a given surface with known properties (soil, leaf and canopy properties) on a pixel-by-pixel basis from the SLC model (see Section 3.2). In case of rugged terrain, a DEM is applied to derive slope, aspect, sky view factor and terrain view factors to estimate topography-induced effects (see Section 3.1). In fact, for a specified scenario only atmospheric coefficients have to be determined. Six atmospheric coefficients plus extraterrestrial solar irradiance are determined by applying an atmospheric radiative transfer code like MODTRAN [45,46]. In this study MODTRAN4 was used to calculate the atmospheric coefficients. We implemented a modified version of the MODTRAN interrogation technique [13,26] to estimate atmospheric coefficients. We briefly explain MODTRAN and then the modified MIT approach is given in more detail in the next Section.

3.4. Atmospheric Modeling

Modified MODTRAN Interrogation Technique (MIT): MODTRAN (MODerate resolution atmospheric TRANsmittance and radiance code) was initially developed from the LOW resolution TRANsmittance 7 (LOWTRAN 7) radiative transfer model. It models radiation scattering and absorption through the atmosphere. Atmospheric multiple scattering is computed using the scaled DISORT (DIScrete Ordinate Radiative Transfer) algorithm and the correlated k algorithm is used to model the absorption for regions presenting considerable absorption effects. MODTRAN allows flexibility in the input parameters where both user inputs and predefined profiles for different atmospheric scenarios are accepted.

Table 3. Standard MODTRAN outputs with the given headers, column number, formula and relevant descriptions.

MODTRAN (*.tp7 or *.7sc)	Column	Formula	Description
FREQ	1	-	band-center wavelength
TRAN	2	${\tau }_{oo}$	direct transmittance from the ground to the sensor
PTH THRML	3	-	upwelling atmospheric emitted radiance received at the sensor (negligible in the reflective region)
THRML SCT	4	-	thermal scattering of the path radiance
SURF EMIS	5	-	surface emission received at the sensor
SOL SCAT	6	$\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}\left[{\rho }_{so}+\frac{({\tau }_{ss}+{\tau }_{sd})}{1-a{\rho }_{dd}}a{\tau }_{do}\right]$	solar multi-scattered radiance received at the sensor(PATH)
SING SCAT	7	-	solar single-scattered radiance
GRND RFLT	8	$\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}\left[\frac{{\tau }_{sd} + {\tau }_{ss}}{1-a{\rho }_{dd}}a{\tau }_{oo}\right]$	ground reflected radiance(GRT)
DRCT RFLT	9	$\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}\left[{\tau }_{ss}a{\tau }_{oo}\right]$	direct solar reflected radiance received at the sensor(GSUN)
TOTAL RAD	10	TOTAL RAD = SOL SCAT + GRND RFLT	PTH THRML + SURF EMIS + SOL SCAT + GRND RFLT (GTOT)
REF SOL	11	${E_{\text{s}}}^o{\tau }_{ss}{\tau }_{oo}$	SOL@OBS multiplied by ${\tau }_{ss}{\tau }_{oo}$
SOL@OBS	12	${E_{\text{s}}}^o$	extraterrestrial solar irradiance
DEPTH	13	$b=-\ln ({\tau }_{oo})$	negative Log of the transmittance column TRAN (total optical thickness)

Table 3. Standard MODTRAN outputs with the given headers, column number, formula and relevant descriptions.

MODTRAN (*.tp7 or *.7sc)	Column	Formula	Description
FREQ	1	-	band-center wavelength
TRAN	2	${\tau }_{oo}$	direct transmittance from the ground to the sensor
PTH THRML	3	-	upwelling atmospheric emitted radiance received at the sensor (negligible in the reflective region)
THRML SCT	4	-	thermal scattering of the path radiance
SURF EMIS	5	-	surface emission received at the sensor
SOL SCAT	6	$\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}\left[{\rho }_{so}+\frac{({\tau }_{ss}+{\tau }_{sd})}{1-a{\rho }_{dd}}a{\tau }_{do}\right]$	solar multi-scattered radiance received at the sensor(PATH)
SING SCAT	7	-	solar single-scattered radiance
GRND RFLT	8	$\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}\left[\frac{{\tau }_{sd} + {\tau }_{ss}}{1-a{\rho }_{dd}}a{\tau }_{oo}\right]$	ground reflected radiance(GRT)
DRCT RFLT	9	$\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}\left[{\tau }_{ss}a{\tau }_{oo}\right]$	direct solar reflected radiance received at the sensor(GSUN)
TOTAL RAD	10	TOTAL RAD = SOL SCAT + GRND RFLT	PTH THRML + SURF EMIS + SOL SCAT + GRND RFLT (GTOT)
REF SOL	11	${E_{\text{s}}}^o{\tau }_{ss}{\tau }_{oo}$	SOL@OBS multiplied by ${\tau }_{ss}{\tau }_{oo}$
SOL@OBS	12	${E_{\text{s}}}^o$	extraterrestrial solar irradiance
DEPTH	13	$b=-\ln ({\tau }_{oo})$	negative Log of the transmittance column TRAN (total optical thickness)

MODTRAN has been successfully used for both at sensor radiance simulation [11,13,14,26,27,47] and for atmospheric correction [44,48] from the visible to the thermal infrared wavelength region. The computational cost of a MODTRAN run is dependent on the calculation mode, as well as whether using DISORT with various streams. MODTRAN takes an ASCII file (with *.tp5 suffix) as input and delivers some standard output files, while the files with suffix “*.tp7” and “*.7sc” are of interest for use in the simulation procedure. The standard MODTRAN outputs do not directly provide the unknown atmospheric coefficients required in the four streams RT approach and they have to be computed from such outputs. Table 3 lists the standard MODTRAN4 outputs along with formulas and explanation of each column.

Verhoef and Bach [26,27] have shown that the six atmospheric coefficients in Equation (17) can be derived from three MODTRAN runs for surface albedos at 0, 0.5 and 1 using the so-called MODTRAN Interrogation Technique (MIT). These runs are done assuming uniform Lambertian surface reflectance, given atmospheric state and angular configuration. We have modified the algebra used by Verhoef and Bach [26,27] to improve computational efficiency. In the modified MIT, only two MODTRAN runs are required to extract the atmospheric coefficients. The two runs can be with any spectrally flat surface albedos between 0 and 1. It is recommended to avoid using zero albedo since it leads to some miscalculations in the extracted parameters from MODTRAN [33]. For instance, with a black surface assumption, the scaling in the scaled-DISORT option by MODTRAN is not applied over the diffuse fluxes, which results in miscalculation of the spherical albedo parameter. We have used two different spectrally flat surface albedos, i.e., 0.5 and 1. The modified MIT needs the solar-scattered path radiance (SOL SCAT in Table 3) and the ground-reflected radiance contribution (GRND RFLT in Table 3) for surface albedos at 0.5 and 1 plus the radiance contribution due to ground-reflected sunlight (DRCT RFLT in Table 3) for surface albedo at 1.

If we assume a uniform Lambertian surface reflectance with albedo $a$ (in this case $a=r_{so}=r_{do}=\overline{r_{sd}}=\overline{r_{dd}}$), the TOA radiance signal in Equation (17) turns into the simplified equation

$$L_{TOA}=\underset{PATH}{\underbrace{\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}({\rho }_{so}+\frac{\left({\tau }_{ss} + {\tau }_{sd}\right)a}{1-a{\rho }_{dd}}{\tau }_{do})}}+\underset{GSUN}{\underbrace{\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}{\tau }_{ss}a{\tau }_{oo}}}+\underset{GSKY}{\underbrace{\frac{{E_{\text{s}}}^o{\mu }_{\text{s}}}{\text{π}}\frac{{\tau }_{sd} + {\tau }_{ss}a{\rho }_{dd}}{1-a{\rho }_{dd}}a{\tau }_{oo}}}$$

(18)

Here, GSKY stands for the radiance contribution due to sky irradiance in the direction of viewing. We first extract ${\tau }_{ss}$ by substituting ${\tau }_{oo}$ and ${E_{\text{s}}}^o$ (direct outputs from TRAN and SOL@OBS in Table 3) into GSUN for surface albedo at 1,

$${\tau }_{ss}=\frac{GSU{N}_1\times \text{π}}{{E_{\text{s}}}^o{\mu }_{\text{s}}\times {\tau }_{oo}}$$

(19)

It should be noted that, as shown in Verhoef and Bach [27], some equations (in Table 3) are inaccurate in absorption bands of the atmosphere due to finite bandwidth. For accurate applications, one needs to make two runs at TOA and two at BOA, especially if an accurate surface irradiance is needed. Next, the spherical albedo ${\rho }_{dd}$ is computed from the combination of two GRT for surface albedos at 0.5 and 1,

$${\rho }_{dd}=\frac{GR{T}_1-2\times GR{T}_{0.5}}{GR{T}_1-GR{T}_{0.5}}$$

(20)

Once ${\rho }_{dd}$ and ${\tau }_{ss}$ are estimated, ${\tau }_{sd}$ can be derived from

$${\tau }_{sd}=\left[\frac{GR{T}_1\left(1-{\rho }_{dd}\right)}{GSU{N}_1}-1\right]{\tau }_{ss}$$

(21)

The diffuse atmospheric transmittance from the ground to the sensor ${\tau }_{do}$ can also be derived from the combination of two sets of PATH and GRT for surface albedos at 0.5 and 1

$${\tau }_{do}=\left[\frac{PAT{H}_1-PAT{H}_{0.5}}{GR{T}_1-GR{T}_{0.5}}\right]{\tau }_{oo}$$

(22)

Finally, the atmospheric bidirectional reflectance ${\rho }_{so}$ is derived using PATH and GRT for surface albedo at 1,

$${\rho }_{so}=\left[PAT{H}_1-(\frac{GR{T}_1*{\tau }_{do}}{{\tau }_{oo}})\right]\frac{\text{π}}{{E_{\text{s}}}^o{\mu }_{\text{s}}}$$

(23)

Comparison between the modified MIT versus the original MIT [26,27] revealed that both methods provide very similar results with negligible differences. The Root Mean Squared Difference (RMSD) between the two methods, across the entire spectral range, was 3.4 × 10⁻⁵ and 3.9 × 10⁻⁵ in the estimation of ${\rho }_{so}$ and ${\tau }_{do}$, respectively.

Therefore, for a given atmosphere and a known sun–target–sensor geometry, all the six atmospheric coefficients are derived in this way. These atmospheric coefficients are spectrally dependent and as a common practice they can be assumed constant for an image due to the homogeneity of the atmosphere over small areas. It is more realistic, however, to model the atmospheric effects on a pixel wise manner (e.g., in a variable retrieval procedure) because of the high spatial variability of water vapor and aerosols. In that case, instead of direct calls to the atmospheric model, an alternative is to build up a large look up table where all the expected input combinations are covered. The computational cost of creating such look up tables is high and can be reduced using the modified MIT method.

4. Model Implementation

TOA radiance over rugged terrain can be modeled using Equation (17). The unknowns in this equation are defined in Table 2 and we explain here how to determine such unknowns in this section. First, the topography related factors are derived from the DEM including the slope, aspect, sky view factor and the list of pixels that are seen from each pixel. Second, sun and observer zenith and azimuth angles are calculated per pixel using the slope and aspect maps and the binary factor $\mathrm{\zeta }$ is estimated as well. Third, for a given surface with known soil, leaf and canopy properties, the four canopy directional reflectance factors; $r_{so}$, $r_{do}$, $r_{sd}$, $r_{dd}$ are calculated using the SLC per pixel considering the local sun and observer zenith angles and the relative azimuth angle. The spatial filtering kernel of the adjacency effect is then applied to the directional-hemispherical and the bi-hemispherical reflectance factors to obtain $\overline{r_{sd}}$ and $\overline{r_{dd}}$ in Equation (17) considering a radius of 1 km around the target pixel. Fourth, the atmospheric coefficients are derived from two MODTRAN4 runs using the modified MIT. The relevant equations (Equation (19) to Equation (23)) to derive such coefficients are given in Section 3.4. The next step is to estimate the factors $F_{sun}$ and $F_{sky}$ as well as $E_{terr}$ using the equations in Section 3.1. In this study, we calculate all the wavelength dependent factors at intervals of 1 nm and the calculated spectral radiance is later convolved with the sensor spectral response functions to generate sensor-specific radiometric data.

5. Design of the Numerical Experiment

To illustrate the effect of topography on the surface reflectance and TOA radiance, we have compared two scenes generated with SLC: the first one for the actual, flat case, site of Barrax and the second one generated using exactly the same spatial patterns of vegetation, soil and atmospheric properties, but laid over the DEM of a mountainous area. The data set is synthetic in the sense that an actual DEM of one study area (Alpine site in Austria) was combined with ground measurements and a validated radiative transfer model for a different (flat) area (Barrax). The advantage of a synthetic dataset, in this case, is that the uncertainty sources are known and it allows for investigating cases that cannot be done with actual data. The observations and ground truth data collected during the ESA SEN3EXP field campaign were used for this study. The field campaign was carried out at the Barrax (39°3′N, 2°60′W) test site and surrounding area from 20 to 24 June 2009 as part of the development process for ESA’s Sentinel-3 Earth observation mission [49]. Barrax is an agricultural test site located in southeastern Spain on the La Mancha plateau characterized by relatively heterogeneous landscape with large uniform irrigated agricultural fields and rain-fed crops. The area consists of irrigated fields, harvested fields and bare soil, which provide an appropriate environment for simulating heterogeneous anisotropic surface behavior. The DEM was taken from a mountainous region in the Austrian Alps, near Salzburg (47°41′N, 13°01′W). The elevation varies between 457 and 1526 m with an average slope of 24° including both north-facing and south-facing hill slopes. Figure 5 (top row) shows the slope and aspect maps of the rugged surface used for this study. The aspect angle is referenced from the north.

Figure 5. Upper row: The slope (left) and aspect (right) maps used for generating the rugged terrain surface. Lower row: RGB (870, 670 and 540 nm) band combination of the simulated CHRIS-Proba image flat case (left) and over rugged terrain (right) using Equations (12) and (13) for the calculation of sky view factor and terrain reflected radiance, respectively. The letters represent two pixels: “A” (a pixel exposed to the sun) and “B” (a self-shadowed pixel) (see Section 6.1 for explanation).

Figure 5. Upper row: The slope (left) and aspect (right) maps used for generating the rugged terrain surface. Lower row: RGB (870, 670 and 540 nm) band combination of the simulated CHRIS-Proba image flat case (left) and over rugged terrain (right) using Equations (12) and (13) for the calculation of sky view factor and terrain reflected radiance, respectively. The letters represent two pixels: “A” (a pixel exposed to the sun) and “B” (a self-shadowed pixel) (see Section 6.1 for explanation).

This experimental design allows for a correct evaluation of the topographic effects, since we can compare the TOA radiance of the flat case with the rugged case for the same land-cover and the same atmosphere. The spatial information on soil, leaf and canopy properties were derived from an input TOA radiance CHRIS-Proba image of 19 June 2009 over the Barrax area (flat case) as explained in [14]. The surface reflectance factors were calculated from the soil, leaf and canopy properties using the SLC model considering the local illumination and view directions (see Section 3.2). The atmospheric coefficients were derived by applying the modified MIT using two MODTRAN runs. Next, having taken into account the topographic effects on the irradiance, the TOA radiance is calculated for the entire synthetic scene (rugged case). The effects of sky view factor, terrain view factor and solar illumination are then evaluated by comparing the actual scene for the flat case with the synthetic scene for the rugged case.

6. Results and Discussion

6.1. Topography Effects on the TOA Radiance

Figure 5 (lower row) shows the simulated TOA radiance CHRIS-Proba image (flat case, lower left) and the simulated TOA radiance image over the rugged terrain (rugged case, lower right). Visually comparing the images, the effect of higher irradiance on south facing facets is very clear, as well as the lower irradiance on north facing facets. For instance, Figure 6 (left) presents a rugged-case self-shadowed pixel (labeled “B” in Figure 5) with slope of 42° and aspect of 319° compared with the same flat-case pixel, while the right hand figure presents a pixel exposed to the sun (labeled “A” in Figure 5) with slope of 20° and aspect of 111°.

Figure 6. Comparison of the simulated TOA radiances over rugged terrain versus flat-case image for a self-shadowed pixel (left) and a pixel exposed to the sun (right).

Figure 6. Comparison of the simulated TOA radiances over rugged terrain versus flat-case image for a self-shadowed pixel (left) and a pixel exposed to the sun (right).

In the left-hand figure, the self-shadowed pixel (blue) receives most irradiance from diffuse radiation causing lower reflected radiance, while the right-hand figure indicates the extra irradiance coming into the pixel due to the sun position and the topographic effects of neighboring areas. The full flat- and rugged-case scenes were compared and the result is illustrated by Figure 7, where the RMSD values between the simulated flat- and rugged-case scenes are shown versus the slope and aspect angles. A positive correlation is seen between slope, aspect and RMSD. The smaller RMSD are attributed to areas with small slope and aspect characterizing flat areas and the RMSD increase with an increase in the slope and aspect. A considerable difference in RMSD is seen over north facing slopes even for intermediate slopes.

The illustrations above document the impact of topography on TOA radiance. In a previous study [14], we have documented the accuracy of the forward radiative transfer model for the Barrax site. Particularly, we have shown that the radiative transfer model calculates accurate TOA radiance by comparison with actual CHRIS PROBA image data at different view zenith angles. According to the principle of symmetry, this is equivalent to show that the model is accurate for varying sun zenith angle at a fixed view zenith angle, which applies to the observation of rugged terrain.

Apart from the general impact of topography, it is also interesting to check how sky view and terrain reflected radiance influence the TOA radiance. In the following sections we quantify the effect of these factors on the TOA radiance.

Figure 7. RMSD values between the simulated and actual images versus slope and aspect angles of the study area.

Figure 7. RMSD values between the simulated and actual images versus slope and aspect angles of the study area.

6.2. Sky View Factor Effect

In order to investigate the impact of sky view factor on TOA radiance, two different methods were used to calculate the sky view factor. Figure 8 (left and middle) displays the sky view factor calculated using the horizon (i.e., Equation (12)) and the approximate (i.e., Equation (10)) methods. The mean difference between the sky view factors calculated with the two methods is around 0.04. Considerable differences are seen (Figure 8 right), however, in the valleys and steep slopes (up to 0.30 difference). The approximate method generally overestimates the sky view factor due to the neglected effect of the neighboring areas.

Figure 8. Sky view factor calculated using the horizon (left) and approximate (middle) methods and the difference (right) between the two methods.

Figure 8. Sky view factor calculated using the horizon (left) and approximate (middle) methods and the difference (right) between the two methods.

We conducted two experiments to study the effect of sky view factor on the output radiance. We first simulated TOA radiance without applying sky view factor and second, we compared simulated TOA radiances applying two different sky view factors. All the other effects were kept identical for both cases and only in the second case the method to calculate the sky view factor from the DEM was changed.

Figure 9 (left) illustrates the RMSD of all the CHRIS-Proba bands between the simulated radiance with and without applying sky view factor. RMSD values range from zero to 5 (mW/m²·sr·nm), though the increasing discrepancy is attributed to the steep and north-facing slopes. This quantifies the effect of sky view factor on the output radiance and shows that neglecting this factor can cause considerable uncertainty in forward TOA radiance modeling in particular for steep slopes.

Figure 9. RMSD between simulated TOA radiance with and without applying sky view factor (left) and RMSD between simulated TOA radiance applying two different sky view factors (right).

Figure 9. RMSD between simulated TOA radiance with and without applying sky view factor (left) and RMSD between simulated TOA radiance applying two different sky view factors (right).

The RMSD between simulated TOA radiance using different sky view factors is shown in Figure 9 (right). A maximum RMSD of 0.8 (mW/m²·sr·nm) is seen in this case over north facing and steep slopes. Although, the RMSD are small in this case compared to the case when neglecting sky view factor, it still introduces uncertainty in the forward TOA radiance modeling. These results suggest that a precise estimation of the sky view factor is therefore essential to reduce the uncertainty in the TOA radiance modeling over rugged terrain. Therefore, the horizon method is recommended for the calculation of the sky view factor because of its fast calculation and also more accurate results than the approximate method.

6.3. Terrain Reflected Radiance Effect

Terrain reflected radiance from adjacent slopes has been reported to have non-negligible impact on the TOA radiance in rugged terrain. To quantify the influence of this factor, we have calculated this contribution using two different approaches, the trigonometric approximation (i.e., Equation (14)) and the exact calculation of the radiance reflected by adjacent pixels (Equation (13)). We then compared the simulated TOA radiance without applying the terrain reflected radiance (i.e., $V_{terrain}=0$) against the radiance simulated using the two different terrain calculated radiances (from the trigonometric and exact methods). The simulation using the trigonometric approximation was done with the assumption of Lambertian reflection from the adjacent pixels, while in the case of the exact method the anisotropic behavior of the surroundings is taken into account. This will furthermore quantify the impact of non-Lambertian reflection from the neighboring pixels on TOA radiance.

Figure 10 shows the RMSD values between the TOA radiance simulations without considering terrain reflected radiance (i.e., $V_{sky}$ from the horizon method, $V_{terrain}=0$) versus the simulations applying this effect calculated from the trigonometric approximation (left) and calculated from the exact method (right). RMSD values are larger for the trigonometric approximation indicating an overestimation of the terrain reflected radiance effect with this approach. On the other hand, the exact method presents very small RMSD over the entire area except the valleys with steep slope. It has been already reported in the literature that this effect is significant for snow covered areas, deep valleys, low sun elevation images and pixels in shadow and can be neglected for other cases due to its intensive computational cost and low impact on the output radiance [38,42,50]. The differences between the two RMSD values in Figure 10 can also be related to the assumption on Lambertian reflection from the neighboring pixels, but this assumption seems to have a secondary impact since there is a high correlation between the steep slopes and the RMSD values.

Figure 10. RMSD between simulated TOA radiance without applying terrain reflected radiance versus simulated radiance applying terrain reflected radiance calculated with the approximate method (left) and the exact method (right).

Figure 10. RMSD between simulated TOA radiance without applying terrain reflected radiance versus simulated radiance applying terrain reflected radiance calculated with the approximate method (left) and the exact method (right).

The calculation of terrain reflected radiance using the exact method is computationally expensive, since the analysis has to be iterated over all the pixels. In this case study, the computation time for the calculation of the terrain reflected radiance using the exact method (the case in Figure 5 (lower right)) was almost 15 h, while the trigonometric approximation took less than a minute to calculate this effect. This computational cost makes the use of the exact method very limited; in particular, it becomes impractical for variable retrieval applications. On the other hand, the trigonometric approximation considers only the slope, which cannot be representative of the amount of the terrain illuminating a given pixel.

We have further carried out an experiment to quantify how significant may be the impact of neglecting topographic effects in terms of variable retrieval. We simulated TOA radiance for different canopy realizations and compared the RMSD between different cases to those found due to topographic effects. This provides an indication as to whether the topographic effects might be neglected without major impacts on retrieval. The TOA radiances (Figure 11) were simulated for different cases varying two important canopy variables: LAI and fractional vegetation cover (fCover) for a corn pixel representative of the study area. Three LAI values, i.e., 2, 2.5 and 3 (m²/m²), and three fCover, i.e., 0.58, 0.65 and 0.72 were used for simulations. The values were chosen according to the actual ground measurements of a representative corn pixel directly illuminated in the ESA SEN3EXP field campaign. The leaf chlorophyll content (40 µg/cm²), leaf dry matter content (0.005 g/cm²), leaf water thickness (0.019 cm), leaf internal structure variable (1.6) and soil moisture (0.2) were set constant for all simulations (see [14] for further details on actual field conditions). The sun zenith angle was 30° and spherical leaf inclination distribution was assumed.

Figure 11. TOA CHRIS-Proba radiance simulation using different LAI values (left) and fCover values (right) for a representative corn pixel.

Figure 11. TOA CHRIS-Proba radiance simulation using different LAI values (left) and fCover values (right) for a representative corn pixel.

The RMSD values between different TOA radiance simulations are given in Table 4. A change of 0.5 in LAI results in RMSD = 4.4 or 3.9 (mW/m²·sr·nm) depending on the LAI value. A change of 0.07 in fCover also causes the RMSD of 3.2 or 2.7 (mW/m²·sr·nm) depending on the fCover value. The RMSD are significantly higher for larger changes as shown in the table.

Table 4. RMSD values between different TOA radiance simulations.

Between Simulations with	RMSD (mW/m2·sr·nm)
LAI of 2 and LAI of 2.5	4.4
LAI of 2.5 and LAI of 3	3.9
LAI of 2 and LAI of 3	8.2
fCover of 0.58 and fCover of 0.65	3.2
fCover of 0.65 and fCover of 0.72	2.7
fCover of 0.58 and fCover of 0.72	9

Table 4. RMSD values between different TOA radiance simulations.

Between Simulations with	RMSD (mW/m2·sr·nm)
LAI of 2 and LAI of 2.5	4.4
LAI of 2.5 and LAI of 3	3.9
LAI of 2 and LAI of 3	8.2
fCover of 0.58 and fCover of 0.65	3.2
fCover of 0.65 and fCover of 0.72	2.7
fCover of 0.58 and fCover of 0.72	9

Referring to Figure 9, for example, if the sky view factor is not taken into account over rugged terrain, it would yield an uncertainty up to 5 (mW/m²·sr·nm) in the simulated TOA radiance. Accordingly, we can interpret Table 4 as indicating that neglecting this effect can be translated into a change of 0.5 (20% uncertainty) in LAI or similarly 0.07 (10% uncertainty) in fCover in the retrieval of these variables. This means that without considering this effect, provided that all the other effects (atmospheric, radiometric, etc.) are considered, an over/underestimation of more than 0.5 in LAI is expected in the retrieval procedure. To understand the importance of the uncertainty caused by neglecting topography effects on the retrieval, it is interesting to know that, for instance, the Global Monitoring for Environment and Security (GMES) user committee has defined an accuracy goal of 10% in the retrieval of vegetation variables to guarantee that the user requirements are met for different applications [51]. To this end, topography effects have to be taken into account when retrieving vegetation properties over rugged terrain. The extended four-stream approach presented in this study can be used to account for such effects in the forward simulation of TOA radiance.

7. Conclusions

In this study, we aimed at investigating the topographic effect on TOA radiance and variable retrieval by quantifying this effect in terms of considering different topographic components such as sky view factor and terrain reflected radiance. To do so, we presented an extension of “the four-stream radiative transfer theory” previously proposed by Verhoef and Bach [13,26,27]. The extension is mainly focused on two aspects in TOA radiance modeling. First, the study gives a detailed account of the main topography-induced effects on TOA radiance simulation and the relevant equations are given. Second, a new formulation is proposed for the derivation of the six atmospheric coefficients using MODTRAN in order to make it more computationally efficient, as well as to avoid using zero surface albedo, which causes some miscalculations in the results.

To illustrate the topography effects on the TOA radiance, we compared two scenes generated with SLC: the first one for the actual, flat case, site of Barrax and the second one generated using exactly the same spatial patterns of vegetation properties, but laid over the DEM of a mountain area. Surface reflectance factors were calculated taking into account the impact of topography on the sun–target–sensor geometry. The simulated radiances for shadowed pixels and those exposed to the sun showed considerable differences. It was shown that neglecting sky view factor causes uncertainty in the forward TOA radiance modeling up to 5 (mW/m²·sr·nm) specifically over steep slopes. Two different methods (approximate and horizon methods) were applied to calculate sky view factor and it was conclude that the horizon method delivers better results for forward TOA radiance modeling. The effect of reflected terrain irradiance was also modeled using two different approaches (approximate and exact methods), where the approximate method seems to overestimate this effect. Although the exact method calculates the terrain reflected radiance more accurately, the computational cost of this method limits its use for applications such as variable retrieval. Furthermore, it was shown that, for instance, RMSD of <5 (mW/m²·sr·nm) due to neglecting topographic effects can cause an over/underestimation of more than 0.5 in LAI (or 0.07 in fCover) in variable retrieval.

This contribution demonstrated that the topographic effect has to be taken into account in variable retrieval. The approach presented here can be used to model TOA radiance over rugged terrain for a range of applications. For instance, retrieving surface properties from a multi-sensor approach using TOA radiance data would require such an approach to accurately simulate the radiance over mountainous areas.