Real-time retrieval of Leaf Area Index from MODIS time series data

Remote Sensing of Environment 115 (2011) 97–106

Contents lists available at ScienceDirect

Remote Sensing of Environment j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / r s e

Real-time retrieval of Leaf Area Index from MODIS time series data Zhiqiang Xiao a,⁎, Shunlin Liang b, Jindi Wang a, Bo Jiang a, Xijia Li a a b

State Key Laboratory of Remote Sensing Science, School of Geography, Beijing Normal University, Beijing, 100875, China Department of Geography, University of Maryland, College Park, MD, 20742, USA

a r t i c l e

i n f o

Article history: Received 18 January 2010 Received in revised form 7 August 2010 Accepted 18 August 2010 Keywords: Real-time inversion Leaf Area Index SARIMA model Time series analysis Ensemble Kalman Filter MODIS

a b s t r a c t Real-time/near real-time inversion of land surface biogeophysical variables from satellite observations is required to monitor rapid land surface changes, and provide the necessary input for numerical weather forecasting models and decision support systems. This paper develops a new inversion method for the realtime estimation of the Leaf Area Index (LAI) of land surfaces from MODIS time series reflectance data (MOD09A1). It consists of a series of procedures, including time series data smoothing, data quality control and real-time estimation of LAI. After the historical LAI time series is smoothed by a multi-step Savitzky–Golay filter to determine the upper LAI envelope, a Seasonal Auto-Regressive Integrated Moving Average (SARIMA) model is used to derive the LAI climatology. Based on the climatology from the SARIMA model to evolve LAI in time, a dynamic model is then constructed and used to provide the short-range forecast of LAI. Predictions from this model are used with Ensemble Kalman Filter (EnKF) techniques to recursively update biophysical variables as new observations arrive. The validation results produced using MODIS surface reflectance data and field-measured LAI data at eight BELMANIP sites show that the real-time inversion method is able to efficiently produce a relatively smooth LAI product. In addition, the accuracy is significantly improved over the MODIS LAI product. © 2010 Elsevier Inc. All rights reserved.

1. Introduction Leaf Area Index (LAI) is defined as one-half the total green leaf area per unit of horizontal ground surface area (Chen & Black, 1992). LAI is an important vegetation biophysical variable that is widely applied for crop growth monitoring and yield estimation, land surface process simulation and global change studies. Currently, many methods estimate LAI from remote sensing data by utilizing the statistical relationship between LAI and spectral vegetation indices (Liang, 2004; Wang et al., 2007), physical model inversion (Qin et al., 2008; Gascon et al., 2004; Xiao et al., 2009), or other nonparametric methods (Smith, 1993; Fang & Liang, 2003; Kalacska et al., 2005). Each of these methods features both advantages and limitations. Since model inversion methods are physically based and can be adjusted for a wide range of circumstances (Kimes et al., 2000), radiative transfer models are increasingly used in the inverse mode to estimate LAI from remotely sensed data (Myneni et al., 1997; Jacquemoud & Baret, 1993). The inversion process is inherently problematic due to the multiple solutions possible and uncertainties within the model and measurements (Combal et al., 2002). Therefore, as much information as possible must be used in the inversion process to improve the accuracy of surface variable estimates. One approach is to use spatial ⁎ Corresponding author. Research Center for Remote Sensing and GIS, School of Geography, Beijing Normal University, Beijing, 100875, China. E-mail address: [email protected] (Z. Xiao). 0034-4257/$ – see front matter © 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2010.08.009

contextual information (Xiao et al., 2008), particularly multi-temporal information. Lauvernet et al. (2008) developed a multi-temporal patch ensemble inversion method to account for the spatial and temporal constraints in the inversion process. Koetz et al. (2005) proposed a method to estimate LAI based on multi-temporal remote sensing observations. The inversion was initially solved exclusively by radiometric information from each observation, and a simple semimechanistic canopy structure dynamic model (CSDM) was used to reconcile the radiometrically estimated LAI. Subsequently, the reconciled LAI was integrated into the inversion as a priori information. Xiao et al. (2009) developed a temporally integrated inversion method to estimate LAI from the time series MODIS reflectance data. Liu et al. (2008) proposed a data assimilation approach that combines the satellite observations of a MODIS albedo with a dynamic leaf model. These methods assume that all observations are available in advance and therefore suitable for analysis of historical data. In practice, the ability to track natural hazards and disasters (such as fires, floods and hurricanes) in real-time/near real-time depends upon a source of accurate and timely data. Remotely sensed data can generate biophysical variables that are highly accurate and frequent, which can used to carry out real-time/near real-time monitoring and warning. This process requires new algorithms that can exploit temporal observations for the real-time retrieval of these biophysical variables. Biophysical variables possess inherent temporal change rules and many dynamic models have been developed to characterize the

98

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

behavior of these variables over time. Using LAI as an example, many mechanical or semi-mechanical models have been proposed to describe LAI dynamics (Jamieson et al., 1998). However, these models contain many parameters and running them often requires extensive details of meteorological and/or agricultural knowledge. So, these models are often not suitable for real-time inversion on a regional or global scale. Over the past several decades, many time series analysis techniques have been developed to extract patterns of variability from the time series of satellite data sets for various applications (Hannachi, 2007; Chuvieco et al., 2008; Alhamad et al., 2007). These techniques are roughly divided into two groups: frequency domain approaches (Alhamad et al., 2007; Zhang et al., 2008) and time domain approaches (Alhamad et al., 2007; Röder et al., 2008). Although the methods described above feature certain advantages, they are primarily suitable for extracting phenological information and smoothing the time series, rather than offering effective predictive capabilities. More realistic details about trends and the seasonality of the time series would improve predictive capabilities. In this study, we developed a real-time inversion method to estimate LAI from MODIS data using coupled dynamic and radiative transfer models. In the following section, we introduce the radiative transfer model, followed by a discussion about how to apply multistep Savitzky–Golay filtering to eliminate contaminated pixels. We then use the Seasonal Auto-Regressive Integrated Moving Average (SARIMA) method to model the MODIS LAI time series and forecast future LAI values that are used to construct the dynamic model. Next, we outline how to use the Ensemble Kalman Filter (EnKF) technique to estimate real-time LAI from time series MODIS reflectance data. Section 3 shows the performance of this new method by using the actual values of LAI at several BELMANIP (Benchmark Land Multisite Analysis and Intercomparison of Products) sites (Garrigues et al., 2008). The final section provides brief discussions and conclusions.

The real-time inversion method takes advantage of the EnKF to recursively update LAI by coupling a canopy radiative transfer model with a dynamic model. A flow chart outlining this method is shown in

MODIS reflectance data (MOD09A1)

Quality control

Data with high quality

Radiative transfer models describe the relationship between canopy characteristics and reflectance, and many have been developed to obtain land surface biophysical parameters (Kuusk, 1994; Jacquemoud & Baret, 1990; Liang & Strahler, 1993). In our parameter retrieval, we chose the Markov Chain Reflectance Model (MCRM) developed by Kuusk (1995, 2001) as the forward model to simulate canopy reflectance. This model incorporates the Markov properties of stand geometry into an analytical, multispectral canopy reflectance model, which makes the model more flexible and applicable. The MCRM can calculate the angular distribution of canopy reflectance for a given solar direction in spectral wavelengths ranging from 400 to 2500 nm (Kuusk, 1995). MCRM model has already been used to estimate LAI, and the inversion practice shows that it produces good results (Fang et al., 2003; Fang & Liang, 2003). The parameters needed to run the forward MCRM are summarized in Table 1. Four parameters to which the remote sensing data are sensitive were identified as control variables to invert the MCRM model; they are denoted by an asterisk in Table 1. Values of the other parameters and effective ranges of the control variables are also displayed in Table 1.

LAI products

Simulated reflectance

Dynamic model

Radiative transfer model

……

Predicted LAI

time series analysis

Since 2000, MODIS LAI (MOD15A2) and surface reflectance products (MOD09A1) have been routinely produced. It is well known that the quality of MODIS reflectance data is affected by many factors, such as clouds, aerosols, water vapor and ozone. Though many of these effects are removed through atmospheric correction, the remaining ones are often significant. These effects, along with the MODIS algorithm (main/backup), affect the quality of the MODIS LAI product. Therefore, the MODIS LAI

Table 1 The parameters needed to run the MCRM.

Ensemble Kalman filter

LAI climatology

2.1. Radiative transfer model

2.2. Savitzky–Golay filtering to eliminate contaminated pixels

2. Methodology

……

Fig. 1. The existing multi-year MODIS LAI product is smoothed by a multi-step Savitzky–Golay filter to attain the upper LAI envelope, and a SARIMA model is used to derive the LAI climatology from smoothed LAI time series data. Then, a simple dynamic model is developed to evolve LAI over time and used to provide the short-range forecast of LAI. A radiative transfer model is used to link canopy characteristics (e.g., LAI) with MODIS reflectance. The EnKF technique is used to estimate realtime LAI by combining the predictions from the dynamic model and the MODIS reflectance data with the highest quality data. Important components of the new algorithm are described in details below.

MODIS LAI product

Fig. 1. Flow chart of LAI real-time inversion from MODIS time series data.

Parameters

Symbol

Values

Unit

Solar zenith angle Relative azimuth angle Viewing zenith angle Angstrom turbidity factor Ratio of leaf dimension and canopy height Leaf Area Index⁎ Markov parameter⁎

sza raa vza beta sl ul sz rnc eln thm SLW Cab wat mat vai rsl1 rsl2 rsl3 rsl4

0–90 0–180 0–90 0.1 0.15 0–10.0 0.4–1.0 0.9 0.0 90.0 100 0.4 150 99.6 1.2 0.05–0.4 −0.1–0.1 0.0 0.0

Degree Degree Degree – – m2/m2 – – Degree Degree g/m2 % of SLW % of SLW % of SLW – – – – –

Factor for refraction index Eccentricity of the leaf angle distribution Mean leaf angle of the elliptical LAD Leaf specific weight Chlorophyll AB content Leaf water content Leaf dry matter content Leaf structure parameter Weight of the first Price function⁎ Weight of the second Price function⁎ Weight of the third Price function Weight of the fourth Price function ⁎ Free parameters.

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

and surface reflectance values are fluctuated in a time series, especially in the growing season. One way to smooth data and suppress disturbances is to use a filter that is applied to a series of equally spaced data values: fi, i = 1, 2, …, K. The simplest type of filter replaces each data value fi by a linear combination gi of nearby values in a window: k

g i = ∑ c j fi j = −k

ð1Þ

+ j

where k is window width, which determines the degree of smoothing and is set to be 5 in our study, and, cj, is weight of data value fi + j. For the case of moving averaging: cj = 1/ (2k + 1). The moving average method preserves the area and mean position of a seasonal peak, but alters both the width and height. The idea of Savitzky–Golay filtering is to find filter coefficients, cj, so that the latter properties are preserved (Press et al., 1992). For each data valuefi, a quadratic polynomial function f(x) = c1 + c2x + c3x2 is fitted to all 2k + 1 points in the moving window using the least-squares method, and then gi is set as the value of that quadratic polynomial function at position i. Note that the atmosphere effects on surface reflectances generally cause increases in reflectance in the red band rather than in the NIR band, resulting in decreases in the retrieved LAI (Fernandes et al., 2003; Chen et al., 2006). To account for negatively biased noise, the fitting is done in multiple steps (Jönsson & Eklundh, 2002). In the first step, the coefficients are obtained by minimizing a merit function with the same weight for all data points. Data points below the model function of the first fit are considered less important, thus in the second step, the minimization is redone with the weight of the low data points decreased by a factor of 0.1 in our study. This multi-step procedure leads to a smoothed curve adapted to the upper envelope of the LAI values in the time series (shown in Fig. 2). In our study, the smoothed time series curve is used to derive LAI climatology. In this study, the contaminated reflectance data was filtered according to Normalized Different Vegetation Index (NDVI), and calculated according to the reflectance values in the red and NIR band. For the same reason as LAI, the atmospheric effects on surface reflectances generally cause negatively biased noise within the NDVI values. The multistep Savitzky–Golay filtering procedure is first used to calculate the upper envelope of the NDVI profile, denoted by NDVI_Env(i), i=1, 2,…, 46. If the NDVI profile, denoted by NDVI(i), satisfies the following condition, the reflectance data is deemed to be contaminated by clouds or other factors. Otherwise, it is considered to be cloud-free and of high quality: jNDVI ðiÞ−NDVIX EnvðiÞjN α × NDVIX EnvðiÞ

ð2Þ

where α is a threshold and is set to 0.2 in our study.

Fig. 3 depicts the MODIS reflectance time series in the red and NIR band at the center pixel of the Bondville site (detailed information about the site is described in Section 3.1). There are abrupt rises at days 153, 177, 201 and 233 during the crop growing season and the reflectance values at day 169 are missing. The NDVI temporal profile, calculated according to the reflectance values in the red and NIR band, is also shown in Fig. 3. It is obvious that the abrupt drops in the NDVI temporal profile correlate with the abrupt rises in the reflectance temporal profile. The method described above is applied to check the MODIS time series reflectances at the Bondville site, and the data points assessed as contaminated are marked by circles in Fig. 3. These data are filtered out and not utilized to update biophysical variables. 2.3. Dynamic model 2.3.1. LAI climatology based on the SARIMA model Generally, LAI climatology is determined by the averaging of LAI over a period of time. For example, we can compute LAI climatology from the highest quality multi-year MODIS LAI data (quality control values less than 32) as follows: clim

LAIt

=

1 M ∑ LAIt ðyearÞ M year = 1

ð3Þ

is the average of LAI values from multi-year MODIS LAI where LAIclim t data (after quality control); LAIt(year) represents LAI data for each year. Clearly, this is a very simplified way to obtain LAI climatology. However, using this method to calculate LAI climatology may be problematic with regard to some vegetation types. In some areas, crops are planted in a fixed order of succession, which would result in an alternating LAI maximum. The field at the Bondville site is a good example of this scenario. The MODIS LAI time series for a pixel from years 2000 to 2004 is shown in Fig. 4a. In 2000, 2002 and 2004, the crop was soybeans with a maximum LAI of 3.5, 3.3, and 3.9, respectively. While in 2001 and 2003, the crop was maize with a maximum LAI of 6.5 and 6.8. In forested areas, LAI values may increase year to year as the forest grows. Fig. 4b shows the monthly LAI calculated from monthly litter fall collections (Martin & Jokela, 2004) at the Donaldson Tract Mid-Rotation site (29.755° N, 82.163° W) from years 1999 to 2003. Following clear-cutting from 1988 to 1989, the site was replanted early in 1990 and the LAI shows seasonal patterns fluctuating with an increasing yearly trend. To solve this problem, we use Seasonal Auto-Regression Integrated Moving Average (SARIMA) models to forecast LAI one year ahead of time, due to the strong seasonality inherent to the vegetation LAI time series.

1.0

1.0

6.0

Red Band NIR Band NDVI Contaminated Obs NDVI Envelope

MODIS LAI upper envelope of LAI 5.0

0.8

Reflectance

4.0

3.0

2.0

1.0

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0.0 2002

2003

2004

Year Fig. 2. MODIS LAI time series and the upper envelope obtained by multi-step Savitzky– Golay filtering.

0.8

NDVI

LAI (m2/m2)

99

0.0 1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2001 Fig. 3. MODIS time series reflectances in the red and NIR band and NDVI temporal profiles with circles to mark the contaminated data at the Bondville site.

100

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

a

According to the Box–Jenkins method (Box & Jenkins, 1976), fitting a SARIMA model to time series data includes three stages: identification, estimation and diagnostic checking. The identification stage serves to determine the estimated orders of the SARIMA model. During the estimation stage, the coefficients of the model are calculated using the method of moments, least square methods, or maximum likelihood methods. Diagnostic checks of the model are performed to reveal possible model inadequacies, and to assist in the selection of the best model. In this study, all analyses are performed using R statistical software, and the selection of SARIMA model was conducted using Akaike's information criterion (AIC) (Akaike, 1974). After the best model is chosen, it is used to forecast future outcomes based on the known data. Since noise signals caused by clouds and other poor atmospheric conditions are negatively biased, an upper LAI envelope smoothing strategy is required. Therefore, the multi-step Savitzky–Golay filter (as described in the previous section) is applied to provide relatively smoothed time series LAI data, using the best models from the smoothed LAI historical data that was utilized to forecast a 1-year lead-time. As shown in Fig. 5a, data from 2000 to 2003 was used to develop the SARIMA model for MODIS LAI, and the dataset in 2004 is used for model validation. Of all the models tested, an ARIMA(2, 0, 0) × (1, 0,0)46 model was found to best fit the data. The equation used is presented in Eq. (9):

8.0 7.0

LAI (m2/m2)

6.0 5.0 4.0 3.0 2.0 1.0 0.0 2000

2001

2002

2003

2004

Year

b

8.0

LAI (m2/m2)

7.0

6.0

5.0

4.0

Zt = 1:358Zt−1 −0:587Zt−2 −0:616Zt−46 + 0:836Zt−47 −0:361Zt−48 + ωt

ð9Þ

3.0

2.0 2001

2002

a

2003

Year Fig. 4. (a) MODIS LAI for the lower right corner pixel of the 7 × 7 km area surrounding the tower at the Bondville site from years 2000 to 2004; and (b) the monthly LAI of evergreen needleleaf forest at the Donaldson Tract Mid-Rotation site from years 1999 to 2003.

Over the last three decades, SARIMA models have been used for economic forecasting and have achieved great success in the commercial and industrial fields (Helfenstein, 1996; Box & Jenkins, 1976). Today, SARIMA models are the most popular linear model for forecasting a seasonal time series. The SARIMA model of Box and Jenkins, generally denoted as ARIMA(p, d, q) × (P, D, Q)s is given by: s
s
D d s
ΦP B ϕp ðBÞ 1−B ð1−BÞ yt = ΘQ B θq ðBÞωt

2

2

p

θq ðBÞ = 1 + θ1 B + θ2 B + ⋯ + θq B

ð5Þ q

ð6Þ

The seasonal auto-regressive and moving average components are represented by polynomials ΦP(Bs) and ΘQ(Bs) of orders P and Q: s
s 2s Ps ΦP B = 1−Φ1 B −Φ2 B −⋯−ΦP B s


ΘQ B

s

2s

= 1 + Θ1 B + Θ2 B

ð7Þ Qs

+ ⋯ + ΘQ B

ð8Þ

SARIMA PRE LAI MODIS LAI AVE SG filtering LAI MODIS LAI

6.0 5.0 4.0 3.0 2.0 1.0 0.0 2002

ð4Þ

where yt denotes the time series data at time t, t = 1, 2,…, k, ωt is the estimated residual at time t, and is the usual Gaussian white noise process, s is the length of the seasonal period, d is the number of regular differences, D is the number of seasonal differences, and B represents the backward shift operators. The ordinary auto-regressive and moving average components are represented by polynomials ϕp (B) and θq(B) of orders p and q, respectively: ϕp ðBÞ = 1−ϕ1 B−ϕ2 B −⋯−ϕp B

8.0 7.0

LAI (m2/m2)

2000

2003

2004

Year

b

8.0 SARIMA PRE LAI LAI AVE Field LAI

7.0

LAI (m2/m2)

1999

6.0

5.0

4.0

3.0

2.0 2001

2002

2003

Year Fig. 5. Comparison of the observed data with forecasts based on the best SARIMA models: (a) for the year 2004 at the Bondville site; and (b) for the year 2003 at the Donaldson Tract Mid-Rotation site.

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

where Zt denotes the LAI value at time t, t = 1, 2,…, 230. Forecasts based on the optimal model for 2004, and the smoothed MODIS LAI data created from the application of the multi-step Savitzky–Golay filter, are shown in Fig. 5a. To better compare and assess the predicted data, the averages of MODIS LAI data from 2000 to 2003 are also shown in Fig. 5a. The predicted LAI profile follows the smoothed LAI data very closely, though the LAI values are slightly larger than the smoothed LAI in the growing season, while the average of MODIS LAI shows abrupt rises or drops in the growing season. It is clear that the accuracy of the predicted LAI has been improved over the average of MODIS LAI, especially during the growing season. Similar results are shown in Fig. 5b for the Donaldson Tract MidRotation site. When compared to the monthly LAI in 2003, forecasts based on the selected best model for the year 2003 are significantly improved over the average of the monthly LAI data from 1999 to 2002. 2.3.2. Dynamic model based on LAI climatology The LAI climatology based on SARIMA models describes the general change tendency of LAI for a 1-year lead-time extracted from historical time series data. In fact, there are always some differences between real LAI and the climatology in a certain year, which could be caused by temperature and precipitation changes, or other disturbances. Therefore, when new observations are available, the biophysical variables must be updated each time. When the Ensemble Kalman Filter (EnKF) is used, the propagation of the state vector requires the use of time-evolving models for the biophysical variables. In this study, the following dynamic model is constructed based on the climatology to evolve LAI in time and used to provide the shortrange forecast of LAI: LAIt = Ft × LAIt−1

ð10Þ

where the term LAIt − 1 represents LAI data at the previous time-step, and the term LAIt represents LAI data at the current time-step; the term Ft is a linear state-transition operator and is formulated as: Ft = 1 +

1 dZt × Zt + ε dt

ð11Þ

where Zt is the temporal profile from LAI climatology, and ε = 10 − 3 prevents a null denominator. The dynamic model in Eqs. (9) and (10) means that LAI is updated by first using a prior estimate of the previous time-step, and then adding to it the climatological derivative of LAI, with time weighted by the ratio of the prior estimate and background estimate. Furthermore, the ratio is used to preserve units and is not necessarily a priori based on physical reasoning. The model incorporates the canopy LAI climatology. Predicated LAI from the dynamic model was input into the radiative transfer model, to simulate reflectance of different bands, which was then compared with MODIS reflectance to correct the LAI predicted using the EnKF techniques. A similar dynamic model was used by Samain et al. (2008) to describe BRDF seasonal evolution for each ECOCLIMAP land cover classes. Using LAI as an example, this approach outlines a real-time evolving model. This method may be applied to other variables that possess inherent temporal change rules as a function of canopy phenology. However, in our study the variables (except LAI) denoted by an asterisk in Table 1 are evolved by a stationary model, which means that the operator Ft in Eq. (10) has been set to the constant value of 1. These parameters and LAI are recursively updated with EnKF techniques as new observations arrive. 2.4. The Ensemble Kalman Filter The EnKF is a sophisticated sequential data assimilation method and was designed to facilitate data assimilation into non-linear

101

models within a Kalman gain scheme. The EnKF was originally proposed by Evensen (1994), and has since been developed further and evaluated in a large number of studies (Burgers et al., 1998; Evensen, 2003). The EnKF applies an ensemble of model states to represent the statistical errors of the model estimate and uses ensemble integrations to predict these errors forward in time. EnKF uses an analysis scheme that operates directly on the ensemble of model states when the observations are assimilated. It has proven capable of efficiently handling strongly non-linear dynamics and large state spaces, and has gained popularity due to its simple conceptual formulation and relative ease of implementation. Let x be an n-dimensional model state vector, and define a matrix A∈ℜn×N containing N ensemble members of the model state vector, A = ðx1 ; x2 ; …; xN Þ. The ensemble mean is stored in each column of a matrix A∈ℜn×N , and the ensemble perturbation matrix A′ ∈ℜn×N can be defined as A′ = A−A. Then, the standard analysis equation is: a

T

A = A + A′ A′ H

T

 −1 T T HA′ A′ H + R ðD−HAÞ

ð12Þ

where H∈ℜm×n (with m being the number of measurements) is the measurement operator relating the model state to the observations, R∈ℜm×m represents the measurement error covariance matrix, and D∈ℜm×N is the perturbed observations ensemble matrix. The superscript ‘a’ denotes the analyzed state vector ensemble and the superscript ‘T’ denotes a matrix transpose. Eq. (12) imposes the condition that the measurement operator H is linear. Therefore, it becomes invalid when the MODIS reflectance data is assimilated into a dynamic model, because the canopy radiative transfer model hð⋅Þ represents the non-linear function of the dynamic model states. It is possible to overcome this limitation by augmenting the dynamic model state vector with a diagnostic variable, which is the model prediction of the measurement (Evensen, 2003). Thus, a new model state vector can be defined for each   ^T = xT ; hT ðxÞ . The analysis equation then can ensemble member as x ˆ n×N ^ as the new ensemble matrix with N be written by defining A∈ℜ ensemble members, where nˆ is the sum of the size of the model state vector n and the number of measurement equivalents added to the ˆ ^′ ∈ℜ n×N original model state vector, and A as the new ensemble perturbation matrix:  −1   a ^′ T H ^′ T H ^ ^′ A ^A ^T H ^A ^T + R A = A + A′ A D−H

ð13Þ

m×nˆ ^ where H∈ℜ is a linear observation operator that transforms between the augmented state vector and the observation. In our study, the model states include the variables denoted by an asterisk in Table 1, so the size of the model state vector n is equal to 4. ^ is formed by augmenting the matrix A with simulated The matrix A reflectances corresponding to the MODIS bands used to inverse LAI. Specifically, the simulated reflectances become part of the state variables of the dynamic model at the observation times. In this paper, the number of MODIS bands used to retrieve LAI is equal to 6, so nˆ is equal to 10. At the same time, we emphasize the importance of the method by which the observation error covarianceR is specified, because it has great influence on the Aa . In our experiments, we use a diagonal error covariance matrix.

3. Results analysis 3.1. Data sets We chose MODIS data over eight BELMANIP sites to test this new algorithm. The eight BELMANIP sites we chose are: Sud-Ouest, Bondville, Mead, Konza, Tonzi, Larose, Harvard Forest and Walker Branch. The attributes of these BELMANIP sites are given in Table 2.

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

The Sud-Ouest site features croplands according to the IGBP classification. The land cover is composed mainly of crops (corn, soybean, and sunflower) and grasslands. Bondville is an agricultural site in the midwestern United States near Champaign, Illinois. The field was continuous no-till with alternating years of soybean and maize crops (Meyers & Hollinger, 2004). The Mead site is located at the University of Nebraska Agricultural Research and Development Center near Mead, Nebraska. Since 2001, the field in this site was under no-till with alternating years of maize and soybean crops, and the crop was soybean in 2004. The Konza site is located in the Flint Hills region of northeastern Kansas; vegetation at this site was tallgrass prairie. In addition to grassland, there are areas of gallery forest and some croplands in the northern part of the study site. The Tonzi Ranch site is located on a privately owned California ranch subject to grazing. The vegetation cover possesses two distinct layers, with blue oak trees in the overstory and annual grasses in the understory. The oak trees are deciduous and dormant in the winter wet season, with an average canopy height of 7.1 m and maximum LAI of 0.6. The grass understory is active only during the winter/spring wet season—from the end of October to the middle of the following May. The maximum LAI of the grass is around 1.0 (Baldocchi et al., 2004). The Larose site is quite flat, and vegetation at the site is forest according to the IGBP classification. The land cover is composed mainly of boreal forest (conifers and deciduous trees) and wetland (grass and shrub). The Harvard Forest site is located in western Massachusetts and is predominantly a closed hardwood forest system with some small patches of conifer and mixed hardwood-conifer stands. Additionally, there are a few marshy lowland areas, pastures and a rural residential development. The Walker Branch site is located on the U. S. Department of Energy's Oak Ridge Reservation near Oak Ridge in Anderson County, Tennessee, and has relatively uniform surface coverage by a deciduous broadleaf forest. The real-time inversion method described above is applied to the MODIS data of these BELMANIP sites. The MODIS collection 5, 8-day, 1-km LAI products (MOD15A2) are used to compute the LAI climatology and construct the dynamic model to evolve LAI in time. The MODIS collection 5, 8-day, 500-m surface reflectance product (MOD09A1) is used to retrieve LAI profiles. One or 2-year data was processed at each site, and the year chosen varies somewhat based on the availability of LAI high resolution reference maps derived from ground measurements. At five of the eight sites, ground LAI data was collected (Table 3) using both direct and indirect methods. At the Bondville site, LAI was measured from destructively sampled corn and soybean plants. At the Konza and Harvard Forest sites, LAI was estimated indirectly using a Li-Cor LAI 2000 canopy analyzer with corrections for clumping and similar factors. At these sites, Landsat Enhanced Thematic Mapper Plus (ETM+) imagery was used to develop high resolution LAI estimates, which are directly linked to field measurements using the methods described by Gower et al. (1999). At the Sud-Ouest and Larose sites, hemispherical images were processed using the method proposed by Lang and Yueqin (1986) to derive true LAI measurements. Additionally, high resolution LAI maps were derived from the determination of the transfer function between the reflectance values Table 2 Selected BELMANIP site information. Site name

Location

Latitude (°)

Longitude (°)

Land cover types

Sud-Ouest Bondville Mead Konza Tonzi Larose Harvard Walker

France USA USA USA USA Canada USA USA

43.506 40.006 41.180 39.080 38.432 45.381 42.532 35.959

−1.237 −88.292 −96.440 −96.620 −120.966 −75.217 −72.188 −84.287

Croplands Croplands Croplands Grassland Grassland Boreal forest Mixed forest Deciduous broadleaf forest

Table 3 LAI field measurement information for the BELMANIP sites. Site name

Year

DOY

True LAI

Effective LAI

Measurement method

Sud-Ouest Bondville

2002 2000 2000 2000 2000 2001 2001 2003 2000

201 186 224 158 238 169 228 231 236

1.9 2.5 3.2 – – – – 5.6 4.2

1.2 – – 2.1 2.1 3.8 2.8 3.5 –

Hemispherical photographs Destructive harvest

Konza

Larose Harvard

LAI 2000

Hemispherical photographs LAI 2000

of the SPOT image acquired during, or around, the ground campaign and LAI measurements (Baret et al., submitted). The General Cartographic Transformation Package (GCTP) is used to transform the high resolution LAI data in UTM WGS84 to the Integerized Sinusoidal Projection (ISIN) used in MODIS data. The high resolution LAI data are up-scaled to a 500 m resolution with a spatial average sampling method. To assess the quality of the retrieved LAI values, we carried out direct comparisons with the mean LAI data. 3.2. Validation results The real-time inversion method is performed to retrieve LAI over the BELMANIP sites with different vegetation types at the pixel level. At the Sud-Ouest site, the LAI climatology for the year 2002 is forecasted by the SARIMA model from the upper envelope of 2000 to 2001 MODIS LAI data, which describes the general change tendency of LAI in 2002. Based on the climatology, a dynamic model is constructed to evolve LAI in time, and then used to provide the short-range forecast of LAI. Forecasts from the dynamic model are used with the EnKF techniques to estimate real-time LAI from the time series MODIS reflectance data. Fig. 6 shows the time series of LAI values with their standard deviation using our new method, MODIS LAI values, LAI climatology, and the BELMANIP mean LAI for July 2, 2002—for a better evaluation of the quality of the retrieved LAI values. The real-time retrieved and MODIS LAI profiles are in very good agreement at the beginning and the end of the growing season, while the climatology is much higher from Julian days 73 to 161 in 2002. During the peak of the growing season, the MODIS LAI profile shows dramatic fluctuations, while the retrieved LAI profile from our new method is relatively smooth and in good agreement with the climatology. Comparatively speaking, the real-time retrieved LAI estimate is closer to BELMANIP mean LAI, than MODIS LAI values. Fig. 7 shows the retrieved LAI profiles at the center pixel of the Bondville site for the years 2000 and 2001. Because the MODIS LAI 3.0 Retrieved LAI MODIS LAI SARIMA Prediction Field LAI

LAI (m2/m2)

102

2.0

1.0

0.0 1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2002 Fig. 6. Retrieved LAI temporal profiles by the real-time inversion method at the center pixel of the Sud-Ouest site.

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

product has only been routinely produced since 2000, there is not enough data available to build a SARIMA model to forecast LAI in 2000 or 2001. At the same time, we noted that at the Bondville site the field was continuous no-till with alternating years of soybean and maize crops. Therefore, the LAI climatology is derived from LAI averages of the 7 × 7 km pixels surrounding the tower in 2000 or 2001, because these pixels belong to the same vegetation type according to the MODIS Land cover product (MOD12Q1). The profiles of LAI averages fluctuate, especially in the growing season. The multi-step Savitzky– Golay filtering is then used to derive the LAI envelopes and a dynamic model of LAI evolution in time is constructed based on these LAI envelopes. Given the initial value of the model state variables, the dynamic model is used to predict the state variable values for the next time. The state variable values are updated when the reflectance data is available. Fig. 7a shows the results retrieved by the real-time inversion method in 2000. The time series of LAI values are displayed with their standard deviation. For the comparison, the MODIS LAI values at the center pixel of a 7 × 7 km area are also shown in Fig. 7a. Both the real-time retrieved and MODIS LAI profiles match reasonably well during the entire growing season. The differences between the profiles are minimal, although the MODIS LAI profile delivers a few fluctuations during the peak of the growing season. In addition, both LAI profiles closely agree with the BELMANIP mean LAI on July 4, 2000. We also noted that the MODIS LAI values on days 217 and 225 are missing, because there is not high quality reflectance data on these days. However, based on the dynamic model, the real-time inversion method can produce reasonable LAI values that effectively agree with the BELMANIP mean LAI on August 11, 2000. The LAI profile retrieved by the real-time inversion method, and the corresponding MODIS LAI

a

5.0 Retrieved LAI MODIS LAI Field LAI

LAI (m2/m2)

4.0

3.0

2.0

1.0

0.0 1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2000

b

values in 2001, are shown in Fig. 7b. Due to cloud contamination and variable atmospheric conditions, fluctuations of the MODIS profile at this pixel are evident—especially during the growing season. Furthermore the MODIS LAI values at days 169 and 177 are missing. Although there is not high quality reflectance data available on these days, the real-time inversion method can remove the abrupt rises or drops caused by noise and provide reasonable LAI values derived from the dynamic model. At the Mead site, multi-step Savitzky–Golay filtering is used to derive the upper envelope of MODIS LAI data. Using the LAI envelope from 2000 to 2003, a SARIMA model is constructed to forecast the LAI climatology for 2004. Based on the climatology, a dynamic model is developed to evolve LAI over time and used to provide a short-range forecast of LAI. Forecasts from the dynamic model are used with the EnKF technique to estimate real-time LAI from the time series MODIS reflectance data. The real-time retrieved LAI values with their standard deviation are shown in Fig. 8. For comparison, the MODIS LAI and LAI climatology are also shown in Fig. 8. The MODIS LAI value on day 185 is missing, but the real-time inversion method can be used to attain reasonable LAI values. At the same time, some abrupt dips within the real-time retrieved LAI profile are found during the crop growing season. These dips result primarily from the inaccuracy of the SARIMA model's forecasted climatology at this pixel. In Fig. 8, the peak value of the climatology occurs around day 193. Before day 193, the dynamic model evolves steadily in a quick upward tendency, and then takes a subsequent rapid and steady downturn. However, in 2004 the real peak value of LAI at this pixel is on day 225. Therefore, in the absence of new observations during the peak of the crop growing season, the LAI values are propagated in a quick downward trend using the dynamic model. Furthermore, as new observations arrive, the updated LAI values are much larger than the corresponding LAI values predicted by the dynamic model, which results in the abrupt dips. At the Konza site, the envelope of the LAI averages for the 7 × 7 km pixels surrounding the tower is used to construct a dynamic model for the same reason as the Bondville site. Fig. 9 depicts the retrieved grassland LAI profiles from our new method, along with MODIS LAI values, at the center pixel of the Konza site for years 2000 and 2001. Overall, the real-time retrieved and MODIS LAI profiles are in very good agreement at the beginning and end of the growing season, while during the peak of the growing season, MODIS LAI values are significantly overestimated in comparison with the retrieved LAI values derived from our new method. Fig. 9a shows the real-time retrieved LAI values with their standard deviation and MODIS LAI values in 2000. It is observed that MODIS LAI values are much higher than the BELMANIP mean LAI for July 4, 2000; except for fluctuations of the MODIS LAI values during 5.0

5.0

Retrieved LAI MODIS LAI SARIMA Prediction

Retrieved LAI MODIS LAI 4.0

LAI (m2/m2)

4.0

LAI (m2/m2)

103

3.0

2.0

3.0

2.0

1.0

1.0

0.0

0.0 1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2001 Fig. 7. Retrieved LAI temporal profile by the real-time inversion method at the center pixel of the Bondville site in 2000 (a) and 2001 (b).

1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2004 Fig. 8. Retrieved LAI temporal profile by the real-time inversion method at the center pixel of the Mead site.

104

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

a

3.0 Retrieved LAI MODIS LAI Field LAI

LAI (m2/m2)

4.0

LAI (m2/m2)

Retrieved LAI MODIS LAI

5.0

3.0

2.0

2.0

1.0

1.0 0.0 0.0

241 265 289 313 337 361 17 41 65 89 113 137 161 185 209 233 1

DOY/2001-2002

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2000

5.0 Retrieved LAI MODIS LAI Field LAI

LAI (m2/m2)

4.0

3.0

2.0

1.0

0.0 1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2001 Fig. 9. Retrieved LAI profile by the real-time inversion method at the center pixel of the Konza site for the years 2000 (a) and 2001 (b).

the peak of the growing season. Additionally, the MODIS LAI values at days 209 and 217 are missing. Compared to the MODIS LAI product, the LAI profile retrieved by the real-time method is much closer to the BELMANIP mean LAI. This is illustrated by the real-time retrieved LAI profile and MODIS LAI values in 2001, which are displayed in Fig. 9b. The MODIS LAI values at days 169 and 177 are missing, but the realtime inversion method can produce reasonable LAI values using the dynamic model. By comparison, the accuracy of the retrieved LAI values from our new method are significantly improved over the MODIS LAI product, when compared to the BELMANIP mean LAI. Fig. 10 shows the retrieved grassland LAI profile at the center pixel of the Tonzi Ranch site from Julian day 241, 2001 to Julian day 233, 2002. Following the same procedure used for the Bondville site, the envelope of the LAI averages for the 7 × 7 km pixels surrounding the tower is used to construct a dynamic model; the real-time retrieved LAI values with their standard deviation and MODIS LAI values are shown in Fig. 10. The MODIS LAI values fluctuated in the grass growing season, with the exception of a missing MODIS LAI value, and abrupt drops on day 361, 2001 and day 1, 2002. At this pixel, the MODIS LAI values are overestimated in comparison with the real-time retrieved LAI values during the peak of the growing season. Fig. 11 shows the real-time retrieved LAI profile at the center pixel of the Larose site in 2003. Because the MODIS LAI values contain dramatic fluctuations, the means of the 7 × 7 km pixels surrounding the tower at each time point from 2000 to 2002 were used to derive the upper envelope of MODIS LAI data. The means were also used to construct a SARIMA model to derive LAI climatology for the year 2003. The real-time retrieved LAI values with their standard deviation, LAI climatology, MODIS LAI values, and the BELMANIP mean LAI for

August 19, 2003 are shown in Fig. 11. The MODIS LAI values contain significant fluctuations throughout the growing season, but the retrieved LAI profile from our new method is relatively smooth. In other words, the accuracy of the real-time retrieved LAI outperforms the MODIS LAI product when compared to the BELMANIP mean LAI data at this site. Fig. 12 depicts the retrieved LAI profile for mixed forest at the center pixel of the Harvard Forest site in 2002 and includes the MODIS LAI values. Because of aerosol contamination and other factors, the MODIS LAI profile depicts dramatic fluctuations during the growing season, while the seasonal trajectory of the real-time retrieved LAI is relatively smooth. For a better comparison, the BELMANIP mean LAI for August 24, 2002, with the standard deviation, were plotted in Fig. 12. At the Harvard Forest site, the field measurements were acquired using the LICOR-2000 instrument, but the BELMANIP LAI data for this site were not corrected for foliage clumping, because the clumping index values were not obtained for the site. Arguably, the BELMANIP LAI data are underestimated in comparison with the true LAI values (Pisek & Chen, 2007). However, without using a clumping factor there is less than a 10 percent error associated with eastern deciduous forests. Further, the BELMANIP LAI data for this site also were not corrected for woody area, so these two effects may cancel partly. Because MODIS LAI values at the Walker Branch site contain dramatic fluctuations, the means of the 7 × 7 km pixels surrounding the tower at each time point from 2000 to 2004 were used to derive the upper envelope of MODIS LAI data and construct a SARIMA model. This was done to determine the LAI climatology for the year 2005. Fig. 13 shows the real-time retrieved LAI values with their standard 8.0 Retrieved LAI MODIS LAI SARIMA Prediction Field LAI

7.0 6.0

LAI (m2/m2)

b

Fig. 10. Retrieved LAI profile by the real-time inversion method at the center pixel of the Tonzi Ranch site.

5.0 4.0 3.0 2.0 1.0 0.0 1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2003 Fig. 11. Real-time retrieved LAI profile at the center pixel of the Larose site in 2003.

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106 8.0 Retrieved LAI MODIS LAI Field LAI

7.0

LAI (m2/m2)

6.0 5.0 4.0 3.0 2.0 1.0 0.0 1

25 49 73 97 121 145 169 193 217 241 265 289 313 337 361

DOY/2002 Fig. 12. Retrieved LAI temporal profiles by the real-time inversion method at the center pixel of the Harvard Forest site in 2002.

deviation at the center pixel of the Walker Branch site in 2005. The LAI climatology and MODIS LAI values are also shown in Fig. 13, which shows that the MODIS LAI values contain a significant amount of fluctuations throughout the growing season. Additionally, some fluctuations within the real-time retrieved LAI profile are also found during the growing season. Comparatively, the retrieved LAI profile from our new method is accurate in the growing season and in good agreement with the LAI climatology extracted from historical time series data.

4. Conclusions and discussions The urgent need for greater hazard/disaster monitoring and warning capabilities demands the real-time, or near-real-time, generation of highly accurate biophysical variables from remotely sensed data. In this study, a new inversion method was developed to estimate real-time LAI from time series MODIS data. The real-time inversion method recursively updates LAI by combining predictions from dynamic models with MODIS reflectance data, and offers the advantage of producing LAI by propagating through a dynamic mode in order to overcome the cloud- and aerosol-related deficiencies of MODIS reflectance data. The results described in this paper have shown that the real-time inversion method is able to produce relatively smooth LAI profiles, and that the accuracy of the real-time retrieved LAI is significantly improved over the MODIS LAI product; as compared to field-measured LAI data. Another advantage of the 8.0

Retrieved LAI MODIS LAI SARIMA Prediction

7.0

LAI (m2/m2)

6.0 5.0 4.0 3.0 2.0

105

method is that it allows for uncertainty of the real-time retrieved LAI values by using the EnKF framework, which may be critical information for LAI product applications. The approach developed in this paper can be generalized to update the biophysical variables of seasonality models for various applications. This includes the recent work by Quaife et al. (2008) and Stockli et al. (2008) who have also demonstrated that the assimilation of satellite data products significantly improves the estimates of model states and parameters. Nevertheless, we have noted that the real-time retrieved LAI profiles are affected by the accuracy of the climatology forecasted by the SARIMA model. In Fig. 8, the peak value of the SARIMA model's forecasted climatology occurs around day 193, while the real peak value of LAI at this pixel is on day 225, which causes abrupt drops or rises within the real-time retrieved LAI profile during the vegetation growing season. In Fig. 11, the LAI trends at the Larose site show unreasonably low values in winter, because the MODIS LAI is inherently biased over snow. However, the real-time method does not account for such large changes in the understory reflectance, even with the variable price coefficient. Fortunately, the accuracy of the climatology forecasted by the SARIMA model can be improved with a longer time series of MODIS LAI data. An alternative approach to solving these problems is to improve the performance of the dynamic model, so that it better characterizes the intrinsic change rules of biophysical variables. The dynamic model used in this study is limited, because it is assumed to be a time-invariant system in which all quantities governing the model's behavior remain constant with time. In fact, there is no guarantee that model behavior does not change over time. In other words, model adjustment over time may be required. Therefore, one potential focal point for future research is to construct more complex models based on the physical and biological process of ecosystems, and develop an integrated and algorithmic framework for dual state-parameter estimation. This framework would use EnKF to estimate both model state variables and parameters with simultaneously given observations. Lastly, based on the existing multi-year MODIS LAI product, a SARIMA is constructed to forecast general change tendency of LAI for a 1-year lead-time. It can be used for long-term monitoring and warning for hazard/disaster prevention. However, there are always some differences between real LAI and the climatology within a given year. Therefore, it is necessary to feedback the updated LAI, before forecasting general change tendency of LAI for the next 1-year leadtime, after satellite data becomes available to update the parameters of SARIMA model. Therefore, further considerations of this feedback loop after their results should also be explored. Acknowledgements This research was supported by the Chinese 973 Program under grant 2007CB714407, NASA under grant NNX09AN36G, the National Natural Science Foundation of China under grants 40871163 and 40701102, the EU Seventh Framework Programme (CEOP-AEGIS), and the National Science and Technology Support Project under grant 2008BAC34B03. We are grateful to Dr. Carol Russell for editing the earlier version of this manuscript. We would also like to thank the CEOS LPV group for providing ground validation data on their Web site (http://lpvs.gsfc.nasa.gov/LPV_LAI10km.html). We are also most grateful to the anonymous reviewers for their helpful comments.

1.0

References 0.0 1

25

49

73

97

121 145 169 193 217 241 265 289 313 337 361

DOY/2005 Fig. 13. Retrieved LAI profiles by the real-time inversion method at the center pixel of the Walker Branch site in 2005.

Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716−723. Alhamad, M. N., Stuth, J., & Vannucci, M. (2007). Biophysical modelling and NDVI time series to project near-term forage supply: Spectral analysis aided by wavelet denoising and ARIMA modeling. International Journal of Remote Sensing, 28(11–12), 2513−2548.

106

Z. Xiao et al. / Remote Sensing of Environment 115 (2011) 97–106

Baldocchi, D. D., Xu, L. K., & Kiang, N. (2004). How plant functional-type, weather, seasonal drought, and soil physical properties alter water and energy fluxes of an oak–grass savanna and an annual grassland. Agricultural and Forest Meteorology, 123(1–2), 13−39. Baret, F., Weiss, M., Allard, D., Garrigues, S., Leroy, M., Jeanjean, H., Fernandes, R., Myneni, R. B., Privette, J.L., Morisette, J., Bohbot, H., Bosseno, R., Dedieu, G., Di Bella, C., Duchemin, B., Espana, M., Gond, V., Gu, X. F., Guyon, D., Lelong, C., Maisongrande, P., Mougin, E., Nilson, T., Veroustraete, F., Vintilla, R. (submitted for publication). VALERI: a network of sites and a methodology for the validation of medium spatial resolution land satellite products, Remote Sensing of Environment. Box, G. E. P., & Jenkins, G. M. (1976). Time series analysis: Forecasting and control. San Francisco: Holden-Day. Burgers, G., van Leeuwen, P. J., & Evensen, G. (1998). Analysis scheme in the ensemble Kalman filter,. Monthly Weather Review, 126(6), 1719−1724. Chen, J. M., & Black, T. A. (1992). Defining leaf area index for non-flat leaves. Plant, Cell & Environment, 15, 421−429. Chen, J. M., Deng, F., & Chen, M. Z. (2006). Locally adjusted cubic-spline capping for reconstructing seasonal trajectories of a satellite-derived surface parameter. IEEE Transactions on Geoscience and Remote Sensing, 44(8), 2230−2238. Chuvieco, E., Peter, E., & Trishchenko, A. P. (2008). Generation of long time series of burn area maps of the boreal forest from NOAA-AVHRR composite data. Remote Sensing of Environment, 112(5), 2381−2396. Combal, B., Baret, F., Weiss, M., Trubuil, A., Mace, D., Pragnere, A., Myneni, R., Knyazikhin, Y., & Wang, L. (2002). Retrieval of canopy biophysical variables from bidirectional reflectance using prior information to solve the ill-posed inverse problem. Remote Sensing of Environment, 84(1), 1−15. Evensen, G. (1994). Sequential data assimilation with a non-linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research, 99, 10,143−10,162. Evensen, G. (2003). The Ensemble Kalman Filter: Theoretical formulation and practical implementation,. Ocean Dynamics, 53(4), 343−367. Fang, H., & Liang, S. (2003). Retrieve LAI from Landsat 7 ETM+ data with a neural network method: Simulation and validation study. IEEE Transactions on Geoscience and Remote Sensing, 41(9), 2052−2062. Fang, H., Liang, S., & Kuusk, A. (2003). Retrieving leaf area index using a genetic algorithm with a canopy radiative transfer model. Remote Sensing of Environment, 85, 257−270. Fernandes, R., Butson, C., Leblanc, S. G., & Latifovic, R. (2003). Landsat-5 TM and Landsat-7 ETM + based accuracy assessment of leaf area index products for Canada derived from SPOT-4 VEGETATION data. Canadian Journal of Remote Sensing, 29(2), 241−258. Garrigues, S., Lacaze, R., Baret, F., Morisette, J., Weiss, M., Nickeson, J., Fernandes, R., Plummer, S., Shabanov, N. V., Myneni, R., et al. (2008). Validation and intercomparison of global Leaf Area Index products derived from remote sensing data. Journal of Geophysical Research, 113(G02028). Gascon, F., Gastellu-Etchegorry, J. P., Lefevre-Fonollosa, M. J., & Dufrene, E. (2004). Retrieval of forest biophysical variables by inverting a 3-D radiative transfer model and using high and very high resolution imagery. International Journal of Remote Sensing, 25(4), 5601−5616. Gower, S., Kucharik, C., & Norman, J. (1999). Direct and indirect estimation of leaf area index, fAPAR, and net primary production of terrestrial ecosystems. Remote Sensing of Environment, 70, 29−51. Hannachi, A. (2007). Pattern hunting in climate: A new method for finding trends in gridded climate data. International Journal of Climatology, 27(1), 1−15. Helfenstein, U. (1996). Box–Jenkins modelling in medical research. Statistical Methods in Medical Research, 5(1), 3−22. Jacquemoud, S., & Baret, F. (1990). PROSPECT: A model of leaf optical properties spectra. Remote Sensing of Environment, 34(2), 75−91. Jacquemoud, S., & Baret, F. (1993). Estimating vegetation biophysical parameters by inversion of a reflectance model on high spectral resolution data.. Crop structure and light microclimate: Characterization and applications. Paris, France: INRA (pp. 339−350).. Jamieson, P. D., Semenov, M. A., Brooking, I. R., & Francis, G. S. (1998). Sirius: a mechanistic model of wheat response to environment variation,. European Journal of Agronomy, 8(3–4), 161−179. Jönsson, P., & Eklundh, L. (2002). Seasonality extraction by function fitting to timeseries of satellite sensor data. IEEE Transactions on Geoscience and Remote Sensing, 40(8), 1824−1832. Kalacska, M., Sanchez-Azofeifa, G. A., Caelli, T., Rivard, B., & Boerlage, B. (2005). Estimating leaf area index from satellite imagery using Bayesian networks. IEEE Transactions on Geoscience and Remote Sensing, 43(8), 1866−1873.

Kimes, D. S., Knyazikhin, Y., Privette, J. L., Abuelgasim, A. A., & Gao, F. (2000). Inversion methods for physically-based models. Remote Sensing Review, 18, 381−440. Koetz, B., Baret, F., Poilve, H., & Hill, J. (2005). Use of coupled canopy structure dynamic and radiative transfer models to estimate biophysical canopy characteristics. Remote Sensing of Environment, 95(1), 115−124. Kuusk, A. (1994). A multispectral canopy reflectance model. Remote Sensing of Environment, 50(2), 75−82. Kuusk, A. (1995). A Markov chain model of canopy reflectance. Agricultural and Forest Meteorology, 76(3), 221−236. Kuusk, A. (2001). A two-layer canopy, reflectance model. Journal of Quantitative Spectroscopy and Radiative Transfer, 71(1), 1−9. Lang, A. R. G., & Yueqin, X. (1986). Estimation of leaf area index from transmission of direct sunlight in discontinuous canopies. Agricultural and Forest Meteorology, 37, 229−243. Lauvernet, C., Baret, F., Haucoet, L., Buis, S., F.X., & Dimet, L. (2008). Multitemporalpatch ensemble inversion of coupled surface-atmosphere radiative transfer models for land surface characterization. Remote Sensing of Environment, 112(3), 851−861. Liang, S. (2004). Quantitative remote sensing of land surfaces. New York: John Wiley and Sons, Inc.. Liang, S., & Strahler, A. H. (1993). An analytic BRDF model of canopy radiative transfer and its inversion. IEEE Transactions on Geoscience and Remote Sensing, 31(5), 1081−1092. Liu, Q., Gu, L., Dickinson, R. E., Tian, Y., Zhou, L., & Post, W. M. (2008). Assimilation of satellite reflectance data into a dynamical leaf model to infer seasonally varying leaf areas for climate and carbon models. Journal of Geophysical Research, 113, D19113. Martin, T. A., & Jokela, E. J. (2004). Developmental patterns and nutrition impact radiation use efficiency components in southern pine stands. Ecological Applications, 14(6), 1839−1854. Meyers, T. P., & Hollinger, S. E. (2004). An assessment of storage terms in the surface energy balance of maize and soybean. Agricultural and Forest Meteorology, 125(1– 2), 105−115. Myneni, R. B., Ramakrishna, R., Nemani, R., & Running, S. W. (1997). Estimation of global leaf area index and absorbed par using radiative transfer models. IEEE Transactions on Geoscience and Remote Sensing, 35(6), 1380−1393. Pisek, J., & Chen, J. M. (2007). Comparison and validation of MODIS and VEGETATION global LAI products over four BigFoot sites in North America. Remote Sensing of Environment, 109, 81−94. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (1992). Numerical recipes in FORTRAN. The art of scientific computing, 2nd Edition. Cambridge: Cambridge University Press. Qin, J., Liang, S., Li, X., & Wang, J. (2008). Development of the adjoint model of a canopy radiative transfer model for sensitivity study and inversion of leaf area index. IEEE Transactions on Geoscience and Remote Sensing, 46(7), 2028−2037. Quaife, T., Lewis, P., Kauwe, M. D., Williams, M., Law, B. E., Disney, M., & Bowyer, P. (2008). Assimilating canopy reflectance data into an ecosystem model with an Ensemble Kalman Filter. Remote Sensing of Environment, 112, 1347−1364. Röder, A., Udelhoven, T., Hill, J., del Barrio, G., & Tsiourlis, G. (2008). Trend analysis of Landsat-TM and -ETM + imagery to monitor grazing impact in a rangeland ecosystem in Northern Greece. Remote Sensing of Environment, 112(6), 2863−2875. Samain, O., Roujeana, J. L., & Geiger, B. (2008). Use of a Kalman filter for the retrieval of surface BRDF coefficients with a time-evolving model based on the ECOCLIMAP land cover classification. Remote Sensing of Environment, 112, 1337−1346. Smith, J. A. (1993). LAI inversion using a back-propagation neural network trained with a multiple scattering model. IEEE Transactions on Geoscience and Remote Sensing, 31 (5), 1102−1106. Stockli, R., Rutishauser, T., Dragoni, D., O'Keefe, J., Thornton, P. E., Jolly, M., Lu, L., & Denning, A. S. (2008). Remote sensing data assimilation for a prognostic phenology model. Journal of Geophysical Research, 113, G04021. Wang, F., Huang, J., Tang, Y., & Wang, X. (2007). New vegetation index and its application in estimating leaf area index of rice. Rice Science, 14(3), 195−203. Xiao, Z., Wang, J., Liang, S., Qu, Y., & Wan, H. (2008). Retrieval of canopy biophysical variables from remote sensing data using contextual information. Journal of Central South University of Technology, 15(6), 877−881. Xiao, Z., Liang, S., Wang, J., Song, J., & Wu, X. (2009). A temporally integrated inversion method for estimating leaf area index from MODIS data. IEEE Transactions on Geoscience and Remote Sensing, 47(8), 2536−2545. Zhang, X., Sun, R., Zhang, B., & Tong, Q. (2008). Land cover classification of the North China Plain using MODIS_EVI time series. International Journal of Photogrammetry and Remote Sensing, 63(4), 476−484.