Ground Reﬂectance Retrieval on Horizontal and Inclined Terrains Using the Software Package REFLECT

: This paper presents the software package REFLECT for the retrieval of ground reﬂectance from high and very-high resolution multispectral satellite images. The computation of atmospheric parameters is based on the 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) routines. Aerosol optical properties are calculated using the OPAC (Optical Properties of Aerosols and Clouds) model, while aerosol optical depth is estimated using the dark target method. A new approach is proposed for adjacency effect correction. Topographic effects were also taken into account, and a new model was developed for forest canopies. Validation has shown that ground reﬂectance estimation with REFLECT is performed with an accuracy of approximately ± 0.01 in reﬂectance units (for the visible, near-infrared, and mid-infrared spectral bands), even for surfaces with varying topography. The validation of the software was performed through many tests. These tests involve the correction of the effects that are associated with sensor calibration, irradiance, and viewing conditions, atmospheric conditions (aerosol optical depth AOD and water vapour)


Introduction
Multispectral satellite imagery, especially at high spatial resolution (30 m or finer), represents an invaluable source of information for decision making [1] in various domains in connection with natural resources management, environment preservation, or urban planning and management.Common applications of remotely sensed imagery concern vegetation cover (species, biomass, damages, crop yields, etc.), surface conditions (land use, types of materials, soil physics, and chemistry, etc.) or the inventory of natural resources [2].This involves the transition from a relative scale (gray scale) of measured physical variables (radiances) to a scale of biophysical variables (biomass, leaf cover, etc.) or a system of classes (e.g., type of land use).This transition is generally accomplished by combining a series of models, algorithms, equations, and/or classification methods.Depending on the spatial, spectral, and radiometric resolution of the sensors, it is possible to identify various classes and conditions of objects.The basic assumption is that the signal recorded by the sensor is directly related to the state and/or nature of the target objects on the ground.
Variation in object reflectance in the solar spectrum is the key parameter for many optical remote sensing applications.However, in addition to object reflectance values, sensor measurements are affected by many factors.Some of them are interrelated, and can be summarised as follows: (1) the atmosphere: composition and optical properties of gases and aerosols; (2) the sun: its position in the sky during data acquisition and its spectral irradiance; (3) the terrain: elevation and topography; (4) ground reflectivity: reflectance anisotropy and reflectivity of the surrounding area; and (5) the sensor: calibration, spectral sensitivity, viewing angle, and altitude.In this article, the REFLECT software package is presented, which considers all of these aspects for the radiometric correction of multispectral satellite images and ground reflectance retrieval.
There are many methods that can be used for ground reflectance retrieval [3].These methods could be grouped into two main approaches.The first is empirical, and the second is based on the physical modelling of the signal.Empirical methods require a series of recognisable targets on the images of invariable and known reflectance.Therefore, it is therefore by regression to establish the relationship between the signal measured by the sensor and the ground reflectance.Using this relationship, the reflectance of any pixel on the images is then estimated [4][5][6].The underlying hypothesis here is that the coefficients of the equations, which link the numerical values to the reflectance, are constants for the entire image.However, these coefficients are variable over time, and at best can be used only for one date.The physical modelling approach is by far the most widely used.The atmospheric effects (additive and multiplicative) and ground irradiances (direct and diffuse) are estimated by applying approximate solutions of the radiative transfer equation in the earth-atmosphere system.These solutions are incorporated in atmospheric codes simulating the desired parameters as a function of a series of inputs.
The purpose of this paper is to present the theoretical and practical bases of the REFLECT software package for ground reflectance retrieval based on physical modelling (6S atmospheric code).The first version of this software was developed in the remote sensing laboratory of the Department of Geography at the University of Montreal in order to estimate the reflectance of water in coastal environments based on Landsat-7 ETM+ images [7].A second version was developed in collaboration with Agriculture and Agri-Food Canada's Horticulture Research and Development Centre (HRDC), St-Jean-sur-Richelieu, Quebec.The aim of calculating the ground reflectance of crops is based on Landsat-7 ETM+ images [8,9].The current version that is presented in this paper includes several new features concerning: (a) terrain with variable topography; (b) adjacency effects; (c) atmospheric conditions; and (d) observation conditions and sensor properties.In order to better present the theoretical and practical bases of REFLECT, this article is subdivided into eight sections.The "Material and Methods" are described in Sections 2-4, including a brief overview of the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) code, basic elements of REFLECT, and the methodology that was used for validation.The "Results and Discussions" are presented in Sections 5-7, where validation tests of the ground reflectance retrieval are shown.These tests involve the correction of effects that are associated with the sensor calibration, irradiance, and viewing conditions, atmospheric conditions (aerosol optical depth AOD and water vapour), adjacency, and topographic conditions.Section 8 concludes the paper.

The 6S Atmospheric Code
The version of the 6S code that was used in REFLECT is the one published in 2006 by Vermote, E.F. et al. [10].Conversely, 6S can estimate the reflectance of a ground object based on spectral radiance measured by a sensor under atmospheric conditions that are also defined by the user.The 6S code assumes a cloud-free atmosphere stratified in parallel layers, in which each one is characterised by its content in gas molecules and aerosol particles.The atmospheric constituents vary in space (water vapour and aerosols).A way of addressing this variability is proposed in REFLECT.
The multi-layer structure of the atmosphere in the 6S code divides the atmospheric column into 34 levels between the ground and the top of the atmosphere, allowing the elevation of the target and the altitude of the sensor to be considered.Thus, the code accounts for the reduction in the number of scattering particles (molecules and aerosols) and the absorbing gases that are above a target at a certain altitude.The atmospheric pressure and temperature profiles are used to adjust the optical depths of the atmospheric constituents, which are more concentrated in the lower layers.Standard models for the vertical distribution of atmospheric gases (US62, MIDSUM, etc.) are used therein.
Gaseous absorption is computed using the Goody model for H 2 O and the Malkmus model for O 2 , CO 2 , O 3 , N 2 O, and CH 4 .Scattering by atmospheric gases is modelled by the Rayleigh theory, while aerosol scattering is described by the Lorenz-Mie theory [11].The single scattering hypothesis is valid solely if the optical path is very thin (τ << 1) and the scattering albedo is very low (SC << 1).Since these conditions do not always hold [12], multiple scattering is considered using the successive-orders-of-scattering method as an approximation of the radiative transfer equation solution [13].
The 6S code assumes that the atmosphere is limited by a horizontal ground surface.Three cases are then accounted for: (a) the target object is part of a Lambertian surface that extends to infinity; (b) the target object is part of an anisotropic surface (bidirectional properties) that extends to infinity; and (c) the target object occupies the central portion of a regular surface (a circle) of limited extent with Lambertian properties, and is surrounded by a Lambertian surface that extends to infinity with a level of reflectance that is different from that of the target object.Among these three cases, the latter one most closely approaches the conditions that are usually observed in remote sensing.This configuration was selected for REFLECT.However, significant modifications were introduced, as described in Section 3.2.
The code formulates the problem of radiative transfer in terms of reflectance factors, and not in terms of radiance.All of the radiances that are generated in the Earth-atmosphere system are expressed relative to the Lambertian radiance (L sat Lamb ) of a target with unit reflectance, which was illuminated and observed under the same conditions as the target, and placed at the top of the atmosphere.This radiance can be expressed as follows (see symbols in Abbreviations): Hence, the radiance L sat , which is measured by the sensor at a spectral band designated by its central wavelength λ, is described in terms of apparent reflectance factor ρ sat by: For the sake of readability, the references to wavelength and the viewing and illumination geometry are sometimes omitted in the equations.
Although atmospheric scattering may be considered independent from gaseous absorption, water vapour presents a problem, since it is often concentrated in the lower atmospheric layers where the aerosols are also concentrated.To take this into account, the code proposes the following formulation (the case of a Lambertian surface): The atmospheric reflectance is divided into two components that were introduced by the gas molecules (first term on the right) and aerosols (second term on the right).The index OG means "Other Gases" than water vapour.Equation (3) can lead to three different formulations (index j = 1, 2, 3).The first corresponds to an extreme situation where water vapour has a minimal effect on aerosol scattering, since it is concentrated at a layer located below the aerosol layer (j = 1).In this case, we find the formulation usually cited in the literature: The formulation corresponding to j = 3 in Equation ( 3) also describes an extreme case where the water vapour has a maximal effect, since it is mainly concentrated at a layer above the aerosol layer.Finally, the case j = 2 assumes that half of the water vapour is concentrated at the same layer as the aerosols, and therefore, the water vapour has an effect that is halfway between the two extreme cases.
In Equations ( 3) and ( 4), the term T ↓ (θ s )t ↑ dir (θ v )ρ tar contains the useful signal ρ tar (L tar in Figure 1) that must be estimated from satellite measurements by eliminating the multiplicative effects of the atmosphere (total transmittance in the descending path and direct transmittance in the ascending path).The term T ↓ (θ s )t ↑ di f f (θ v )ρ env is due to the radiance reflected by the environment surrounding the observed target (adjacency effect) that was introduced into the field of view of the sensor by atmospheric scattering (L env in Figure 1).This radiance is also considered as a noise (additive), since it does not contain information about the target viewed by the sensor; it is significant for low-reflectance targets.The term in the denominator is added in order to account for multiple reflections between the ground and the atmosphere, which is characterised by its spherical albedo, S alb .This component can be significant in case of a very hazy atmosphere and a highly reflective environment.

REEFLECT: Principles and Main Developments
This section presents the basic elements of the proposed software package: atmospheric properties, surface properties, topographic terrain with variable topography, adjacency effect, observation conditions, and sensor properties.

Atmospheric Properties
This section shows the modifications made to the atmospheric properties used in radiative transfer calculations of the 6S code for both gaseous absorption and aerosol scattering.

Gaseous Absorption
Except for water vapour and ozone, the content of gases in the atmosphere is fixed and known.It is also possible to better represent actual atmospheric conditions by using meteorological data for the location and date of acquisition [14].Therefore we have introduced this option in REFLECT by incorporating data on ozone content variation according to latitude and month of the year, as well as the formula of Leckner [14] (Equations (5a) and (5b)) for water vapour content estimation [15].
where w is the total water vapour content along the length of the atmospheric column in g cm −2 , H r is the relative humidity as a fraction of 1, T a is the ambient temperature at ground level in Kelvin, and P s is the partial pressure of water vapour in saturated air given by the equation:

Aerosol Scattering
The knowledge of aerosol properties is essential in remote sensing, since their effect on solar radiation is present throughout the optical spectrum.These properties are calculated using the Mie theory based on the microphysical properties of the particles [16].However, the chemical composition of aerosols is quite variable, and depends both on the geographic distribution of the sources and on atmospheric dynamics.Determining the proportion of the various types of aerosols at a given location and at a particular time is therefore not easy.In most cases, standard models are used that assume the presence of a mixture of typical aerosol particles (dust, water-soluble, soot, etc.) with known and stable properties.As used in 6S, typical atmospheres (continental, urban, maritime, etc.) are constructed with different mixtures of these particles.
The atmosphere includes a type of aerosol that is called water-soluble, whose particle size increases when they are in contact with water vapour molecules.This contributes to modifying the refractive index and the effective section of the particles, and hence all of the optical properties of this type of aerosol, including the AOD [17][18][19].Therefore, for the atmosphere types that are rich in water-soluble aerosols, as continental and urban mixtures that are found in the regions most frequently targeted by remote sensing, the AOD is expected to be higher for high values of air relative humidity.
We examined the relationship between the AOD measured by the Microtops II Sun photometer (Solar Light Company, Philadelphia, PA, USA.) at two experimental sites (Saint-Jean-sur-Richelieu and St-Valentin, Quebec, Canada), and the corresponding relative humidity measured at the closest weather station (L'Acadie, Quebec, Canada). Figure 2 shows that there is a proportional relationship (r = 0.62) between the AOD at 550 nm and the relative humidity.We can assume that this relationship is due to the proportion of water-soluble aerosols that are present in the continental atmosphere whose optical depth increases with humidity.Thus, we updated the aerosol model of the 6S code (published in 1983 by the World Meteorological Organization) with a more recent model (OPAC: Optical Properties of Aerosols and Clouds), which was presented by Hess, M. et al. [18], and takes into account the effect of relative humidity on the properties of aerosols and thus on the values of the AOD.

The OPAC model
The OPAC model [18] describes the optical and microphysical properties of 10 aerosols that are considered the most typical, as well as their mixtures that are usually found in the atmosphere.The optical properties are: the extinction coefficient (EX), the scattering coefficient (SC) (or scattering albedo), the asymmetry parameter (ASY), and the phase function (PH).The data are available for 61 wavelengths between 0.25-40 µm, and eight relative humidity values (for the water-soluble type).OPAC proposes a larger number of atmosphere models (aerosol mixtures) than 6S.
The particle types that are considered in the OPAC model are: (1) water-insoluble (INSO), which is equivalent to dust-like in 6S; (2) several mineral-type aerosols corresponding to the desert dust particles in 6S; (3) water-soluble (WASO), which is composed of various particles such as sulfates, nitrates, and other organic matter, whose size and optical properties are a function of the air relative humidity; the properties of this type of aerosol are considered stable in 6S; and (4) soot, which is identical to that in 6S, which is made up of very absorbent (black) particles composed of carbon.OPAC proposes the following mixtures: continental clean, continental average, continental polluted, urban, desert, maritime clean, maritime polluted, maritime tropical, Arctic, and Antarctic.The properties of a mixture are the averages weighted by the percentages of each aerosol in each mixture.A comparison of the OPAC dust-type aerosols (INSO and mineral) with the 6S dust aerosol (present in the continental and urban mixtures) has shown that the optical properties (EX, SC, and ASY) of the OPAC particles are more variable than those in the 6S [20].The WASO aerosols in OPAC are characterised by higher extinction and scattering coefficients than the same type of aerosol that are used in the 6S code (Figure 3).Also, the higher the relative humidity, the higher these coefficients are as well, since the particles are larger (and capture more photons by their larger effective section).The asymmetry parameter of WASO in OPAC is higher for short wavelengths and lower for wavelengths greater than 1.5 µm.The phase functions of INSO in OPAC and of dust in the 6S code are fairly close and very directional [20].For most wavelengths, the phase function curves of the WASO aerosols in the 6S code are close to those of OPAC for moderate humidity percentages (Figure 4).

Calculation of relative humidity profiles
Using the WASO aerosol from the OPAC model requires knowing the relative humidity in the layer of the atmosphere where the aerosol-scattering properties are calculated.The relative humidity profile is reconstructed using the formula of Leckner, B. [14] (Equations (5a) and (5b)), which is based on the temperature and pressure profiles according to the atmospheric model used (e.g., MIDSUM).The parameters of the water-soluble aerosols for the relative humidity values calculated at various heights are determined by means of a linear interpolation.
Methods for estimating the AOD and choice of an approach for REFLECT Three approaches are commonly for measuring and characterising aerosols: intensive in situ measurement campaigns; ground-based networks; and indirect estimates by remote sensing [21].The AOD can also be estimated from the visibility measured in weather stations.Measurement campaigns allow to characterise the physical, chemical, and optical properties of aerosols [22,23].However, since the photoradiometers that are used are very expensive, these measurements are limited to a certain number of stations.The largest ground-based network is AERONET (AErosol RObotic NETwork), which was established by NASA [24].This network currently includes more than 450 stations that perform spectral radiometric measurements that are capable of determining the most important properties of aerosols (single scattering albedo, size distribution, phase function, asymmetry factor, AOD).Other local networks also exist, such as the Canadian AEROCAN network (which includes about 10 stations and is part of the AERONET network), the EARLINET network (European Aerosol Research Lidar Network), and the AD-Net network (Asian Dust Network).These measurement networks are of insufficient spatial density for specific earth observation applications using high and very high-resolution satellite images.
Global indirect measurements of the AOD are regularly performed from satellite observations above dark targets with low and invariable reflectance values [21].The satellite approach is greatly appreciated, since it allows monitoring the great spatial and temporal variability of the AOD.For example, the AVHRR (Advanced Very High-Resolution Radiometer) sensor provides AOD estimates in the visible and near infrared bands with a resolution that is close to 1 km [25][26][27][28].The POLDER (Polarisation and Directionality of the Earth Reflectance) sensor is used to measure the properties of the aerosols above the oceans and continents with a resolution of approximately 6 km [29,30].Since 2000, the MODIS (Moderate Resolution Imaging Spectroradiometer) sensor has been used to deduce the AOD over the oceans and portions of the ground with daily coverage and for seven multispectral bands between 0.47-2.13µm, including a blue band for which the ground reflectance is low and aerosol scattering is high [31][32][33][34][35][36].The AOD estimated from MODIS and MISR (Multi-angle Imaging Spectroradiometer) sensors has a relative error of about 20% compared to AERONET measurements [37,38].The AOD data obtained from the MODIS sensor are provided at a daily scale for a latitude/longitude grid of 1 • [37], with a spatial resolution of 10 km.This resolution is not adequate for the correction of the spatial variability of the AOD in a HR or VHR image (pixels of 30 m or less).Models to estimate AOD at 550 nm from visibility are proposed in the literature, but the results are very approximate [39,40].
We performed a series of in situ measurements of AOD by a Microtops II Sun photometer in the vicinity of Montreal, Quebec, Canada.These measurements were compared to data from the nearest AEROCAN station located in Sherbrooke, Quebec, which was more than 100 km from the measurement sites.In addition, the data from this station were not available for all of the dates of interest.The AOD measured by the Microtops was also compared to estimates performed based on visibility measured by the weather station in St-Hubert, Quebec.For estimating the AOD based on visibility, we used the empirical formula: AOD 550,VIS = 3/(VIS − 3), yielding values intermediate between the model of So [39] and Qiu [40].The comparison showed that the data obtained from the AEROCAN network (Sherbrooke, Quebec, station), when available, are fairly close to the Microtops measurements.The AOD estimated based on visibility (AOD 550,VIS ) is often far from the measured values, which is probably because the visibility data are not precise measurements (e.g.,: 5 km, 15 km, 23 km . . .).
Given the state of the art concerning the sources of AOD data, the approach adopted in REFLECT is the dark target method applied to the images to be corrected.This approach is explained in Section 3.4.

New AOD spectral extrapolation model
For each wavelength of the spectral bands of each sensor, the scattering parameters ) are evaluated by interpolation of their values that were calculated in the 10 reference wavelengths (400 nm, 488 nm, 515 nm, 550 nm, 633 nm, 694 nm, 860 nm, 1536 nm, 2250 nm, and 3750 nm).In the 6S code, the AOD is estimated at these 10 wavelengths by extrapolation based on the AOD at 550 nm (AOD 550 ) according to the power law: The coefficient a varies between 0 and 4, with the standard value being 2. For each spectral band, these parameters are calculated by integrating the interpolated values for all of the wavelengths covered by the spectral band, taking into account the spectral sensitivity weighting of the sensor and the exoatmospheric solar spectral irradiance [41].
Bouroubi, M.Y. et al. [9] pointed out that the extrapolation law of AOD 550 that was used in the 6S code was not adapted to our study area (Montérégie, Quebec), which was possibly because of a dominance of small particles.The study of this issue led us to propose another law based on measurements (at 380 nm, 870 nm, 936 nm, and 1020 nm) that were taken by the Microtops in agricultural (Montérégie, Quebec), urban (Montreal, Quebec), and forest (Laurentians, Quebec) areas.Figure 5 shows that the power law that was used in the 6S code cannot follow the measured values of AOD across wavelengths.For a = 1, the law is suitable for short wavelengths (around 400 nm), but not for long wavelengths (800 nm and more).Conversely, the curve with a = 2 is close to the measurements for long wavelengths, but deviates for short wavelengths.We therefore propose a Gaussian function (Equation ( 6)) that closely matches the average measurements for the five wavelengths offered by the Microtops.
The values of the coefficients are K 1 = 1.91 and K 2 = 0.65.The type of site (rural or urban) does not affect this relationship [20].

Surface Properties
The reflectance of the environment (adjacent pixels) affects the digital value of the viewed pixel due to scattering within the viewing angle of the sensor.This secondary source of radiance becomes important if the target object is surrounded by a highly reflective object [15].In this case, the adjacency effect can be similar in magnitude to the reflectance of the target [42].According to several authors, the reflectance of the environment is the sum of the reflectances of all of the surfaces surrounding the target, and is weighted by a function f(r), which depends on their distance r to the target and the diffusivity of the atmosphere [42][43][44]: The environment function f(r) is calculated by taking into account the scattering properties of gases and aerosols.This formulation was included by Tanré, D. et al. [42] in the 5S (Simulation of the Satellite Signal in the Solar Spectrum) code and by Vermote, E.F. et al. [43] in the 6S code.The latter offers an anisotropic model to better account for the scattering properties of the atmospheric constituents in a specific direction.
Calculation of the reflectance of the environment according to the rings method proposed by Vermote, E.F. et al. [45] was introduced by Cavayas, F. et al. [8] in the first version of the REFLECT program.This method assumes that the target occupies the central portion of a circular area of radius r.The tests performed with the rings method led to the conclusion that the phenomenon of adjacency should be more significant than estimated by Vermote's formulation [45].An approach with the individual contribution of each pixel of the environment represents more accurately the actual conditions of a scene with a heterogeneous spatial reflectivity.Also, the specular reflection of certain pixels of the environment must be accounted for.These ideas were behind the development of the method used in REFLECT, which was based on: (1) the anisotropic environment function model proposed in the 6S code [45], (2) the independent contribution of surrounding pixels to the adjacency effect when calculating the weighted average of Equation (7), and (3) the presence of both diffuse (Lambertian) and specular components in the reflectance of the environment (Figure 6); hence: where ρ e is reflectance of each surrounding pixel (individual contribution) at a distance r, and f is the anisotropic environment function.The function g, which was inspired from the Phong model [46], takes into account the specularity of each pixel (Figure 7), and it is given by:  The coefficients k d and k s are, respectively, the proportions of diffusely and specularly reflected radiances, with k d + k s = 1.The calculation window must respect a compromise between two conditions: (i) accounting for all of the pixels of the environment for which f (r) has a significant value; and (ii) computation time savings.
In the case of a nadir view (Figure 6), the specular direction γ sp is given by: cos γ sp = cos θ s cos θ e (i, j) + sin θ s sin θ e (i, j) cos(ϕ sp − ϕ e (i, j)) The angles θ s and ϕ s are respectively solar zenith and azimuth, and ϕ sp = ϕ s + 180 • .For a pixel (i,j): ).sign(j) where r(i,j) is the distance between the environment pixel and the central pixel, and H sat is the height of the satellite.However, since H sat >> r, it is also possible to consider a single value θ e (i, j) = θ v .
Observation of the specular effect and estimation of the parameters k d , k s , and n Figure 8 shows a near-infrared (NIR) Landsat-7 ETM+ image of the Hertel Lake on Mount Saint-Hilaire, Montérégie, Quebec (45 • 32 40"N; 73 • 09 00"W), which is in the middle of a forest.This is an example of a very dark surface surrounded by very bright surfaces.We can observe that the profiles of the digital numbers of the NIR band follow a certain slope when they are extracted in the direction parallel to the solar azimuth.The profiles of the pixels extracted in the direction orthogonal to the solar azimuth did not show this trend.This phenomenon was observed for other dark surfaces surrounded by bright surfaces and for other sensors (data not shown) suggesting the presence of a specular effect.The parameters (k d , k s ) of Equation ( 9) can be deduced from these observations.In Figure 8, the apparent reflectances (simplified formula) for the dark pixels 1 and 2 are calculated, considering only the effect of the bright pixels on the dark pixels by: For point 1: Subtracting the apparent reflectances of pixels 1 and 2 gives: By calculating the apparent reflectances for the example of Figure 9 and using the atmospheric parameters given by REFLECT, we obtain k s = 0.14 and k d = 0.86.These values are in accordance with those found by the Columbia-Utrecht Reflectance and Texture Database (CUReT) [47].The exponent "n" is more difficult to estimate; we took the value n = 50, which is intermediate between the cases illustrated in Figure 6 where n = 20 and n = 100.

New Model for Topographic Correction
The formulation used in REFLECT is a generalisation of Equations ( 3) and ( 4) to include terrain with variable topography, as well as non-Lambertian reflectances.We explain here this formulation with expressions involving irradiances and radiances, as a first step, and then translate all of the values into reflectance terms, as used in the 6S code.The various sources of irradiance whose intensity is modulated by an inclined surface are: (1) solar irradiance transmitted directly by the atmosphere; (2) solar irradiance scattered by the atmosphere, which reaches the surface from any direction of the sky (hemispherical source); and (3) total solar irradiance (direct and diffuse), which is reflected by the surrounding area, and which illuminates the surface viewed after multiple reflections between the Earth and the atmosphere (hemispherical source).These three types of irradiance are given by the following equations: -Direct solar irradiance: E dir = E 0 t ↓ dir cos i s , in which i s is the angle of incidence of the sun's rays on an inclined terrain; it is given by: cos i s = cos θ s cos β + sin θ s sin β cos(ϕ s − α).The angles β and α are, respectively, the slope and orientation of the inclined surface.
-Diffuse sky irradiance: E di f f ,sky = E Horiz di f f ,sky g sky , where the factor g sky takes into account the losses in intensity of the sky irradiance relative to the irradiance incident to a horizontal terrain.These losses are incurred because parts of the sky are masked by the surface itself, given its inclination.
-Diffuse irradiance due to multiple Earth-atmosphere reflections: g env ; the factor g env considers the losses because parts of the hemisphere of the sky are masked due to the inclination of the surface itself and of the surrounding area (environment).
These irradiances are reflected in the direction of the sensor by the surface, which has its own factors of bidirectional reflectance and hemispherical-directional reflectance.Hence, the radiances can be written according to their source as follows: -Ground radiance due to direct radiation: Ground radiance due to diffuse sky radiation: L di f f ,sky = (1/π)E di f f ,sky ρ HD -Ground radiance due to multiple Earth-atmosphere reflections: Normalising these radiances in accordance with the 6S code (Equations ( 1) and ( 2)) gives: ρ B and ρ HD are respectively the bidirectional and hemispherical-directional components of the target reflectance ρ tar .All of these reflectances are transmitted directly to the sensor with losses due to the atmospheric scattering depending on the transmittance t ↑ dir .In addition, the adjacency effect can be written as a first approximation: ρ adj = T ↓ t ↑ di f f ρ env .By homogenising all of these reflectances for the denominator 1 − ρ env S and ignoring (in keeping with the 6S code) all of the terms involving the products of two reflectances and the spherical albedo, because their value approaches zero, we will finally get the equation: Equation ( 12) is the REFLECT form of Equation ( 4), which was used in the 6S code.It allows separately modulating the direct and diffuse components of solar radiation.For the direct component, the cosine correction is applied.The definition of the correction factor g sky of the diffuse component requires the use of a diffuse irradiance on the inclined plane model.Loutzenhiser, P.G. et al. [48] presented a review of these models, some of which have already been used in remote sensing studies [49][50][51][52][53].The one developed by Temps and Coulson [54] is the most elaborate according to Mefti, A. et al. [55].It is given by Equation ( 13):  The red-infrared relationship is used to differentiate vegetation from bare soil and water surfaces in the extraction of dense forest pixels (Figure 11).Among these pixels, only those that have an appropriate reflectance in the red band are kept.The condition to select a pixel (with the digital number DN) as dense forest is based on three thresholds: (DN N IR − DN Red > T 1 ) and (T 2 ≤ DN Red < T 3 ).To find the T 1 threshold, the digital numbers corresponding to reflectances of 0.02 in the red band and 0.15 in the NIR band are simulated for the conditions under which the image was acquired.The atmospheric parameters are calculated using REFLECT atmospheric routines by assuming the presence of a clear atmosphere (AOD = 0.05).The T 3 threshold is the first mode of the red band histogram, which was calculated for the pixels that meet the first criterion.The T 2 threshold was added to avoid selecting the border pixels of lakes and the edges of cloud-shaded areas.It was set at T 2 = T 3 − 3.In order to give the user more options to choose the "acceptance" thresholds of dark targets, we added to the algorithm a manual selection of the "degree of severity" of the darkness conditions.This option enables the user, for water-type dark targets, to vary the thresholds from MIN NIR to MOD NIR for the NIR and some digital numbers around MOD VIS for the visible bands.The same principle is applied for threshold T 1 when looking for dense forest pixels.
Another improvement was made to the search process for dark targets due to the following observation: in the case of images with a heterogeneous spatial distribution of the AOD, the histogram thresholding technique selects only dark targets that are located in the parts of the image where the atmosphere is the clearest.Thus, REFLECT allows dividing the scene into NxM sub-images and conducting the search in each of them.This local search enables detecting the spatial variation of the AOD and correcting it.It should be noted that Ouaidrari, H. and Vermote, E.F. [56] also proposed dividing the image into 4 × 4 quadrants to search for the dark targets in an image with a non-homogeneous atmosphere.Ahern, F.J. et al. [57] and Teillet, P.M. et al. [58] also addressed this issue.

Sensors Integration
We have incorporated in REFLECT most of the HR and VHR sensors that are currently used for Earth observation applications.These sensors are Landsat-5 TM, Landsat-7 ETM+, Landsat-8 OLI, Terra ASTER, SPOT-1 to SPOT-7, Sentinel-2, Ikonos, QuickBird, Pleiades, and WorldView-2 and WorldView-3.New sensors will be incorporated as they come on line.Incorporating a sensor requires: (1) the calibration coefficients "gain" and "offset" for each spectral band (provided with the image as metadata) that are used to convert the digital numbers (given in eight or 16 bits) to apparent radiances in W.m -2 .sr - .µm - according to sensor technical documentation; and (2) the spectral sensitivities SS k (λ) used for the averaging of spectral properties (L sat k , T gas , t ↓ dir , t ↓ di f f , t ↑ dir , t ↑ di f f , ρ atm , or S alb ) given at the wavelengths λ (with a step of ∆λ = 2.5 nm in the 6S code) within the limits of the spectral band k.Hence, for any parameter p, this integration is performed by the equation: E 0 (λ) is the exoatmospheric solar spectral irradiance.The calibration coefficients are determined during pre-launch operations, but their values drift over the years with the age of the optical system.Certain methods are proposed for making post-launch corrections [59,60], and these coefficients are regularly updated by the satellite operators.Since the satellites cover almost the entire planet, the areas observed often have very different reflectances; consequently, several gains are available on certain sensors to optimise the accuracy of the images (better contrast without saturation).For example, REFLECT allows the user to choose predefined calibration coefficients when it is possible (examples: ETM+ and ASTER) [61].Otherwise, values from the image metadata file should be used.

Software Implementation
The REFLECT software was programmed in C++ language.It was provided with a graphic interface and organised in modular form to reduce the computational burden on the main program, which calculates ground reflectance.A fast computation option based on the generation of internal look-up tables of the atmospheric parameters for a combination of AOD and terrain altitude values (which distinguish the various pixels of the image) considerably reduces the computation time.Another tool is provided in the interface for the calculation of atmospheric parameters.

Methodology and Data Used for Validation
The accuracy of ground reflectance retrieval by REFLECT was evaluated by means of several tests and using various types of images and related data.These tests involved: a.
Calculation of atmospheric parameters: Temperature and relative humidity data for the image acquisition dates were obtained from the Environment Canada website [62] to estimate the water vapour content in the atmosphere as well as gaseous transmittances, as explained in Section 3.1.
The AOD that was used for the calculation of the diffusion parameters was calculated from images using the dark targets method.b.
Correction of sensor gain effect and atmospheric effects: we used five Landsat-5 TM images, two Landsat-7 ETM+ images, one SPOT-1 HRV image, and one SPOT-5 HRG image.Images were acquired over several years between June and August for an area (intersection of the scenes of these images) of approximately 80 km × 65 km around Montreal, Quebec, Canada (Figure 12; Table 1).By comparing the fifth and 95th percentiles of the DNs and ground reflectance values calculated by REFLECT, we assessed the correction of sensor gain, atmospheric effects (additive and multiplicative), and conditions of illumination and observation.c.
Comparison of spectral signatures of selected materials: The ground reflectance values of three materials (water, asphalt, and vegetation) identified on the images of the scene in Figure 12 were compared with their simulated values from the USGS ASTER spectral library.The reflectance values of these materials for the spectral bands of the sensors used are simulated using following Equation (similar to Equation ( 14)): R band,material = ( where ρ material is the spectral reflectance from the ASTER library for a given material, SS band is the spectral sensitivity of the sensor by band, and E 0 is the exoatmospheric spectral solar irradiance.d. Comparison with ASD measurements: The ground reflectance values calculated for four ETM+ images acquired between 2000-2002 and one WorldView-2 image from 2011 were compared with the spectroradiometer (ASD FieldSpec HH Pro) measurements taken in agricultural fields in the Montérégie region, Quebec, Canada.The ASD reflectance values, which were measured concurrently with image acquisition, were incorporated by the spectral band of the images that was used according to the principle of Equation ( 16).Validation tests (a) to (c)-where all of the surfaces are considered horizontal-are discussed together in Section 5. Tests (d) and (e) are presented in Section 6, which is devoted to vegetation targets.Tests (f) and (g) deal mainly with topographic corrections, and are presented in Section 7.

Reflectance Retrieval for Global Scenes
This section deals with ground reflectance retrieval for a region comprising Montreal, Quebec, Canada, which is common to nine Landsat and SPOT images (Table 1; Figure 11).This operation validates the radiometric corrections that are associated with sensor gain, illumination and observation conditions, and additive and multiplicative atmospheric effects.After an atmospheric parameter calculation step, analysis of the digital numbers and ground reflectance values was carried out as follows: (1) by comparing the fifth and 95th percentiles of the digital numbers and ground reflectance (ρ ground , which is ρ tar of Equation ( 3), and includes values of all of the available images; and (2) by comparing the spectral signatures from raw (DNs) and corrected images (ρ ground ) of selected materials from the USGS ASTER library (water, asphalt, and vegetation).

Gaseous Transmittances
Gaseous transmittances are calculated by estimating the water vapour content in the air from the air relative humidity (H r ) and ambient temperature (T a ) data obtained from Canada's National Climate Archive for the Pierre Elliott Trudeau Airport station, which is located approximately at the centre of the scene delimited around Montreal (Table 1).The effect of H r (lower T gas for higher H r , and vice versa) can be seen in the NIR band because of the water vapour absorption lines between the 800-840 nm wavelengths [20].Gaseous transmittances are close to one in the short wavelengths (blue band), and decline with increasing proximity to the MIR (mid-infrared) wavelengths, mainly owing to water vapour absorption.

AOD from Dark Targets Method and Diffusion Parameters
The AOD at 550 nm is estimated using the dark targets method for all of the images (Table 2).An average value was taken because the AOD spatial distribution did not show any significant differences.We can see that most of the images have a clear atmosphere.The atmospheric parameters also show the magnitude of the additive effect in the short wavelengths (blue band).However, owing to the new interpolation rule AOD(λ) = f (AOD550, λ) given by Equation ( 6) of Section 3.1.2,the diffusion parameters (t dir ↓ , t dir ↑ , t diff ↓ , t diff ↑ , ρ atm and S alb ) are no longer overestimated in the blue band compared to tests carried out with the original rule of the 6S code (data not shown).

Comparison between Low and High Values of DN and Ground Reflectance
The images listed in Table 1, which were all acquired in the summer, were used to demonstrate the homogenisation of dynamic ranges by the radiometric corrections.We examined the minima (fifth percentiles) and maxima (95th percentiles) of the DNs and ground reflectance values (ρ ground ) calculated by REFLECT.The minima provide an indication of the correction of the additive atmospheric effects and sensor offsets, while the maxima show the correction of the multiplicative atmospheric effects and sensor gain.Table 3 shows that the differences between the DNs are significant for the different dates (atmospheric conditions) and different sensors (gains).
The DNs of the SPOT-1 HRV image of 1 August 1987 are particularly low, and those of the SPOT-5 HRG image of 29 July 2006 are particularly high.The reason for this is that SPOT has variable gain settings for each spectral band.The gain values in the green, red, and NIR bands for the HRV image of 1 August 1987 were 0.852, 0.794, and 0.888, respectively, versus 2.33, 2.8, and 1.31 in the same bands for the HRG image of 29 July 2006.The fifth and 95th percentiles of ρ ground are much less variable; differences are around 0.01 for the low reflectance values (visible bands), and around 0.02 for the high reflectance values (NIR and MIR).

Comparison between Values of DN and Ground Reflectance for Selected Materials
Pixels corresponding to clear water (lakes), asphalt (highways), and vegetation (deciduous trees) surfaces were selected (Table 4) on all of the previously used Landsat and SPOT images.The purpose of this test was to compare the ground reflectance values calculated by REFLECT for these materials from the various images and compare these reflectance values with those obtained from the USGS ASTER library.The latter values are calculated as indicated by Equation ( 16).
As we can see from Table 5, there is significant variability in the DNs by type of sensor.The variability by date is relatively less significant for the TM and ETM+ sensors.The reason for this is that the images were acquired under the same conditions of: (i) sensors gain per spectral band (unique for TM and fixed to high in visible bands and low in infrared bands for ETM+); (ii) season and time (similar illumination and observation conditions); and (iii) clear skies.However, the levels of the DNs from the SPOT HR images are not only different from those of the TM and ETM+ sensors; they also differ among themselves.As mentioned above, this is due to the difference in gains for the two acquisitions.Hence, the variation in DNs is affected more by the sensor than by the atmospheric conditions corresponding to clear skies in our case.The ground reflectance values derived from the various images show significantly lower variability compared to the DNs (Table 6).
For each material and each spectral band, the ground reflectance values are very similar for all of the sensors and all of the acquisition dates.In most cases, the variation does not exceed ±0.01 in the reflectance units.In the case of vegetation (deciduous trees), the reflectance values derived from the images are lower than those given in the ASTER library, which was probably because of mixed pixels.
Note that for the different sensors, we assigned unique values for the reflectance values of the materials from the ASTER library (Table 6).The reason for this is that the differences between the integrated reflectance values by spectral band of the sensors used (TM, ETM+, HRV, and HRG) were very small (less than 0.01).15).These measurements were compared with the ground reflectance values that were obtained with REFLECT.
-Landsat-7 ETM+ images The ASD measurements corresponding to the Landsat-7 ETM+ images showed significant variability by date and measurement point (Figure 13), indicating a certain variability in vegetation density.The comparison of the measured reflectance values and the reflectance values that are estimated using REFLECT (Figure 14) shows that the estimates are fairly dispersed relative to the measurements with mean absolute errors (|estimate-measurement|) of 0.012 in the blue band, 0.014 in the green band, 0.019 in the red band, and 0.046 in the NIR, which corresponds to a mean relative absolute error (|estimate-measurement|/measurement) of approximately 10%.These deviations seem to be acceptable, since the images were obtained under variable atmospheric conditions and have spatial resolutions (30 m) that are very different from those of the ASD measurements (1 circular metre).In the case of the image of 5 June 2000, a systematic overestimation is observed, which was probably due to the presence of a fine cloud layer that led to a high additive effect in all of the the spectral bands (non-selective diffusion).If the AOD is increased to correct the additive effect of the blue and green bands, the NIR reflectance is increased, since the transmittances (T↓ and T↑) are reduced without increasing the atmospheric reflectance (ρ atm ).
Differences between the estimated reflectance from the WorldView-2 image and ASD measurements are slightly less than those found for ETM+ images (mean relative absolute error of less that 10%) (Figure 15).The good correction of the coastal blue band indicated that the AOD spectral extrapolation formulae that was proposed in REFLECT gave correct AOD values for the lower wavelengths (close to 400 nm).The original power law used in 6S gave over estimation of AOD and thus overcorrection in this band (data not shown).The reflectances that were estimated for NIR-1 and NIR-2 (particularly affected by water vapour absorption) bands were very close, indicating a good correction of water vapour effect in the NIR domain.

Formosat-2 Images of a Vegetal Surface: Atmospheric and Adjacency Effects
This test shows the atmospheric corrections-including the adjacency effect-for nine Formosat-2 images that were provided in apparent reflectance values (16-bit format).The AOD was determined on dense vegetation dark targets only, since water was too affected by specular reflection (Figure 16a).In order to show the effect of atmospheric conditions as well as the adjacency effect on agricultural applications, we selected a vegetated surface that was located in the Boucherville area, Quebec, Canada, with very bright edges in the red band (Figure 16b).For the nine Formosat-2 images used (acquired on 5 June, 19 June, 21 June, 22 June, 25 June, 27 June, 28 June, 2 July, and 3 July 2005), we calculated the apparent NDVIs and extracted the profiles by field width (Figure 17a).Figure 16b shows that the apparent NDVIs are very different, even for dates separated by only one or two days; this is due to the difference in atmospheric conditions.The adjacency effects can also be seen on the NDVI profiles: except for the mixed pixels (containing soil and vegetation), the pixels at the edge of the field have lower NDVIs than the pixels in the centre (Figure 17b).Correction of the atmospheric effects reduces the NDVI underestimation obtained with apparent reflectances [64].In addition, with correction of the adjacency effect, the edge pixels have values closer to those of the centre.By excluding the ground pixels (numbered 1 to 4 and 14 to 18) and the mixed pixels (numbered 5, 6, and 13), the comparison between the apparent and corrected NDVIs of the field pixels (7 to 12) with the centre pixels (9 and 10) shows that the NDVIs that were calculated from the corrected images are less affected by the adjacency effect than the apparent NDVI.In addition, the adjacency effect is greater for pixels 11 and 12, which were located on the side of the field where the reflectance values in the red band from the neighbouring surface are higher (northeast direction, pixels 14 to 18).Correction of the adjacency effect for pixels 11 and 12 is clearly expressed by the reduction in the difference between their NDVIs and those of the pixels in the centre of the field.

Reflectance Retrieval for Inclined Surfaces
We conducted two validation tests for the topographic corrections that were used in REFLECT: the first for flat inclined planes (building roofs) in a high northern latitude, and the second for a high-elevation mountainous terrain covered by forests.In the first case, the choice of high latitudes results in a high solar zenith angle and very low atmospheric effects because of the clear sky.The choice of high elevation in the second case also minimises the incidence of aerosol particles, which are fairly rare at these elevations.

Flat Inclined Surfaces on an Ikonos Image
On an Ikonos image dating from 29 August 2002 covering an area of the Paulatuk region in the Northwest Territories, Canada, we identified large buildings with inclined symmetrical roofs.Since solar elevation was very low, the effect of the difference in irradiance on the two sides of the roof is quite significant (Figure 20).The roof slopes can be approximately estimated (ignoring the difference in azimuth between the sensor and the surface) from the ratio of the apparent lengths of each of the half roofs by means of a trigonometric calculation (Figure 21).The ratio between the half roof with a large apparent width and the half roof with a small apparent width can be expressed as follows: A B = L+d L−d = k; therefore: d = L k−1 k+1 .From the view incidence angle, we obtain: tan(i v ) = d h ; therefore: h = d tan(i v ) .The slope of the roof can be expressed as follows: tan Finally: β = arctan( (k−1) (k+1) tan(i v ) ) Angle i v represents the view incidence angle, and is calculated by taking into account the curvature of the earth; the result is i v = 28.44 • .For buildings 1, 2, and 3 (Figure 19), the k values are 8/6, 7/5, and 4/2.5 pixels, which gives β slopes of 15 • , 18 • , and 24 • , respectively.
By considering these slopes with the buildings' orientations, the illumination (solar zenith of 60.06 • and solar azimuth of 182.75 • ) and viewing parameters (incidence angle of 28.44 • ), the REFLECT radiometric corrections that were used with and without accounting for the topographic effects gives the results provided in Table 7. From this table, we can see that there are significant differences between the DNs of the light and dark sides of the three buildings.The same is true for the ground reflectance values without topographic corrections.We can also see the effectiveness of the topographic corrections in finding reflectance values that are very close for the two sides of the same building.Hence, the proposed topographic correction method works well for flat surfaces.

Forest Canopy in Inclined Terrain on SPOT-1 Images
For the validation (and adaptation) of the topographic corrections in forested environments, we used three SPOT-1 images (5 September 1988, 24 September 1989, and 3 September 1990) of the Tarn region, France (43 • 39 N; 02 • 42 E), with a highly varied relief (elevations between 300-1200 m).The result shows that the distribution of the DNs of the NIR band for low slopes (β < 8 • ) is not much affected by the orientation of the pixels α relative to the solar azimuth ϕ s .On the other hand, for slopes of more than 18 • , surfaces oriented towards the sun (|α − ϕ s | <20 • ) have DNs that are clearly higher than those of surfaces with low slopes.However, surfaces with slopes of more than 18 • and oriented in the direction opposite the sun (|α − ϕ s | > 160 • ) have DNs that are barely lower than the surfaces with a low slope.Hence, in a forest environment, the magnitude of the slope effect varies depending on the orientation of the surfaces relative to the sun.This observation has also been reported by Ekstrand,S. [65].
For atmospheric corrections, AOD values were estimated using the dark targets method on the lakes in the scene.Values obtained were 0.05, 0.05, and 0.07 for the 1988, 1989, and 1990 images, respectively.First, we used REFLECT to calculate the ground reflectance values without considering the topographic effects.The overestimation of the reflectance values in inclined terrains oriented toward the sun is observed.However, surfaces oriented in the direction opposite to the sun are not significantly underestimated relative to the flat terrains.This could be explained by the geotropic nature of trees and the amount of shading within the canopy that varies depending on the slope exposition to the sun [66].Hence, by applying the topographic correction of the model that we developed for flat inclined surfaces, a very high overcorrection is obtained for the surfaces that are oriented in the direction opposite the solar azimuth.The model that was developed for inclined planes is not applicable to a forest canopy.Therefore, one would expect that the equivalent slope of a forest canopy would be gentler than that of the terrain in which it is located.We therefore propose a slope correction that depends on the orientation of the pixel relative to the solar azimuth given by: The parameters β and β' are actual and corrected slopes, respectively.The function κ(ϕ s − α), which is less than or equal to 1, is determined by analysing the data derived from the SPOT images.
In order to find the equivalent slope that adapts our topographic correction model to the forest canopy, we looked for the multiplicative factor κ (by varying it from 1 to 0), which minimises the RMSE between the ground reflectance of inclined terrains (three slope ranges: 6 • < β < 12 • , 12 • < β < 18 • , and β > 18 • ) and of flat terrains (β < 6 • ) for all of the forested area of the three SPOT-1 images used, and all of the forest stands combined.The obtained values of κ can be smoothed based on the difference between the solar azimuth ϕ s and the orientation of the surface α, which was noted as δϕ, by the function : where A = 0.1, B = 0.9, and C is an exponent that defines the shape of the curve.All that remains to determine is whether these values are of the same order of magnitude under other conditions (other images, other forested environments, etc.).The application of topographic correction with the equivalent slopes principle yields very good results for the NIR band (the most affected by topographic effects in our example) of the SPOT-1 images used; the mean error is 0.02 reflectance units.This error is very small compared to that obtained with the Minnaert model [67], which is frequently used for this type of application [68,69], and yields an error in the order of 0.1 reflectance units.For example, we can see on the NIR band of the image of 5 September 1988 (Figure 22), in the circled areas, that the ground reflectance values of the inclined surfaces (steep slope) oriented toward the sun (low |orientation-azimuth|) have been homogenised relative to the digital numbers.

Conclusions
The objective of this research was to develop a tool (REFLECT software package) for ground reflectance restitution, which considers all of the radiometric effects affecting multispectral images, including illumination and viewing conditions, sensor properties, atmospheric conditions, and topography.REFLECT uses the formulation and the routines of the 6S code, and introduces several new features such as: the incorporation of a more recent aerosol model (OPAC); the proposal of a law of extrapolation of the AOD at 550 nm to the other wavelengths that are better adapted to the studied area (Quebec, Canada); accounting for the adjacency effect by means of a model inspired by the Phong model; improvement of the dark targets method that is used for estimating the AOD with the possibility of subdividing the image and the semi-automatic choice of the target "darkness" criterion; topographic corrections for smooth surfaces by means of a model that separates the direct and diffuse components of solar radiation and for forest canopy by means of the "equivalent slopes" principle; and the reduction of processing time in the case of large images by generating internal look-up tables (LUTs) of atmospheric parameters.
Validation tests showed that the ground reflectance is retrieved with a good accuracy: the variability of this reflectance for the same surfaces was approximately ±0.01 among the atmospheric conditions and sensors for all of the spectral bands.The topographic correction that was proposed in REFLECT allowed reducing the error of ground reflectance retrieval for flat inclined surfaces (roofs) observed by an Ikonos image from 0.1 without topographic corrections to 0.01.This topographic correction model adapted to forest canopy reduced the error of NIR reflectance retrieval from three SPOT images from 0.08 (without topographic correction) to less than 0.02 (with the "equivalent slopes" principle).These performances are very satisfactory compared to ATCOR (used in PCI Geomatics) or FLAASH (used in ENVI), according to the comparative study of Fuyi, T. et al. [70].

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

Figure 1 .
Figure 1.Solar radiation components received by a remote sensing sensor according to the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) code.

Figure 2 .
Figure 2. Relation between aerosol optical depth (AOD) at 550 nm measured by Microtops-II instrument in Montérégie, Québec, Canada, and relative humidity given by the closest meteorological station.

Figure 3 .
Figure 3.Comparison between aerosol properties extinction coefficient (EX), scattering coefficient (SC), and asymmetry parameter (ASY) of the Optical Properties of Aerosols and Clouds (OPAC) and 6S models.

Figure 5 .
Figure 5.Comparison of average relative AOD to its value at 550 nm with 6S and proposed models.

Figure 6 .
Figure 6.Adjacency effect with a specular component of the reflectance of the environment.

Figure 7 .
Figure 7. Function g(θ s , θ v , ∆ϕ) with a Lambertian and a specular component of the reflectance of the environment for different values of θ s , k d , k s , and n.

Figure 8 .
Figure 8. Specular effect of the reflectance of the environment in the NIR band of a Landsat-7 ETM+ image acquired on 8 June 2001 and covering Hertel lake in Mont St-Hilaire, Quebec, Canada.

Figure 9 .
Figure 9. Specular adjacency effect of a bright pixel on a dark pixel according to the position relative to sun azimuth.

) 3 . 4 .
Dark Targets Method for AOD EstimationThe dark targets method developed in REFLECT for AOD estimation is based on clear deep water and very dense forest pixels identification by an automatic (or semi-automatic) histogram thresholding in the NIR and visible spectral bands.For clear deep-water pixels, REFLECT uses the NIR band to distinguish water from land.By keeping among the water pixels those having the lowest values in the available visible bands, turbid or shallow water pixels are avoided.The thresholds in the NIR and in the visible bands (examples: blue, green, and red) are defined, respectively, as the first local minimum (MIN NIR ) and the first local modes (MOD VIS = MOD blue , MOD green , and MOD red ) of the histograms of the image digital numbers (Figure10).However, for an image contaminated by clouds, and to avoid selecting cloud-shadowed pixels, a stricter condition can be imposed in the NIR band by defining a lower threshold than MIN NIR ; for example, the first mode of the histogram MOD NIR .

Figure 11 .
Figure 11.Histogram thresholding of digital number (DN) differences of the NIR and red bands (e.g., image ETM+ of 8 June 2001) for the detection of dense forest dark targets.
e. Correction of adjacency effects: A series of nine Formosat-2 images (Taiwanese experimental sensor operated by SPOT Image, a spatial resolution of 8 m, fast revisit time [63]) acquired between 5 June and 3 July 2005 near Montreal were used to show correction of the adjacency effect (in addition to correction of sensor gain and atmospheric effects) for a vegetated surface surrounded by highly reflective surfaces in the red band.f.Topographic corrections for flat surfaces: One Ikonos image acquired on 29 August 2002 over the Paulatuk region, Northwest Territories, Canada, was used to test REFLECT's topographic effects correction model on flat inclined surfaces (roofs of large buildings).g.Topographic corrections for a forest canopy: Three SPOT-1 HR images acquired during the same period in 1988, 1989, and 1990 over the mountainous region of Tarn, France, as well as a DEM (digital elevation model) of the region enabled us to validate the adaptation of the topographic correction model to a forest canopy.

Figure 12 .
Figure 12.Scene used to compare DNs and ground-based reflectances for five Landsat-5 TM images, two Landsat-7 ETM + images, 1 HRV SPOT image and 1 SPOT 5 HRG images.

6. Reflectance Retrieval on Vegetation Targets 6 . 1 .
Comparison with ASD Measurements for Agricultural Targets A series of reflectance measurements were made with the ASD FieldSpec HH Pro spectroradiometer in agricultural fields (corn at several growth stages) located in the Montérégie, Quebec, Canada, concurrently with the acquisition of the Landsat-7 ETM+ images (5 June 2000, 26 July 2001, 11 August 2001, and 14 August 2002) and WorldView-2 images (5 July 2011) (Figure

Figure 13 .
Figure 13.ASD spectroradiometer measurements in agricultural fields located in Montérégie, Quebec, Canada.The blue, green, red, and NIR spectral bands of ETM+ are indicated by blue, green, red, and purple lines, respectively.

Figure 15 .
Figure 15.Comparison of reflectances measured (ASD) and estimated (REFLECT) for the World-View-2 image of 5 July 2011.

Figure 17 .
Figure 17.Profile of apparent NDVI (normalized difference vegetation index) extracted over the width of the plant surface selected on the Formosat-2 images that were used.

Figure 18
Figure18shows the red and NIR apparent reflectance values of the study area, as well as the environment reflectance values that were calculated as explained previously.

Figure 18 .
Figure 18.Example of apparent reflectance and environmental reflectance in the red and NIR bands for the selected farm field on the Formosat-2 images used.

Figure 19
Figure 19 illustrates the NDVIs calculated from ground reflectance values retrieved by REFLECT.

Figure 19 .
Figure 19.Corrected NDVI profile extracted on the width of the selected field on the Formosat-2 images.

Figure 21 .
Figure 21.Calculation of the slope of the rooftop parts from viewing geometry of the Ikonos image.

Figure 22 .
Figure 22.Example of topographic effects reduction by the equivalent slopes model applied to an NIR band of the SPOT 1 image acquired on 5 September 1988 in the Tarn region, France.

Table 1 .
Gaseous transmittances calculated for Landsat and SPOT images used for validation.Temperature and relative humidity were used for atmospheric water vapour estimation.MIR: mid-infrared.

Table 2 .
AOD values calculated with REFLECT dark objects method and corresponding atmospheric parameters (t dir ↓ , t diff ↓ , t dir ↑ , t diff ↑ , ρ atm , and S alb ) for Landsat and SPOT images used for validation.

Table 3 .
The fifth and 95th percentiles of DNs and ground reflectances for Landsat and SPOT images.

Table 4 .
Sites used for comparing ground reflectances of selected materials with their values given by the USGS ASTER library.All of these sites are around Quebec, Canada.

Table 5 .
Averages of DNs of the selected materials on Landsat and SPOT images.

Table 6 .
Averages of ground reflectances of the selected materials on Landsat and SPOT images.

Table 7 .
DNs and ground reflectances without and with topographic corrections for bright and dark sides of the rooftops in Ikonos Paulatuk.
of the direct solar radiation in the ascending and descending paths t ↑ diff , t ↓ diff diffusion transmittances of the diffuse solar radiation in the ascending and descending paths w total water vapour content in the atmospheric column dir diffusion transmittances