Sequential Method with Incremental Analysis Update to Retrieve Leaf Area Index from Time Series MODIS Reflectance Data

High-quality leaf area index (LAI) products retrieved from satellite observations are urgently needed for crop growth monitoring and yield estimation, land-surface process simulation and global change studies. In recent years, sequential assimilation methods have been increasingly used to retrieve LAI from time series remote-sensing data. However, the inherent characteristics of these sequential assimilation methods result in temporal discontinuities in the retrieved LAI profiles. In this study, a sequential assimilation method with incremental analysis update (IAU) was developed to jointly update model states and parameters and to retrieve temporally continuous LAI profiles from time series Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance data. Based on the existing multi-year Global Land Surface Satellite (GLASS) LAI product, a dynamic model was constructed to evolve LAI anomalies over time. The sequential assimilation method with an IAU technique takes advantage of the Kalman filter (KF) technique to update model parameters, uses the ensemble Kalman filter (EnKF) technique to update LAI anomalies recursively from time series MODIS reflectance data and then calculates the temporally continuous LAI values by combining the LAI climatology data. The method was tested over eight Committee on Earth Observing Satellites-Benchmark Land Multisite Analysis and Intercomparison of Products (CEOS-BELMANIP) sites with different vegetation types. The results indicate that the sequential method with IAU can precisely reconstruct the seasonal variation patterns of LAI and that the LAI profiles derived from the sequential method with IAU are smooth and continuous. OPEN ACCESS Remote Sens. 2014, 6 9195


Introduction
Leaf area index (LAI), which is defined as one-half of the total green leaf area per unit of horizontal ground surface area [1], is an important biophysical parameter and is related to ecological processes, such as photosynthesis, respiration and transpiration.It is now widely used in climate and ecology models, including ecological system function models, crop growth models and net primary productivity models.Many methods have been developed to retrieve LAI from remote-sensing data.In general, methods to estimate LAI from satellite-sensor data can be classified into empirical methods [2,3] and physical methods [4][5][6].Empirical methods are based mainly on constructing statistical relationships between LAI and a variety of vegetation indices.The main limitation of these empirical approaches is their dependence on measurement conditions, sites and vegetation types.Physical methods are based on the inversion of canopy radiative-transfer models that simulate spectral top-of-canopy bidirectional reflectance by means of physical principles and provide an explicit connection between canopy biophysical variables and the resulting canopy reflectances [7].Because physical methods can be adjusted for a wide range of situations, an increasing number of studies tend to estimate biophysical variables using canopy radiative-transfer models.However, currently available empirical and physical methods generally use only satellite observations acquired at a specific time to retrieve LAI, which may cause temporally discontinuous and inaccurate LAI estimation.The retrieved LAI profiles exhibit many time series fluctuations, and some LAI values are missing, especially during vegetation growing seasons, because of persistent cloud cover.
To improve the quality of LAI retrievals, increasingly more attention is being paid to retrieving LAI using data assimilation techniques that are based on a reasonable consideration of the dynamic change rules of LAI and the time series satellite observations [8,9].Currently, two types of algorithms suitable for approaching the problem of data assimilation have been developed: one is non-sequential assimilation, such as four-dimensional variational (4D-VAR) data assimilation algorithms, and the other is sequential assimilation, such as the Kalman filter (KF) and ensemble Kalman filter (EnKF) [10].Compared with 4D-VAR data assimilation algorithms, sequential assimilation is a type of real-time assimilation system [11] and has the advantage of utilizing time series satellite observations for the real-time/near-real-time retrieval of LAI or other land surface biophysical variables.Many studies have used the sequential assimilation method to estimate real-time biogeophysical variables from satellite observations for rapid land surface change monitoring.Samain et al. [12] developed a time-evolving model based on ECOCLIMAP land cover to describe bidirectional reflectance distribution function (BRDF) seasonal evolution and used a Kalman filter to retrieve surface BRDF coefficients recursively each time that fresh observations became available.This approach ensured continuous production of surface BRDF parameters and improved the time consistency of the albedo product.Xiao et al. [13] developed a sequential assimilation method for the real-time estimation of LAI from MODIS time series reflectance data.The LAI values were updated recursively by combining predictions from dynamic models and MODIS reflectance data using an EnKF technique.The results were significantly improved over MODIS LAI when compared to field-measured LAI data.Nagarajan et al. [14] assimilated L-band microwave observations into a coupled soil-vegetation-atmosphere transfer model using an EnKF technique to improve root-zone soil moisture estimates.The assimilation method offered significant reductions in the root mean square error and the average standard deviation.
Although sequential assimilation methods can assimilate satellite observations into dynamic models sequentially during a forward integration, major drawbacks of these assimilation methods are the discontinuity between the model forecasts and the analysis estimates around the observation times and the introduced shock at the model restart stage [15][16][17].This type of temporal discontinuity in the retrieved LAI profiles is inconsistent with the process of vegetation succession, during which the LAI values vary continuously with time, except in the case of disturbances (such as fires, floods and hurricanes).One method to solve this problem is called the incremental analysis update (IAU), which was first designed for assimilation systems in atmospheric circulation models to reduce initial imbalances in the model forecast [18][19][20].The IAU methods incorporate analysis increments into the model hindcast in the given assimilation interval in a smooth manner [21,22] and have been widely used in atmospheric [23,24] and oceanic general-circulation models [16,20].
The LAI profiles retrieved by the sequential assimilation methods are also affected by the accuracy of the prediction from the dynamic models.In many previous studies, the dynamic models were assumed to be time-variant systems in which all quantities governing the model's behavior remain constant with time, which affects the accuracy of the prediction from the dynamic models [13].In fact, the dynamic models can be subject to many sources of uncertainty and are usually time-variant in the process of evolution [25].To improve the performance of the dynamic model, an alternative approach is to develop dual state-parameter estimation methods in which the model states and parameters are estimated simultaneously from given observations [26].In recent years, studies in many fields, such as hydrology [25,26] and aerography [27,28], have applied these dual state-parameter estimation methods to improve the forecasts of state variables.
Therefore, this paper aims at developing a sequential assimilation method with incremental analysis update to jointly update model states and parameters and retrieve temporally continuous LAI profiles from time series satellite observations.In this method, the dynamic model is used to forecast LAI based on the LAI climatology to overcome the problem of data missing.Meanwhile, the IAU method is chosen to smooth the fluctuation caused by the sequential assimilation.The assimilation method is tested using MODIS data at some flux sites and validated using ground LAI measurements.
The organization of this paper is as follows.Section 2 introduces several details of the implementation of the new assimilation method.We introduce the radiative transfer model, followed by a description about how to construct a dynamic model using an autoregressive (AR) model from Global Land Surface Satellite (GLASS) LAI time series.We then outline the techniques to update model parameters and states from time series MODIS reflectance data.Next, the IAU technique is described to smooth out discontinuities between the forecasts and the analysis estimates around the observation times.The experimental data are also described in this section.Section 3 demonstrates the performance of this new assimilation method over eight sites.Discussions are presented in Section 4, and the final section provides brief conclusions.

Methodology and Data
LAI climatology, defined as the average of multi-year LAI data, describes the general change tendency of LAI extracted from historical time series data.However, due to differences in temperature, precipitation and various other factors, there is always a difference between the real LAI and the LAI climatology in a given year.Let us define this difference as the LAI anomaly.In other words, the LAI anomaly represents the amount by which the real LAI is increased relative to the LAI climatology.Obviously, if the LAI anomaly is accurately known, the LAI value can be calculated on each day of the year by adding together the accurately estimated LAI anomaly and the LAI climatology.In this study, the sequential assimilation method with an IAU technique takes advantage of the KF technique to update model parameters, uses the EnKF technique to update LAI anomalies recursively from time series MODIS reflectance data and then calculates the temporally continuous LAI values by combining the LAI climatology.A flow chart outlining this method is shown in Figure 1.The existing multi-year GLASS LAI product was used to calculate the LAI climatology and LAI anomalies.Then, a dynamic model was constructed using the LAI anomaly time series to evolve the LAI anomalies over time and was used to provide a short-range forecast of LAI anomalies.EnKF was then used to estimate the analysis of LAI anomalies by combining the predictions and the MODIS reflectance data after quality control.After this, an IAU method was developed to smooth out the discontinuities between the forecasts and the analysis estimates in the adjacent observation times.In addition, the newly obtained LAI anomaly was added to the LAI anomaly time series, and the first of the time series was dropped simultaneously.Then, the parameters of the dynamic model were updated with the new LAI anomaly time series using the KF technique before the next data-assimilation cycle to generate accurate forecasts.Detailed descriptions of the new algorithm are given in the following subsections.

Radiative-Transfer Model
The radiative-transfer model is a physical model that links canopy reflectance and canopy properties.In this study, a two-layer Analytical Canopy Reflectance Model (ACRM) developed by Kuusk [29,30] was used for canopy reflectance modelling.Table 1 lists the parameters needed to run the ACRM.Four parameters of the ACRM, which are sensitive to remote sensing data, were chosen as control variables in this study and are marked by asterisks in Table 1.

Dynamic Model
Many dynamic models have been developed to characterize the behavior of biophysical variables over time.Generally, these dynamic models can be classified into two groups: mechanical or semi-mechanical models and statistical models [26].
Mechanical or semi-mechanical models explain the course of growth on the basis of underlying physiological processes in relation to the environment [31].Many mechanical models have been developed for different purposes [32][33][34][35].However, most of these models have many input parameters and are driven by meteorological data.In most cases, it is hard to acquire all of these parameters accurately or the driving data at regional or global scales to run the models, which will have an impact on the accuracy of the model simulation.
Statistical models incorporate inherent change rules to evolve biophysical variables over time through statistical optimization methods.One type of commonly used statistical model, including the AR model, the linear transfer function model and the neural-network updating model, aims to predict state variables through time series analysis techniques [36].Typically, statistical models have the advantage of short-term prediction, because of their ability to represent the relationships between input and output data flexibly.In this study, the AR model was used to describe the evolution of LAI anomalies over time based on LAI climatology.
The standard linear AR model can be described by the following equation [37]: Although the parameters of the AR model can be determined as long as the input LAI anomaly time series is given, there is no guarantee that the model behavior will not change over time.Therefore, an adaptive AR model was used in this study to enable the model parameters to vary as the state estimates were updated.
The adaptive AR model can be described by the following equation: where , i t a are the time-varying parameters of the AR model.At the initial stage, the adaptive AR parameters , i t a were calculated directly from a time series of LAI anomalies.After one data-assimilation cycle, a new LAI anomaly was calculated.The first LAI anomaly was dropped from the time series, and the new LAI anomaly was added to the time series of LAI anomalies to update the parameters of the adaptive AR model using a KF technique described in Section 2.3.Then, the adaptive AR model with the updated parameters was used to predict the LAI anomalies for the next data assimilation cycle.The above procedure was repeated, and the parameters of the adaptive AR model were recursively updated to generate accurate forecasts.

Parameter Estimation of the Dynamic Model Using a Kalman Filter
After each data assimilation cycle, the parameters of the adaptive AR model were estimated.Many algorithms have been developed to estimate the time-varying parameters of adaptive AR models [38,39].In this study, a Kalman filter [40,41] was used to estimate the time-varying parameters recursively each time a new set of LAI anomalies became available.
To use the Kalman filter to calculate the parameters of the adaptive AR model, Equation (4) can be written in the form of Equations ( 3) and (4): where F I is the transition matrix, which is assumed to be time-invariant identity matrix for convenience.Therefore, the AR parameter can be described as: Additionally, the analysis equation of the adaptive AR model parameters can be given by: where the Kalman gain t K and the a-priori state-error correlation matrix 1, are defined as: The superscript " a " denotes the analyzed results, and the superscript "T " denotes a matrix transpose.The one-step prediction error e t is calculated as: where The remaining unknown parameters are z t and t w .Schlögl [37] summarized different methods to utilize parameters of the dynamic model to estimate t z and t w .Here, the estimation methods developed by Jazwinski [42] and Penny and Roberts [38] are used: where and UC is short for the update coefficient, which determines the speed of adaptation.

State Estimation of the Dynamic Model Using an Ensemble Kalman Filter
After parameter estimation of the dynamic model, a short-range forecast of LAI anomalies was provided.When observational data were available, the state estimation was conducted with the predictions from the dynamic model.In the linear case, the Kalman filter was proposed as an optimal recursive data-processing algorithm.For nonlinear problems, the extended Kalman filter (EKF) can be used after linearization of the model using a first-order Taylor series approximation.However, the EKF can lead to unstable results when the nonlinearity in the system is strong [43].To overcome the drawbacks of EKF, Evensen [43] first proposed the EnKF as a sequential data-assimilation method to handle the nonlinearity problem.
The EnKF uses an ensemble of model states to represent the statistical errors in the model estimate and uses ensemble integrations to forecast these errors forward in time.Suppose that ξ is an n-dimension model state vector and that being the new ensemble matrix.n ˆ is the sum of the size of the model state vector n and the number of measurement equivalents, and is the new ensemble perturbation matrix.The new analysis equation can then be written as follows: where is a linear measurement operator relating the augmented state vector to the observations.
In this study, the model states included the variables denoted by asterisks in Table 1.Therefore, the size n of the model state vector was 4. The new ensemble matrix Ŷ combined the four variables of the model state vector and the surface reflectance for six MODIS spectral bands, meaning that n was equal to 10. Measurement error covariance matrix R was defined according to the typical theoretical accuracy of six bands of MODIS reflectance data.

Incremental Analysis Updates
During the state estimation process, the predicted LAI anomalies from the adaptive AR model and the MODIS reflectance data were combined to produce the analysis of LAI anomalies using the EnKF technique.However, the inherent characteristics of sequential assimilation methods result in temporal discontinuities between the forecast and the analysis of LAI anomalies.In this study, an incremental analysis update technique was used to solve this problem.
The IAU technique incorporated analysis increments into model hindcasts and forecasts in a smooth manner [21].Suppose that times 1 − j t , j t and 1 + j t are three adjacent observation times, and time from 1 − j t to 1 + j t is defined as an assimilation interval.The prediction from the dynamic model is calculated from 1 − j t to j t .At observation time j t , the EnKF is used to provide the analytical data in the way of the sequential assimilation method.At this time, the increment between the predicted result and the analysis result is computed.Then, the increment is added proportionally to the time series from 1 The dynamic model is run again from time 1 − j t to 1 + j t , and the IAU process is repeated until the next observation time 1 During IAU integration, the increment of LAI is described as the difference between the analytical result from the EnKF, a t x , and the forecast from the dynamic model, f t x .It is applied as shown in Equation ( 16): where t+1 x is the updated LAI anomaly at time t; t f(X ) is the forecast of the adaptive AR model, where ( , , , ) λ t is a weighting factor for time, which means that the further away the time is from the central time j t , the smaller the added part is.As the time period changes with the time of the observed data, ( ) λ t is calculated according to Equation ( 17): where: The updated LAI anomaly is added with the LAI climatology on that day as the final result of retrieval.In addition, it is also added to the time series of LAI anomalies, and the first of the time series was dropped simultaneously.Then, the parameters of the dynamic model were updated utilizing the updated LAI anomaly time series to forecast the LAI anomaly on the next day.

Data and Preprocessing
This new assimilation method was tested over eight BELMANIP sites [45] with different biome types.The eight sites chosen were Bondville, Demmin, Utiel, Konza, Donga, the Harvard Forest, Larose and the Sonian Forest.Detailed information on these sites is shown in Table 2.The Bondville site is located in the Midwestern United States near Champaign, Illinois.This is an agricultural district with alternating years of soybean and maize crops [46].The Demmin site is located in the Mecklenburgische Seenplatte district, Mecklenburg-Western Pomerania, Germany.The land cover is composed mainly of crops, including barley, wheat, corn, grassland, rapeseed and potatoes [47].The Utiel site is located in the eastern part of Spain.The main land-cover type is vineyards and orchards, followed by natural Mediterranean vegetation and smaller fractions of industrial and urban areas.The topography is generally flat, with slightly hilly regions, mainly covered by natural Mediterranean vegetation [48].The Konza site is located in the Flint Hills region of Kansas, USA.The main vegetation type here is tallgrass prairie, with small areas of gallery forest and cropland in the northern part of the study site.Donga is one of the twelve departments of Benin in West Africa.The Donga site consists of a mosaic of agricultural fields and natural vegetation, including patches of wooded savannah and gallery forests.The dominant vegetation is wooded savannah, including trees, such as Berlinia grandiflora, Daniela olivieri and Isoberlinia.The continuous herbaceous cover is made up of Andropogon, Hyparrhenia and indigo [47].The Harvard Forest site is located in western Massachusetts, USA.It is a closed hardwood forest system with conifer and mixed hardwood-conifer areas.There are also small areas of marshy lowland, pastures and rural residential.The Larose site is located in the southeastern part of Canada.The vegetation type is forest, including boreal forest (conifers and deciduous trees) and wetland (grasses and shrubs) [47].The Sonian Forest lies in the southeastern part of Brussels, Belgium.The land-cover type is mainly forests, such as beeches, oaks and chestnut trees.The undergrowth is sometimes very dense, including ferns and grasses [47].At the Bondville, Konza and Harvard Forest sites, high-resolution LAI reference maps were generated from Landsat Enhanced Thematic Mapper Plus (ETM+) imagery with LAI ground measurements using the methods proposed by Gower et al. [49].At the Demmin, Utiel, Donga, Larose and Sonian sites, LAI ground measurements were derived from hemispherical images as proposed by Lang and Yueqin [50], using SPOT images and the LAI ground measurements to calculate high-resolution LAI reference maps.To validate the accuracy of the retrieved results, the high-resolution LAI reference maps were transferred from the UTM WGS84 projection into the sinusoidal projection used in the MODIS data and were aggregated to 1-km resolution using the spatial-average method.The mean LAI reference map values for all eight sites are also given in Table 2.
The 8-day, 1-km GLASS LAI product derived from the MODIS surface reflectance data was used to calculate the LAI climatology and the LAI anomaly time series.The GLASS LAI product was retrieved using general regression neural networks, which were trained by the fused time series LAI values from MODIS and CYCLOPES LAI products and the reprocessed time series MODIS reflectance values at the BELMANIP sites [51].The GLASS LAI product is spatially complete and temporally continuous, and the temporal smoothness of the GLASS LAI product is the best amongst the existing global LAI products [51].
For each site, the MODIS Collection 5, 8-day, 500-m land-surface reflectance product (MOD09A1) was used as the observed data to retrieve one year's LAI profile using the new method.The specific year chosen for each site was based on available field measurements.The quality of the MODIS surface reflectance data was affected by many factors, including clouds, aerosols, water vapor and ozone.Although many of these effects were removed through atmospheric corrections, the remaining effects could sometimes be very large, requiring further processing.In this study, the contaminated reflectance data were filtered using the NDVI profile calculated from the reflectance values in the red and NIR bands and the upper envelope of the NDVI profile [13].Only high-quality reflectance data were used to estimate LAI in this study.The MODIS reflectance data were aggregated into 1-km resolution to maintain a spatial resolution consistent with the MODIS and CYCLOPES LAI products using an average method.

Analysis of Results
The EnKF with IAU was used to retrieve LAI from time series MODIS reflectance data at the selected sites with different biome types.For comparison, EnKF without the IAU method was also performed over the selected sites.For better comparison and analysis of the retrieved results, both methods used the adaptive AR model with EnKF to retrieve LAI profiles from time series MODIS reflectance data.For clarification, in the following figures, the LAI values derived by the EnKF with IAU are denoted by "LAI_EnKF_IAU", whereas the LAI values derived by the EnKF without IAU are denoted by "LAI_EnKF".For comparison, the main-algorithm retrievals for the MODIS LAI product and the mean LAI reference map values, which are denoted by "Reference LAI", are also shown in the following figures.
At each site, the LAI climatology was calculated from 2000 to 2010 using GLASS LAI data.The LAI anomaly time series were calculated by subtracting the LAI climatology from the GLASS LAI value.Based on the time series of LAI anomaly as the input data, an adaptive AR model was constructed to predict the LAI anomaly for the next time step.In addition, the predictions from the dynamic models and the MODIS reflectance data were combined to update the LAI using the EnKF coupled with the ACRM.Then, the discontinuities between the analysis from the EnKF and the prediction from the dynamic model were removed using the IAU technique.
At the Utiel site (Figure 2a), the LAI anomaly time series was utilized from 2006 to 2007, providing the input data of the adaptive AR model.As seen, the LAI profile retrieved by EnKF with IAU is in agreement with the MODIS LAI profile at the beginning and the end of the crop growing season.From Day 211, the retrieved LAI values were slightly higher than the MODIS LAI values.The retrieved LAI values were closer to the reference LAI compared with the MODIS LAI at this site.By comparison with the LAI profile retrieved by EnKF without IAU, the LAI profile retrieved by EnKF with IAU tackled the discontinuities at observation times.
Figure 2b shows the retrieved results at the Demmin site in 2004.The LAI anomaly used represents 2002 to 2003.In 2004, the MODIS LAI displayed spurious temporal variations, such as sudden drops in LAI and missing data, which were possibly due to residual cloud or atmospheric effects.Good-quality observation data were chosen for the assimilation process.The LAI values retrieved by EnKF with IAU were higher than the MODIS LAI during the growing season.However, the retrieved LAI value was much closer to the reference LAI on Day 164.The difference between the retrieved LAI value and the reference LAI was 1.0, while the difference between the MODIS LAI and the reference LAI was 2.95.In the whole curve of LAI values retrieved by EnKF without IAU, unreasonable discontinuities always exist at the observation time.The IAU technique can effectively handle this type of fluctuation, which is caused by the sequential assimilation method, and make the LAI profile conform to the process of plant growth.At the Konza site (Figure 3b), for the same reason as at the Bondville site, the LAI anomaly was calculated in 2000.At the beginning and end of the growing season, the curves matched well.On Day 169, the IAU results were close to the reference LAI in the absence of the MODIS LAI.From Days 185 to 240, they were slightly higher than the MODIS LAI.On Day 228, the difference between the reference LAI and the retrieved result was 0.6, while the difference between the reference LAI and the MODIS LAI was 1.1.Some stair-step jumps existed in the curve of LAI values retrieved from EnKF without IAU during the whole growing process, but when the IAU method was used, these discontinuities were removed completely.
At the Donga site (Figure 4a), the LAI anomaly was calculated from 2003 to 2004.During the whole growing season, the MODIS LAI was erratic and displayed some low values, which were not realistic for this site.Excluding those observation data that exhibit poor quality, the LAI values retrieved from EnKF without IAU still exhibited large-amplitude discontinuities reproduced from the fluctuating values of the observed data, but exhibited relatively reasonable changing tendencies, especially during the peak growing season compared with the MODIS LAI.However, the IAU method damped the spurious discontinuities, included in the LAI values without IAU and presented a continuous result.The reference LAI was on Day 172.Although the MODIS LAI indicated an abrupt dip in value on that day and the difference with the reference LAI was 1.15, the IAU result still remained close to the reference, with the difference being 0.08. Figure 4b shows the changes in LAI at the Harvard site.Because of atmospheric conditions and other reasons, there were some missing values and dramatic fluctuations in MODIS LAI.Some instability also existed in the LAI values retrieved from EnKF without IAU at the beginning of the growing season, but it was completely damped in the retrieved results from EnKF with IAU.Compared with the reference LAI on Day 236, the retrieved results from EnKF with IAU were higher, and the difference with the reference LAI was 0.82.However, because clumping and undergrowth were not taken into account in the LAI measurements at this site [52], the reference LAI here was actually an effective LAI.This indicates that some underestimation may have occurred in comparison with the true LAI values [53].
For the Larose site, the real-time retrieved LAI profile in 2003 is depicted in Figure 5a.As the input data of the adaptive AR model, the time series of LAI anomaly ran from 2001 to 2002.Only a few numbers from the MODIS LAI retrieved from the main algorithm are shown in the profile.Utilizing the assimilation with reflectance data, the retrieved results from EnKF without IAU were relatively smooth, except for a few jumps in the profile.When added to the IAU process, this method not only simulated the process of plant growth according to the dynamic model, but also updated the curve realistically according to the growth rhythm of vegetation.The reference LAI was obtained on Day 219.Although the quality of reflectance data was poor at this time, the retrieved results from EnKF with IAU still presented continuous results.
At the Sonian site (Figure 5b), the time series of the LAI anomaly was from 2002 to 2003.Dramatic fluctuations and data missing existed in MODIS LAI during the whole growing season.The retrieved LAI without IAU removed the abrupt rises and drops and provided reasonable LAI values, despite some stair-step shapes at observation times.However, the IAU method successfully handled this problem.Even though the MODIS LAI displayed unstable changes during the whole growing season, the retrieved result from EnKF with IAU was closer to the reference LAI.On Day 174, the difference between the MODIS LAI and the reference LAI was about 3.8, while the difference of the retrieved result was 0.83.

Discussions
A sequential assimilation method with incremental analysis updates was developed to retrieve temporally continuous LAI profiles from time series satellite observations.The method estimated states and parameters of dynamic models simultaneously from the given observations.The model behavior is adjusted after each data assimilation cycle, which leads to better characterization of the intrinsic change rules of biophysical variables.The results described in this study have shown that the adjusted dynamic model can be corrected according to the new added observation data after each data assimilation cycle.
To overcome the discontinuities between the updated LAI from the EnKF and the predicted LAI from the dynamic model, which result from the inherent characteristics of the sequential assimilation methods [13], an IAU technique was introduced to produce continuous and smooth LAI profiles that are consistent with the process of vegetation succession in this study.In fact, the ensemble Kalman smoother (EnKS) is also an effective method to smooth profiles from time series observations.It applies future observations backwards in time using the ensemble covariances [9].However, every time a new dataset is available during the forward integration, an analysis is computed by EnKS for all previous times up to this time.Thus, the EnKS is much more CPU-intensive and also requires greater storage proportional to the total number of analysis times over which the statistics are to be applied [16].In contrast, IAU is computationally more efficient.On the downside, the inversion method can remove abrupt spikes and dips, which may preclude a description of more local (in time) impacts of disturbances, such as fire, disease and insect damage.
In addition, the accuracy of the retrieved LAI profiles depends on the observation error and the forecast error from the dynamic model.The observation error is related to instabilities in observation, such as atmospheric, cloud contamination and calibration errors.Though quality control is used, there may still exist some errors, because the quality flags cannot ensure being absolutely right.The forecast error is mainly from the prediction process of the dynamic model.It had been shown on each day using the straight grey line.After each update, the forecast error was small, and it would gradually increase over time.

Conclusions
A new sequential assimilation method with IAU was developed in this study to retrieve LAI from time series MODIS reflectance data.The LAI climatology was calculated by averaging multi-year GLASS LAI data, and the LAI anomaly was obtained by subtracting the LAI climatology from the GLASS LAI value for each point in time.With the LAI anomalies as the input data, an adaptive AR model was constructed to evolve the LAI anomalies over time, and the parameters of the AR model were updated after each data assimilation cycle.The predictions from the adaptive AR model and the MODIS reflectance data were combined to update the LAI anomalies using the EnKF technique.The IAU technique was used to remove the discontinuities between the predicted LAI anomalies from the adaptive AR model and the updated LAI anomalies.The retrieved result was finally obtained from the updated LAI anomalies using the IAU technique with the added LAI climatology.The method jointly updated model states and parameters from given observations, and an IAU technique was introduced to eliminate the discontinuities between the updated LAI from the EnKF and the predicted LAI from the dynamic model.The results have demonstrated that the new method can produce temporally continuous and smooth LAI profiles from time series satellite observations.Compared with the MODIS LAI product, the retrieved LAI values are closer to the ground-based reference LAI values.In a forthcoming study, more extensive validation for all cover types will be carried out over those sites with multiple field measurements in a single year.

Figure
Figure 1.Flow chart of the LAI inversion method with incremental analysis updates.AR, autoregressive; ensemble Kalman filter (EnKF), ensemble Kalman filter.

Figure 2 .
Figure 2. Retrieved LAI profiles by the inversion method with incremental analysis update (IAU) at the center pixel of (a) the Utiel site in 2008 and (b) the Demmin site in 2004.

Figure 3 .
Figure 3. Retrieved LAI profiles by the inversion method with IAU at the center pixel of (a) the Bondville site in 2000 and (b) the Konza site in 2001.

Figure 4 .
Figure 4. Retrieved LAI profiles by the inversion method with IAU at the center pixel of (a) the Donga site in 2005 and (b) the Harvard site in 2002.

Figure 5 .
Figure 5. Retrieved LAI profile by the inversion method with IAU at the center pixel of (a) the Larose site in 2003 and (b) the Sonian site in 2004.

1 .
Flow chart of the LAI inversion method with incremental analysis updates.AR, autoregressive; ensemble Kalman filter (EnKF), ensemble Kalman filter.

Table 1 .
Parameters needed to run the two-layer Analytical Canopy Reflectance Model (ACRM).
previous LAI anomalies, and p is the order of the AR model.The input LAI anomalies from time 1 − t to time p t − were used to calculate the parameters of the AR model, and the predicted LAI anomalies t x can be obtained from Equation (1).
observation ensemble matrix.The superscript " a " denotes the analyzed state vector ensemble.Equation (14) imposes the condition that the measurement operator H is linear.It is invalid when the MODIS reflectance data are assimilated into a dynamic model, because the canopy reflectance [44]epresents the non-linear function of the dynamic model states.This limitation can be overcome by augmenting the dynamic model state vector with a diagnostic variable, which is the model prediction of the measurement[44].Therefore, a new model state vector can be defined as ˆ[ , ( )]

Table 2 .
Characteristics of the selected sites.