Extended Data-Based Mechanistic Method for Improving Leaf Area Index Time Series Estimation with Satellite Data

Leaf area index (LAI) is one of the key parameters in crop growth monitoring and global change studies. Multiple LAI products have been generated from satellite observations, many of which suffer from data discontinuities due to persistent cloud contamination and retrieval algorithm inaccuracies. This study proposes an extended data-based mechanistic method (EDBM) for estimating LAI time series from Moderate Resolution Imaging Spectroradiometer (MODIS) data. The data-based mechanistic model is universalized to supply the LAI background information, and then the vegetation canopy radiative-transfer model (PROSAIL) is coupled to calculate reflectances with the same observation geometry as MODIS reflectance data. The ensemble Kalman filter (ENKF) is introduced to improve LAI estimation based on the difference between simulated and observed reflectances. Field measurements from seven Benchmark Land Multisite Analysis and Intercomparison of Products (BELMANIP) sites and reference maps from the Imagine-S project La Albufera, Spain site were used to validate the model. The results demonstrate that when compared with field measurements, the LAI time-series estimates obtained using this approach were superior to those obtained with the MODIS 500 m resolution LAI product. The root mean square errors (RMSE) of the MODIS LAI product and of the LAI estimated with the proposed method were 1.26 and 0.5, respectively. When compared with reference LAI maps, the results indicate that the estimated LAI is spatially and temporally consistent with LAI reference maps. The average differences between EDBM and the LAI reference map on the selected four days was 0.32.


Introduction
Leaf area index (LAI) is a key parameter in crop growth monitoring, crop yield estimation, land surface process simulation, and global change studies.LAI time series are critical for global land and climate change research.Currently, global LAI products are routinely produced from data acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) [1], the Multiangle Imaging SpectroRadiometer (MISR) [2], the Advanced Very High Resolution Radiometer (AVHRR) [3], and a low-resolution instrument onboard the Système Probatoire d'Observation de la Terre (SPOT/VEGETATION) [4].However, due to persistent cloud contamination, instrument problems, and retrieval algorithm inaccuracies, the continuity of the LAI product still needs to be improved.

Methodology
The EDBM method consists of three components used to estimate the LAI time series.First, the UDBM model is constructed using historical MODIS 500 m LAI and reflectance data.The reflectance data for the date to be estimated are then input to the UDBM model to obtain primary estimates of LAI as an estimation background (LAI UDBM ).Second, the radiative-transfer model (RROSAIL) is used to link LAI with remotely sensed bidirectional reflectance.The LAI UDBM is input together with other state variables into PROSAIL to calculate the simulated bidirectional reflectance, which contains the same angular information as the current observations.The state variables are determined according to a global sensitivity analysis method.Third, the ENKF is implemented to update the LAI UDBM according to the difference between the observed reflectance and the PROSAIL simulated reflectance to obtain the final estimate (LAI EDBM ), and the state variables are also updated for further use in the next LAI UDBM calculation.A flowchart of this process is shown in Figure 1.Flowchart of the proposed estimation approach.The UDBM model is constructed using time-series reflectance and LAI data, the PROSAIL model is executed to simulate reflectance, and an ensemble Kalman filter is used to update LAIUDBM according to the difference between observed and simulated reflectances to obtain the final estimation.

Theory of the Data-Based Mechanistic Modeling
The term "data-based mechanistic (DBM) modeling" means that the model is inferred statistically from time-series data in an inductive manner.It involves a scientific understanding of complex parameterized models [21].The DBM approach is a modeling method based on time-series data, and objective inference from data is the primary driving force.Only the basic form of the model structure is known before the DBM model is constructed; the model order and parameter values must be identified from training data.This method has been used to model water quality in rivers [26] and many different systems in diverse areas of application [20,27].The basic form of the model is presented in Equation (1): where is the model output, , is the model input data, and and are time-series sequencerelated polynomials that are related to the model input.is the backward shift operator, which has the form = , and is the model estimation error.The forms of and are shown in Equation ( 2): where and are model polynomial operator parameters which are determined during model construction and , , , are model order numbers that are related to the polynomial operator and Flowchart of the proposed estimation approach.The UDBM model is constructed using time-series reflectance and LAI data, the PROSAIL model is executed to simulate reflectance, and an ensemble Kalman filter is used to update LAI UDBM according to the difference between observed and simulated reflectances to obtain the final estimation.

Theory of the Data-Based Mechanistic Modeling
The term "data-based mechanistic (DBM) modeling" means that the model is inferred statistically from time-series data in an inductive manner.It involves a scientific understanding of complex parameterized models [21].The DBM approach is a modeling method based on time-series data, and objective inference from data is the primary driving force.Only the basic form of the model structure is known before the DBM model is constructed; the model order and parameter values must be identified from training data.This method has been used to model water quality in rivers [26] and many different systems in diverse areas of application [20,27].The basic form of the model is presented in Equation (1): where y t is the model output, u i,t is the model input data, and A and B are time-series sequence-related polynomials that are related to the model input.L is the backward shift operator, which has the form L n y t = y t−n , and ε t is the model estimation error.The forms of A and B are shown in Equation (2): where a and b are model polynomial operator parameters which are determined during model construction and n, i, j, k are model order numbers that are related to the polynomial operator and must be identified during the modeling procedure.To construct a DBM model for LAI estimation, historical long-term observations of surface reflectance and the LAI product are used for training.The reflectances of the red (first), near infrared (second), and shortwave infrared (seventh) MODIS bands are often used to estimate LAI [28,29].Hence, when using the DBM modeling method to estimate LAI, the same reflectance data bands were selected.The DBM model used for LAI estimation can be described as in Equation (3): where the LAI t is the estimated LAI at time t and R 1 t , R 2 t , R 7 t are the MODIS observed reflectances of bands 1, 2, and 7 respectively.
Assuming that L n LAI t = LAI t−n , Equation (3) can be rewritten as Equation ( 4): which means that a current LAI value can be deduced from its current and previous time-node reflectance observations and from former time-node LAI values.The simplified refined instrumental variable (SRIV) method [30] is used to calculate the model polynomial operator parameter values (a and b).The modeling accuracy is evaluated by the corresponding coefficients of determination (R 2 ) and the Young information criterion (YIC) [27].R 2 indicates how well the model explains the data, and the YIC represents the complexity of the model structure.The lower the YIC, the simpler the model structure is.The standard for selecting the model structure is to choose the structure with the highest R 2 and the lowest YIC.Once the DBM model has been constructed, it can predict future LAI values based on known historical time-series data.

Universal Data-Based Mechanistic Model Construction
To avoid pixel-by-pixel construction of the DBM model, which requires a large amount of historical data and calculations and may fail when the LAI variation pattern changes due to climate change and human activities, the DBM is first revised for regional use to obtain regional LAI climatology information.For this purpose, the UDBM model is introduced.The primary goal of the UDBM model is to provide a common model structure for pixels with similar land cover types or LAI dynamics in the same region.The basic requirement of the model is that it should represent the integrated LAI variation for the same type of vegetation, or for different types of vegetation but similar dynamic LAI levels.The results obtained using the UDBM model are used as a priori background information to obtain the final estimates.
When constructing the UDBM dynamic model, pixels covered with vegetation are divided into two types: those with a relatively high LAI level (here forest pixels are defined as having a high LAI level), and those with a relatively low LAI level, such as brush, crops, and grass.Data from 20 forest pixels at BELMANIP sites were collected to obtain the first series of the reflectance-LAI dataset.Another 20 pixels of data from non-forest vegetation types were collected to generate the second series of the reflectance-LAI dataset.Considering that the MODIS LAI product may be affected by cloud contamination during the growing season [9], Savitzky-Golay (SG) filtering [31,32] was used subsequently to smooth the LAI time series.The filter window size was set to 4, and LAI quality control (QC) data were used as the weighting file.For high-quality data (QC < 31), the weight was set to 1; for acceptable LAI data (32 < QC < 63), the weight was set to 0.5; and for poor-quality data (QC > 64), the weight was set to 0. Filtered datasets were then used to construct the two UDBM models, namely UDBM 1 for forest and UDBM 2 for other land cover types.The structures of the two UDBM models are given in Equations ( 5) and (6).

PROSAIL Model
The vegetation canopy radiative-transfer model establishes the relationship among surface reflectance, canopy structure parameters, and canopy spectral parameters by forward simulation of the radiation transfer in and out of the canopy [33][34][35].In the parameter retrieval strategy used here, the PROSAIL model, described in detail by Jacquemoud [33] and Kuusk [36], was used as the forward model to simulate canopy directional reflectance.PROSAIL is a combination of a leaf optical properties model (PROSPECT) and a canopy bidirectional reflectance model (SAIL).The outputs of the PROSPECT model (leaf reflectance and transmittance) are input directly into the SAIL model.The input parameters for PROSAIL are listed in Table 1.Six parameters that are globally sensitive to the output [37,38] were defined as the state variables in the ENKF scheme described in Section 2.3.They are marked with an asterisk in Table 1.The initial LAI values were determined by the LAI background from the UDBM model.The spectral resolution of the PROSAIL output reflectance is 1 nm, which is not consistent with the MODIS reflectance bandwidth.Therefore, PROSIAL model simulate reflectance was first converted to match the MODIS bands using the MODIS band 1, 2, and 7 spectral response functions (the band ranges and spectral response functions are shown in Table 2 and Figure 2).The conversion formula is given in Equation (7): where R λ is the reflectance calculated using PROSAIL with wavelength λ and f (λ)

Ensemble Kalman Filter
The ENKF [39,40] is an extension of the Kalman filter (KF) [41,42].In general, the parameters of the processing model are stable, which might introduce increased errors over time.The ENKF was designed to facilitate data assimilation into nonlinear processing models within a Kalman gain scheme.The formulation is mainly according to [39,40], which has the same basic form as the KF.However, instead of propagating a model error covariance matrix, the ENKF uses an ensemble of model states to represent the error statistics [43].
Let x be an n-dimensional model state variable, and define a matrix A ∈ × that contains N ensemble members of the model state variable set; ∈ × is the matrix of ensemble perturbations (the ensemble minus its mean); the ensemble covariance matrix ∈ × is defined as Equation ( 8): For each measurement, we define the perturbed observation as Equation ( 9): Then the observation ensemble is given as Equation ( 10): The observation error covariance matrix derived from the observation perturbation can be expressed as Equation ( 11): where = ( , , … ) ∈ × and ∈ × is the observation operator.The standard analytical equation for ENKF is as Equation ( 12): where the superscript a denotes the analyzed set of state variables and the superscript T denotes the matrix transpose.Equation ( 12) imposes the condition that the observation operator H is linear, and therefore it becomes invalid when the PROSAIL canopy radiative transfer model represents a nonlinear function of the model states.This situation can be overcome by forming augmented state variables [40].Equation ( 12) subsequently becomes Equation ( 13):

Ensemble Kalman Filter
The ENKF [39,40] is an extension of the Kalman filter (KF) [41,42].In general, the parameters of the processing model are stable, which might introduce increased errors over time.The ENKF was designed to facilitate data assimilation into nonlinear processing models within a Kalman gain scheme.The formulation is mainly according to [39,40], which has the same basic form as the KF.However, instead of propagating a model error covariance matrix, the ENKF uses an ensemble of model states to represent the error statistics [43].
Let x be an n-dimensional model state variable, and define a matrix A ∈ R n×N that contains N ensemble members of the model state variable set; A ∈ R n×N is the matrix of ensemble perturbations (the ensemble minus its mean); the ensemble covariance matrix P e ∈ R n×n is defined as Equation ( 8): For each measurement, we define the perturbed observation as Equation ( 9): Then the observation ensemble is given as Equation ( 10): The observation error covariance matrix derived from the observation perturbation can be expressed as Equation (11): where γ = ( 1, 2, . . .N ) ∈ R m×N and H ∈ R m×n is the observation operator.The standard analytical equation for ENKF is as Equation ( 12): where the superscript a denotes the analyzed set of state variables and the superscript T denotes the matrix transpose.Equation ( 12) imposes the condition that the observation operator H is linear, and therefore it becomes invalid when the PROSAIL canopy radiative transfer model represents a nonlinear function of the model states.This situation can be overcome by forming augmented state variables [40].Equation ( 12) subsequently becomes Equation ( 13): where Â ∈ R n×N and Â ∈ R n×N are formed by augmenting the state variable ensemble by the predicted reflectances and Ĥ ∈ R m× n is the observation operator that transforms between the augmented state variables and the observations.The main ENKF procedure includes four steps: first, for each measurement, virtual observations are generated by adding a normal disturbance.Here 100 virtual observations are generated according to the stability of the estimated LAI.Then, forecasting is performed with the dynamic model; (in this case the UDBM model).Third, estimation is updated with remote-sensing observations according to the estimation error and covariance, and finally the updated state variables and error variance are entered into the next assimilation cycle.In the ENKF assimilation scheme, the bidirectional reflectance with the current solar and observation angular information is simulated using the PROSAIL model.The PROSAIL model is the observation operator in Equation ( 13).The LAI calculated from UDBM is the initial input.The state variables are the six PROSAIL sensitive parameters listed in Table 1.The reflectances corresponding to MODIS bands 1, 2, and 7 are used in this procedure.Hence, the number of augmented state variables ( n) is 9.The observation error covariance R e is specified using a diagonal error covariance matrix.The initial values and variances of the state variables are listed in Table 3.

Satellite Datasets
The MODIS sensor contains seven spectral bands in the shortwave range for land use, and bands 1, 2, and 7 are always used to estimate LAI [28,29].The reflectances of these three bands were also used to estimate LAI in the present study.MOD09A1 provides the best possible set of observations during an 8-day period with 500 m resolution [44], including directional reflectance values, quality assessment, and the day of the year for each pixel, as well as solar, view zenith, and relative azimuth angles.The MODIS global LAI product (MOD15A2H) is produced every 8 days at 500-m resolution, including the LAI value, quality rating, and the standard deviation [45].In the present study, 20 BELMANIP forest (Table 4) and non-forest sites (Table 5) were selected separately, and MOD09A1 reflectance and MOD15A2H LAI data for 2008 at these sites were downloaded from the Land Processes Distributed Active Archive Center.An SG filter was first applied to the LAI data, and then the paired reflectance and LAI dataset was used for UDBM model construction with the method introduced in Section 2.1.
When estimating LAI with the EDBM model, only reflectance (MOD09A1) and land cover type data (MCD12Q1) are needed; these were also downloaded from the Land Processes Distributed Active Archive Center.Garrigues, Lacaze [24] listed 41 BELMANIP sites with 81 LAI reference values obtained through a coordinated international effort towards direct validation conducted by the Committee on Earth Observation Satellites Working Group on Calibration & Validation Land Product Validation Sub-Group (CEOS-WGCV LPV).Generally, field optical measurement methods obtain effective estimates of LAI [46], and satellite products give an approximation of the true LAI [47].So effective LAI must be converted to true LAI when validating remote-sensing products using field measurements.Therefore, in this study, sites for EDBM evaluation were selected where the field measurements took into account the clumping effect [48] and where the land cover types included broadleaf forest, mixed forest, savanna, and crops.With the clumping index (C index ) available, we can easily transfer field observed effective LAI (LAI e f f ective ) to true LAI (LAI true ) with Equation (14).
Sites located in the Southern Hemisphere were excluded because the UDBM models are constructed for Northern Hemisphere use.Bigfoot sites with typical seasonal features were also included.Site information is listed in Table 6.

LAI Reference Maps
Since 2013, the Image-S project has been conducting field campaigns at various sites in different countries to collect ground data to validate satellite-derived biophysical products, including LAI, FAPAR, FCover, and others.Field measurements over the main vegetation types have also been upscaled using high-resolution satellite imagery to generate reference maps at 5 km × 5 km and 20 km × 20 km.In this study, the 20 km × 20 km reference map at La Albufera, Spain, was chosen to validate the proposed method.The site is located in the La Albufera Natrral Park, in the east of coast of Spain.The central coordinate is 39.2743 • N−3164 • E. The study area has a subtropical Mediterranean climate and a rich variety of flora.The location land cover types of the research area is shown Figure 3, the classification is based on the GlobeLand30 dataset [49], most of the land area is cultivated land in this area.Field campaigns were carried out on 17 June, 15 July, 7 August, and 22 August.Field observations were first used to generate the LAI regression function, after which 30 m LAI reference maps covering the four days were produced.This study used the ground-derived true LAI map as the reference map; the 30 m reference map was first re-projected to a sinusoidal projection and then resampled to 500 m using a cubic resampling strategy.The MODIS product that was closest in time to the reference map was then selected.The data used in this study are listed in Table 7. was closest in time to the reference map was then selected.The data used in this study are listed in Table 7.

Estimated LAI Temporal Profile
Figure 4 shows the estimated LAI (LAIEDBM) time series compared with the MODIS LAI product (LAIMODIS), the LAI background (LAIUDBM), and the LAI reference (LAIREF).For each site, the forest or non-forest UDBM model was selected corresponding to the land cover type.UDBM models were used as the dynamic equation and provided the LAI estimation background.Reflectance and its observation geometry were the inputs to the LAI estimation.
At the HARV mixed forest site, LAIMODIS fluctuated severely, especially during the growing season, and data were missing for days 97 and 217 in 2000, 169 and 177 in 2001, and 81 in 2002.LAIUDBM depicted the LAI variation pattern, but showed some discrepancies with the LAIMODIS and LAIREF.Compared with LAIREF, LAIEDBM had the highest accuracy, with the mean difference between LAIEDBM and LAIREF being 0.02, whereas the difference between LAIMODIS and LAIREF was −0.3.
At the CHEQ mixed forest site, there was a large difference between the LAIMODIS (6.2) and LAIREF (3.05).LAIEDBM updated estimates based on the LAI background supplied by the UDBM model were more consistent with LAIREF.
At the KNOZ and Sud-Ouest sites, the UDBM non-forest model was used to obtain the estimation background.For both sites, the background peak value was between 1 and 2. At the KONZ site, LAIUDBM was underestimated compared with LAIREF and LAIMODIS.However, the proposed method improved estimation accuracy.The final estimation fitted LAIREF very well, with a mean difference between LAIEDBM and LAIREF of 0.11.At the Sud-Ouest site, LAIMODIS was 1.1, which is very small.LAIUDBM is also very small and consistent with the MODIS data, but much lower than LAIREF.The LAIEDBM upgraded the estimates, and fitted LAIREF well, but the LAI profile is very different with LAIMODIS.

Estimated LAI Temporal Profile
Figure 4 shows the estimated LAI (LAI EDBM ) time series compared with the MODIS LAI product (LAI MODIS ), the LAI background (LAI UDBM ), and the LAI reference (LAI REF ).For each site, the forest or non-forest UDBM model was selected corresponding to the land cover type.UDBM models were used as the dynamic equation and provided the LAI estimation background.Reflectance and its observation geometry were the inputs to the LAI estimation.
At the HARV mixed forest site, LAI MODIS fluctuated severely, especially during the growing season, and data were missing for days 97 and 217 in 2000, 169 and 177 in 2001, and 81 in 2002.LAI UDBM depicted the LAI variation pattern, but showed some discrepancies with the LAI MODIS and LAI REF .Compared with LAI REF , LAI EDBM had the highest accuracy, with the mean difference between LAI EDBM and LAI REF being 0.02, whereas the difference between LAI MODIS and LAI REF was −0.3.
At the CHEQ mixed forest site, there was a large difference between the LAI MODIS (6.2) and LAI REF (3.05).LAI EDBM updated estimates based on the LAI background supplied by the UDBM model were more consistent with LAI REF .
At the KNOZ and Sud-Ouest sites, the UDBM non-forest model was used to obtain the estimation background.For both sites, the background peak value was between 1 and 2. At the KONZ site, LAI UDBM was underestimated compared with LAI REF and LAI MODIS .However, the proposed method improved estimation accuracy.At the Alpilles broadleaf crop site, LAI UDBM , which presented the common state of low-LAI plants, had the same variation pattern as LAI MODIS , but at a slightly higher level.The EDBM model results were the best compared with LAI REF , with a difference between LAI EDBM and LAI REF of only −0.07.
At the Larzac savanna site, LAI MODIS was much higher than LAI REF .The UDBM model estimates were lower than LAI MODIS , and LAI EDBM had the highest accuracy compared with LAI REF .8. The root mean square error (RMSE), bias, and mean absolute error (MAE) of LAI EDBM were 0.5, 0.12, and 0.30, respectively, which indicates that high estimation accuracy can be achieved using the proposed EDBM method.

Comparison with LAI Reference
At the ARGO broadleaf crop site, MODIS LAI values were missing for days 217 and 225 in 2000.The EDBM model upgraded the estimates to fit LAIREF better, but gave a higher peak value during the growing season than the LAIMODIS.
At the Alpilles broadleaf crop site, LAIUDBM, which presented the common state of low-LAI plants, had the same variation pattern as LAIMODIS, but at a slightly higher level.The EDBM model results were the best compared with LAIREF, with a difference between LAIEDBM and LAIREF of only −0.07.
At the Larzac savanna site, LAIMODIS was much higher than LAIREF.The UDBM model estimates were lower than LAIMODIS, and LAIEDBM had the highest accuracy compared with LAIREF.

Comparison with LAI Reference
Figure 5 shows a scatter diagram of the LAIREF, LAIMODIS, and LAIEDBM.The results show that LAIMODIS was the most discrete.LAIEDBM fitted LAIREF very well because most of the points are very near to the y = x line.Statistical results for the LAI estimation errors are listed in Table 8.The root mean square error (RMSE), bias, and mean absolute error (MAE) of LAIEDBM were 0.5, 0.12, and 0.30, respectively, which indicates that high estimation accuracy can be achieved using the proposed EDBM method.

Comparison with the Reference Maps
We then used the proposed EDBM model to generate the LAI map at La Albufera site.UDBM models constructed in Section 2.1 were used to generate the estimation background.Land cover type data were used to determine which UDBM model should be chosen for each pixel.Four days' reference LAI maps during the growing season in 2014 were generated based on field observations.LAIEDBM was directly compared with the reference LAI map. Figure 6 shows a comparison of LAIMODIS, LAIEDBM, and reference LAI maps. Figure 7 shows the difference between estimated LAIs and reference LAI maps.Table 9 shows the mean value of LAIMODIS, LAIEDBM, and the reference LAI maps for the four days.Figures 6 and 7 and Table 9 show that, LAIEDBM properly depicted the spatial pattern of surface LAI and captured the temporal variation.On 17 June, the reference LAI was very low, but LAIMODIS and LAIEDBM were higher than the reference LAI.

Comparison with the Reference Maps
We then used the proposed EDBM model to generate the LAI map at La Albufera site.UDBM models constructed in Section 2.1 were used to generate the estimation background.Land cover type data were used to determine which UDBM model should be chosen for each pixel.Four days' reference LAI maps during the growing season in 2014 were generated based on field observations.LAI EDBM was directly compared with the reference LAI map. Figure 6 shows a comparison of LAI MODIS , LAI EDBM , and reference LAI maps. Figure 7 shows the difference between estimated LAIs and reference LAI maps.Table 9 shows the mean value of LAI MODIS , LAI EDBM , and the reference LAI maps for the four days.Figures 6 and 7 and Table 9 show that, LAI EDBM properly depicted the spatial pattern of surface LAI and captured the temporal variation.On 17 June, the reference LAI was very low, but LAI MODIS and LAI EDBM were higher than the reference LAI.On 15 July, MODIS and EDBM model overestimated obviously compared to the reference LAI.The mean values of LAIMODIS, LAIEDBM, and reference LAI maps were 1.52, 1.5, and 1.12, respectively.However, the spatial pattern of LAIEDBM was better than that of LAIMODIS compared to the reference LAI maps.The mean values of LAI MODIS , LAI EDBM , and reference LAI maps were 1.52, 1.5, and 1.12, respectively.However, the spatial pattern of LAI EDBM was better than that of LAI MODIS compared to the reference LAI maps.On 7 August, MODIS slightly undrestimated LAI, the mean difference between LAI MODIS and reference LAI map was −0.19, while EDBM model overestimated a little.Difference between LAI EDBM and LAI MODIS was 0.33.
On 22 August, MODIS underestimated over large areas; the mean LAI MODIS of the map was 1.27, whereas the mean values of LAI EDBM and the reference map were 2.19 and 1.91, respectively.
For all four days, the differences between LAI EDBM and the reference LAI maps were spatially stable, with no substantial overestimation or marked underestimation.Difference between LAI MODIS and the reference LAI maps were more acute than difference between LAI EDBM and the reference LAI maps.Especially on 15 July, MODIS has an obvious overestimation of 0.3 and on 22 August, there is an obvious underestimation of 0.64.
Remote Sens. 2017, 9, x FOR PEER REVIEW 14 of 19 On 7 August, MODIS slightly undrestimated LAI, the mean difference between LAIMODIS and reference LAI map was −0.19, while EDBM model overestimated a little.Difference between LAIEDBM and LAIMODIS was 0.33.
On 22 August, MODIS underestimated over large areas; the mean LAIMODIS of the map was 1.27, whereas the mean values of LAIEDBM and the reference map were 2.19 and 1.91, respectively.
For all four days, the differences between LAIEDBM and the reference LAI maps were spatially stable, with no substantial overestimation or marked underestimation.Difference between LAIMODIS and the reference LAI maps were more acute than difference between LAIEDBM and the reference LAI maps.Especially on 15 July, MODIS has an obvious overestimation of 0.3 and on 22 August, there is an obvious underestimation of 0.64. Figure 8 shows a histogram of LAI differences among LAIMODIS, LAIEDBM, and the reference LAI maps.The first row is the differences between LAIEDBM and the reference LAI, and the second row is the differences between LAIMODIS and the reference LAI maps.Each column shows the difference for Figure 8 shows a histogram of LAI differences among LAI MODIS , LAI EDBM , and the reference LAI maps.The first row is the differences between LAI EDBM and the reference LAI, and the second row is the differences between LAI MODIS and the reference LAI maps.Each column shows the difference for one day.The histogram illustrates that on 17 June, and 15 July EDBM performed similar with the MODIS algorithm.On 7 and 22 August, the error distribution of EDBM model tends to be normal, while the MODIS has obvious underestimation on these two days.
Table 9 shows the mean value of LAI MODIS , LAI EDBM , and the reference LAI map on four dates.The differences between LAI EDBM and the reference LAI maps on 17 June, 15 July, 7 August, and 22 August were 0.35, 0.21, 0.3, and 0.24, respectively, which meets the 0.5 accuracy requirement for LAI use.
Remote Sens. 2017, 9, x FOR PEER REVIEW 15 of 19 one day.The histogram illustrates that on 17 June, and 15 July EDBM performed similar with the MODIS algorithm.On 7 and 22 August, the error distribution of EDBM model tends to be normal, while the MODIS has obvious underestimation on these two days.Table 9 shows the mean value of LAIMODIS, LAIEDBM, and the reference LAI map on four dates.The differences between LAIEDBM and the reference LAI maps on 17 June, 15 July, 7 August, and 22 August were 0.35, 0.21, 0.3, and 0.24, respectively, which meets the 0.5 accuracy requirement for LAI use.

Impact of Surface Inhomogeneity
The DBM was constructed at pixel scale.Theoretically, it can avoid the problem of surface inhomogeneity as all the data are the pixel scale value.The proposed method combines sites with the same LAI level to construct the UDBM model and then uses this model for sites with the same land cover type or the same LAI level.Although homogeneity varies from site (pixel) to site (pixel), which may produce some inaccuracy during LAIUDBM estimation, the UDBM model can still supply the LAI estimation background.
Furthermore, the inherent spatial resolution of the MODIS 500 m reflectance data varies according to the viewing geometry [50,51].When reflectance data for multiple days are composed to generate LAI product, the spatial characteristics are even more complicated.Considering that the UDBM model requires training data that can represent the situation at the pixel level, the spatial discrepancy between reflectance data from different days and the spatial discrepancy between MODIS reflectances and LAI data were not taken into account.The results indicate that the UDBM model can provide the LAI estimation background.

Impact of Errors Associated with the MODIS LAI Product and SG Filter on UDBM Modeling
The MODIS LAI was produced using a 3D radiative-transfer model at an 8-day interval [45].If no candidate biome/canopy models passed the comparison test for a given pixel, a backup algorithm was triggered to estimate the LAI Normalized Difference Vegetation Index (NDVI), but the estimation procedure failed occasionally, and therefore the MODIS LAI fluctuated frequently for most sites [24].In general, the LAI variation pattern in the time series was continuous according to the crop growth rule.This LAI fluctuation, which can be considered as a retrieval error to some degree, had to be eliminated to obtain the LAI time series.
In the present study, an SG filter was implemented to obtain a smooth LAI profile.Only highquality MODIS LAI data were filtered to minimize the error.The smoothed LAI profile adequately represented the LAI variation pattern and was then used to construct the UDBM model.Use of high-

Impact of Surface Inhomogeneity
The DBM was constructed at pixel scale.Theoretically, it can avoid the problem of surface inhomogeneity as all the data are the pixel scale value.The proposed method combines sites with the same LAI level to construct the UDBM model and then uses this model for sites with the same land cover type or the same LAI level.Although homogeneity varies from site (pixel) to site (pixel), which may produce some inaccuracy during LAI UDBM estimation, the UDBM model can still supply the LAI estimation background.
Furthermore, the inherent spatial resolution of the MODIS 500 m reflectance data varies according to the viewing geometry [50,51].When reflectance data for multiple days are composed to generate LAI product, the spatial characteristics are even more complicated.Considering that the UDBM model requires training data that can represent the situation at the pixel level, the spatial discrepancy between reflectance data from different days and the spatial discrepancy between MODIS reflectances and LAI data were not taken into account.The results indicate that the UDBM model can provide the LAI estimation background.

Impact of Errors Associated with the MODIS LAI Product and SG Filter on UDBM Modeling
The MODIS LAI was produced using a 3D radiative-transfer model at an 8-day interval [45].If no candidate biome/canopy models passed the comparison test for a given pixel, a backup algorithm was triggered to estimate the LAI Normalized Difference Vegetation Index (NDVI), but the estimation procedure failed occasionally, and therefore the MODIS LAI fluctuated frequently for most sites [24].In general, the LAI variation pattern in the time series was continuous according to the crop growth rule.This LAI fluctuation, which can be considered as a retrieval error to some degree, had to be eliminated to obtain the LAI time series.
In the present study, an SG filter was implemented to obtain a smooth LAI profile.Only high-quality MODIS LAI data were filtered to minimize the error.The smoothed LAI profile adequately represented the LAI variation pattern and was then used to construct the UDBM model.Use of high-quality LAI data and smoothing guaranteed continuity of the training data (the data for UDBM training) and modeling accuracy.

Errors Introduced by the Radiative-Transfer Model
When estimating LAI for the non-growing period, the LAI value may fluctuate, and it is difficult to obtain a smooth LAI profile, as had been predicted (there are no leaves left, and therefore the LAI value should be stable).This problem may be caused by the use of the PROSAIL model, which is mainly related to leaf biochemical content, but during the non-growth period, few leaves are left.Yin, Li [52] found that the canopy structure representation of the radiative transfer (RT) model has substantial implications for LAI.

Conclusions
The MODIS LAI suffers from cloud contamination, instrument problems, and retrieval algorithm inaccuracies and therefore has some missing data and fluctuations.Studies based on the database mechanistic modeling approach have demonstrated the superiority of LAI time-series estimation, especially for avoiding the problem of surface inhomogeneity.However, modeling failures where the LAI variation pattern changes due to climate change and other human activities may stymie the use of this approach in large-area LAI estimation.This study has proposed a UDBM model to supply the LAI dynamic and estimation background.Time-series MODIS observations of typical land cover types are combined to construct a specific UDBM model representing the main LAI dynamics for pixels with the same land cover type or the same LAI level.In this manner, regional LAI background information can be obtained.Once bidirectional information has been introduced, the ENKF method can be used to update the LAI to obtain the final estimates, which are very close to the LAI reference at sites with various land covers.Compared with field reference observations, the proposed method had high accuracy, with RMSE = 0.5, bias = 0.12, and MAE = 0.30 (the corresponding values for the MODIS product were 1.26, −0.22, and 0.98, respectively).Four LAI maps at the La Albufera, Spain site were generated using the proposed EDBM method, and the results were compared with the MODIS data and the LAI reference map.The EDBM estimated LAI properly depicted the spatial pattern of surface LAI and captured the temporal variation.The average difference between EDBM and the LAI reference map on the selected four days was 0.32 and meet the LAI accuracy requirement of 0.5 for general use.
Nevertheless, this method still has drawbacks.Only two UDBM models (UDBM 1 for forest land cover type and UDBM 2 for the others) were constructed in this study.Obtaining more accurate background estimates will require development of a more reasonable method of modeling the LAI background by constructing a UDBM model for every land cover type based on a detailed LAI level classification.

Figure 1 .
Figure1.Flowchart of the proposed estimation approach.The UDBM model is constructed using time-series reflectance and LAI data, the PROSAIL model is executed to simulate reflectance, and an ensemble Kalman filter is used to update LAIUDBM according to the difference between observed and simulated reflectances to obtain the final estimation.

Figure 1 .
Figure1.Flowchart of the proposed estimation approach.The UDBM model is constructed using time-series reflectance and LAI data, the PROSAIL model is executed to simulate reflectance, and an ensemble Kalman filter is used to update LAI UDBM according to the difference between observed and simulated reflectances to obtain the final estimation.

Figure 3 .
Figure 3. Location and land cover types of La Albufera site.The left figure is the location of the research area, the middle figure is the the 20 km × 20 km research area, the right figure is the land cover map according to the GlobeLand30 classification system.The ocean water is not included in the classification system, and was added by the authors.

Figure 3 .
Figure 3. Location and land cover types of La Albufera site.The left figure is the location of the research area, the middle figure is the the 20 km × 20 km research area, the right figure is the land cover map according to the GlobeLand30 classification system.The ocean water is not included in the classification system, and was added by the authors.

Figure 4 .
Figure 4. (a-g) LAI temporal profile comparison of MODIS LAI (LAIMODIS), the LAI estimation background (LAIUDBM), and the proposed EDBM method (LAIEDBM).The site names and land cover types of sites are listed below each figure.The green dashed line is LAIMODIS, the magenta line is LAIUDBM, the red line is LAIEDBM, and the gray vertical lines are the estimation error of the EDBM model.

Figure 4 .
Figure 4. (a-g) LAI temporal profile comparison of MODIS LAI (LAI MODIS ), the LAI estimation background (LAI UDBM ), and the proposed EDBM method (LAI EDBM ).The site names and land cover types of sites are listed below each figure.The green dashed line is LAI MODIS , the magenta line is LAI UDBM , the red line is LAI EDBM , and the gray vertical lines are the estimation error of the EDBM model.

Figure 5
Figure 5 shows a scatter diagram of the LAI REF , LAI MODIS , and LAI EDBM .The results show that LAI MODIS was the most discrete.LAI EDBM fitted LAI REF very well because most of the points are very near to the y = x line.Statistical results for the LAI estimation errors are listed in Table8.The root mean square error (RMSE), bias, and mean absolute error (MAE) of LAI EDBM were 0.5, 0.12, and 0.30, respectively, which indicates that high estimation accuracy can be achieved using the proposed EDBM method.

2 )Figure 5 .
Figure 5. Scatter diagram of the LAI reference and estimated LAIs.The green dots are LAIMODIS, and the red triangles are LAIEDBM.

Figure 5 .
Figure 5. Scatter diagram of the LAI reference and estimated LAIs.The green dots are LAI MODIS , and the red triangles are LAI EDBM .

Figure 6 .
Figure 6.Map comparison of LAIEDBM, LAIMODIS, and the reference LAI map.In each row, the leftmost map is LAIMODIS, and the second and third maps are LAIEDBM and reference LAI map.Each row is one day's data.

Figure 6 .
Figure 6.Map comparison of LAI EDBM , LAI MODIS , and the reference LAI map.In each row, the leftmost map is LAI MODIS , and the second and third maps are LAI EDBM and reference LAI map.Each row is one day's data.

Figure 7 .
Figure 7. Map of differences among LAIMODIS, LAIEDBM, and the reference LAI maps.In each row, the left map shows the difference between LAIMODIS and reference LAI map, and the right map shows the difference between LAIEDBM and reference LAI map.Each row contains one day's data.

Figure 7 .
Figure 7. Map of differences among LAI MODIS , LAI EDBM , and the reference LAI maps.In each row, the left map shows the difference between LAI MODIS and reference LAI map, and the right map shows the difference between LAI EDBM and reference LAI map.Each row contains one day's data.

Figure 8 .
Figure 8. Histogram of LAI differences among LAIMODIS, LAIEDBM, and the reference LAI map.

Figure 8 .
Figure 8. Histogram of LAI differences among LAI MODIS , LAI EDBM , and the reference LAI map.

Table 1 .
PROSAIL model parameters and their ranges.

Table 2 .
is the spectral response function on wavelength λ of each MODIS band.Reflectance with the same geometry of MODIS reflectance data is simulated.Spectral response functions of MODIS bands 1, 2, and 7. Band range of MODIS bands 1, 2, and 7.

Table 3 .
Initial values and variances of the ENKF state variables.

Table 4 .
Information on the forest sites selected for UDBM 1 modeling.

Table 5 .
Information on the non-forest sites selected for UDBM 2 modeling.

Table 6 .
Characteristics of the sites selected for model validation.

Table 7 .
Date matching of the reference and MODIS LAI maps.

Table 7 .
Date matching of the reference and MODIS LAI maps.
The final estimation fitted LAI REF very well, with a mean difference between LAI EDBM and LAI REF of 0.11.At the Sud-Ouest site, LAI MODIS was 1.1, which is very small.LAI UDBM is also very small and consistent with the MODIS data, but much lower than LAI REF .The LAI EDBM upgraded the estimates, and fitted LAI REF well, but the LAI profile is very different with LAI MODIS .Remote Sens. 2017, 9, x FOR PEER REVIEW

Table 8 .
Statistical results for LAIMODIS and LAIEDBM compared with LAIREF.

Table 8 .
Statistical results for LAI MODIS and LAI EDBM compared with LAI REF .

Table 9 .
Mean values of LAI MODIS , LAI EDBM , and LAI REF at La Albufera site.