Скачать 218.86 Kb.
|ASSESSMENT OF RADIOMETRIC CORRECTION TECHNIQUES IN ANALYZING VEGETATION VARIABILITY AND CHANGE USING TIME SERIES OF LANDSAT IMAGES|
Sergio M. Vicente-Serrano1*, Fernando Pérez-Cabello2 and Teodoro Lasanta1
1 Instituto Pirenaico de Ecología, CSIC (Spanish Research Council), Campus de Aula Dei, P.O. Box 202, Zaragoza 50080, Spain
2 Departamento de Geografía. Universidad de Zaragoza. C/ Pedro Cerbuna 12. 50009. Zaragoza. Spain.
Abstract. The homogeneity of time series of satellite images is crucial when studying abrupt or gradual changes in vegetation cover via remote sensing data. Various sources of noise affect the information received by satellites, making it difficult to differentiate the surface signal from noise and complicates attempts to obtain homogeneous time series. We compare different procedures developed to create homogeneous time series of Landsat images, including sensor calibration, atmospheric and topographic correction, and radiometric normalization. Two seasonal time series of Landsat images were created for the middle Ebro Valley (NE Spain) covering the period 1984–2007. Different processing steps were tested and the best option selected according to quantitative statistics obtained from invariant areas, simultaneous medium-resolution images, and field measurements. The optimum procedure includes cross-calibration between Landsat sensors, atmospheric correction using complex radiative transfer models, a non-lambertian topographic correction, and a relative radiometric normalization using an automatic procedure. Finally, three case studies are presented to illustrate the role of the different radiometric correction procedures when analyzing and explaining gradual and abrupt temporal changes in vegetation cover, as well as temporal variability. We have shown that to analyse different vegetation processes with Landsat data, it is necessary to accurately ensure the homogeneity of the multitemporal datasets by means of complex radiometric correction procedures. Failure to follow such a procedure may mean that the analyzed processes are non-recognizable and that the obtained results are invalid.
Key-words. Landsat time series, TM-ETM+ cross-calibration, atmospheric correction, relative normalization, vegetation change, Ebro Valley, Spain
The Landsat program for Earth observation has provided invaluable information on the Earth’s surface characteristics over the past four decades (Cohen and Goward, 2004). Landsat images have been widely used for land cover mapping and the creation of vegetation inventories at different spatial scales (e.g., Bossard et al., 2000). Moreover, the systematic archiving of Landsat data makes this information highly valuable for retrospective analyses of land surface characteristics. Together with other types of satellite images such as NOAA-AVHRR (National Oceanic and Atmospheric Administration-Advanced Very High Resolution Radiometer) and MODIS (Moderate Resolution Imaging Spectroradiometer), multitemporal Landsat data have been widely used in recent decades to identify changes in land cover (e.g., Lenney et al., 1996). The spatial and spectral resolution of Landsat images makes these data highly suitable in analyzing both abrupt and gradual changes in vegetation cover and monitoring environmental processes such as degradation and desertification (Almeida-Filho and Shimabukuro, 2002), deforestation (Cohen et al., 2002; Huang et al., 2007), habitat fragmentation (Millington et al., 2003), forest succession (Song and Woodcock, 2003; Song et al., 2007), overgrazing (Jano et al., 1998; Pickup and Chewing, 1994), rangeland monitoring (Hostert et al., 2003), and vegetation recovery after natural disturbance such as volcanism (Lawrence and Ripple, 1999) and forest fires (Viedma et al., 1997; Lozano et al., 2007).
Nevertheless, there exist constraints in using Landsat data for multitemporal studies because of problems in obtaining homogeneous time series, not affected by non-surface noise, and in which the images are comparable between different dates since the data only report on surface conditions. Multitemporal satellite-image datasets are affected by different sources of noise related to the stability of sensors, changes in satellite responsivity, changes in illumination, atmospheric effects, etc. These problems are not unique to Landsat images, but they are more difficult to overcome in Landsat images compared with images compiled by other satellites because the low temporal resolution of Landsat images makes it impossible to apply simple homogenization procedures such as composite creation (Holben, 1986) and temporal filtering (Van Dijk et al., 1987), which markedly reduce the non-surface noise.
Notable efforts have been made to reduce non-surface noise in Landsat images, including attempts to calibrate the sensor to correct lifetime radiometric trends (Teillet et al., 2004; de Vries et al., 2007), cross-calibrate the images obtained from different sensors (Teillet et al., 2006; Röder et al., 2005), minimize atmospheric noise (Chavez, 1989; Ouaidrari and Vermote, 1999; Liang et al., 2001), and reduce the influence of topography (Civco, 1989; Gu and Gillespie, 1998; Pons and Solé, 1994); however, other studies have shown that the application of accurate sensor calibrations and complex atmospheric corrections does not guarantee the multitemporal homogeneity of Landsat datasets (Schroeder et al., 2006) because complete atmospheric properties are difficult to quantify and simplifications are commonly assumed. This problem has led to the development of relative radiometric normalization techniques based on adjustments to the radiometric properties of an image time series to match that of a single reference image. In recent years, efforts have been made to develop methods to select invariant pixels for a reliable application of relative normalization techniques (Du et al., 2002; Chen et al., 2005; Canty et al., 2004).
Protocols have been proposed in processing multitemporal Landsat datasets (e.g., Hall et al., 1991; Hill and Strum, 1991; Han et al., 2007), comprising the following steps: i) geometric correction, ii) calibration of the satellite signal to obtain Top of the Atmosphere Radiance, iii) atmospheric correction to estimate surface reflectance, iv) topographic correction, and v) relative radiometric normalization between images obtained on different dates. Radiometric processing is recommended to be done prior to geometric processing since this resampling step generally smoothness the data set. Nevertheless, Landsat users in Europe commonly do not have access to data that have not already been geometrically corrected (the Level 1 System Corrected -1G- from the European Space Agency (ESA) is considered the standard product for most users since previous levels require extensive processing). When these protocols are applied, several decisions must be made due to the different procedures available at each step.
The high degree of interest in using multitemporal Landsat datasets is in contrast with the small number of studies that have tested different procedures of radiometric correction with the aim of obtaining temporally homogeneous images. The majority of studies have tested the influence of atmospheric correction (Song and Woodcock, 2003; Norjamäki and Tokola, 2007) or relative normalization (Yuan and Elvidge, 1996; Tokola et al., 1998; Heo and Fitzhugh, 2000; Olthof et al., 2005) on the temporal stability of time series of Landsat images. In contrast, few studies have tested the influence of each of the above steps in obtaining temporally stable time series of images and few have assessed the relative importance of using simple or complex techniques at each step. Existing studies have shown that relative normalization is the most critical step in obtaining temporal stability in a series of images (Schroeder et al., 2006; Janzen et al., 2006), although the calibration procedure also has a noticeable effect on the final results (Paolini et al., 2006) and atmospheric correction is important in obtaining accurately estimate surface reflectance's and accurate magnitudes of vegetation indices over time (Song and Woodcock, 2003).
Few studies have analyzed the role of complete radiometric correction protocols in processing multitemporal Landsat data when a number of different vegetation processes are of interest. Nevertheless, this is a crucial topic when analysing vegetation multitemporal dynamic using Landsat data. Thus different results have been found for land classification (Song et al., 2001; Paolini et al., 2006; Norjamäki and Tokola, 2007) and forest succession (Song and Woodcock, 2003; Schroeder et al., 2006) as a function of the radiometric correction applied.
The present study involves a complete evaluation of various radiometric correction processes required to obtain accurate time series of Landsat imagery, including calibration, atmospheric correction, topographic correction, and relative normalization. The objective is to identify the influence of different processing steps on multitemporal spectral reflectance trajectories developed with Landsat data and also to detect their performance to analyse different vegetation processes. Three case studies related to different processes of vegetation change and variability are provided to illustrate how different results can be obtained as a function of the employed radiometric correction protocol and the ecological process of interest: land cover change, forest regeneration after fire and ecosystem response to climate variability. The case studies were carried out in a complex ecological region in which several natural and human-induced processes of changes in vegetation cover have been recorded in recent decades.
2. Study area
The study area is the central Ebro Valley, Spain (Landsat Path 199 Row 31), located in the northernmost semiarid region in Europe. Figure 1 shows the location of the study area, including a detailed land cover map. Surrounded by mountain chains, the Ebro Valley has a Mediterranean climate with continental characteristics, with marked spatial and seasonal variations in precipitation; the dry season occurs during the summer months. The elevation ranges from less than 300 m above sea level in the middle of the Ebro Depression to more than 2000 m a.s.l. in the Prepyrenees (North) and 2000 m a.s.l. in the highest peaks of the Iberian Range (South). In the central part of the valley, mean annual precipitation is 326 mm, with marked seasonality. The study area shows a strongly negative water balance (precipitation minus potential evapotranspiration), greater than 900 mm. The dominant land covers include steppes and herbaceous cultivation in dry farming areas. Coniferous forests (mainly Pinus halepensis) cover the few slopes of the region. Vegetation distribution is strongly controlled by aridity (Vicente-Serrano et al., 2006), and drought conditions have a marked influence on vegetation cover and activity (Vicente-Serrano, 2007). The lithology of the area is characterised by limestones and gypsums (Peña et al., 2002) that contribute to its aridity, since the soils are unable to retain the water as a consequence of the high hydraulic conductivity.
3.1. Satellite imagery database
To determine the patterns of vegetation variability and change, it is important to take into account seasonal variations in the distributions of different vegetation types. In the middle Ebro Valley, herbaceous species, shrubs, and forests show contrasting seasonal variations in vegetation activity (Vicente-Serrano et al., 2006). The peak activity for herbaceous species and shrubs occurs in spring; for forests in summer. These seasonal fluctuations make it difficult to capture the full range of vegetation processes with just a single season of imagery.
Abrupt changes in vegetation activity commonly occur as a consequence of climate seasonality; consequently, it is important that the capture dates of images are similar in different years. It is not possible to combine images taken in different months, as this would have a strong influence on the temporal homogeneity of the dataset.
We reviewed all of the available Landsat-Thematic Mapper (TM) and -Enhanced Thematic Mapper (ETM+) images in the archives of the European Space Agency (ESA). Since most studies based on Landsat data consider ETM+ and TM radiometry to be comparable, We used TM data from 1984 and, then, switched to the ETM+ data when it became available in 1999 due to the better calibration of this sensor (Teillet et al., 2001); we subsequently switched back to using TM date due to the ETM+ SLC failure in 2003. Frequent cloud cover in spring means that only the month of March has a sufficient number of images to enable an analysis of vegetation activity. Months in summer pose fewer problems in terms of obtaining reliable clouds-free images; therefore, we selected the month of August to create a second time series. A total of 28 Landsat-TM and -ETM+ images taken between 1984 and 2007 were acquired from ESA, 16 corresponding to the summer season and 12 to spring. Table 1 lists the dates and types (TM or ETM+) of the selected images.
3.2. Geometric correction
The images were orthorectified using control points and a 1-m digital elevation model obtained from stereo pairs of aerial photographs, and resampled to 30 m to match the TM and ETM+ resolutions. The image taken on 8 August 2000 with good visibility and being free of clouds, was registered using orthorectified digital aerial photographs as a reference. The rest of the images were co-registered to this image using control points. The images were orthorectified following the method of Palà and Pons (1995), which includes elevation data in performing a polynomial geometric correction. The X and Y root mean square error was less than 15 m (0.5 pixels) for all images, guaranteeing a precise geometric match among them. After geometric correction, cloud-covered and cloud shadows were manually digitized and eliminated.
A precise calibration is required to convert DN (Digital Numbers) to satellite radiances in W m–2 sr–1 m–1. This is highly problematic for Landsat 5 imagery because calibration coefficients are not constant over time. A large variation, on the order of 20%, has been observed between the prelaunch gain coefficients and postlaunch cross-calibration (Chander et al., 2004; Chander et al., 2007). Moreover, the different spectral responses of TM and ETM+ sensors introduces comparability problems between Landsat 5 and Landsat 7 images (Teillet et al., 2001).
For Landsat 5 imagery, the European Spatial Agency (ESA) has used constant calibration dynamic ranges since 1984 (http://earth.esa.int/pub/ESA_DOC/landsat_FAQ/#_Toc69122047), embedded within the product format. Chander et al. (2004) and Teillet et al. (2004) demonstrated an exponential decay of the solar reflective bands of Landsat 5-TM since 1984, with some differences between bands. We corrected the calibration coefficients embedded in the ESA products with reference to the time after launch of Landsat 5 in 1984, applying the equations formulated by Teillet et al. (2004). This procedure is recommended by the ESA in recalibrating Landsat products (Saunier and Rodríguez, 2006). The coefficients of calibration for Landsat 7-ETM+ were obtained according to upper and lower at-satellite radiance indicated in the Landsat-7 Science Data Users Handbook (http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_toc.html) using the values corresponding to the date (before or after July 1, 2000) and the type of gain (high or low, as embedded within the product format).
Although the Landsat 5-TM and Landsat 7-ETM+ bands are commonly considered to be comparable, previous studies have demonstrated important differences between the two that may affect comparability, suggesting the need to apply cross-calibration procedures to both sensors because for the majority of bands (2, 3, 4, and 7) the TM sensor underestimates the radiance values regarding ETM+ (Teillet et al., 2001; Vogelmann et al., 2001).
Vogelmann et al. (2001) developed a simple procedure to cross-calibrate the Landsat-TM and ETM+ images, applicable to Level 1G formats. These authors proposed empirically derived slope and intercept values to convert Landsat 7 ETM+ DNs to Landsat 5 TM DNs based on two simultaneous images taken on June 2, 1999. We used these slope and intercept values in the present study to adjust the Landsat 7-ETM+ DNs to Landsat 5-TM DNs.
Satellite radiance was obtained for Landsat 5-TM quantized calibrated pixel values in DNs, Landsat 7-ETM+ quantized calibrated pixel values in DNs, and the cross-calibrated Landsat 7 ETM+ quantized calibrated pixel values in DNs to Landsat 5-TM quantized calibrated pixel values in DNs according to:
where Lsat is the satellite radiance in W m–2 sr–1 μm–1 for band . Grescale and Brescale are the calculated band-specific rescaling factors. Satellite radiances were converted to Top Of the Atmosphere (TOA) reflectances according to
where is the TOA reflectance for band , d is the earth–sun distance in astronomical units, ESUN is the mean solar exoatmospheric irradiance for band , and s is the solar zenith angle in degrees. ESUN values were obtained from Chander and Markham (2003) for TM images and from the Landsat-7 Science Data User Handbook for ETM+ images (http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_toc.html).
3.4. Atmospheric correction
Although TOA reflectance values are widely used in inventory and ecosystem studies, they do not take into account signal attenuation by the atmosphere, which strongly affects the intercomparability of satellite images taken on different dates. Upward and downward irradiance is modified by two atmospheric processes: absorption by gases and scattering by aerosols and water molecules (Vermote et al., 1997a; Lu et al., 2002). There are noticeable daily variations in atmospheric gas concentrations and aerosol volumes; these must be taken into account in multitemporal studies. Moreover, the effect of atmospheric processes is spectral-dependent (Teillet and Holben, 1993; Cachorro et al., 2000), affecting the magnitude of band ratio operations such as vegetation indices (Myneni and Asrar, 1994).
Different models have been developed to minimize the noise introduced by atmospheric processes on the signal received by the satellite, ranging from simple methods based on information contained in the image [e.g., Dark Object Subtraction (DOS)-based methods (Chavez, 1996)] to complex radiative transference models such as SMAC (Rahman and Dedieu, 1994), 6S (Vermote et al., 1997a), MODTRAN (Berk et al., 1999), and ATCOR (Richter, 1996) that simulate the atmosphere/light interactions between the sun surface and surface-sensor trajectories.
In this paper we tested two methods developed to minimize the atmospheric effects on Landsat time series: i) one based on a modification of the DOS method that includes some atmospheric information (Dark Target Approach -DTA-), and ii) one based on a complex radiative transfer code that includes some atmospheric information commonly available in global climate databases and some parameters obtained from the images themselves.
3.4.1. DTA method
Dark Target Approach (DTA)-based methods assume that some areas of the images have near-zero reflectance. Although there are various DTA-based methods (Chavez, 1996), Song et al. (2001) reported better results when including in the method atmospheric transmittance and Rayleigh scattering. We therefore followed this approach in the present study. Surface reflectance can be obtained as follows (Song et al., 2001):
where is the radiance (W m–2 sr–1 m–1) in areas with zero or near-zero reflectance, 0 is the optical depth of the atmosphere, u0 is the cosine of the solar zenith angle, uv is the cosine of the zenithal view angle, and Edown is the diffuse irradiance at the surface due to the scattered solar flux in the atmosphere (W m–2 m–1). The most critical step is related to the calculation of Lhaze, for which there are a number of alternatives (Chavez, 1989, 1996). In this study, we calculated Lhaze according to the method of Song et al. (2001) and Schroeder et al. (2006):
where Ldark is the lowest radiance in the image whose value is taken from at least 1000 pixels (Teillet and Fedosejevs, 1995). The estimation of 0 is complex as it requires in situ measurements, which are unavailable for the dates of the selected images. We therefore used constant values of 0 for the entire set of images, based on USA standard values (Dozier, 1989) that are in agreement with the values reported by Yang and Vidal (1990) for NE Spain. 0 values can be consulted in Pons and Solé-Sugrañes (1994). Edown for a Rayleigh atmosphere was estimated as zero aerosol optical depth at 550 nm using the 6S radiative transfer code, following Schroeder et al. (2006).
3.4.2. Radiative transfer model
Among the available codes, we used the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) code (Vermote et al., 1997a) to invert the TOA radiance and obtain surface reflectance values from TM and ETM+ data assuming lambertian targets. 6S code takes into account gaseous transmission and Rayleigh and aerosol transmission and intrinsic reflectance. Although 6S enables the inclusion of detailed atmospheric conditions in several layers, as obtained from radiosonde data, such data are commonly unavailable. 6S code also enables the user to restrict the number of inputs and to assume certain constants. In this study, we used ozone concentrations and column water vapor as gaseous inputs for 6S calculations because they are the most critical parameters for atmospheric correction (Arino et al., 1997). Daily ozone concentrations (in atm cm–2) can be obtained from Total Ozone Mapping Spectrometer (TOMS) data, at a resolution of 1.25 1.00, from the NASA GSFC Data Active Archive Center (http://toms.gsfc.nasa.gov/ozone/ozone_v8.html). For the period for which ozone data were unavailable (between 1994 and 1996), the data were calculated following Van Heuklon (1979) as a function of the day of year:
where lat is the latitude in degrees, long is the longitude in degrees, and d is the day of year.
The column water vapor (in g cm–2) is obtained from NCEP reanalysis daily data at a resolution of 2.5 2.5 (http://dss.ucar.edu/datasets/ds090.0/). Elevation was also included by means of the DEM described in Section 3.2.
The most critical parameter for atmospheric correction is the Aerosol Optical Thickness (AOT). Although it is ideal is to obtain simultaneous measurements corresponding to the overpass of the satellite, in the study area there are no systematic measurements of this parameter, this being the case in most regions of the world. To solve this problem, various methods have been developed to obtain reliable estimations of AOT (Liang et al., 1997; Wen et al., 1999). Most such methods are based on the relatively strong aerosol radiance effect in low-surface-reflectance areas compared with bright areas, as it is then possible to apply an inversion procedure in assessing AOT. In this study, we chose to use the Dense Dark Vegetation (DDV) method (Liang et al., 1997). This approach is based on the empirical relationship observed between the visible and infrared (IR) bands in areas with dark and dense vegetation cover. The method is based on the weak influence of aerosols on the mid-IR signal, as most aerosol particles sizes are smaller than the wavelength in this part of the spectrum. Kaufman et al. (1996) proposed the following relationship for dark vegetation canopies:
where , , and are the surface reflectances for bands 1, 3, and 7 of Landsat-TM images, respectively.
We assumed that each Landsat scene has uniform AOT, and, following Song et al. (2001), dense and dark vegetation areas were identified over the whole image (and NDVI (Normalized Difference Vegetation Index) > 0.1). Average TOA reflectance values were obtained from bands 1, 3, and 7 in these regions, and band 7 was used to predict surface reflectance values for bands 1 and 3 according to the empirical relationship stated above. We iteratively ran the 6S code with the ozone and water vapor values corresponding to each image and considering a continental aerosol model. Following Song et al. (2001), AOT for 550 nm was defined in a range from 0.01 to 2.0 with a step size of 0.01 for each iteration. AOT was set at 550 nm when modeled reflectance in bands 1 and 3 matched the reflectance obtained via the empirical relationship indicated in the equation. Values of AOT for images with different dates ranged from 0.048 to 0.96, with generally higher values for the August time series. AOT estimations, together with ozone and water vapor values, were included as input for the 6S code to atmospherically correct the 28 selected images.
3.5. Topographic correction
Another source of artificial noise is the modification of illumination conditions by topography. Although the effect of topography on illumination conditions is complex, even including reflections from adjacent slopes (Proy et al., 1989), it is usually simplified for the purpose of analysis: shaded areas show less than expected reflectance, whereas sunny areas show the opposite pattern.
Various methods are available to topographically correct satellite imagery (Civco, 1989; Pons and Solé, 1994; Gu and Gillespie, 1998). In this paper, we tested two of the most commonly used methods: the first assumes a lambertian behavior of the surface and the second considers non-lambertian effects.
The illumination conditions can be modeled following the cosine law of spherical geometry (Civco, 1989):
where is the cosine angle between the solar incident angle and the local surface normal, is the solar zenith angle, is the zenith angle of the normal to the surface, is the topographic aspect angle, and is the solar azimuth angle.
Illumination conditions corresponding to the date and timing of satellite overpass were determined for each pixel of the Landsat images using the DEM described in Section 3.2 and by applying the formulation stated above. The DEM was used at a spatial resolution of 15 m to guarantee the reliability of the derived DTMs (Digital Terrain Models) in the topographic correction procedure. To calculate the spatially distributed values of and , we employed the Geographic Information System (GIS) MiraMon.
Based on a lambertian assumption, the reflectance of a horizontal surface can be determined by Eq. 8.
Lambertian models usually overestimate radiance for the high incidence angles of slopes facing away from the sun (Vincini and Frazzi, 2003). Simple models have been proposed to take into account the non-lambertian properties of the surface. These models are based on Minnaert’s theory, in which a constant K, derived empirically for each band, enables the characterization of the non-isotropic conditions of the surface. Minnaert’s original proposal has been modified in a number of studies (Colby, 1991; Teillet et al., 1982; Vinzini and Frazzi, 2003); among these methods, the C-correction (Teillet et al., 1982) has shown superior performance because it successfully retains the spectral characteristics of each band and significantly reduces reflectance variability for homogeneous vegetation classes (Riaño et al., 2003). Following C-correction, the reflectance corresponding to a horizontal surface can be obtained as follows:
where is obtained empirically from the entire image following , where . This procedure removes the common dependence of to more efficiently than other methods (Vinzini and Frazzi, 2003).
3.6. Relative radiometric normalization
The commonly assumed simplifications employed in atmospheric and topographic corrections mean that they usually fail to completely remove non-surface noise. To obtain improved temporal homogeneity of satellite imagery, it is common to apply relative normalization between images (e.g., Yuan and Elvidge, 1996; Coppin and Bauer, 1994; Tokola et al., 1998). Relative radiometric normalization is sometimes used as the sole correction procedure (Caselles and López-García, 1989), and is usually preferred to atmospheric correction for the purpose of change detection (Olsson, 1995; Janzen et al., 2006). Thus, relative normalization can be applied directly to DNs, radiance, TOA reflectance, or surface reflectance values.
Relative normalization is generally based on a linear comparison of image statistics for images obtained on different dates. One of the images, commonly the most recent or least affected by atmospheric effects, is considered as the reference to which the rest of the images are adjusted. Among the methods proposed for relative normalization, linear regression is the most commonly used and widely recommended (Yuan and Elvidge, 1996). The validity of the assumption of a linear relationship between the reference image and the image to be normalized has been confirmed in several studies (Schott et al., 1988; Caselles and López, 1989; Hall et al., 1991; Heo and FitzHugh, 2000). This assumption greatly simplifies the normalization process, based on matching the reflectance values of the image to be normalized to the reference image as follows:
where a and b are the linear regression parameters.
The most critical decision to be made is the sampling of targets/pixels in estimating the regression parameters, as it is necessary to identify constant reflectors between dates and to assume that reflectance differences in these targets are due to non-surface noise. Several procedures can be followed in selecting the constant targets, of which the most widely used is based on a visual selection of non-variant targets [Pseudo-Invariant Features (PIFs)] (Schott et al., 1988). These areas should contain minimum amounts of vegetation, be located in relatively flat areas to minimize the effects of illumination differences, and cover a wide range of reflectance values (from dark to bright areas) (Eckhardt et al., 1990). Targets such as sand, asphalt, and water are commonly selected as PIFs (Caselles and García, 1989; Coppin and Bauer, 1994). In this paper, we followed a PIFs normalization (Schott et al., 1988) based on the selection of 280 invariant pixels to be used for the calculation of a and b in each TM and ETM+ band and image. PIFs were identified visually in urban areas and areas of asphalt, water, dark vegetation and sand, for which an invariant reflectance can be assumed over the study period. The last images for the series of March (13/03/2007) and August (01/08/2006) were considered as references, and the rest were matched to these images.
We also used a simple radiometric correction method termed Temporally Invariant Cluster (TIC) (Chen et al., 2005), based on the visual identification of centers of high frequency in density scatterplots. The method assumes that these centers correspond to areas in which there are no changes between the reference image and the image to be normalized (Chen et al., 2005). This method is less time-consuming than PIFs because although identification of the TIC centers is not automatic, it is faster than the identification of PIFs. Using this procedure, the values of only two points are needed to obtain the regression line and the a and b coefficients that intersect the centers of high density.
Finally, we used an automatic procedure for relative normalization. Automatic techniques have the advantage of being less time-consuming than other techniques and maintaining criteria for target selection. Pseudo-Invariant Feature Regression (PIFR) (Du et al., 2001, 2002) is an automatic method that has provided excellent results in the relative normalization of Landsat images (Janzen et al., 2006; Olthof et al., 2005; Paolini et al., 2006). The PIFR method applies Principal Component Analysis (PCA) between each pair of images to obtain the PIFs. The first component represents a least-square regression between overlap areas and the second residual variation. It is recommended to remove outliers prior to model calculation, as they may affect model performance (Heo and FitzHugh, 2000). Although the problem of outliers is minimized in our dataset (as cloud cover was removed previously), the PIFR method removes outliers iteratively on the second principal component to discard changed pixels. We used an iterative process in which the standardized residuals obtained in the first run where used to remove those pixels that exceeded a certain threshold. An initial threshold of +/–1.28 (10% probability according to the normal distribution) was set to remove outliers after the first run. Thresholds corresponding to +/–1% probability were added in subsequent runs to remove pixels to be excluded from the model. The percentage of variance explained by the model was assessed for each run and a minimum of 95% was chosen in selecting the model and the a and b parameters.
3.7.1. Calibration and relative normalization
To validate the calibration and relative normalization procedures, we identified Pseudo-Invariant (PI) validation pixels for which it is assumed that reflectivity did not vary over time. The identification of these pixels was independent regarding the PIFs used for relative normalization. These pixels were selected from topographically corrected reflectance values and prior relative radiometric normalization using a combined automatic and manual technique. First, the Coefficients of Variation (CV) in each of the six Landsat reflective bands was calculated for each pixel in the time series for March and August. Average CV values were obtained in the different bands of the two series, with the aim of obtaining a unique image that represents the average variability of each pixel independently of the spectral band and season. We considered as PI validation pixels those with average CV values below 0.05. We sampled 300 such pixels, including those in urban areas, dark coniferous forests, water bodies, areas of human infrastructure, and areas of bright sand. None of the selected pixels was subjected to modifications in land cover over the period of interest; reflectance values can be considered to have been stable over time. Temporal differences in reflectance values may be attributed to non-surface noise related to factors such as calibration and atmospheric influence.
The performance of TM and ETM+ calibration and the cross-calibration between ETM+ and TM images was assessed by comparing the TOA reflectance obtained for the PI validation pixels. The Mean Absolute Error (MAE) was calculated for each spectral band, considering as a reference the last images in the time series of August (08/01/2006) and March (03/13/2007). PI validation pixels were also used to assess the relative radiometric normalization methods.
3.7.2. Atmospheric correction
Atmospheric correction methods described in section 3.4 were applied to geometrically corrected and calibrated/cross-calibrated images. Atmospheric correction applied to these images was evaluated following two different procedures. The first approach involved a comparison of atmospherically corrected surface reflectance values with near-simultaneous MODIS reflectances. Although MODIS includes spectral bands with a number of differences compared with Landsat bands (in terms of bandpasses and spectral response), they are generally considered to be equivalent (Masek et al., 2006). Atmospheric correction of MODIS images can be performed with a greater degree of robustness than Landsat images due to the improved onboard capabilities of MODIS (Vermote et al., 1997b, 2002). MODIS surface reflectances were obtained from the MOD09GAV5 daily reflectance dataset (http://lpdaac.usgs.gov/modis/mod09ghkv4.asp) for the same region as that covered by the Landsat dataset. The Landsat–MODIS comparison period is 2001–2007, for which 5 MODIS images were near nadir. Landsat TOA and surface reflectance values obtained from DTA and 6S were aggregated to 500 m spatial resolution to match the MODIS reflectance dataset.
Atmospheric correction was also evaluated by comparing Landsat surface reflectance with field surface reflectances obtained simultaneous with the overpass of Landsat 5. For this purpose, we selected a homogeneous area characterized by xeric vegetation and very low degrees of herbaceous and shrub cover. The study area (0°40'W, 41°27'N) covers 28.2 ha and is a plane platformal area unaffected by topographic influence on illumination conditions. On 08/01/2006 and 03/13/2007 we used an Ocean Optics USB2000 field spectroradiometer to obtain a random sample of field reflectance measures (~200) simultaneous with the overpass of the Landsat-5 satellite (+/– 0.5 h). This spectroradiometer makes 10 IFOV observations over the 400–900 nm bandwidth with a 0.3 nm sampling interval. A spectralon reference panel was used to obtain reflectance values. To ensure a high signal-to-noise ratio, integration time was adjusted according to illumination conditions, and each measure was calculated as the mean of 20 individual spectra. The spectral curves were integrated to simulate the TM bands, taking into account the relative spectral response of each TM band (Teillet et al., 2001) and compared with the TOA reflectance and DTA and 6S surface reflectances obtained from Landsat imagery. We applied a one-way analysis of variance to determine if DTA and 6S surface reflectance are significantly different from the field surface reflectances. Tamhane’s pairwise comparisons test was used to produce the multiple comparisons where the variances are unequal. The significance threshold was set at 0.05.
3.7.3. Topographic correction
Topographic corrections were assessed via temporal comparisons of the reflectance trajectories before and after topographic corrections in an homogeneous Pinus halepensis forest located on slopes of 20 in the central Ebro Valley (semi-arid conditions, average annual precipitation = 320 mm, 41°43'N,0°28'W). The forest has the same density in south-facing and north-facing slopes. Ten pixels were selected in a south-facing slope of the forest, besides additional ten pixels in a north-slope facing to assess the effects of topographic correction procedures on the derived trajectories.