Multi-sensor NDVI data continuity: Uncertainties and implications for vegetation monitoring applications

Remote Sensing of Environment 100 (2006) 67 – 81 www.elsevier.com/locate/rse

Multi-sensor NDVI data continuity: Uncertainties and implications for vegetation monitoring applications Willem J.D. van Leeuwen *, Barron J. Orr, Stuart E. Marsh, Stefanie M. Herrmann The University of Arizona, Office of Arid Lands Studies, Arizona Remote Sensing Center, 1955 E. Sixth Street, Tucson, AZ 85719, United States Received 3 September 2004; received in revised form 22 July 2005; accepted 1 October 2005

Abstract Consistent NDVI time series are paramount in monitoring ecological resources that are being altered by climate and human impacts. An increasing number of natural resource managers use web-based geospatial decision support tools that integrate time series of both historical and current NDVI data derived from multiple sensors to make better informed planning and management decisions. Representative canopy reflectance and NDVI data were simulated for historical, current and future AVHRR, MODIS and VIIRS land surface monitoring satellites to quantify the differences due to sensor-specific characteristics. Cross-sensor NDVI translation equations were developed for surface conditions. The effect of a range of atmospheric conditions (Rayleigh scattering, ozone, aerosol optical thickness, and water vapor content) on the sensor-specific reflectance and NDVI values were evaluated to quantify the uncertainty in the apparent NDVI for each sensor. MODIS and VIIRS NDVI data are minimally affected by the atmospheric water vapor, while AVHRR NDVI data are substantially reduced by water vapor. Although multi-sensor NDVI continuity can be obtained by using the developed cross-sensor translation equations, the interactions between the spectral characteristics of surface vegetation and soil components, sensor-specific spectral band characteristics and atmospheric scattering and absorption windows will introduce uncertainty due to insufficient knowledge about the atmospheric conditions that affect the signal of the Earth’s pixels at the time of data acquisitions. Processing strategies and algorithm preferences among data streams are also hindering cross-sensor NDVI continuity. D 2005 Elsevier Inc. All rights reserved. Keywords: NDVI continuity; MODIS; AVHRR; Vegetation monitoring; Uncertainty

1. Introduction Spectral vegetation index data have been used to investigate the interactions between climate and landscape ecosystems, monitor the effects of floods, drought, fire and desertification, aid with land management and sustainability, investigate climate change and carbon sequestration, and assess natural resources, agricultural production and food aid (Myneni et al., 1997; Nemani et al., 2003; Seelan et al., 2003; Yang et al., 1998). The time series data provide a powerful tool to learn from past events, monitor current conditions (van Leeuwen et al., 2004), and prepare for future change. Comparison of current vegetation data records with historic long-term averages have

* Corresponding author. E-mail address: [email protected] (W.J.D. van Leeuwen). 0034-4257/$ - see front matter D 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2005.10.002

been used to support ecosystem monitoring (Orr et al., 2004), and help evaluate the impact of rising global temperature and CO2 levels and provide evidence of the impact of the 1989 and 1998 El Nin˜o events around the world (Nemani et al., 2003). Global, regional and local natural resource survey and assessment strategies are increasingly incorporating remotely sensed imagery to monitor current and historical vegetation dynamics and often rely on the combined use of multi-sensor vegetation data. A rising number of national, regional and local users and applications are employing geospatial tools that incorporate time series of spectral vegetation index data and other reference data such as roads, rivers and soil information for spatially and temporally explicit natural resource and agricultural monitoring. Although a variety of satellite sensor options are now available, practical considerations (i.e. data and processing costs, free distribution, the inherent tradeoff between spatial and temporal resolution, and the influence of cloud cover)

68

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

favor platforms that provide frequent images that are systematically processed into products useful for the assessment of vegetation. Two sensors among those that currently meet these criteria are the NOAA Advanced Very High Resolution Radiometer (AVHRR; Cracknell, 1997) and NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS; Justice et al., 1998). Since SPOT VEGETATION and Sea-viewing Wide Field-of-view Sensor (SeaWiFS) NDVI data are not freely distributed, these data are not included in the analysis. Among other products, both AVHRR and MODIS reflectance data are transformed into the Normalized Difference Vegetation Index (NDVI; Myneni et al., 1995; Tucker, 1979), the most widely used vegetation index. The following examples of web-based decision support systems make use of MODIS and AVHRR data to address operational needs of stakeholders and illustrate the importance of data continuity to the applications user community. RangeView (http://rangeview.arizona.edu/) serves the needs of natural resource managers in the southwest who need to assimilate disparate data sets and visualize time-varying spatial relationships. This system serves to enhance users’ understanding and management of the environment and natural resources as they vary across the landscape and change through time. The application has the flexibility to handle raster, vector and tabular data at varying spatial and temporal resolutions, such as biweekly vegetation indices derived from MODIS (Huete et al., 1999, 2002) and AVHRR (Eidenshink, 1992, in press) obtained from the Land Processes Distributed Active Archive Center (LP DAAC). Ranchers, natural resource managers and biologists use the tool to enhance ground-based monitoring efforts, identify anomalous trends, and compare current climate conditions with analogous conditions in past years for insight into potential future vegetation dynamics. The Production Estimates and Crop Assessment Division (PECAD) in the USDA Foreign Agricultural Service (FAS) is currently using some AVHRR and SPOT-VEGETATION, but mostly MODIS vegetation time series data in a decision support system (DSS) called CADRE (Crop Condition Data Retrieval and Evaluation) to provide global agricultural production outlooks and conditions that affect global food security commodity prices (Reynolds, 2001). PECAD’s webbased geospatial crop monitoring tool is called ‘‘CropExplorer’’ (http://www.pecad.fas.usda.gov/cropexplorer/). CropExplorer permits assessment of agricultural regions with AVHRR, MODIS and associated climate data. Analysts use these tools to produce yield and production estimates for key commodities. Both RangeView and PECAD provide analytical tools to assess vegetation trends using the data behind the AVHRR and MODIS images. As users become more dependent on these data for decision making, temporal continuity becomes more and more important. Sometimes unknown to the user, the temporal sequence of data provided comes from more than one satellite sensor. Since the lifetime of satellites like AVHRR and MODIS are expected to be about 5 years, new satellites are being put into polar orbits to maintain Earth resource observing capabilities. The scope of this research will be limited to 1-km

resolution data (or better) from AVHRR and AVHRR-like sensors that are free of charge and assist the previously described agricultural and natural resource monitoring applications. Since the 8-km NOAA-7 and NOAA-9 data available for 1982– 1989 are too coarse for the described user applications, the presented research results pertain mostly to the 1-km satellite data distributed by the LP DAAC. However, the presented results for AVHRR-11 are to a certain degree representative for AVHRR-7 and AVHRR-9 because of the similarities in their band spectral response functions. From 1989 to 1994, the NOAA-11 satellite was the source for 1-km AVHRR data. From 1995 to 2000, the NOAA-14 satellite provided AVHRR sensor coverage. From 2001 until the beginning of 2004, the NOAA-16/AVHRR data was used for vegetation monitoring. When that system became unstable, NOAA-17/AVHRR took over. All AVHRR data (1989 to present) obtained from the LP DAAC have been corrected for the effects of atmospheric water vapor but not for aerosols. Since 2000, the NDVI data derived from the Terra/MODIS satellite sensors are being regularly used because they provide higher spatial resolution, enhanced atmospheric corrections and more precise geo-registration. Clearly, issues of continuity across the AVHRR – MODIS time series have implications for the analysis of NDVI data and the monitoring of vegetation dynamics (Cihlar et al., 2004; Steven et al., 2003). Satellite-derived seasonal greenness/NDVI data have the potential to provide temporal indicators of the onset, end, peak and duration of vegetation greenness as well as the rate of growth, senescence and periodicity of photosynthetic activity (Reed et al., 1994; Yang et al., 1998). Through long-term data analysis, trends in such ecological indicators can also be assessed (Orr et al., 2004). The ultimate goal of the time series analysis of the historical biophysical data will be to improve our forecasting capability and make inferences about climate and drought conditions while improving critical vegetation mapping capabilities associated with critical needs such as production estimates, habitat assessments and the dynamics of fuels. The continuity of the NDVI values simulated for different satellites and sensors are being investigated to provide insight into the uncertainty of multisensor time series data. Our goal, therefore, was to quantify the uncertainty in a socalled ‘‘seamless’’ 1-km time series of the NDVI data for both AVHRR-11/14/16/17 and MODIS. The Visible/Infrared Imager Radiometer Suite (VIIRS) data characteristics are included in the simulation as much as possible to prepare for the transition to the VIIRS data stream (¨ 2009) at which time both AVHRR and MODIS data acquisitions will discontinue. Although applicable to other disciplines, this research has specifically focused on the use of different sensors for vegetation monitoring (AVHRR-14, AVHRR-16, MODIS, and VIIRS) in the context of their operational use by natural resource managers and those assessing agricultural production. Practical uses of these data to aid understanding our changing environment must be based on a quantitative appreciation of the uncertainty between different data sources and sensors. Although the NDVI data derived from the AVHRR sensors are

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

the heritage which forms our baseline, an increasing number of data sets derived from a range of sensors and diverse processing streams need to be evaluated by means of accuracy assessments to make judicious use of them. The results of this and other ongoing research related to data uncertainty derived from multiple discontinuous sensors should be of significant value if we are going to be able to make the best informed management decisions. The objectives of this research were to focus on: (a) quantifying the effect of multi-sensor band spectral response characteristics on spectral reflectance and NDVI values for a range of leaf area index (LAI) values; (b) examining the impact of atmospheric correction strategies applied to the AVHRR and MODIS NDVI data provided by the LP DAAC and used by a significant user community; and (c) providing translation equations to facilitate a seamless NDVI time series using successive sensors (AVHRR-14, AVHRR-16, MODIS, VIIRS). A detailed overview of the impact of Rayleigh scattering, ozone, aerosol and water vapor on the NDVI values for each sensor are presented to better inform the increasing number of users about the continuity issues associated with these data. This research is expanding on the work of Tanre´ et al. (1992) which described an atmospheric correction algorithm for NOAA-9 –AVHRR and the resulting uncertainties in reflectance and NDVI products. Although important, the scope of this research does not include the effects of calibration accuracy (Rao & Chen, 1995, 1999; Slater et al., 2001; Vermote & Kaufman, 1995) sensor degradation (Los, 1998), NDVI compositing techniques (Cihlar et al., 2004; Holben, 1986; van Leeuwen et al., 1999) and data sampling/filtering (Chen et al., 2004; Sellers et al., 1994), (sub)pixel-based cloud (Derrien et al., 1993), and snow mask accuracy (Hall et al., 2002) and quality assurance (Roy et al., 2002), geo-registration (Wolfe et al., 2002), differences in spatial resolution with view angle (van Leeuwen et al., 1997), effects of topography (Proy et al., 1989; Richter, 1997) and the bidirectional reflectance distribution function (BRDF; Schaaf et al., 2002) on the NDVI values (Roujean et al., 1992). BRDF and topography are coupled and complex problems that affect the NDVI. The

69

interpretation of the NDVI needs to take these effects into account. Fig. 1 illustrates that the non-water-vapor-corrected AVHRR-16 data exhibits much lower NDVI values than both MODIS- and water-vapor-corrected AVHRR data. Until 2002, LP DAAC only provided non-water-vapor-corrected AVHRR data. Since the water vapor correction is based on estimates of water vapor derived from the Total Ozone Mapping Spectrometer (TOMS) satellite or climatology, the question was asked if these corrections would provide a more accurate NDVI product or if they would introduce significant errors. It is noted that the accuracy of water vapor estimates is changing with the season and location and are hard to quantify. Climatology data will deviate from the actual water vapor concentrations dependent on dry and humid air weather patterns. Based on a visual inspection (Fig. 1), the water-vapor-corrected AVHRR data are closer to the MODIS data than the non-water-vapor-corrected AVHRR data. It should be noted that no aerosol correction was applied to the AVHRR data. The MODIS NDVI data seem to be slightly higher in the southwest of Arizona, but lower in many other regions compared to the water-vapor-corrected AVHRR data. Some clouds seem to have obstructed some of the northern regions in the MODIS image. These kind of images (Fig. 1), supplemented with images displaying the difference from long term average images and difference from previous period images, are now being provided on an operational basis (e.g. RangeView.arizona.edu). The question that must therefore be asked is if we can use the long term AVHRR NDVI average as a baseline to compute a difference from average for the MODIS NDVI time series. This paper addresses some of the issues that need to be solved to be able to do this within a defined set of confidence limits. 2. Data and methods Numerical experiments were conducted to evaluate the atmospheric uncertainty among NDVI data derived from multiple sensors and the possibility to derive cross-sensor translation equations to create a continuous NDVI time series no data water -0.2 0

State of AZ

0.1

0.2 0.3 0.4 0.5 0.6

MODIS 0

100 km

AVHRR Water vapor corr.

AVHRR Non-Water vapor corr.

0.7 0.8 0.9 1.0

Fig. 1. NDVI comparison between atmospherically corrected MODIS data (left; Terra; composite period: Nov 19 – Dec 2, 2001), water-vapor-corrected AVHRR-16 data (middle; NOAA; composite period: Nov 16 – 29, 2001), and non-water-vapor-corrected AVHRR-16 data (right; NOAA; composite period: Nov 16 – 29, 2001).

70

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

data set. The baseline data used in this research are the actual AVHRR and MODIS NDVI data archives from LP DAAC that are being used to create value added decision support tools for natural resource managers (van Leeuwen et al., 2004). Natural resource managers use images that represent the difference from average NDVI (1989 to present for AVHRR or 2000 to present for MODIS), difference from previous period NDVI values or difference from last year NDVI values in combination with pasture boundaries and allotments and other reference layers as ‘‘value added decision support tools’’ to make decisions about, for example grazing, cattle rotation, and wildfire management. Since it is not clear how cross-sensor NDVI data are best used in time series analysis, some insight into the uncertainties in multi-sensor NDVI data can be provided by a simulation experiment. A direct comparison of the data from multiple sensors cannot be done due to the lack

of overlapping data acquisitions before the year 2000. Although the MODIS and AVHRR data sets overlap (2000 and forward), they will be compared in more detail in a separate research activity with due attention to the different compositing intervals, algorithms and BRDF effects. This research used a canopy radiative transfer model in combination with an atmospheric radiative transfer model to generate the necessary data sets to meet our objectives. The following eight modeling and analysis steps were executed to derive the multi-sensor NDVI differences and translation equations (Fig. 2). 1. The sensor-specific band spectral response functions and representative hyperspectral soil and leaf reflectance data were used to generate red and near-infrared (NIR) soil, snow and leaf component spectral reflectance data.

Fig. 2. Data input and modeling overview to derive multi-sensor NDVI translation equations and uncertainty due to atmospheric correction inaccuracies.

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81 Table 1 Full width of the spectral bands at half of the maximum spectral transmissivity for the red and NIR bands of the AVHRR-11, AVHRR-14, AVHRR-16, AVHRR-17, MODIS and future VIIRS sensors Sensor

FWHM bandwidth (nm) Red

NIR

AVHRR-11 AVHRR-14 AVHRR-16 AVHRR-17 MODIS VIIRS

572 – 698 575 – 705 587 – 687 589 – 680 620 – 670 600 – 680

715 – 985 720 – 1000 733 – 986 734 – 988 841 – 876 846 – 885

2. The component spectra were used as an input to the canopy model which generates red and NIR canopy spectra for a range of LAI values and canopy backgrounds data (soil and snow backgrounds) for a range of LAI values (simulating different vegetation cover types). 3. The canopy reflectance data were used to compute surface NDVI values. 4. Analysis of multi-sensor reflectance and NDVI differences and development of cross-sensor NDVI translation equations. 5. Sensor-specific apparent reflectance data (using the red and NIR spectral response functions) were simulated for a range of surface reflectance values and atmospheric parameters (Rayleigh, ozone, water vapor, and aerosol optical thickness). 6. Polynomial multiple regression models were developed for the red and NIR apparent reflectance data as a function of the surface reflectance and atmospheric parameters for interpolation purposes. 7. The canopy apparent reflectance data were used to compute NDVI values affected by atmospheric effects.

SPECTRAL RESPONSE, REFLECTANCE

RED

71

8. Multi-sensor NDVI uncertainty are analyzed and crosssensor translation equations are developed. 2.1. Canopy reflectance and NDVI simulations The Scattering by Arbitrarily Inclined Leaves (SAIL) model (Verhoef, 1984, 1998) was encoded into Interactive Data Language (IDL) and used to simulate the spectral reflectance values for a range of LAI values (LAI = 0, 0.1, 0.3, 0.5, 1., 1.5, 2., 3., 4.0, 5.0) and spectral bands (red, NIR) at a constant nadir view angle and 40- sun angle. LAI is defined as the ratio of the one sided area of vegetation elements to the ground area. The surface NDVI was calculated from these red and NIR reflectance values. These SAIL-derived reflectance data were then used to quantify the effect of sensor-specific spectral response functions and atmosphere on the NDVI continuity. A spherical leaf angle distribution was used for all SAIL simulations. The sensor-specific spectral bandwidths are listed in Table 1. Spectral inputs to SAIL are simulated as realistically as possible using several hyperspectral bright and dark soil and snow data from the ASTER spectral library (http://speclib.jpl. nasa.gov/), and herbaceous and deciduous leaf spectra from Asner et al. (1998). These component spectra were adapted according to the specific spectral response functions of AVHRR, MODIS, and VIIRS (Fig. 3). The output from the SAIL model allows us to estimate the effects of spectral bandwidths of the different sensors on the NDVI without the disturbing effects of the atmosphere. The sensor-specific normalized spectral responses, the leaf and soil reflectance and leaf transmittance data were weighted and matched by interpolating the reflectance and transmittance data to the same wavelength as the spectral response functions.

NIR

MODIS-red MODIS-NIR AVHRR-16-red AVHRR-16-NIR AVHRR-14-red AVHRR-14-NI VIIRS-red VIIRS-NIR Alfisol Herbaceous leaf

1.0

0.8

0.6

0.4

0.2

0.0 0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

WAVELENGTH (µm)

Fig. 3. Spectral bandwidths (top) associated with the spectral response functions for the red and NIR bands of the AVHRR-14, AVHRR-16, MODIS, and VIIRS instruments. A soil and vegetation spectral signature are shown to illustrate the importance of the position of these band passes as they affect the reflectance response of a surface consisting of a mixture of contributing components e.g. soil, vegetation, water, litter. Alfisols are soils developed under temperate forests of the humid midlatitudes.

72

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

The effect of the spectral band response for MODIS and AVHRR were compared in terms of NDVI for herbaceous vegetation, deciduous vegetation, dark and bright soils and snow targets. Linear MODIS-to-AVHRR NDVI transfer functions (which can improve the synergistic use of the MODIS and AVHRR data) were developed based on the generation of simulated spectral reflectance data for a variety of surface targets including pure soils and snow, mixed vegetation – soil and mixed vegetation – snow. All comprised of LAI values between 0 and 5. Fig. 3 and Table 1 also show a comparison of the spectral band response functions and bandwidths using the full-widthhalf-maximum (FWHM) criterion. For VIIRS, currently only the approximate FWHM is known. 2.2. Atmosphere radiative transfer simulations The Second Simulation of the Satellite Signal in the Solar Spectrum (6S) atmospheric radiative transfer code (Vermote et al., 1997a, 1997b) was used to examine the magnitude of the atmospheric effects and to derive the uncertainties in reflectance values and associated NDVI for surfaces at sea level. Apparent at-sensor reflectance values were simulated for specified surface reflectance values using the 6S atmospheric radiative transfer model for four satellite sensors commonly used to produce NDVI values for vegetation monitoring. The atmospheric effects were simulated for the MODIS, AVHRR14/16 and VIIRS spectral bands using the sensor-specific normalized spectral response functions (Fig. 3). AVHRR-11 and AVHRR-14 as well as AVHRR-16 and AVHRR-17 have very similar spectral response curves (Fig. 3). Therefore, the spectral response curves from AVHRR-14 and AVHRR-16 were used as representatives in this analysis. The atmospheric radiative transfer simulations were carried out utilizing the spectral response functions of the red and NIR bands of the respective sensors. The top of the atmosphere or apparent reflectance (q app) can be expressed in terms of surface reflectance and atmospheric correction parameters. For a Lambertian, homogeneous and flat surface the following equation has been developed (Vermote et al., 1997a, 1997b):  qapp ðhv ; hs ; uÞ ¼ sg ðhv ; hs Þ qi ðhv ; hs ; uÞ þ stw ðhv ; hs Þqs =ð1  Sqs Þ where q app qi

ð1Þ

Apparent reflectance Intrinsic atmospheric reflectance due to Rayleigh and aerosol scattering qs Surface reflectance of a ground target sg Gaseous transmittance s tw Two-way atmospheric transmittance due to Rayleigh and aerosol scattering S Spherical albedo of the atmosphere due to Rayleigh and aerosol scattering. The following inputs were used to examine the atmospherically affected apparent reflectance data with respect to the

corrected surface reflectance values. The uncertainties in the atmospherically corrected reflectance data were estimated by simulating a range of realistic values for the 6S input parameter values. Four atmospheric parameters that are reported to significantly affect the top of atmosphere (TOA) readings were varied one at a time: & Rayleigh effect (atmospheric elevation at 0, 1, 2, and 3 km above sea level), & Ozone content (0, 0.1, 0.2, 0.3, 0.4 and 0.5 cm-atm), & Aerosol optical thickness (0, 0.1 and 0.5) using the continental aerosol model, and & Water vapor content (0, 1, 2, 3, 4, and 5 g/cm2). Reflectance and NDVI uncertainty estimates can then be derived by assigning measures of uncertainty in the elevation (Rayleigh), ozone content, water vapor content and aerosol optical thickness. Aerosol type, as a potential element of uncertainty (Miura et al., 2001), was not included in this analysis. A discussion of the effect of stratospheric aerosols of volcanic origin (e.g. Mt Pinatubo volcano eruption in 1991) can be found in Vermote et al. (1997a, 1997b). It should be noted that more absorption would take place if a biomass burning aerosol model was used. The levels selected for the aforementioned parameters are in the same range as the U.S. standard atmosphere profile (H2O = 1.42 g/cm3, O3 = 0.344 cm) and tropical atmospheric profile (H2O = 4.12 g/cm3, O3 = 0.247 cm). Also, since the water vapor content for a tropical environment will generally be higher than the U.S. standard water vapor parameters, these kinds of atmospheric corrections can easily introduce more errors rather than better accuracy. The ‘‘6S’’ surface reflectance data were taken to be representative of dark and bright soils, vegetation and snow cover. The red and NIR surface reflectance levels used were q k = 0.05, q k = 0.1, q k = 0.4, and q k = 0.9. 2.3. Simulating the effects of atmosphere on surface reflectance and NDVI Polynomial regression models were developed using 5th(Rayleigh), 10th- (aerosols) and 11th- (ozone and water vapor) order polynomial equations that allow for the computation of the apparent spectral reflectance (dependent variable) values for any canopy surface spectral reflectance value (independent variable) and any of the four specified atmospheric condition/parameter values (independent variable). The residual errors in the reflectance values were less than 0.0005, with correlation coefficients between 1 and 0.9999. The polynomial model is based on the following general functional relationship:
s qapp ð2Þ k ¼ f qk ; sv where q kapp

Are the apparent canopy reflectance values for the red or NIR wavebands and sensor-specific spectral response function.

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

q ks

73

Are the surface reflectance values for the red or NIR wavebands and sensor-specific spectral response function. sv Are the atmospheric parameter values for Rayleigh scattering based on elevation (km), Ozone concentration (in 0.001 Dobson units or cm), Water vapor column (g H2O/cm2-atm) and aerosol optical thickness at 550 nm. Thirty-two models were produced to represent the four sensors, two spectral bands, and four atmospheric parameters. The red and NIR surface reflectance data obtained from the SAIL model were than transformed into the corresponding apparent atsensor reflectance values, using these models. The apparent atsensor NDVI values were computed for the four sensors and range of LAI and vegetation and canopy backgrounds.

Estimations of the uncertainty in NDVI, dependent on sensor differences and partial atmospheric transmissivity corrections (s), were made by quantifying the difference (d) between the NDVIapparent and NDVIsurface. Uncertainties in surface NDVI values are derived using the following equation:

2.4. Cross-sensor surface NDVI translation equations

d2 ¼ NDVIapparent

The simulated surface reflectance data from the two vegetation and three canopy backgrounds are used to compute sensor-specific NDVI data. The relationships between these NDVIsensor data could then be evaluated to derive translation equations that could allow us to create a time series of NDVI data that are independent of the sensor-specific spectral band response functions. The simulated NDVIMODIS data is used as a baseline because the current MODIS data are more accurate for vegetation monitoring than the AVHRR data. The differences (d) between the NDVI values derived from each sensor are computed by:

NDVI uncertainty due to differences between different sensors and atmospheric characteristics can be estimated similarly:

d0 ¼ NDVIMODIS  NDVIsensor :

ð3Þ

2.5. Cross-sensor uncertainties in the apparent NDVI The sensitivity of the surface NDVI to atmospheric effects (Rayleigh, ozone, aerosols and water vapor) are evaluated and compared for each of the four selected sensors to provide quantitative information about the importance of the atmospheric corrections that are currently applied on an operational basis. Uncertainties in apparent at-sensor reflectance data due to inaccuracies in the estimation of atmospheric parameters were assessed using the simulated data. Thus the sensitivities of the NDVI data to uncertainties in atmospheric parameters were determined for each of the four sensors that are the focus of this research.

d1 ¼ NDVIsurface  NDVIapparent;s

ð4Þ

The uncertainty in apparent at-sensor reflectance and associated NDVI values are also simulated and analyzed based on the introduction of some realistic levels of error in the atmospheric parameters. The uncertainty in the NDVIsensor derived from each sensor due to an error in an atmospheric parameter can be estimated using the following equation: s1

 NDVIapparent

ð5Þ

s2 :

d3 ¼ NDVIapparent þ sensor 1  NDVIapparent þ sensor 2 :

ð6Þ

3. Results 3.1. Multi-sensor component spectra Red and NIR background reflectance values for different targets (herbaceous and deciduous leaves, two soil types, and snow) for each sensor resulting from the convolution analysis are summarized in Table 2. Similarly, hyperspectral reflectance data of representative backgrounds, a mollisol (an organic matter-rich dark soil), an alfisol (bright soil), and snow (Table 2) were used to derive red and NIR background reflectance values for each sensor. These data are the basic input to the canopy simulation model (SAIL). Substantial differences exist among the sensors for each of the component spectra. The largest differences (2% to 3%) occur most consistently between the MODIS and AVHRR-14 red leaf spectra. The smallest differences among the sensors can be observed for the soils and snow background spectra. Among the sensors, each of the NIR component spectra is within about 1% of each other.

Table 2 Multi-sensor red and NIR component spectra MODIS

AVHRR-16

AVHRR-14

VIIRS

Red Herbaceous reflectance Herbaceous transmittance Deciduous reflectance Deciduous transmittance Mollisol reflectance Alfisol reflectance Snow reflectance

0.0972 0.0365 0.0623 0.0344 0.0815 0.2672 0.9485

MODIS

AVHRR-16

AVHRR-14

VIIRS

0.3814 0.3637 0.4541 0.4208 0.1673 0.3572 0.8156

0.3785 0.3637 0.4492 0.4165 0.1675 0.3576 0.8093

0.3869 0.3716 0.4654 0.4320 0.1739 0.3623 0.8168

NIR 0.1072 0.0472 0.0695 0.0440 0.0773 0.2602 0.9521

0.1205 0.0610 0.0881 0.0631 0.0808 0.2642 0.9495

0.1073 0.0453 0.0743 0.0493 0.0792 0.2633 0.9518

0.3870 0.3710 0.4652 0.4312 0.1712 0.3587 0.8300

74

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

a

b AVHRR-16 0.25

0.4

LAI=0

NIRsensor surface reflectance

REDsensor surface reflectance

0.3 AVHRR-14 VIIRS

0.2 0.15 0.1 0.05 LAI=5 0 0

0.05

0.1

0.15

0.2

0.25

0.3

AVHRR-16 AVHRR-14

LAI=5

VIIRS

0.35

0.3

LAI=0 0.25 0.25

REDMODIS surface reflectance

0.3

0.35

0.4

NIRMODIS surface reflectance

Fig. 4. Cross-sensor red (a) and NIR(b) reflectance comparisons for a simulated canopy with LAI ranging from 0 to 5, a bright soil background and herbaceous leaf spectra.

3.2. Cross-sensor translation equations for surface NDVI values An example of the multi-sensor SAIL simulated reflectance and NDVI values are presented in Figs. 4 and 5. Fig. 4a shows that the relationship between the MODIS-based red reflectance values and the AVHRR-14/16 and VIIRS derived red reflectance values and are very similar. The lower range (around 0.05 or higher LAI values) of red surface reflectance values are slightly higher for the VIIRS and AVHRR-16 sensors and especially for the AVHRR-14 sensor. The position of the red MODIS band spectral response function is centered on the vegetation absorption spectra and increasing soil reflectance spectra, which, combined with its narrow width, cause the minor difference in red reflectance values. Fig. 4b shows that the relationship between the MODIS-based NIR reflectance values and VIIRS-derived NIR reflectance values are very similar. The AVHRR-14/16-derived NIR reflectance values are generally lower than the MODIS based NIR reflectance values. The position of the narrow NIR MODIS and VIIRS band spectral response functions are located at the

plateau of the vegetation NIR reflectance spectra and increasing soil reflectance spectra, while the wider NIR AVHRR-14 and AVHRR-16 band spectral response functions swerve into the visible portion of the vegetation spectra. This will cause lower NIR reflectance values for the AVHRR sensors than for the MODIS and VIIRS sensors. The combined canopy red and NIR reflectance behavior would indicate that the NDVIMODIS would be higher than the NDVIAVHRR-14, NDVIAVHRR-16 and NDVIVIIRS values as can be observed in Fig. 5. Fig. 5 also shows that the dynamic range of MODIS is larger than that of the other three sensors, AVHRR-14 having the lowest dynamic range. The NDVI –LAI plot for each sensor is divergent because the spectral response functions between sensors (e.g. AVHRR-14 and AVHRR-16) are different. When we simulate the red and NIR reflectance values based on the soil and vegetation spectral signatures (SAIL) and spectral band response functions, the NDVI become sensor-specific and cause a difference in NDVI for the same LAI values. Fig. 6 shows the multi-sensor NDVI simulations and the associated translation equations as a function the MODIS NDVI values. Although some NDVI 0.9

0.8

0.8

0.7

0.7 0.6

0.6 NDVIsensor

NDVI

0.5 0.4 0.3 0.2 0.1

AVHRR-16 AVHR AVHRR-14 VIIRS

MODIS AVHRR-16 AVHRR-14 VIIRS

0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 LAI

Fig. 5. Effect of four different (MODIS, AVHRR-16, AVHRR-14, VIIRS) red and NIR band spectral response functions on the NDVI values for herbaceous vegetation and a bright soil and a range of LAI values.

0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.1 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 NDVIMODIS

Fig. 6. Multi-sensor NDVI data and translation equations (NDVIsensor = aNDVIMODIS  b) taking into account the band spectral response functions of the selected sensors.

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

variability exists, a linear regression function was found to be adequate in translating the AVHRR and VIIRS NDVI data to MODIS equivalent NDVI data. The average difference between NDVI MODIS and NDVIAVHRR-16 is about 1.6% T 1.4%. The highest average difference exists between NDVIMODIS and NDVIAVHRR-14 and is about 4.3% T 2.3%. The lowest average difference exists between NDVIMODIS and NDVIVIIRS and is about 1.1% T 1.3%. If the translation to NDVIMODIS is not performed the differences among multi-sensor NDVI values are about 1% for low NDVI values for the AVHRR-14/16 and are almost negligible for VIIRS (Table 3). For the higher NDVI, values the difference between NDVIMODIS and NDVIVIIRS is 2%, between NDVIMODIS and NDVIAVHRR-16 about 3%, and between NDVIMODIS and NDVIAVHRR-14 about 7% (Table 3). 3.3. Multi-sensor NDVI — correction equations and atmospheric uncertainties The next sections evaluate the effects of Rayleigh, ozone, aerosols and water vapor on the NDVI. Oxygen absorption is incorporated in the calculations as it affects mostly the NIR bands. Since the Rayleigh correction of the NDVI is the only correction that is generally accurate on an operational basis, the ozone, aerosol and water vapor corrections are quantified using the apparent Rayleigh NDVI values as a baseline. 3.3.1. The effect of Rayleigh scattering on multi-sensor NDVI Rayleigh scattering mostly affects shorter wavelengths and thus the red waveband reflectance values will increase with respect to the surface reflectance values. However, both the red and NIR reflectance data for each sensor will be affected differently due to the sensor-specific spectral waveband characteristics. The average and standard deviation of red and NIR reflectance differences between the surface and apparent Rayleigh reflectance values are presented in Table 4. The differences for the red reflectance values are much larger than the NIR differences. Since the magnitudes of the red reflectance values are much smaller for vegetated areas, the effect of the increase in red reflectance values will substantially decrease the NDVI unless the atmospheric correction is performed. Table 5 shows the NDVIsurface to NDVIRayleigh correction equations, based on surface and apparent Rayleigh reflectance values for MODIS, AVHRR-16, AVHRR-14 and VIIRS, to be very similar. An example of the relationship between the simulated NDVIsurface and NDVIRayleigh is shown in Fig. 7. The average difference between the surface and apparent Rayleigh NDVIsensor values are about 6% T 3%. Rayleigh scatter is

75

Table 4 Average and standard deviation of the difference between surface reflectance and apparent Rayleigh reflectance values for red and NIR bands ; ; Sensor dred S d red dNIR S d NIR MODIS AVHRR-16 AVHRR-14 VIIRS

0.0134 0.0147 0.0125 0.0140

0.0051 0.0058 0.0073 0.0054

0.0019 0.0055 0.0062 0.0018

0.0008 0.0046 0.0050 0.0008

substantially reducing the NDVI values by about 2% in the lower dynamic range of the NDVI and by about 10% in the higher part of the NDVI dynamic range. An error in elevation of 1 km will cause an error in the NDVI of ¨ 0.5% T 0.3%. A 2km error in elevation would about double the error in the NDVI. The NDVIMODIS values are the least affected by Rayleigh scattering as the red spectral bandwidth is narrower than the other sensors. It should be noted that the spectral response function of VIIRS is based on preliminary specifications, and that this might still change. The use of the sensor correction equations are intended to approximate the magnitude of the Rayleigh scattering on the NDVI and allows us to translate the NDVIsensor data to the NDVIMODIS baseline. 3.3.2. The effect of ozone absorption on multi-sensor NDVI The red reflectance values for all sensors are affected by ozone absorption. The NIR reflectance values acquired by MODIS and VIIRS are not affected by ozone, but the NIR bands of the AVHRR sensors are, as they tend to include a small part of the visible spectrum. Typical seasonal ozone concentrations for low and high latitudes vary between 0.1 and 0.5 cm-atm. Opposite to Rayleigh scattering, ozone absorption in the red wavelength will slightly increase the apparent NDVI values. The relationship between the simulated NDVIOzone + Rayleigh and NDVIRayleigh is shown in Fig. 8 for the four sensors and an ozone concentration of 0.5 cm-atm. Table 5 shows the NDVIRayleigh + Ozone to NDVIRayleigh correction equations, based on the associated apparent reflectance values for MODIS, AVHRR-16, AVHRR-14 and VIIRS, to be very similar. Since the NDVI increases fairly linearly with higher ozone concentration, NDVI uncertainty estimates can be derived by estimating the uncertainty in the ozone concentration and interpolating the difference between ozone-corrected and -uncorrected NDVI values based on the correction equations in Table 5. The average increase of the apparent ozone-affected NDVIsensor (Ozone concentration of 0.5 cm-atm) values is about 3.5% T 1%. Ozone absorption is increasing the NDVI values by about 4.5% in the lower dynamic range of the NDVI and by about 1.6% in the higher part of the NDVI dynamic range. The average NDVI error

Table 3 NDVI translation equations based on surface observation taking into account the band spectral response functions of AVHRR-14, AVHRR-16, MODIS and VIIRS ; NDVIsensor = a NDVIMODIS  b R2 dNDVI S d NDVI d NDVI = 0.1 d NDVI = 0.9 NDVIAVHRR-16 = 0.9733 NDVIMODIS  0.002 0.9981 0.0159 0.0144 0.0047 0.0260 NDVIAVHRR-14 = 0.9318 NDVIMODIS  0.0071 0.9976 0.0426 0.0234 0.0014 0.0685 NDVIVIIRS = 0.975 NDVIMODIS + 0.0022 0.9985 0.0108 0.0130 0.0003 0.0203 ; The average difference (dNDVI ) and standard deviation of the differences (S d NDVI) are presented as well as examples of the absolute differences [d NDVI] between these sensors and the NDVIMODIS.

76

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

Table 5 NDVIRayleigh, NDVIOzone, NDVIAerosol, and NDVIH2O correction equations based on surface and apparent reflectance values for MODIS, AVHRR-16, AVHRR-14 and VIIRS ; ; R2 d1;NDVI S d 1, NDVI d NDVI = 0.1 d NDVI = 0.9 d2;NDVI S d 2, NDVI NDVI sensor, Rayleigh = a NDVI sensor, surface  b NDVIMODIS, R = 0.9046 NDVIMODIS, s  0.0034 NDVIAVHRR-16, R = 0.8958 NDVIAVHRR-16, s  0.0139 NDVIAVHRR-14, R = 0.8967 NDVIAVHRR-14, s  0.0112 NDVIVIIRS, R = 0.9003 NDVIVIIRS, s  0.004

0.9941 0.9925 0.9934 0.9933

0.0529 0.0665 0.0605 0.0547

0.0337 0.0362 0.0340 0.0346

0.0129 0.0243 0.0215 0.0140

0.0893 0.1077 0.1042 0.0937

0.0046 0.0058 0.0053 0.0048

0.0027 0.0029 0.0027 0.0027

NDVI sensor, Ozone = a NDVI sensor, Rayleigh  b NDVIMODIS, R+O3 = 0.9659 NDVIMODIS, R  0.0451 NDVIAVHRR-16, R+O3 = 0.9628 NDVIAVHRR-16, R  0.052 NDVIAVHRR-14, R + O3 = 0.9694 NDVIAVHRR-14, R  0.046 NDVIVIIRS, R + O3 = 0.9624 NDVIVIIRS, R  0.051

0.9999 0.9999 0.9999 0.9999

0.0293 0.0358 0.0332 0.0341

0.0091 0.0097 0.0077 0.0098

0.0417 0.0483 0.0429 0.0472

0.0144 0.0185 0.0185 0.0172

0.0060 0.0073 0.0067 0.0070

0.0018 0.0019 0.0016 0.0019

NDVI sensor, Rayleigh+Aerosol = a NDVI sensor, Rayleigh  b NDVIMODIS, R + A = 0.8366 NDVIMODIS, R NDVIAVHRR-16, R + A = 0.8216 NDVIAVHRR-16, R NDVIAVHRR-14, R + A = 0.8342 NDVIAVHRR-14, R  0.0008 NDVIVIIRS, R + A = 0.8354 NDVIVIIRS, R  0.0001

0.9828 0.9815 0.9904 0.9819

0.0775 0.0797 0.0699 0.0765

0.0480 0.0492 0.0437 0.0468

0.0173 0.0178 0.0167 0.0165

0.1471 0.1606 0.1500 0.1483

0.0155 0.0159 0.0140 0.0153

0.0096 0.0098 0.0087 0.0094

NDVI sensor, Rayleigh+H 2 O = a NDVI sensor, Rayleigh  b NDVIMODIS, R+H2O = 0.9991 NDVIMODIS, R  0.0019 1 0.0019 0.0003 0.0020 0.0027 0.0006 0.0001 NDVIAVHRR-16, R+H2O = 1.067 NDVIAVHRR-16, R  0.106 0.9995 0.0620 0.0149 0.0993 0.0457 0.0202 0.0045 0.9996 0.0566 0.0118 0.0873 0.0433 0.0184 0.0036 NDVIAVHRR-14, R+H2O = 1.055 NDVIAVHRR-14, R  0.0928 NDVIVIIRS, R+H2O = 0.9967 NDVIVIIRS, R + 0.0027 1 0.0008 0.0005  0.0024 0.0003 0.0003 0.0002 ; The average difference (d1;NDVI ) and standard deviation of the differences (S d 1, NDVI) are presented as well as examples of the absolute differences [d NDVI] between the ; surface and apparent Rayleigh, Ozone, Aerosol and water vapor affected NDVIsensor values. The average difference (d2;NDVI ) and standard deviation of the difference between apparent Rayleigh, Ozone, Aerosol and water vapor NDVI values (S d 2, NDVI) are shown for an elevation uncertainty of 1 km (sea level as reference), an ozone uncertainty of 0.1 cm-atm, an aerosol optical thickness uncertainty of 0.1, and when the uncertainty in atmospheric water vapor content is 1 g/cm2 (assuming an average water vapor content of 2 g/cm2).

due to an uncertainty of 0.1 cm-atm ozone concentration is 0.7% T 0.2% (shown in the last two columns of Table 5). 3.3.3. The effect of aerosol optical thickness on multi-sensor NDVI Both red and NIR reflectance values are affected by atmospheric aerosols. The red reflectance values generally increase in the lower range (smaller than 0.2) and decrease for the higher reflectance ranges (i.e. bright soils and snow). The

NIR reflectance values generally decrease across the range of reflectance values (Fig. 9). These results indicate that the NDVI will generally decrease due to atmospheric aerosols, which is shown in Fig. 10. The reduction in NDVI values due to aerosols is larger for higher aerosol optical thickness and NDVI values (Fig. 10). Similar patterns were found for the other sensors. Low NDVI values are much less affected by aerosols than high NDVI values. Since both red and NIR are independently affected by

0.9

0.9

MODIS 0.8

0.8

0.7

0.7

AVHRR-16 AVHRR-144

0.6

0.6

0.5

NDVIsensor, R+O3=0.5

NDVIRayleigh

MODIS M

0.4 0.3 0.2

VIIRS 0.5 0.4 0.3 0.2

0.1

0.1 0 -0.1 -0.1

0 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 NDVIsurface

Fig. 7. Relationship between the simulated NDVIsurface and NDVIRayleigh with error bars representing the range of values if the Rayleigh corrections were not applied.

-0. 1 -0.1

0

0.1 0. 2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 NDVIsensor,Rayleigh

Fig. 8. Multi-sensor NDVI data and ozone correction equations (Table 5) for ozone concentration equal to 0.5.

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81 0.9

(the last two columns of Table 5). NDVI dynamic ranges are lower for the AVHRR sensors than for MODIS and VIIRS, making the impact of aerosols even more pronounced among NDVIAVHRR sensors.

red

0.8 REFLECTANCE AOT=0.5+R

NIR 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.1 0.2

0.3 0.4 0.5 0.6 0.7

0.8 0.9

REFLECTANCERayleigh

Fig. 9. Red and NIR reflectance response for an aerosol optical thickness (AOT) of 0.5 and a range of LAI and canopy types for the MODIS sensor.

aerosol optical thickness (AOT), more variability among the representative NDVI values is observed resulting in lower correlation coefficients for the correction equations compared to the other atmospheric correction equations. Table 5 shows the NDVIRayleigh + Aerosols to NDVIRayleigh correction equations, based on the associated apparent reflectance values for MODIS, AVHRR-16, AVHRR-14 and VIIRS, to be very similar among the four sensors. The NDVI decreases fairly linearly with higher aerosol concentrations, and thus NDVI uncertainty estimates can be derived by estimating the uncertainty in the aerosol optical thickness concentration and interpolating the difference between AOT and uncorrected NDVI values based on the correction equations in Table 5. The average increase of the apparent AOT-affected NDVIsensor (AOT = 0.5) values is about 7.5% T 5%. AOT is decreasing the NDVI values by about 1.7% in the lower dynamic range of the NDVI and by about 15% in the higher part of the NDVI dynamic range. The average NDVI error due to an uncertainty of 0.1 AOT (Kaufman and Tanre´, 1996) is about 1.5% T 0.9% 0.9 0.8

AOT=0.1

0.7

AOT=0.5

1.0

Surface

REFLECTANCER+H2O=4

NDVIAOT+R

3.3.4. Cross-sensor NDVI uncertainties due to water vapor corrections Water absorption bands are avoided by the MODIS and VIIRS NIR spectral bands, but affect the AVHRR NIR spectral bands. Fig. 11 shows a comparison of the AVHRR-16 with the MODIS red and NIR bands for an atmospheric water vapor content of 4 g/cm2. The NIRAVHRR-16 reflectance values are substantially reduced compared to the NIRMODIS reflectance values. The red reflectance values tend to slightly increase due to water vapor in the atmosphere. AVHRR-14 behaves similarly to AVHRR-16, while VIIRS will behave similarly to MODIS. The apparent NDVIsensor values for MODIS and VIIRS are virtually unaffected by atmospheric water vapor, while for both AVHRR-14 and AVHRR-16 water vapor substantially decreases the surface NDVI values (Fig. 12). The NDVIAVHRR values decrease with increasing atmospheric water vapor content (Fig. 13), while NDVIMODIS and NDVIVIIRS decrease very little. Fig. 13 also shows that lower NDVIAVHRR values are more affected by atmospheric water vapor than the higher NDVIAVHRR values. The NDVI decreases non-linearly with increasing water vapor concentrations, making NDVI uncertainty estimates due to uncertainty in atmospheric water content a challenge. The correction equations in Table 5 provide a general overview of the NDVIsensor response to atmospheric water vapor. The average increase of the apparent water vapor (WV)-affected NDVIMODIS and NDVIVIIRS (WV = 4 g/cm2) values are less than 0.2% (Table 5). The average increase of the apparent WV-affected NDVIAVHRR-16 and NDVIAVHRR-14 (WV = 4 g/ cm2) values are 6.2% T 1.5% and 5.7% T 1.2%. WV decreases the NDVIAVHRR-16/14 values by about 9% to 10% in the lower dynamic range of the NDVIAVHRR-16/14 and by about 4.5% in the higher part of the NDVIAVHRR-16/14 dynamic range. The average NDVIAVHRR-16/14 error due to an uncer-

0.9

0.6

77

0.5 0.4 0.3 0.2 0.1 0

0.8 0.7

MODIS red MODIS NIR AVHRR-16 red AVHRR-16 NIR

0.6 0.5 0.4 0.3 0.2 0.1

-0.1 -0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

NDVIRayleigh

Fig. 10. Example of the effect of aerosol on the NDVIMODIS. The Surface NDVI values are plotted as a reference. The apparent NDVI (Rayleigh and AOT = 0.1) values are just below the 1:1 line.

0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 REFLECTANCERayleigh

Fig. 11. MODIS and AVHRR-16 red and NIR reflectance response to a water vapor column concentration of 4 g/cm2-atm (H2O = 4) and a range of LAI values.

78

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81 0.9

0.9

0.8

MODIS

0.7

AVHRR-16

MODIS AVHRR-16 AVHRR-14 VIIRS 1:1 Linear (VIIRS)

0.8 0.7

AVHRR-14 VIIRS

NDVI sensor, R+WV+AOT

NDVI R+H2O=4

0.6 0.5 0.4 0.3 0.2 0.1

0.4 0.3 0.2

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

NDVIRayleigh

Fig. 12. NDVIsensor response to a water vapor column concentration of 4 g/cm2atm (H2O = 4) and a range of LAI values.

tainty of 1 g/cm2 is about 2% T 0.5% (shown in the last two columns of Table 5). 3.4. Cross-sensor NDVI uncertainties due to water vapor and aerosols The previous sections have shown that the impact of AOT and water vapor cause the largest differences between the NDVIsensor values. The AVHRR sensors are similar to each other, but differ substantially with the MODIS and VIIRS. Fig. 14 provides some insight into the discontinuity between the AVHRR and MODIS-like sensors when the atmospheric correction is inaccurate or not applied. Aerosol corrections are applied to MODIS data, but not to the long-term time series of AVHRR data. When water vapor correction was applied to AVHRR data, accuracy was improved, but not to the level that can be derived from the MODIS bands. Fig. 14 shows that the differences between NDVIsensor are quite large, but are fairly linear, indicating that correction equations can be applied based 0.7 0.6 LAI=2

0.5 NDVI R+H2O

0.5

0.1

0 -0.1 -0.1

0.6

MODIS AVHRR-16 AVHRR-14 VIIRS

0.4 0.3

0

0.1

0.2

0.3 0.4 0.5 NDVIMODIS,R

0.6

0.7

0.8

0.9

Fig. 14. The Rayleigh NDVIMODIS as a function of the multi-sensor NDVI with Rayleigh, WV = 2 g/cm2 and AOT = 0.2.

on atmospheric reference data. The trend lines of these corrections or uncertainties are presented in Table 6. Under these moderately low atmospheric effects, the average NDVIMODIS, VIIRS uncertainty is approximately 3% T 1.8% while the NDVIAVHRR-16/14 uncertainties are at 10% T 1.5% (Table 6). The difference for the lower range of NDVIMODIS, VIIRS = 0.1 has uncertainties less than 1% and are 5% and 7% for the higher NDVI ranges, where NDVIMODIS = 0.9 and NDVIVIIRS = 0.9 respectively. The difference for the lower range of NDVIAVHRR-16/14 = 0.1 has uncertainties of about 10% and are 14% and 17% when the NDVIAVHRR-16 and NDVIVIIRS are equal to 0.9. 3.5. Overview of atmospheric effects on the NDVI for a range of LAI values Multiple atmospheric constituents are impacting the surface NDVI values from each sensor differently. Fig. 15 gives an example of the alterations that each of the major atmospheric variables make on the surface NDVI for a range of LAI values. Without atmospheric corrections, the NDVI values and derived biophysical parameters have large uncertainties. Seasonal atmospheric variability makes the consistency of NDVI time series data less stable and makes inter- and intra-annual comparisons of NDVI data less accurate. 3.6. Other factors affecting NDVI continuity/uncertainty

0.2 0.1

LAI=0

0 0

-0.1 -0.1

1 2 3 4 5 ATMOSPHERIC WATER VAPOR

6

(g/cm2)

Fig. 13. Example of NDVIsensor response to variable water vapor column concentrations of 0 – 4 g/cm2-atm (H2O = 4) and LAI values of 0 and 2 (m2/m2).

The atmospheric state and correction methods (e.g. aerosol optical thickness, water vapor content; Tanre´ et al., 1992), spectral signature of surface types (snow, soil, vegetation cover) and satellite spectral band response functions (Teillet et al., 1997; Trishchenko et al., 2002) are of great importance when comparing NDVI values derived from multiple sensors. The reported uncertainty in these factors can shift NDVI values by more or less than 20%. The MODIS atmospheric correction uncertainty is estimated to be between 1.5% and 5% of the

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81 Table 6 NDVI uncertainty or continuity equations when atmospheric water vapor content is 2 g/cm2 and an AOT of 0.2 ; NDVIsensor, R+AOT+H2O = a NDVIMODIS, R  b R2 d3;NDVI S d 3, NDVIMODIS, R+AOT+H2O = 0.9427 NDVIMODIS, R  0.0014 NDVIAVHRR-16, R+AOT+H2O = 0.9512 NDVIMODIS, R  0.0981 NDVIAVHRR-14, R+AOT+H2O = 0.9049 NDVIMODIS, R  0.0895 NDVIVIIRS, R+AOT+H2O = 0.914 NDVIMODIS, R + 0.0025

0.998 0.9965 0.9996 0.9975

0.0284 0.0913 0.0951 0.0300

NDVI

0.0179 0.0131 0.0174 0.0189

79

d NDVI = 0.1

d NDVI = 0.9

0.0071 0.1030 0.0990 0.0061

0.0530 0.1420 0.1751 0.0749

The correction equations are based apparent NDVI values for MODIS, AVHRR-16, AVHRR-14 and VIIRS, using NDVIMODIS as a reference. The average difference ; (d1;NDVI ) and standard deviation of the differences (S d 3, NDVI) are presented as well as examples of the absolute differences [d NDVI] between the surface and apparent water vapor and AOT affected NDVIsensor values.

least five cloud-free observations that are needed to invert a BRDF model (van Leeuwen et al., 1999). The MODIS compositing algorithm has some rules that constrain the selection of the NDVI to the view angle closest to nadir. Since NDVIAVHRR and NDVIMODIS are derived through different compositing techniques, some systematic bias in the NDVI can be expected; the maximum value NDVI composite technique will be biased towards higher NDVI composited data. When multi-sensor NDVI data are compared directly, the compositing time-step between processing streams varies, introducing timelag-related discontinuities in the NDVI which in turn can be exacerbated by the vegetation dynamics during the growing season. Accurate pixel based cloud mask and data quality will enhance the interpretation of the NDVI values. MODIS data provides a data quality field with information regarding the compositing technique, atmospheric parameters, snow and cloud cover, while AVHRR data just comes with a cloud mask.

reflectance factors (Liang et al., 2002; Petitcolin and Vermote, 2002) can result in large NDVI uncertainties especially for dense vegetation covers. Co-registration is important for pixel to pixel comparisons through time. MODIS geolocation is very accurate, ¨ 50 m at nadir (Wolfe et al., 2002), while AVHRR data are less accurate (¨ 1 km; Rosborough et al., 1994). The finer spatial resolution of MODIS is also contributing to improved co-registration (Wolfe et al., 2002). Calibration (spectral, radiometric) accuracy has been improved over time, benefiting the current MODIS data streams. The calibration uncertainty in the MODIS reflectance factors is estimated to be ¨ 1.8% (Guenther et al., 2002). The calibration results of the historical AVHRR records show much more degradation over time, making some of the intra-sensor NDVI data less stable through time (Cracknell, 1997). The equatorial time of overpass of Terra MODIS is about 10:30 h, while AVHRR16/14 times were 13:30 h. The AVHRR-17 overpass occurs in the morning around 10:30 h. Consequently, directional effects due to sampling (e.g. bi-weekly and seasonal) of surface reflectance values at different view and sun geometry will affect the NDVI (Cihlar et al., 2004). A compositing technique like the maximum value NDVI during a 10-day composite period will often choose the higher off-nadir NDVI values. BRDF corrections that can be applied to consistently use nadir reflectance observations are feasible (Schaaf et al., 2002) with an uncertainty of about 5% (Liang et al., 2002), but might introduce gaps in the data due to the limited availability of at

4. Discussion and summary Canopy and atmospheric radiative transfer models were used to quantify the effect of multi-sensor band spectral response characteristics on spectral reflectance and NDVI data for a range of LAI values. The development of multi-sensor translation equations was based on canopy simulations and this could facilitate a seamless NDVI time series from multiple sensors under ideal and comparable conditions. The impact of 0.8

0.8

AVHRR-16

0.7

0.7

0.6

0.6

0.5

0.5

0.4

surface Rayleigh O3 aerosol=0.5 WV=2 WV=2 AOT=0.2

0.3 0.2 0.1 0

NDVI

NDVI

MODIS

0.4 surface Rayleigh O3 aerosol=0.5 WV=2 WV=2 AOT=0.2

0.3 0.2 0.1 0

0

1

2

3 LAI

4

5

0

1

2

3

4

5

LAI

Fig. 15. The different impacts of atmospheric components on the NDVI for MODIS and AVHRR-16. Vertical arrows show the change in NDVI due to different atmospheric parameters, while the horizontal arrow illustrates the magnitude of the corresponding change in LAI values.

80

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81

atmospheric correction strategies applied to the AVHRR and MODIS NDVI data provided by the LP DAAC and used by a wide range of scientists and natural resource managers were quantified with assumptions regarding the accuracy or uncertainty in the atmospheric parameter retrievals. The different processing streams and associated uncertainties in the atmospheric characterization and corrections introduce uncertainties in multi-sensor NDVI continuity in attempts to translate NDVI values from one sensor to another. However, the translation equations and enumeration of the atmospheric accuracy present a realistic view of the magnitude of these NDVI uncertainties and we hope will lead to an improved and more judicious use of these tools in a range of ecological applications. The results indicate that current multi-sensor NDVI applications will benefit most if atmospheric corrections are adequately addressed and translation equations applied. Some multi-sensor discrepancies might be more complex, but can be overcome by using error reducing analysis techniques like data smoothing, and normalization (z-score). The MODIS equation can be used as a baseline to translate other NDVIsensor time series data. Use of the presented equations can help estimate uncertainty among multi-sensor NDVI. Some of the equations will need to be amended depending on which sensor one wants to compare. It should be noted that if MODIS data are taken as a baseline, MODIS overcorrection in the red wavebands can substantially increase the NDVI values. MODIS aerosol retrievals over bright targets are often less accurate and the derived NDVI values might thus be less certain (Hsu et al., 2004). Aerosols affect the red and NIR differently, while water vapor mostly affects the NIR band of AVHRR. Future NDVI continuity research will entail validation of this research by comparing coinciding actual reflectance and NDVI data sets acquired on the basis of daily and 16-day composited timesteps of MODIS and AVHRR using similar geometry and a range of vegetation types. Validation of these results needs to be conducted using temporally and spatially coinciding data. Acknowledgements This research was sponsored by a Synergy grant which was monitored by Raytheon Systems Company. The authors would like to thank Dr. Wout Verhoef for providing the SAIL code, Dr. Greg Asner for the use of hyperspectral component spectra and Dr. Eric Vermote for making the 6S code available. These data are distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at the U.S. Geological Survey (USGS) Center for Earth Resources Observation and Science (EROS) http://LPDAAC.usgs.gov. The valuable comments by the reviewers were much appreciated. References Asner, G. P., Wessman, C. A., Schimel, D. S., & Archer, S. (1998). Variability in leaf and litter optical properties: Implications for canopy BRDF model inversions using AVHRR, MODIS, and MISR. Remote Sensing of Environment, 63, 200 – 215.

Chen, J., Jo¨nsson, P., Tamura, M., Gu, Z., Matsushita, B., & Eklundh, L. (2004). A simple method for reconstructing a high-quality NDVI timeseries data set based on the Savitzky – Golay filter. Remote Sensing of Environment, 91, 332 – 344. Cihlar, J., Latifovic, R., Chen, J., Trishchenko, A., Du, Y., Fedosejevs, G., et al. (2004). Systematic corrections of AVHRR image composites for temporal studies. Remote Sensing of Environment, 89(2), 217 – 233. Cracknell, A. P. (1997). The Advanced Very High Resolution Radiometer. London’ Taylor and Francis 534 pp. Derrien, M., Farki, B., Harang, L., Le Gle´au, H., Noyalet, A., Pochic, D., et al. (1993). Automatic cloud detection applied to NOAA-11/AVHRR imagery. Remote Sensing of Environment, 46, 246 – 267. Eidenshink, J. C. (1992). The 1990 conterminous U.S. AVHRR data set. Photogrammetric Engineering and Remote Sensing, vol. 58(6). (pp. 809 – 813). American Society for Photogrammetry and Remote Sensing. Eidenshink, J. C. (in press). A 16-year time series of 1-km AVHRR satellite data of the conterminous United States and Alaska. Photogrammetric Engineering and Remote Sensing. American Society for Photogrammetry and Remote Sensing. Guenther, B., Xiong, X., Salomonson, V. V., Barnes, W. L., & James, Y. (2002). On-orbit performance of the Earth Observing System Moderate Resolution Imaging Spectroradiometer; first year of data. Remote Sensing of Environment, 83(1 – 2), 16 – 30. Hall, D. K., Riggs, G. A., Salomonson, V. V., DiGirolamo, N. E., & Bayr, K. J. (2002). MODIS snow-cover products. Remote Sensing of Environment, 83(12), 181 – 194. Holben, B. N. (1986). Characteristics of maximum-value composite images from temporal AVHRR data. International Journal of Remote Sensing, 7, 1417 – 1434. Hsu, C. N., Tsay, Si-Chee, King, M. D., & Herman, J. R. (2004). Aerosol properties over bright-reflecting source regions. IEEE Transactions on Geoscience and Remote Sensing, 42(3), 557 – 569. Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83(1 – 2), 195 – 213. Huete, A. R., Justice, C., & van Leeuwen, W. J. D. (1999) MODIS vegetation index algorithm theoretical basis document, Version 3. http://modis.gsfc. nasa.gov/data/atbd/land_atbd.html. Justice, C., Hall, D., Salomonson, V., Privette, J., Riggs, G., Strahler, A., et al. (1998). The Moderate Resolution Imaging Spectroradiometer (MODIS): Land remote sensing for global change research. IEEE Transactions on Geoscience and Remote Sensing, 36(4), 1228 – 1249. Kaufman, Y. J., & Tanre´, D. (1996). Strategy for Direct and indirect methods for correcting the aerosol effect on remote sensing: From AVHRR to EOS – MODIS. Remote Sensing of Environment, 55, 65 – 79. Liang, S., Fang, H., Chen, M., Shuey, C. J., Walthall, C., Daughtry, C., et al. (2002). Validating MODIS land surface reflectance and albedo products: Methods and preliminary results. Remote Sensing of Environment, 83, 149 – 162. Los, S. N. (1998). Estimation of the ratio of sensor degradation between NOAA AVHRR channels 1 and 2 from monthly NDVI composites. IEEE Transactions on Geoscience and Remote Sensing, 36(1), 206 – 213. Miura, T., Huete, A. R., Yoshioka, H., & Holben, B. N. (2001). An error and sensitivity analysis of atmospheric resistant vegetation indices derived from dark target-based atmospheric correction. Remote Sensing of Environment, 78(3), 284 – 298. Myneni, R. B., Hall, F. B., Sellers, P. J., & Marshak, A. L. (1995). The interpretation of spectral vegetation indices. IEEE Transactions on Geoscience and Remote Sensing, 33, 481 – 486. Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G., & Nemani, R. R. (1997). Increased plant growth in the northern high latitudes from 1981 – 1991. Nature, 386, 698 – 701. Nemani, R. R., Keeling, C. D., Hashimoto, H., Jolly, W. M., Piper, S. C., Tucker, C. J., et al. (2003, June 6). Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science, 300, 1560 – 1563.

W.J.D. van Leeuwen et al. / Remote Sensing of Environment 100 (2006) 67 – 81 Orr, B. J., Casady, G. M., Tuttle, D. G., van Leeuwen, W. J. D., Baker, L. E., McDonald, C. L., et al. (2004). Phenology and trend indictors derived from spatially dynamic bi-weekly satellite imagery to support ecosystem monitoring. Proceedings of the 5th Conference on Research and Resource Management in Southwestern Deserts — Biodiversity and Management of the Madrean Archipelago II: Connecting Mountain Islands and Desert Sea. May 11 – 15, 2004. Tucson (AZ). Petitcolin, F., & Vermote, E. (2002). Land surface reflectance, emissivity and temperature from MODIS middle and thermal infrared data. Remote Sensing of Environment, 83, 112 – 134. Proy, C., Tanre´, D., & Deschamps, P. Y. (1989). Evaluation of topographic effects in remotely sensed data. Remote Sensing of Environment, 30, 21 – 32. Rao, C. R. N., & Chen, J. (1995). Inter-satellite calibration linkages for the visible and near-infrared channels of the Advanced Very High Resolution Radiometer on NOAA-7,-9, and -11 spacecraft. International Journal of Remote Sensing, 16, 1931 – 1942. Rao, C. R. N., & Chen, J. (1999). Revised post-launch calibration of the visible and near-infrared channels of the Advanced Very High Resolution Radiometer (AVHRR) on the NOAA-14 spacecraft. International Journal of Remote Sensing, 20, 3485 – 3491. Reed, B. C., Brown, J. F., VanderZee, D., Loveland, T. R., Merchant, J. W., & Ohlen, D. O. (1994). Measuring phenological variability from satellite imagery. Journal of Vegetation Science, 5, 703 – 714. Reynolds, C. (2001). Input data sources, climate normals, crop models, and data extraction routines utilized by PECAD. Third International Conference on Geospatial Information in Agriculture and Forestry, Denver, Colorado, November 5 – 7. Richter, R. (1997). Correction of atmospheric and topographic effects for high spatial resolution satellite imagery. International Journal of Remote Sensing, 18(5), 1099 – 1111. Rosborough, G. W., Baldwin, D. G., & Emery, W. J. (1994). Precise AVHRR image navigation. IEEE Transactions on Geoscience and Remote Sensing, 32(3), 644 – 657. Roujean, J. -L., Leroy, M., & Dechamps, P. -Y. (1992). A bidirectional reflectance model of the earth’s surface for the correction of remote sensing data. Journal of Geophysical Research, 97(D18), 20455 – 20468. Roy, D. P., Borak, J. S., Devadiga, S., Wolfe, R. E., Zheng, M., & Descloitres, J. (2002). The MODIS land product quality assessment approach. Remote Sensing of Environment, 83(1 – 2), 62 – 76. Schaaf, C. B., Gao, F., Strahler, A. H., Lucht, W., Li, X., Tsang, T., et al. (2002). First operational BRDF, Albedo and Nadir reflectance products from MODIS. Remote Sensing of Environment, 83, 135 – 148. Seelan, S. K., Laguette, S., Casady, G. M., & Seielstad, G. A. (2003). Remote sensing applications for precision agriculture: A learning community approach. Remote Sensing of Environment, 88(12, 30), 157 – 169. Sellers, P. J., Tucker, C. J., Collatz, G. J., Los, S. O., Justice, C. O., Dazlich, D. A., et al. (1994). A global 1- by 1- NDVI data set for climate studies: Part II. The generation of global fields of terrestrial biophysical parameters from the NDVI. International Journal of Remote Sensing, 15(17), 3519 – 3545. Slater, P. N., Biggar, S. F., Palmer, J. M., & Thome, K. J. (2001). Unified approach to absolute radiometric calibration in the solar-reflective range. Remote Sensing of Environment, 77, 293 – 303.

81

Steven, M. D., Malthus, T. J., Baret, F., Xu, H., & Chopping, M. J. (2003). Intercalibration of vegetation indices from different sensor systems. Remote Sensing of Environment, 88(4), 412 – 422. Tanre´, D., Holben, B. N., & Kaufman, Y. J. (1992). Atmospheric correction algorithm for NOAA – AVHRR products: Theory and applications. IEEE Transactions on Geoscience and Remote Sensing, 30(2), 231 – 248. Teillet, P. M., Staenz, K., & William, D. J. (1997, July). Effects of spectral, spatial, and radiometric characteristics on remote sensing vegetation indices of forested regions. Remote Sensing of Environment, 61(1), 139 – 149. Trishchenko, A. P., Cihlar, J., & Li, Z. (2002). Effects of spectral response function on surface reflectance and NDVI measured with moderate resolution satellite sensors. Remote Sensing of Environment, 81(1), 1 – 18. Tucker, C. J. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment, 8, 127 – 150. van Leeuwen, W. J. D., Huete, A. R., Didan, K., & Laing, T. (1997). Modeling bi-directional reflectance factors for different land cover types and surface components to standardize vegetation indices. 7th Int. Symp. Phys. Measurements and Signatures in Remote Sensing, Courcheval (pp. 373 – 380). van Leeuwen, W. J. D., Huete, A. R., & Laing, T. W. (1999). MODIS vegetation index compositing approach: A prototype with AVHRR data. Remote Sensing of Environment, 69, 264 – 280. van Leeuwen, W. J. D., Orr, B., Marsh, S. E., Casady, G., Tuttle, D., Baker, L., et al. (2004). The use of MODIS products in a geospatial decision support tool for rangeland monitoring. In J. D. Greer (Ed.), Remote Sensing for Field Users Proceedings of the Tenth Forest Service Remote Sensing Applications Conference, Salt Lake City, Utah, April 5 – 9. American Society for Photogrammetry and Remote Sensing. Verhoef, W. (1984). Light scattering by leaf layers with application to canopy reflectance modeling, the SAIL model. Remote Sensing of Environment, 16, 125 – 141. Verhoef, W. (1998). Theory of radiative transfer models applied in optical remote sensing of vegetation canopies, Ph.D. dissertation, pp. 1 – 310. Vermote, E. F., & Kaufman, Y. J. (1995). Absolute calibration of AVHRR visible and near-infrared channels using ocean and cloud views. International Journal of Remote Sensing, 16, 2317 – 2340. Vermote, E., Saleous, N. E., Kaufman, Y. J., & Dutton, E. (1997a). Data preprocessing: Stratospheric aerosol perturbing effect on the remote sensing of vegetation: Correction method for the composite NDVI after the Pinatubo eruption. Remote Sensing Reviews, 15, 7 – 21. Vermote, E., Tanre´, D., Deuze´, J. L., Herman, M., & Morcrette, J. J. (1997b). Second simulation of the satellite signal in the solar spectrum: An overview. IEEE Transactions on Geoscience and Remote Sensing, 35-3, 675 – 686. Wolfe, R. E., Masahiro, N., Fleig, A. J., Kuyper, J. A., Roy, D. P., Storey, J. C., et al. (2002). Achieving sub-pixel geolocation accuracy in support of MODIS land science. Remote Sensing of Environment, 83, 31 – 49. Yang, L., Wylie, B. K., Thieszen, L. L., & Reed, B. C. (1998). An analysis of relationships among climate forcing and time-integrated NDVI of grasslands over the U.S. northern and central great plains. Remote Sensing of Environment, 65, 25 – 37.