898
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009
Towards a Generalized Approach for Correction of the BRDF Effect in MODIS Directional Reflectances Eric Vermote, Member, IEEE, Christopher O. Justice, and François-Marie Bréon
Abstract—We propose a new method for the estimation and the correction of directional effect on satellite derived surface reflectance time series. The correction assumes that the shape of the Bidirectional Reflectance Distribution Function (BRDF) varies much more slowly than the reflectance magnitude. The BRDF shape is quantified by two parameters R and V which account primarily for surface scattering and volume (canopy) scattering processes. Assuming a constant shape throughout the year permits an easy inversion of the two parameters that are used for a normalization of the reflectance time series to a standard observation geometry. Using this method, the corrected time series of data from the Moderate Resolution Imaging Spectroradiometer show much less high frequency variability than the original input data. A refinement of the method accounts for a slow variation of R and V , linearly with the Normalized Difference Vegetation Index (NDVI) accounting for vegetation dynamics, and produces an even better correction of directional effects and noise reduction in the time series. The nonlinear short-term variations in the reflectances, which are interpreted as “noise,” are strongly reduced, by a factor of two to four. The reduction in noise in the NDVI time series is less, as the index already implicitly corrects for the directional effects. Nevertheless, the noise in the NDVI time series is reduced by typically 30%–40%. Although the main output of processing using the method is the corrected reflectance, estimates of R and V are also generated. Global and regional map of these parameters show coherent patterns consistent with landcover. Index Terms—Directional reflectance, geophysical measurements, Moderate Resolution Imaging Spectroradiometer (MODIS), remote sensing, vegetation.
I. I NTRODUCTION
S
PACEBORNE remote sensing has been used for more than 20 years for monitoring the land surface. In particular, the combination of visible and near-infrared reflectances allows a quantitative estimate of vegetation cover and photosynthetic activity [1]. There are several difficulties and sources of uncertainty when using reflectance series to monitor vegetation dynamics. A first prerequisite is the identification of clouds in the field of view, followed by correction for atmospheric absorption
Manuscript received June 13, 2008; revised August 22, 2008. First published December 12, 2008; current version published February 19, 2009. E. Vermote and C. O. Justice are with the Department of Geography, University of Maryland, College Park, MD 20742 USA (e-mail:
[email protected];
[email protected]). F.-M. Bréon is with the Laboratoire des Sciences du Climat et de l’Environnement, Unité Mixte de Recherche Commissariat à l’Energie Atomique/Centre National de la Recherche Scientifique/Université de Versailles Saint Quentin-en-Yvelines, 91191 Gif-sur-Yvette, France (e-mail:
[email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2008.2005977
and scattering effects on the top-of-atmosphere reflectances. Finally, directional effects on the surface reflectances must be accounted for. Cloud detection and atmospheric correction have been the subject of intense research and are now used in operational processing (e.g., [2]). The quality of these atmospheric corrections is such that other effects, namely, the directional signatures of land surface reflectances, become dominant in the reflectance and vegetation index time series. This paper focuses on the correction of directional effects. Surface reflectance is highly anisotropic. For a given sun angle, the surface reflectance may vary by a factor of two in the near infrared and even larger in the visible part of the electromagnetic spectrum [3]. The largest reflectance variations are observed close to the backscattering direction, which is seldom observed by cross-track sensors [4], [5]. Nevertheless, even when this particular viewing geometry is avoided, the reflectance variations due to the observation geometry can be of the same order as the temporal variations that we are trying to detect [6]. The use of vegetation indices such as the Normalized Difference Vegetation Index (NDVI) partially avoids this difficulty, because the directional effects are similar in the visible and near infrared bands, and therefore the effect is reduced over the individual reflectances by ratioing the two bands [7]. The effect of surface anisotropy on remote sensed satellite data has been the subject of much research over the past 20 years, based both on theory using radiative transfer models, and on measurements either at the surface, or from aircraft or satellite ([8]–[13]). The surface reflectance is described by the Bidirectional Reflectance Distribution Function (BRDF), which is a function of the sun zenith angle θs , the view zenith angle θv , and both azimuths φs and φv with respect to a reference direction. In practice, for most applications, the azimuth variations only depend on the relative azimuth φ = φs − φv . For a simple modeling and correction of the reflectance directional effects, analytical models are needed. Several such models, both linear and nonlinear, have been published. Using multidirectional Parasol data at coarse resolution (6 km) over a large set of representative targets [14] show that simple models with only three free parameters permit an accurate representation of the BRDFs. The best results (low rms residual) were obtained with the linear Ross–Li–Maignan model, a version of the Ross–Li model that accounts for the Hot-Spot process [5]. The ability of simple linear models to reproduce natural target BRDFs opens the way for the correction of directional effects on reflectance time series. However, the question remains as to the choice of the BRDF model, i.e., the determination of its free parameters.
0196-2892/$25.00 © 2009 IEEE
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
VERMOTE et al.: TOWARDS A GENERALIZED APPROACH FOR CORRECTION OF THE BRDF EFFECT
899
The simpler approach is to assume a model, either globally, or based on simple criteria such as biome type or the vegetation index. Recently, analysis of the Polarization and Directionality of the Earth’s Reflectances (POLDER) BRDF database measurements, have shown that it is possible to assume a typical BRDF signature on a biome basis ([15]) and, therefore, apply a priori correction of the BRDF effect. This approach has been applied successfully on wide-swath data from polar-orbiting satellite systems ([16]–[18]). Another approach is to invert the BRDF model parameters over a period when the surface can be assumed to be stable. Such a method is currently used to process measurements from the Moderate Resolution Imaging Spectroradiometer (MODIS) with a composite period of 16 days [13]. One limitation of the method is the assumption of a stable target over the temporal compositing period. Another, which is more problematic in many cases, is that it requires several cloudfree measurements of the target over the compositing period. In addition, the observation geometry of these measurements may not be suitable to properly constrain the BRDF model. It is therefore desirable to design a method for the inversion of the BRDF model without the assumption of a stable target. In this paper, we offer such method. We do assume that the BRDF shape varies slowly in time but account for the fact that the reflectance magnitude may change rapidly. We show that the directional correction provides normalized reflectance time series that are much cleaner than the original input data, and that the BRDF model parameters have a geophysical meaning. II. D ATA AND M ETHODS A. Satellite Data In the following analysis, we use surface reflectance measurements derived from MODIS top-of-the-atmosphere observations. The measurements have been calibrated; corrected for atmospheric absorption and scattering; and filtered for cloud, cloud shadow, and snow and formed into the MODIS surface reflectance Collection 5 (C5) data set [2]. In order to process a manageable data set, we have used the Climate Modeling Grid (CMG) data at a resolution of 0.05◦ (≈5 km at the equator). The study is based on five years of morning observations from the MODIS sensor operating onboard the NASA Terra platform, from early 2000 to the end of 2004. As a specific example of the measurement time series as well as the normalized values, we have selected a site in southern Africa with a large amplitude seasonal variation of the reflectance. The site is located around Kaoma, Zambia, located at [14.79◦ S; 24.79◦ E] and is composed of savanna vegetation, consisting of a mixture of grasses and trees. The time-series signal was extracted from a single pixel (∼ 25 km2 ) from the MODIS CMG product. Fig. 1(a)–(c) shows the reflectance in the visible and near infrared, as well as the NDVI, for this particular site. B. BRDF Model and Notation The analysis of Parasol multidirectional data has shown that, among analytical BRDF models, the Ross–Li–Maignan model
Fig. 1. Time series of surface reflectance and NDVI derived from Terra/ MODIS reflectances at Kaoma (Zambia). For each subplot, the four groups of point correspond to different levels of processing. (Black) Bottom group is for the surface reflectance without directional correction. Next toward (red) the top are the corrected value where V and R have been estimated with the classical method (Section II-D). (Green) Next are the values corrected with the method described in Section II-E. Finally, the upper time series shows the values obtained when accounting for the variation of V and R with the NDVI. To generate this figure, the time series have been shifted by multiples of 0.1, 0.15, and 0.2 for band 1, band 2, and NDVI, respectively.
provides the best fit to the measurements. This model computes the reflectance as the sum of three terms ρ(θs , θv , φ) = k0 + k1 F1 (θs , θv , φ) + k2 F2 (θs , θv , φ) k2 k1 = k0 1 + F1 (θs , θv , φ) + F2 (θs , θv , φ) k0 k0
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
(1)
900
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009
where F1 is the volume scattering kernel, based on the Rossthick function, but corrected for the Hot-Spot process, and F2 is the geometric kernel, based on the Li-sparse reciprocal function. F1 and F2 are fixed functions of the observation geometry, but k0 , k1 , and k2 are free parameters. In the following, we will use V as k1 /k0 and R for k2 /k0 . The F1 and F2 functions are given in [5]. We recall them here for convenience F1 (θs , θv , ϕ) = cos t = m=
F2 (θs , θv , ϕ) =
1 + cos ξ m (t − sin t cos t − π) + π 2 cos θs cos θv
(2)
π 1 4 − ξ cos ξ + sin ξ 3π cos θs + cos θv 2
−1 ξ 1 × 1+ 1+ (3) − ξ0 3
where ξ is the scattering angle. A correction of the directional effect consists in putting the measurement in a standard observation geometry. In the following, the standard observation geometry is for the Sun at 45◦ from zenith, and the observation at nadir. The normalized reflectance is therefore computed as ρN (45, 0, 0) = ρ(θs , θv , φ) ×
In the following text, for simplicity and efficiency of writing, the N superscript always refers to the reference viewing geometry, and we drop the explicit notation to θs , θv , and φ. Equation (4) then is written as 1 + V F1N + RF2N . 1 + V F1 + RF2
Obviously, the linear variation hypothesis is valid only when the period is short. This is why, in practice, we compute the noise only based on the triplets when the two extreme dates are less than 20 days apart. Similarly, to avoid the noise estimate being driven by outliers, we first remove 5% of the triplets that show the largest difference before computing the sum in (6). Fig. 1(a)–(c) shows the reflectance and the NDVI time series for the Kaoma site. The absolute noise is 0.022, 0.044, and 0.019 for the visible, near IR, and NDVI, respectively. In relative terms, the noise is 20.7%, 18.4%, and 4.9%. Note that, for the reflectance time series, use of the term “noise” is incorrect, as there was no attempt to correct for directional effects. The values given above quantify the short-term variability of the observations acquired with varying geometries. For the NDVI, use of the term “noise” is appropriate, as the users expect a value that is indicative of the vegetation cover and state, without consideration of the observation geometry. D. Classic BRDF Inversion
1 + V F1 (45, 0, 0) + RF2 (45, 0, 0) . (4) 1 + V F1 (θs , θv , φ) + RF2 (θs , θv , φ)
ρN = ρ
1 NoiseR(y) = √ N −2
⎞2 ⎛ yi+2−yi − (day −day )−y y i ⎟ i+1 i n−2⎜ i+1 dayi+2−dayi ⎜ ⎟ . (6b) × ⎠ ⎝ yi+1 i=1
2 2 Δ + (tan θs tan θv sin ϕ)2 m 1 1 + cos θs cos θv
The aforementioned equation provides an absolute evaluation of the noise. Alternatively, one may compute a relative noise
(5)
The BRDF inversion technique used for the MODIS BRDF/Albedo product [13] consists in minimizing the difference between observed and modeled directional reflectances over a given time period, assuming no temporal variation in the surface reflectance during the observation period. Given a set of N reflectance measurements (ρi , i = 1, N ) the BRDF parameters (k0 , k1 , k2 ) are inverted by finding the minimum of the function N
k0 + k1 F1i + k2 F2i − ρi
2
.
(7)
i=1
C. Estimation of the Noise To quantify the impact of the correction on either the normalized reflectances or the vegetation indices, we compute an estimate of the noise. This estimate assumes a linear variation between two dates, a few days apart. Given three successive measurements, a triplet i, i + 1, and i + 2, the statistical difference between the center measurement and the linear interpolation between the two extremes quantifies the “noise” Noise(y)
2 n−2 i=1 yi+1 − yi+2−yi (dayi+1 −dayi )−yi dayi+2−dayi . = N −2
(6a)
The main hypothesis of this simple method is that the surface is stable during the compositing period, or rather, that the directional variations are larger than the temporal signal. In the standard MODIS processing, the compositing period is 16 days (orbit cycle of the Terra satellite). Although the surface is unlikely to be stable during the five years analyzed here, we have attempted an inversion of the BRDF parameters for the Kaoma site. The parameters are then used to estimate a normalized reflectance as in (4) which is shown in Fig. 1. Interestingly, the short-term variability in the reflectance is very much reduced, almost by a factor of four. This is an unexpected result which demonstrates that a simple correction of the directional effect can be very efficient. On the other hand, the NDVI does not show an improvement (Fig. 1, bottom). In fact, the normalized reflectances result in an
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
VERMOTE et al.: TOWARDS A GENERALIZED APPROACH FOR CORRECTION OF THE BRDF EFFECT
901
NDVI whose noise is larger than that of the original time series. As a ratio, the NDVI implicitly contains a directional correction because the effects are similar in the two bands. It is therefore more difficult to achieve an improvement to the NDVI than for the individual reflectances. E. Alternative BRDF Inversion The main limitation of the standard method described and applied above is the hypothesis that the target is stable during the period of synthesis. We therefore propose an alternative method that accounts for the fact that the target reflectance changes during the year, but assumes that the BRDF shape variations are limited. Another way of presenting the hypothesis is that k0 , k1 , and k2 vary in time, but k1 and k2 stay proportional to k0 . This hypothesis can be written as ρ(θs , θv , φ, t) = k0 (t) [1 + V F1 (θs , θv , φ) + RF2 (θs , θv , φ)] . (8) Having a set of N observations of surface reflectance acquired on successive days, we assume that the difference between the successive observations is mainly due to directional effects while the variations of k0 (t) are small. One can therefore estimate V and R through the minimization of the day-to-day variations of k0 (t). This minimization criterion can be derived from ρ(ti ) ρ(ti+1 ) ≈ . i i 1 + V F1 + RF2 1 + V F1i+1 + RF2i+1
(9)
One could therefore search for V and R that minimize the sum over time of the squared difference of the two terms above. This leads, however, to a system of equations that can only be solved through iteration. It then appears more efficient to note that (9) is equivalent to ρ(ti ) 1 + V F1i+1 + RF2i+1 ≈ ρ(ti+1 ) 1 + V F1i + RF2i . (10) One can therefore minimize the merit function as follows: 2 N −1 ρi+1 1+V F1i +RF2i −ρi 1+V F1i+1 +RF2i+1 M= . dayi+1 −dayi +1 i=1 (11) Which leads to a linear solution as shown below. Note that the denominator is a simple weighting that gives less weight to the observation couples that are distant in time. As such, the parameter inversion is mostly based on observation couples that are very close in time and which are, therefore, likely to exhibit directional effects, rather than changes in the surface. For simplicity of notation, we define the operators Δi d = dayi+1 − dayi + 1 √ Δi ρ = (ρi+1 − ρi )/ Δi d √ Δi ρF1 = ρi+1 F1i − ρi F1i+1 / Δi d √ Δi ρF2 = ρi+1 F2i − ρi F2i+1 / Δi d.
(12)
Fig. 2. V and R for five classes of NDVI for the Kaoma data set (see text, Section II-F). The error bars have been computed by running the inversion ten times, each time removing 10% of the data set at random. Squares are for band 1 while diamonds are for band 2.
We solve for R and V by the classic derivation of the merit function that is dM dM = =0 dV dR which leads to ⎛ N −1 i i ⎜ i=1 Δ ρF1 Δ ρF1 ⎜ −1 ⎝ N Δi ρF1 Δi ρF2 i=1
N −1 i=1 N −1 i=1
(13)
⎞
Δ ρF1 Δ ρF2 ⎟ ⎟⊗ V ⎠ R Δi ρF2 Δi ρF2 i
i
⎛
N −1
⎞
⎜ − i=1 Δ ρΔ ρF1 ⎟ ⎟. =⎜ −1 ⎝ N ⎠ i i − Δ ρΔ ρF2 i
i
(14)
i=1
An estimate of V and R can therefore be obtained through the computation of sums and a simple matrix inversion. Applied to the Kaoma data set, the method described above led to (V, R) estimates of (0.21, 0.78) in band 1, and (0.11, 1.48) for band 2. These values are close to the typical values for savanna-type surfaces reported by [15]. These values are then used for a further correction of the reflectance time series, as shown in Fig. 1. Compared to the Classic Method described in Section II-D, the new estimate of (V, R) shows a further decrease of the noise in both bands. Recall that the Classic Method leads to a large decrease in the short-term variability of the reflectance
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
902
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009
TABLE I STATISTICAL RESULTS FOR THE FIVE SITES THAT ARE ANALYZED IN DETAIL IN THE PAPER. THE FIRST LINES SHOW NOISE (ABSOLUTE AND RELATIVE) FOR THE ORIGINAL TIME SERIES AND THOSE OF CORRECTED TIME SERIES (SECTION II-E). THE FOLLOWING LINES SHOW THE D IRECTIONAL P ARAMETERS R ETRIEVED BY THE M ETHODS D ESCRIBED IN S ECTIONS II-D AND E
but that the impact was detrimental for the NDVI. The new method described in this section also has a very positive impact on the noise in the NDVI as well, as it decreases from 0.019 for the original time series to 0.011. Clearly, for the Kaoma data set, the (V, R) parameters estimated with the new method allow a much better correction of the directional effects. In the next section, we propose a further refinement accounting for a change of (V, R) with vegetation cover.
are used to define five classes of equal population. We then make an inversion of (V, R) for each of the five populations independently, and make the best linear fit on the results. This leads to a function (two coefficients) for each of the two bands and each of the two parameters, to estimate V and R as a function of the NDVI.
F. Final BRDF Inversion
A. Specific Sites
Although many empirical results demonstrate that the BRDF shape is more stable than its amplitude, one may expect that it varies with the vegetation cover. Indeed, an annual cycle of vegetation leaf up, growth, and senescence changes the surface structure, which has a direct effect on the BRDF. It seems, therefore, likely that V and R will vary seasonally with vegetation cover. In the BRDF formulation that we use, the F1 function is derived from the approximation of radiative transfer processes in an idealized vegetation canopy, while F2 models the BRDF of a rough surface. Therefore, as the vegetation grows, one may expect V to increase and R to decrease. To test this hypothesis, we segmented the Kaoma timeseries data set into five different classes of NDVI with equal population. We then inverted V and R, using the method described in the previous section, for each of these classes. The results are shown in Fig. 2: In band 1 (visible), the trend is not fully significant. On the other hand, the near infrared band clearly shows the expect increase of V and decrease of R with increasing NDVI. From results such as those shown in Fig. 2, one can derive a linear fit among the five retrievals and generate a linear function that gives V and R as a function of the NDVI. The full BRDF inversion process is therefore as follows: We first invert an estimate of (V, R) based on the previous method. This estimate allows the correction of the reflectances from which a normalized NDVI time series is obtained. These NDVI
To further evaluate the new inversion approach, we selected five additional sites representative of contrasting land covers:, an evergreen forest site in Brazil near Belterra [7.20◦ S; 54.40◦ W], a deciduous broadleaf forest in the eastern United States [37.75◦ N; 82.75◦ W], a broadleaf crop site in the Midwest United states [50.00◦ N; 103.20◦ W], croplands in India [25.50◦ N; 84.00◦ E], and the Kaoma savanna site [14.79◦ S; 24.79◦ E]. In all cases, the inversion performed very well with a strong reduction of the reflectance noise in both bands (see Table I). Fig. 3 shows the NDVI time series for both the uncorrected and the corrected data. The noise reduction is not as large as for the Kaoma site, but is nevertheless significant and allows the identification of features in the corrected time series that are difficult to see on the uncorrected data set. Fig. 4 shows for band 2, the R and V parameters for the five sites as a function of the NDVI. On this very limited data set, there seems to be a different behavior for the forest than for other biomes. Interestingly, the deciduous forest has very similar V and R parameters to those of the evergreen forest when the NDVI is at its maximum. The retrieved coefficients are reported in Table I.
III. R ESULTS FOR D IFFERENT L AND C OVERS
B. Global Analysis We applied the same method to the full global data set (i.e., five years of data at a 0.05◦ resolution). For the reflectance
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
VERMOTE et al.: TOWARDS A GENERALIZED APPROACH FOR CORRECTION OF THE BRDF EFFECT
903
Fig. 4. V and R for band 2 for five classes of NDVI for five different targets representative of Savanna, Evergreen Broadleaf forest, deciduous Forest, and crops. The error bars have been computed by running the inversion ten times, each time removing 10% of the data set at random.
Fig. 3. Original NDVI time series and corrected values for four targets of interest. For bottom to top, evergreen forest in Brazil near Belterra [7.20◦ S; 54.40◦ W], deciduous forest in the USA [37.75◦ N; 82.75◦ W], broadleaf crops in the USA [50.00◦ N; 103.20◦ W], and croplands in India [25.50◦ N; 84.00◦ E].
time series, the noise reduction is larger than three almost everywhere (not shown). This confirms the necessity for directional correction of reflectance time series. However, as shown
above, correction of NDVI data is not always beneficial. As most space-borne monitoring of the land surface utilizes the NDVI rather than the reflectance time series, it is necessary to quantify whether directional correction has a positive impact on the quality of the vegetation index time series. This is shown in Fig. 5(a)–(c). Fig. 5(a) shows a global map of the original NDVI noise, as defined in (6a) while Fig. 5(b) shows the same for the corrected time series. Most areas show an NDVI noise lower than 0.03. In the uncorrected NDVI, the smallest noise values are observed over desert areas, which are of little consequence for vegetation monitoring. The largest values are widespread with no specific region or biome standing out. The noise in the corrected time series [Fig. 5(b)] is clearly reduced. Regions with large noise {≈0.03) are limited to equatorial and some high-latitude areas. Most tropical and mid-latitudes regions have a noise of 0.02 or less. Fig. 5(c) shows that the noise reduction is very large (i.e., 40% or more) over widespread areas of the mid-latitude of both hemispheres. The improvement of the NDVI time series is very good over the United States, Mexico, Australia, and South Africa. On the other hand, there is a slight degradation of the time series over some equatorial areas, as well as high latitudes of the Northern hemisphere. The equatorial areas are characterized by large and persistent cloud cover while northern high-latitudes are similarly affected by large cloud cover, snow, and large sun zenith angles. This may explain the poor BRDF correction resulting from an improper inversion of its parameters. The south-eastern part of China also stands out as having a low improvement of the NDVI time series. Here, we believe that the culprit is the large and persistent atmospheric aerosol load:
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
904
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009
Fig. 5. (Top) Global map of the noise [as defined in (6a)] for the original NDVI time series and (middle) the same after directional correction of the reflectance. The lower panel is derived from the ratio of the noise shown in (a) and (b). Areas in white indicate an insufficient number of observations due to persistent cloud or snow cover. For (c), regions in black indicate an increase of the noise induced by the correction. For graphical purposes, the original noise maps at 0.05◦ resolution have been degraded to 0.3◦ resolution by taking the median value in each cell.
The residual of the atmospheric correction is larger than in any other part of the globe and obviously cannot be corrected by the BRDF correction procedure detailed in this paper, which assumes clean surface reflectance data as input. Despite these regions where the correction is questionable, the fact remains that most areas of the globe show a very significant improvement of the time series, which may improve the monitoring of the vegetation dynamics and the accurate determination of phenological parameters [18]. In particular, it
is interesting to note that most crop and grass surfaces, together with savannas and deciduous broadleaf forests, show a large reduction in the noise in the NDVI time series (see Fig. 6). IV. T OWARDS A G EOPHYSICAL I NTERPRETATION OF THE BRDF P ARAMETERS The primary goal of the study presented in this paper is the correction of reflectance and NDVI time series for a better
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
VERMOTE et al.: TOWARDS A GENERALIZED APPROACH FOR CORRECTION OF THE BRDF EFFECT
905
[19] to estimate aerodynamic roughness length Z0 . In this paper, they used a similar BRDF model [20] applied to the POLDER data set. To examine this relationship for the new inversion over desert areas, we use the data set of roughness length collected by Greeley et al. [21] for Namibia, and Death Valley and Lunar Lake USA and the data set collected by Marticorena et al. [20] for arid surface sites in southern Tunisia. Excluding sites with substantial vegetation cover, we compared the R parameter derived from this paper to the aerodynamic roughness length Z0 in Fig. 8. The relationship derived is close to the one derived by Marticorena et al. [19], i.e., R = 0.277 + 0.052 log10 (Z0 ).
Fig. 6. Histogram of the NDVI noise reduction at the original 0.05◦ resolution. The thick plain line is for the full data set. Thinner lines are for main biome classes are shown within the figure.
monitoring of vegetation cover. The BRDF model that is used for that purpose is semiempirical: The F1 and F2 functions are based on simplified modeling of he radiative transfer for an idealized canopy or a very simplified surface. Besides, there is no clear justification for having the reflectance as a linear combination of F1 and F2 together with an isotropic function. Despite these limitations, it is tempting to undertake a physical interpretation of the retrieved V and R parameters. For this objective, we favor the use of band 2 (near infrared) both because the surface reflectance is larger and the atmospheric perturbation smaller than in band 1. As explained in Section II-F, R and V are inverted as a function of the NDVI. In the following, we show the values retrieved for the upper NDVI class, i.e., for each pixel, the 20% of the measurements with the highest NDVI. Fig. 7(a) and (b) shows the V and R parameters, respectively, while Fig. 7(c) shows the mean NDVI for the NDVI upper range class. As expected, there is a good correspondence between the V parameters and the vegetation index: Desert areas have a generally low V parameter around 0.5, while vegetated areas have a larger V parameter between 1 and 2.5, with forested areas between 1 and 1.5, and short vegetation (grassland and crops) having the highest values of up to 2.5. Over desert areas, the topography appears to play an important role on the anisotropy of the surface, as shown by the variations of V over the Sahara Desert, associated with rugged relief. The R parameter shows high values (≈0.2) over both tropical forests, e.g., the Amazon, the Congo Basin, and mountainous regions of the World, e.g., the Himalayas, the Rockies, the Alps, the Rif Mountains, and the Andes. The lowest R values (≈−0.05) can be found for areas of little relief including the Great Plains of the U.S., the Indo Gangetic Plain and the Steppes of Russia, Mongolia and China, the Pampas, and the Canadian Shield. The Sahara Desert shows a wide range of variation in the R parameter. High values can be seen for the Atlas Mountains and the major mountain massifs of the Air, Ahaggar, the Tibesti, and the Adrar des Ifores. The R parameter can be associated with the surface roughness and was used by Marticorena et al.
(15)
V. S UMMARY AND C ONCLUSIONS Land surface reflectances show directional effects that must be corrected to enable proper interpretation of time series. The NDVI implicitly corrects a large fraction of these effects so that NDVI time series are much less affected than the individual visible and near infrared reflectances. Nevertheless, it remains desirable to correct the directional effects in both to ensure an accurate monitoring of vegetation dynamics. It had been shown before that typical directional shapes that may depend on the biome type and/or the vegetation cover, allow a first-order correction of the reflectance time series. This is because the directional signatures are much less variable than the reflectance magnitude. A second-order correction requires the use of a directional signature that depends on the observed target. However, the inversion of the target BRDF from satellite-measured reflectance time series is difficult, as the BRDF measurement requires several observations obtained over a given period, usually 10, 16, or 30 days, and for a large range of observation geometries. During the compositing period, the target may change due to vegetation dynamics, which also introduces uncertainties in the measured BRDF. Besides, the presence of cloud cover may limit the availability of the observations needed to constrain the BRDF model. In this paper, we show a method that allows estimation of the target BRDF even for dynamic targets. The hypothesis is that the directional signature shape varies more slowly than its magnitude. A first version assumes that the BRDF shape is constant over the period of synthesis, which may last over several years. A more complex version assumes that the two parameters that describe the BRDF shape vary linearly with the NDVI. Both methods can be applied to long time series and are, therefore, little affected by the lack of observations during a portion of the year. These methods were tested against five years of surface reflectance data derived from MODIS measurements at 0.05◦ resolution (CMG data set, collection 5). The performance of the correction method is quantitatively evaluated through the apparent noise in the time series, defined as the nonlinear shortterm variability. The results clearly demonstrate the validity of the approach. The noise in the reflectance time series is reduced by a factor three to five and is on the order of the residual of the atmospheric correction [22]. As expected, the improvement is less for the NDVI. Nevertheless, the noise in the NDVI time
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
906
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009
Fig. 7. Global map of the V and R parameters derived by the method shown in this paper applied to the time series of Terra MODIS band 2 CMG5 (2000–2004). V and R are shown for the highest NDVI class of each pixel. V is shown on the top, R on the middle, and the NDVI on the bottom.
series is reduced by typically 30% to 40%, except over desert and equatorial forest areas. The results showed systematic variation of the BRDF parameters related to seasonal vegetation dynamic. These variations are particularly clear in the near-infrared band, presumably because the larger reflectance amplitude and the smaller atmospheric effects than in the visible band, make it easier to measure. The R parameter, related to surface protrusions in the simple modeling, tends to decrease with an increase of the vegetation cover, while the V parameter which accounts for the canopy contribution, tends to increase.
As expected, both parameters show coherent patterns at the global scale. The V parameter appears to be related primarily to the vegetation amount and type, showing higher values for grassland and crops. Conversely, the R parameter shows low values for grassland and crops and higher values for forested areas. Additionally, the R coefficient is higher in desert areas with a high aerodynamics roughness from dunes and rocks and over major cities due to the shadowing effect or high buildings. Additional research is needed to investigate the relationship between the coefficients and other parameters such as leaf area index and canopy structure.
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
VERMOTE et al.: TOWARDS A GENERALIZED APPROACH FOR CORRECTION OF THE BRDF EFFECT
Fig. 8. Comparison of aerodynamic roughness length collected by Marticorena et al. [20] over Tunisia and Greeley et al. [21] over USA and Namibia with the R parameter derived from this paper.
Usual methods of BRDF inversion assume an invariance of the surface during the compositing period. As such, their application is limited by data availability and does not permit the monitoring of the vegetation with a period smaller than the compositing period. Here, we offer a method that explicitly accounts for the surface evolution with time. It can generate a product with the same frequency as the observations, which is particularly useful for monitoring rapid changes of vegetation cover, e.g., for agricultural areas. Note that the BRDF correction can be performed with the (V, R) coefficients that have been derived from previous observations, so that real time processing is possible. The corrected NDVI remains a vegetation index however, and there is no demonstration that it is better related to biophysical parameters. We did not attempt any validation against surface observations. Nevertheless, the reduced short-term variability implies that NDVI anomalies can be more easily and rapidly identified than with their uncorrected counterparts. R EFERENCES [1] C. J. Tucker, “Red and photographic infrared linear combinations for monitoring vegetation,” Remote Sens. Environ., vol. 8, no. 2, pp. 127–150, May 1979. [2] E. F. Vermote, N. Z. El Saleous, and C. O. Justice, “Atmospheric correction of MODIS data in the visible to middle infrared: First results,” Remote Sens. Environ., vol. 83, no. 1/2, pp. 97–111, Nov. 2002. [3] K. T. Kriebel, “Measured spectral bidirectional reflection properties of 4 vegetated surfaces,” Appl. Opt., vol. 17, pp. 253–259, 1978. [4] E. F. Vermote and D. P. Roy, “Land surface hotspot observed by MODIS over central Africa,” Int. J. Remote Sens., vol. 23, pp. 2141–2143, Jun. 2002. [5] F. M. Breon, F. Maignan, M. Leroy, and I. Grant, “Analysis of hot spot directional signatures measured from space,” J. Geophys. Res., vol. 107, no. D16, p. 4282, Aug. 22, 2002. [6] D. P. Roy, P. E. Lewis, and C. O. Justice, “Burned area mapping using multi-temporal moderate spatial resolution data—A bi-directional reflectance model-based expectation approach,” Remote Sens. Environ., vol. 83, no. 1/2, pp. 263–286, Nov. 2002. [7] B. Holben and R. S. Fraser, “Red and near-infrared sensor response to off-nadir viewing,” Int. J. Remote Sens., vol. 5, pp. 145–160, Jan./Feb. 1984. [8] H. Rahman, M. M. Verstraete, and B. Pinty, “Coupled surface-atmosphere reflectance (CSAR) model. 1: Model description and inversion on synthetic data,” J. Geophys. Res., vol. 98, no. D11, pp. 20 779–20 789, Nov. 20, 1993.
907
[9] D. S. Kimes, “Dynamics of directional reflectance factor distributions for vegetation canopies,” Appl. Opt., vol. 22, no. 9, pp. 1364–1372, May 1983. [10] D. S. Kimes, “Modeling the directional reflectance from complete homogeneous vegetation canopies with various leaf-orientation distributions,” J. Opt. Soc. Amer. A, Opt. Image Sci., vol. 1, no. 7, pp. 725–737, Jul. 1984. [11] X. W. Li and A. H. Strahler, “Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: Effect of crown shape and mutual shadowing,” IEEE Trans. Geosci. Remote Sens., vol. 30, no. 2, pp. 276–292, Mar. 1992. [12] M. Schoenermark, B. Geiger, and H. P. Roeser, Reflection Properties of Vegetation and Soil. Berlin, Germany: Wissenschaft und Technik Verlag, 2004. [13] C. B. Schaaf, F. Gao, A. H. Strahler, W. Lucht, X. W. Li, T. Tsang, N. C. Strugnell, X. Y. Zhang, Y. F. Jin, J. P. Muller, P. Lewis, M. Barnsley, P. Hobson, M. Disney, G. Roberts, M. Dunderdale, C. Doll, R. P. d’Entremont, B. X. Hu, S. L. Liang, J. L. Privette, and D. Roy, “First operational BRDF, Albedo nadir reflectance products from MODIS,” Remote Sens. Environ., vol. 83, no. 1, pp. 135–148, Nov. 2002. [14] F. Maignan, F. M. Breon, and R. Lacaze, “Bidirectional reflectance of earth targets: Evaluation of analytical models using a large set of spaceborne measurements with emphasis on the hot spot,” Remote Sens. Environ., vol. 90, no. 2, pp. 210–220, Mar. 30, 2004. [15] C. Bacour and F. M. Breon, “Variability of biome reflectance directional signatures as seen by POLDER,” Remote Sens. Environ., vol. 98, no. 1, pp. 80–95, Sep. 30, 2005. [16] C. Bacour, F. M. Breon, and F. Maignan, “Normalization of the directional effects in NOAA-AVHRR reflectance measurements for an improved monitoring of vegetation cycles,” Remote Sens. Environ., vol. 102, no. 3/4, pp. 402–413, Jun. 15, 2006. [17] F. Maignan, F. M. Breon, C. Bacour, J. Demarty, and A. Poirson, “Interannual vegetation phenology estimates from global AVHRR measurements: Comparison with in situ data and applications,” Remote Sens. Environ., vol. 112, no. 2, pp. 496–505, Feb. 15, 2008. [18] F. Maignan, F. M. Breon, E. Vermote, P. Ciais, and N. Viovy, “Mild winter and spring 2007 over Western Europe led to a widespread early vegetation onset,” Geophys. Res. Lett., vol. 35, no. 2, pp. L02 404.1–L02 404.6, Jan. 24, 2008. [19] B. Marticorena, P. Chazette, G. Bergametti, F. Dulac, and M. Legrand, “Mapping the aerodynamic roughness length of desert surfaces from the POLDER/ADEOS bi-directional reflectance product,” Int. J. Remote Sens., vol. 25, no. 3, pp. 603–626, Feb. 10, 2004. [20] B. Marticorena, M. Kardous, G. Bergametti, Y. Callot, P. Chazette, H. Khatteli, S. Le Hegarat-Mascle, M. Maille, J. L. Rajot, D. Vidal-Madjar, and M. Zribi, “Surface and aerodynamic roughness in arid and semiarid areas and their relation to radar backscatter coefficient,” J. Geophys. Res., vol. 111, no. F3, p. F03 017, Sep. 23, 2006. [21] R. Greeley, D. G. Blumberg, J. F. McHone, A. Dobrovolskis, J. D. Iversen, N. Lancaster, K. R. Rasmussen, S. D. Wall, and B. R. White, “Applications of spaceborne radar laboratory data to the study of aeolian processes,” J. Geophys. Res., vol. 102, no. E5, pp. 10 971–10 983, May 25, 1997. [22] E. F. Vermote and N. Z. Saleous, “Operational atmospheric correction of MODIS visible to middle infrared land surface data in the case of an infinite Lambertian target,” in Earth Science Satellite Remote Sensing, vol. 1. Berlin, Germany: Springer-Verlag, 2006.
Eric Vermote (M’95) received the Ph.D. degree in atmospheric optics from the University of Lille, Lille, France, in 1990. He is currently a Senior Research Scientist with the Department of Geography, University of Maryland, College Park. He is a member of the Moderate Resolution Imaging Spectroradiometer Science Team, the NASA National Polar-orbiting Operational Environmental Satellite System Preparatory Project Science Team, and the Landsat Data Continuity Mission Science Team, and is responsible for the atmospheric correction over land surfaces in the visible-to-middle infrared. His research interests cover radiative transfer modeling, vicarious calibration, atmospheric correction, and aerosol retrieval.
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.
908
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009
Christopher O. Justice received the Ph.D. degree from the University of Reading, Reading, U.K. Since 2001, he has been a Professor and Research Director with the Department of Geography, University of Maryland, College Park. He is currently a Team Member and Land Discipline Chair of the NASA Moderate Resolution Imaging Spectroradiometer (MODIS) Science Team and is responsible for the MODIS Fire Product; a member of the NASA National Polar-orbiting Operational Environmental Satellite System Preparatory Project Science Team; Cochair of the Global Observation of Forest Cover/Global Observation of Land Dynamics Fire Implementation Team, a project of the Global Terrestrial Observing System, a member of the Integrated Global Observation of Land Steering Committee; and a Task Leader for the Global Earth Observation System of Systems Agricultural Monitoring Task. He is Program Scientist for the NASA Land Cover Land Use Change Program. His current research is on land cover and land use change, the extent and impacts of global fire, and global agricultural monitoring.
François-Marie Bréon received the Ph.D. degree in atmospheric sciences and remote sensing from the University of Paris, Paris, France, in 1989. He then held several postdoctoral positions with Scripps Institution of Oceanography, San Diego, CA, the Centre National de la Météorologie, Toulouse, France, and the Meteorological Research Institute, Tsukuba, Japan, before joining the Laboratoire des Sciences du Climat et de l’Environnement, Gif-sur-Yvette, France. He has made specific contributions on the use of directional signatures as well as polarization for Earth monitoring. He is a Coinvestigator on several spaceborne Earth observing missions, including Polarization and Directionality of the Earth’s Reflectances, Parasol, Calipso, and Orbiting Carbon Observatory. He is the Chief Scientist of the Interactions Clouds Aerosols Radiations Etc thematic center for the study of aerosol, cloud, radiation, and their interactions. He has authored or coauthored more than 70 papers in the peer-reviewed literature. His main competences and interests are in clouds and aerosols, radiative transfer, and the remote sensing of the Earth atmosphere and surfaces.
Authorized licensed use limited to: University of Maryland College Park. Downloaded on August 13, 2009 at 16:46 from IEEE Xplore. Restrictions apply.