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

Abstract

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.

1. 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.

Various strategies have been used to improve LAI temporal continuity and estimation accuracy. Researchers have focused on using crop growth models to supply the LAI estimation background, and with the coupled use of radiative transfer models and data assimilation methods, have made great progress in improving LAI estimation accuracy and yield forecasting [5,6,7,8,9,10,11,12]. The crop growth models that supplied the crop growth background and LAI dynamics in their work guaranteed estimation continuity, but were feasible only if the land cover units (crop areas) in the specific regions studied were homogeneous enough that crop-specific biophysical variables could be retrieved and that dynamic changes could be characterized by crop growth models. Bastiaanssen and Ali [13] and Mo, Liu [14] found that the results deteriorated in areas where land cover was fragmented and several types of crops were grown within one pixel.

Attempts have also been made to use data-driven dynamic models instead of crop growth models to overcome pixel-scale heterogeneity. Xiao, Wang [15,16], Jiang, Liang [17], and Li, Mao [18] used temporal observation profiles to model LAI dynamics with statistical methods and further used these data in a data assimilation scheme [9]. Chen, Wang [19] developed a data-based mechanistic (DBM) approach based on the methods of Young [20,21] to estimate LAI time series. Guo, Wang [22] further improved this approach by including bidirectional information in the reflectance time series. Although all these studies used pixel-scale remote-sensing data for LAI dynamic modeling, which can theoretically avoid the problem of surface inhomogeneity, pixel-scale modeling may fail if sudden changes in land-surface conditions occur or if the LAI variation pattern changes due to climate conditions or human activity. A more stable and feasible LAI model is necessary for high-accuracy LAI time-series estimation.

This paper proposes an extended data-based mechanistic (EDBM) method for LAI time-series estimation by using a universal data-based mechanistic (UDBM) model to provide the LAI estimation background. The UDBM model is constructed based on the data time series of the same LAI-level pixels using the methods proposed by Young [20,21] and Chen, Wang [19] and avoids the failure of pixel-scale LAI background estimation. The ENKF assimilation method is then used to update LAI recursively by coupling with a canopy radiative-transfer model (PROSAIL).

The following section presents an introduction to the UDBM model construction. Next, the ENKF data assimilation scheme for estimating LAI using MODIS bidirectional reflectance time-series data is outlined. Section 3 demonstrates the performance of the method using actual LAI values at several Benchmark Land Multisite Analysis and Intercomparison of Products (BELMANIP) sites [23,24]. A comparison with LAI reference maps from the Imagine-S project [25] at the La Albufera, Spain site is also presented. The final section provides a brief discussion and presents conclusions.

2. 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.

2.1. UDBM Model Construction

2.1.1. 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):

$$y_t=u_{1,t}\frac{B_1\left(L\right)}{A\left(L\right)}+u_{2,t}\frac{B_2\left(L\right)}{A\left(L\right)}+\cdots +u_{p,t}\frac{B_p\left(L\right)}{A\left(L\right)}+{\epsilon }_t$$

(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 ${\epsilon }_t$ is the model estimation error. The forms of $A$ and $B$ are shown in Equation (2):

$$\begin{gathered} A\left(L\right)=1+a_{1}L+\cdots +a_n{L}^n \\ {B}_1\left(L\right)=b_{10}+b_{11}L+\cdots +b_{1i}{L}^i \\ {B}_2\left(L\right)=b_{20}+b_{21}L+\cdots +b_{2j}{L}^j \\ {B}_p\left(L\right)=b_{p0}+b_{p1}L+\cdots +b_{pk}{L}^k \end{gathered}$$

(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):

$$LA{I}_t=\frac{b_{10}+b_{11}L+\cdots +b_{1i}{L}^i}{1+a_{1}L+\cdots +a_n{L}^n}{R}_t^1+\frac{b_{20}+b_{21}L+\cdots +b_{2j}{L}^j}{1+a_{1}L+\cdots +a_n{L}^n}{R}_t^2+\frac{b_{30}+b_{31}L+\cdots +b_{3k}{L}^k}{1+a_{1}L+\cdots +a_n{L}^n}{R}_t^7+{\epsilon }_t$$

(3)

where the $LA{I}_t$ is the estimated LAI at time t and $R_t^1$, $R_t^2$, $R_t^7$ are the MODIS observed reflectances of bands 1, 2, and 7 respectively.

Assuming that $L^{n}LA{I}_t=LA{I}_{t-n}$, Equation (3) can be rewritten as Equation (4):

$$\begin{gathered} LA{I}_t=\left[b_{10},b_{11},\cdots ,b_{1i}\right]\left[\begin{array}{c}{R}_t^1\\ R_{t-1}^1\\ ⋮\\ R_{t-i}^1\end{array}\right]+\left[b_{20},b_{21},\cdots ,b_{2j}\right]\left[\begin{array}{c}{R}_t^2\\ R_{t-1}^2\\ ⋮\\ R_{t-j}^2\end{array}\right]+\left[b_{30},b_{31},\cdots ,b_{3k}\right]\left[\begin{array}{c}{R}_t^7\\ R_{t-1}^7\\ ⋮\\ R_{t-k}^7\end{array}\right] \\ -\left[a_1,a_2,\cdots ,a_n\right]\left[\begin{array}{c}LA{I}_{t-1}\\ LA{I}_{t-2}\\ ⋮\\ LA{I}_{t-n}\end{array}\right]+{\epsilon }_t \end{gathered}$$

(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²) and the Young information criterion (YIC) [27]. R² 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² and the lowest YIC. Once the DBM model has been constructed, it can predict future LAI values based on known historical time-series data.

2.1.2. 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₁ for forest and UDBM₂ for other land cover types. The structures of the two UDBM models are given in Equations (5) and (6).

$$\begin{gathered} LA{I}_{UDBM1,t}=-4.13R_t^1+2.969R_{t-1}^1+1.099R_{t-2}^1+4.081R_t^2-4.017R_{t-1}^2-1.272R_t^7 \\ +2.587R_{t-1}^7-0.9419R_{t-2}^7+1.7LA{I}_{t-1}-0.719LA{I}_{t-2} \end{gathered}$$

(5)

$$\begin{gathered} LA{I}_{UDBM2,t}=-2.3453R_t^1+3.0016R_{t-1}^1-0.7932R_t^2+2.9223R_{t-1}^2-3.7622R_{t-2}^2 \\ +1.0405R_t^7-0.713R_{t-1}^7+0.6466R_{t-2}^7+2.0519LA{I}_{t-1} \\ -1.5187LA{I}_{t-2}+0.4308LA{I}_{t-3} \end{gathered}$$

(6)

2.2. 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. PROSAIL model parameters and their ranges.

Model	Parameters	Symbol	Values	Unit
Prospect	Chlorophyll content *	Cab	30	μg cm−2
	Carotenoids content	Car	10	μg cm−2
	Total pigment content	Cbp	0	-
	Water equivalent thickness	Cw	0.015	cm
	Dry matter content *	Cm	0.00125–0.00625	μg cm−2
	Leaf structure index *	$N$	1–2.5	-
Sail	Leaf area index *	LAI	0–8	-
	Mean leaf angle *	ALA	40–85	deg
	Hot spot *	SL	0.01–1	-
	Soil factor	ρs	0.2	-
	Scatter light proportion	SKYL	0.1	-
	Solar zenith angle	θs	0–90	deg
	Observation zenith angle	θv	0–90	deg
	Relative azimuth angle	φsv	0–180	deg

* sensitive parameters.

. 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. Band range of MODIS bands 1, 2, and 7.

Band Number	Band Range (nm)
1	620–670
2	841–876
7	2105–2155

and Figure 2). The conversion formula is given in Equation (7):

$$Ref=\frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{R}_{\lambda }·f\left(\lambda \right)}{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}f\left(\lambda \right)}$$

(7)

where $R_{\lambda }$ is the reflectance calculated using PROSAIL with wavelength λ and $f\left(\lambda \right)$ is the spectral response function on wavelength λ of each MODIS band. Reflectance with the same geometry of MODIS reflectance data is simulated.

2.3. 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 $\mathrm{A}\in R^{n\times N}$ that contains N ensemble members of the model state variable set; $A^{\prime }\in R^{n\times N}$ is the matrix of ensemble perturbations (the ensemble minus its mean); the ensemble covariance matrix $P_e\in R^{n\times n}$ is defined as Equation (8):

$$P_e=\frac{A^{\prime }{\left(A^{\prime }\right)}^T}{N-1}$$

(8)

For each measurement, we define the perturbed observation as Equation (9):

$$d_j=d+{ϵ}_{j}j=1,\dots N$$

(9)

Then the observation ensemble is given as Equation (10):

$$D=\left(d_1,d_2,\dots ,d_N,\right)\in R^{m\times N}$$

(10)

The observation error covariance matrix derived from the observation perturbation can be expressed as Equation (11):

$$R_e=\frac{\gamma {\gamma }^T}{N-1}\in R^{m\times n}$$

(11)

where $\gamma =\left({ϵ}_{1,}{ϵ}_{2,}\dots {ϵ}_N\right)\in R^{m\times N}$ and $H\in R^{m\times n}$ is the observation operator. The standard analytical equation for ENKF is as Equation (12):

$$A^a=A+A^{\prime }{A}^{\prime T}{H}^T{(H{A}^{\prime }{A}^{\prime T}{H}^T+R_e)}^{-1}\left(D-HA\right)$$

(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):

$$A^a=A+A^{\prime }{\widehat{A}}^{\prime T}{\widehat{H}}^T{(\widehat{H}{\widehat{A}}^{\prime }{\widehat{A}}^{\prime T}{\widehat{H}}^T+R_e)}^{-1}\left(D-\widehat{H}\widehat{A}\right)$$

(13)

where $\widehat{A}\in R^{\widehat{n}\times N}$ and ${\widehat{A}}^{\prime }\in R^{\widehat{n}\times N}$ are formed by augmenting the state variable ensemble by the predicted reflectances and $\widehat{H}\in R^{m\times \widehat{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 ($\widehat{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. Initial values and variances of the ENKF state variables.

Parameter	LAI	Cab	Cw	Cm	ALA	Ps
Initial value	1.0	30	0.01	0.001	70	0.2
Variance	0.35	6	0.001	0.0001	9	0.001

.

3. Data

3.1. 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. Information on the forest sites selected for UDBM₁ modeling.

Site Name	Latitude	Longitude	Land Cover Type	Site Name	Latitude	Longitude	Land Cover Type
Lost Creek	46.08	–89.97	Deciduous broadleaf forest	Soignes	50.78	4.42	Deciduous broadleaf forest
Malga Arpaco	46.11	11.70	Mixed forest	Hainich	51.07	10.45	Deciduous broadleaf forest
Sylvania Wilderness Area-Michigan	46.24	–89.34	Mixed forest	Brasschaat (De Inslag Forest)	51.30	4.52	Deciduous broadleaf forest
New Brunswick-Nashwaak Lake 1	46.47	–67.1	Mixed forest	Leinefelde	51.32	10.36	Deciduous broadleaf forest
Neustift/Stubai Valley	47.11	11.31	Mixed forest	Sask-SSA Old Aspen	53.62	–106.19	Mixed forest
Gilching (VALERI)	48.08	11.32	Mixed forest	BOREAS SSA Young Aspen	53.65	–105.32	Mixed forest
Ontario-Groundhog River-Mature Boreal Mixed Wood	48.21	–82.15	Mixed forest	Kannenbruch Forest	53.78	10.6	Mixed forest
Hesse Forest-Sarrebourg	48.67	7.06	Deciduous broadleaf forest	Whitecourt	54.04	–115.79	Mixed forest
Vielsalm	50.30	5.99	Mixed forest	Sask-Fire 1989	54.25	–105.87	Mixed forest
COMPLET 160	50.66	27.89	Mixed forest	Jarvselja (VALERI)	58.29	27.26	Mixed forest

) and non-forest sites (

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

Site Name	Latitude	Longitude	Land Cover Type	Site Name	Latitude	Longitude	Land Cover Type
COMPLET 163 (56.1129°N, −69.3589°W)	56.11	–69.35	savannas	Quebec	54.5	–75.7	shrubs
UCI-1998 burn site	56.63	–99.94	savannas	Audubon Research Ranch-Arizona	31.59	–110.51	shrubs
Quebec Boreal Cutover Site	49.26	–74.03	savannas	Santa Rita Mesquite-Arizona	31.82	–110.86	shrubs
COMPLET	0.58	14.83	croplands	Haouz (VALERI)	31.65	–7.60	shrubs
COMPLET 133	8.57	19.91	savannas	Sky Oaks-Young Stand-California	33.37	–116.62	shrubs
Qianyanzhou	26.73	115.06	savannas	Jornada LTER-New Mexico (JRN1)	32.59	–106.84	shrubs
Ilorin	8.32	4.34	savannas	Maricopa Agricultural Center-Arizona	33.07	–111.97	croplands
Sardinilla Pasture	9.30	–79.63	croplands	Sud-Ouest (VALERI)	43.506	1.23	croplands
Tonzi Ranch-California	38.43	–120.96	savannas	Bondville-Illinois	40.00	–88.29	croplands
Freeman Ranch-Grassland-Texas	29.93	–98.01	savannas	Mead-irrigated maize-soybean rotation site-Nebraska	41.16	–96.47	croplands

) 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.

3.2. LAI Reference Data

3.2.1. Field Data

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 ($LA{I}_{effective}$) to true LAI ($LA{I}_{true}$) with Equation (14).

$$LA{I}_{true}=\frac{LA{I}_{effective}}{C_{index}}$$

(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. Characteristics of the sites selected for model validation.

Site Name	Latitude	Longitude	Land Cover Types	Year	DOY	LAI
HARV	42.53	–72.17	Mixed forest	2000	170	5.1
				2000	217	5.0
				2001	208	5.5
				2002	236	5.4
CHEQ	45.95	90.27	Mixed forest	2002	220	3.05
KONZ	39.09	–96.57	Grasses	2000	159	2.0
				2000	239	2.0
				2001	169	2.9
				2001	228	2.5
Sud-Ouest	43.51	1.24	Grasses and cereal crops	2002	201	2.3
ARGO	40.01	–88.29	Broadleaf crops	2000	186	2.5
				2000	224	3.6
Alpilles	43.81	4.74	Broadleaf crops	2002	201	1.7
Larzac	43.94	3.12	Savannahs	2002	193	0.9

.

3.2.2. 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. Date matching of the reference and MODIS LAI maps.

Map Name	Date Match
MODIS LAI DOY	169	193	217	233
Reference LAI date	17 June	15 July	7 August	22 August

.

4. Results

4.1. 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. 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.

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 LAI_REF better, but gave a higher peak value during the growing season than the LAI_MODIS.

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.

4.2. Comparison with LAI Reference

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

Table 8. Statistical results for LAI_MODIS and LAI_EDBM compared with LAI_REF.

	LAIMODIS (m2/m2)	LAI EDBM (m2/m2)
RMSE	1.26	0.50
BIAS	–0.22	0.12
MAE	0.98	0.30

. 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.

4.3. 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. Mean values of LAI_MODIS, LAI_EDBM, and LAI_REF at La Albufera site.

Date (2014)	LAIMODIS (m2/m2)	LAIEDBM (m2/m2)	LAIREF (m2/m2)
17 June	0.44	0.59	0.30
15 July	1.52	1.50	1.12
7 August	1.5	2.02	1.69
22 August	1.27	2.19	1.91

shows the mean value of LAI_MODIS, LAI_EDBM, and the reference LAI maps for the four days. Figure 6 and Figure 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 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.

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.

5. Discussion

5.1. 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.

5.2. 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.

5.3. 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.

6. 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₁ for forest land cover type and UDBM₂ 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.