Generating a Spatio-Temporal Complete 30 m Leaf Area Index from Field and Remote Sensing Data

Abstract

The leaf area index (LAI) is an important parameter for vegetation monitoring and land surface ecosystem research. Although a variety of LAI products have been generated, the moderate to coarse spatial resolution and low temporal resolution of these products are insufficient for regional-scale analysis. In this study, a modified ensemble Kalman filter model (MEnKF) was proposed to generate spatio-temporal complete 30 m LAI data. High-quality, filtered historical Moderate-resolution Imaging Spectroradiometer (MODIS) LAI data were used to obtain the LAI background, and an LAI temporal dynamic model was constructed based on it. An improved back-propagation (BP) neural network based on a simulated annealing algorithm (SA-BP) was constructed with paired Landsat surface reflectance data and field LAI data to generate a 30 m LAI. The MEnKF was used to estimate the spatio-temporal complete LAI beginning from the LAI peak value position where Landsat observations were available. The spatio-temporal 30 m LAI was estimated in farmland (Pshenichne), grassland (Zhangbei), and woodland (Genhe) sites. The results indicate that the MEnKF-estimated LAI is consistent with the field measurements for all sites (the coefficient of determination (${\mathrm{R}}^2$) = 0.70; root mean squared error (RMSE) = 0.40) and is better than that of the conventional sequence data assimilation algorithm (${\mathrm{R}}^2$ = 0.40; RMSE = 0.78). The regional LAI captures the vegetation growth pattern and is consistent with the Landsat LAI, with an ${\mathrm{R}}^2$ larger than 0.65 and an RMSE less than 0.51. The proposed MEnKF algorithm, which effectively avoids error accumulation in the data assimilation scheme, is an efficient method for spatio-temporal complete 30 m LAI estimation.

1. Introduction

The leaf area index (LAI), defined as the sum of green leaf area per unit of land area, is an important parameter for land surface ecosystems [1]. It provides dynamic information on vegetation growth and provides environment-related policy decision support [2]. High-quality LAI data are important to agricultural, ecological, environmental, and climate change studies [3].

Field LAI measurement using direct and indirect measurement methods [4] are reliable ways to obtain accurate LAI data. However, the time- and labor-consuming nature of field measurements hinders long time-series LAI observations in a large area with traditional optical instruments [5,6]. A more convenient and rapid method is needed to obtain the regional time-series LAI [7]. Remote sensing technology, which has a series of advantages, such as wide coverage, low cost of human and material resources, and fine spatial and temporal resolution, is the most effective means to obtain long time series LAI at regional and global scales [8,9].

There are many widely used remote sensing LAI products, such as Moderate-resolution Imaging Spectroradiometer (MODIS) LAI [10], CYCLOPES LAI [11], Advanced Very High Resolution Radiometer (AVHRR) LAI [12], polarization and directionality of the Earth’s reflectance (POLDER) LAI [13], GEOLAND2 Version 1 (GEOV1) LAI [14], and Global Land Surface Satellite (GLASS) LAI [15,16]. Although these products have been widely used in various fields, their coarse-moderate spatial resolution and missing data have hindered their application at the regional scale. Scholars have made various efforts to improve the spatial resolution and temporal continuity of remote sensing LAI data.

Estimating LAI directly from high spatial resolution remote sensing data is an efficient way to improve the spatial resolution of LAI products. Researchers have used the empirical relationship of the vegetation index (VI) to calculate LAI with different fine resolution multispectral sensor data, such as Landsat, WorldView-2, and Sentinel [5,17,18]. S Chaurasia et al. [19] used an empirical model to estimate LAI with the normalized difference vegetation index (NDVI) and ratio vegetation index (RVI) and found that the NDVI-LAI wheat empirical model is suitable for early wheat growth stages and that RVI-LAI is suitable for later wheat growth stages. Nguy-Robertson et al. [20] tested the performance of fifteen VIs in estimating the green LAI of C3 and C4 crops and found that VIs using the red edge and near infrared (NIR) spectral bands were the most accurate for estimating green LAI in both C3 and C4 crops. The VI-induced empirical model is simple and easy to construct, but it is difficult to extend to other research areas.

To obtain a universal LAI estimation model, researchers have used the radiative transfer (RT) model to estimate LAI with high spatial resolution data [21,22,23]. M Meroni [24] used a 1-D bidirectional canopy reflectance scattering by arbitrarily inclined leaves (SAIL) model coupled with the leaf optical properties spectra model (PROSPECT) to retrieve LAI in a poplar plantation with hyperspectral digital airborne imaging spectrometer (DAIS) data. S Punalekar et al. [25] combined proximal hyperspectral and Sentinel 2 data with a PROSAIL radiative transfer (RT) model to estimate LAI. RT models that are based on the physical theory of light and particle propagation are more accurate and suitable for different situations. However, the model structures are complicated, and multiple input parameters are needed. To make the calculation fast and suitable for large-scale LAI estimation, scholars have used the RT model to perform extensive simulations and construct a lookup table (LUT) for rapid LAI estimation [26,27]. K Richter et al. [28] established an LUT based on the well-established SAILH + PROSPECT RT model for homogeneous canopies to estimate LAI using Sentinel-2 satellite data.

With the rapid development of machine learning algorithms, the neural network method (NN) has become a widely used LAI estimation method [5,27,29,30]. R Prasad [31] retrieved spinach crop parameters, including LAI with back-propagation artificial neural network (BPANN), and the results showed that BPANN is practical for nonlinear datasets. Xue et al. [32] developed a BP neural network improved by a simulated annealing algorithm (SA-BP) model for 30 m LAI estimation, and a comparison of the results with field observations showed that it is an effective method for 30 m resolution LAI estimation. Galal Omer et al. [5] predicted the LAI of six endangered tree species using support vector machines (SVM) and artificial neural networks (ANN) based on WorldView-2 data and demonstrated the dependability of these two methods in predicting LAI using fragmented forest data.

LAI estimated directly from remotely sensed data is always discontinuous due to data missing and retrieval failure [33,34,35]. Combining the high temporal resolution information from coarse-resolution data and fine spatial resolution information is a common way to obtain the spatio-temporal continuous 30 m LAI. Data fusion and data assimilation are two popular methods to improve the spatio-temporal resolution of LAI estimation. The spatial and temporal adaptive reflectance fusion model (STARFM) algorithm is a widely used data fusion method for fine spatio-temporal land surface reflectance, LAI, albedo, net primary productivity (NPP), and NDVI [36,37,38,39] estimation. However, the accuracy was affected by the choosing of input data pairs [40] and land surface heterogeneity [41].

Data assimilation is efficient for time-series LAI estimation [42,43,44]. Ling et al. [45] compared four sequential assimilation algorithms, including the Kalman Filter (KF), the Ensemble Kalman Filter (EnKF), the Ensemble Adjust Kalman Filter (EAKF), and the Particle Filter (PF) in LAI estimation, and found EnKF and EAKF performed best. Xiao et al. [46] used the EnKF assimilation algorithm to get real-time LAI time series estimation from MODIS reflectance data, which efficiently filled the MODIS LAI product data missing. Zhou et al. [47] further improved the estimation algorithm by substituting the crop growth model with an extensive data mechanistic model in the EnKF assimilation scheme. The EnKF assimilation method couples the LAI dynamic/growing model with LAI observation data and is capable of obtaining a reliable time series LAI product with a small number of observations, but the estimation is inaccurate before the first observation is introduced into the data assimilation scheme.

This study aims to propose a Modified Ensemble Kalman Filter (MEnKF) model to improve the accuracy of the 30 m LAI time series estimation. In contrast to the regular EnKF, the MEnKF starts from the date near the peak vegetation growing period where LAI observations exist. Estimation is then performed forward and backward to complete the entire vegetation growing season. The MEnKF avoids a large error in the period when no observation data are introduced. It is efficient in fine spatial resolution LAI time series estimation with high accuracy.

2. Data and Methods

2.1. Research Area

In this study, three research areas with different land surface types, including farmland, woodland, and grassland, were selected according to the availability of field observation data. They are Pshenichne in Ukraine and Zhangbei and Genhe in China.

The Pshenichne study area is located at the Pshenichne site in Ukraine; it is a farmland site and belongs to the 7th Framework Programme (FP7). The center point is 50°04′ N and 30°13′ E. The main land cover types are farmland and artificial land, with little grassland, wetland, and water. The area is mainly plain, with a temperate continental climate. The average temperature is −5.8 °C in January and 19.5 °C in July. The average annual precipitation is 622 mm. Regional crops mainly include winter wheat, spring barley, corn, soybean, and other crops. Due to the variety of crops, there is no typical crop rotation in this area, and the surface distribution is very complicated.

The Zhangbei research area is a grassland site located in Zhangbei County, China, and is a VALidation of European Remote Sensing Instruments (VALERI) site. The majority of the land is grassland and farmland, with a small part of artificial land and water. The site center is 41°16′ N and 114°41′ E. It has a medium-temperate continental monsoon climate. The annual average temperature is 3.2 °C, the annual precipitation is approximately 300 mm, and the average annual number of sunshine hours is 2897.8 h. The terrain of this research area is complex, mainly with hills and undulating plains, and the vegetation is mainly grass.

The Genhe experimental area is a forest site located in the northern section of the Greater Khingan Range in northeastern China. More than 90% of the area is covered by forest; in the forest-free area, grass grows, and a small amount of farmland exists in the northern part of the area. The site center is 51°25′ N and 120°34′ E. The area has a humid temperate forest climate with a continental monsoon climate. It is cold and humid, with long winters and short summers. The rainy season is generally between July and August of each year. The average annual frost-free period is 90 days, the average temperature is −6.3 °C, and the freezing period is over 210 days. The region has a series of mountains, with an average elevation of more than 1000 m and a relative height difference of 100 to 300 m. The main vegetation type in the area is forest and the dominant tree species is Larix gmelinii. The 30 m land cover map used in this study was produced by the National Geomatics Center of China [48]. The land cover maps of the three study areas are shown in Figure 1. In the experiment, only vegetation pixels were included. Water, wetland, and artificial land (lands modified by human activities, including all kinds of habitation, industrial and mining area, transportation facilities, and interior urban green zones and water bodies, etc.) were masked.

2.2. Data Preprocessing

Three kinds of datasets were used for 30 m LAI time series estimation. Landsat data combined with field observations were used for SA-BP model training, and then 30 m LAI was estimated for each Landsat scene based on the SA-BP model. The historical MODIS LAI product was filtered to construct the LAI temporal dynamic model to supply the initial LAI estimation in the data assimilation scheme. A 30 m LAI estimated from Landsat data was used as the LAI observation data.

2.2.1. Landsat Data

The land surface reflectance data from Landsat-7 Enhanced Thematic Mapper (ETM+) and Landsat-8 Operational Land Imager (OLI) images were used in this study. Landsat Level 2 ETM+ and OLI multispectral land surface reflectance data were downloaded from the USGS website. According to the acquisition time of the ground data, the data types downloaded for different study areas were different. High-quality (cloud cover less than 10%) Landsat-8 OLI data of 2014 were downloaded at Pshenichne, Landsat-7 ETM+ data of 2002 were downloaded at Zhangbei, and Landsat-8 OLI data of 2016 were downloaded at Genhe. The Landsat reflectance data used in this study are shown in

Table 1. List of Landsat data used in this study.

Research Area	Sensor	Data
Pshenichne	OLI	06/06/2014
		06/22/2014
		07/24/2014
		08/09/2014
Zhangbei	ETM+	04/27/2002
		05/29/2002
		07/07/2002
		08/17/2002
		09/09/2002
		25/09/2002
Genhe	OLI	05/13/2016
		07/07/2016
		08/01/2016
		09/18/2016

.

Landsat Level-2 surface reflectance was radiometric, geometric, and atmospheric corrected. Due to cloud and cloud shadow contamination, Landsat has some missing data. Missing data were filled based on a modified neighborhood similar pixel interpolator (NSPI) approach [33]. Reference images of the near period were used to fill the target cloud contamination image. According to the spectral temporal information and spectral spatial information, a mapping model between reference images and the target image is constructed; all the polluted pixel values were estimated from clear pixels in reference images.

2.2.2. MODIS LAI Product

The MODIS LAI products (MOD15A2H) have a spatial resolution of 500 m and temporal resolution of 8 days. In this study, MOD15A2H data were downloaded from the EARTHDATA website. For each site, six years of data, including five historical years and the year to be estimated, were downloaded. The data from five historical years were used to construct the LAI background, and the data of the estimated year were used for validation.

Note that the atmosphere effects, such as clouds, aerosols, water vapor, and ozone, cause increases in reflectance in the red band rather than in the near-infrared (NIR) band, resulting in errors in the MODIS LAI product. To account for negatively biased noise, the Multi-step S-G filtering was used to smooth the MODIS LAI data. This multi-step procedure leads to a smoothed curve adapted to the upper envelope of the LAI values in the time series; in this study, the smoothed LAI curve was used to derive the LAI background.

2.2.3. Field Data

Field observation data were collected at each site for 30 m SA-BP model construction and MEnKF estimation result validation.

The field LAI of the Pshenichne research area was collected by the 7th Framework Programme (FP7) project. The FP7 is a long-term major technology research and development program. It is the most invested global technology development plan of the European Union. The site dataset includes LAI, vegetation coverage (fCover), photosynthetically active radiation (FAPAR), and atmospheric top reflectance, which can be used for verifying the authenticity of related remote sensing products. The Pshenichne study area has a total of 141 field LAI data points collected from farmland areas. The data acquisition time was from 2013 to 2015, the acquisition method was the hemispheric digital image method (DHP), 12–15 observations were taken for each sample, and CAN-EYE software was used to calculate effective LAI.

The field LAI of the Zhangbei research area was collected at the grassland area and collated by the European VALERI project. The VALERI project was initiated by the European Aviation Administration in 2000. The site contains datasets such as LAI, vegetation canopy reflectance, and vegetation coverage, which can be used for modeling remote sensing algorithms and verification of corresponding remote sensing data products. There are 40 field LAI data points in the Zhangbei research area. The acquisition time was from July 2002 to August 2002, the acquisition method was DHP, 12 observations were taken for each sample, and CAN-EYE software was used to calculate the effective LAI.

The Genhe research area field LAI data were measured by Beijing Normal University. It is part of the continuation of the Complicate Observations and Multi-Parameter Land Information Constructions on Allied Telemetry Experiment (COMPLICATE) in the Genhe area [49]. The Genhe research area has a total of 71 field forest LAI points. The data acquisition time was 2016, and the acquisition method was a wireless sensor network (LAINet), 15 nodes was installed at one site and observations were averaged to compare with Landsat pixel values.

3. Method

To generate spatio-temporal complete 30 m resolution LAI data, the required process was as follows. (1) The SA-BP neural network model was constructed based on Landsat reflectance data and LAI field observations, 30 m LAI was estimated from Landsat land surface reflectance data using the constructed SA-BP model. The LAI estimated from Landsat data (Landsat LAI) was considered the 30 m LAI observation data in the data assimilation scheme. (2) The high-quality historical MODIS LAI product was averaged and filtered to construct the LAI background. The LAI temporal dynamic model was constructed based on the LAI background. For each site, Landsat observations near the LAI background peak value were set as the estimation beginning point, and forward and backward LAI dynamic models were constructed. (3) The MEnKF model was applied to estimate LAI at time k. The LAI dynamic model provided the initial LAI estimation. If Landsat LAI at time k exists, EnKF was applied to update the initial estimation combined with model and observation error introduced. If Landsat LAI at time k was unavailable, LAI initial estimation was considered as final estimation of time k. (4) LAI estimated at time k was input into the dynamic model for the next time LAI initial estimation. Time-series LAI was generated.

The specific steps are shown in the flowchart in Figure 2.

3.1. Landsat LAI Estimation

Landsat LAI was estimated with the simulated annealing and error back-propagation neural network (SA-BP) [50]. It allowed an improvement of the BP neural network to improve the generalization ability and iteration speed of the model.

In the SA-BP method, during the forward propagation of the neural network, for each neuron node except the input layer node, the relationship between the input and output is

$$\mathrm{net}(i)={\displaystyle {\sum }_{j=1}^n{W}_{ij}{X}_j}-\mu$$

(1)

$$Y_i=f(\mathrm{net}(i))$$

(2)

where $X$ is the signal input, $W$ is the weight of the neuron node, $\mu$ is the bias, $f$ is the excitation function, and $Y$ is the signal output.

The gradient descent method is usually used to continuously adjust the weights until the optimal solution of the weights is obtained. This study uses the simulated annealing algorithm instead of the gradient descent method to adjust the weights. The related parameter settings of this algorithm are shown in

Table 2. List of simulated annealing and error back-propagation neural network (SA-BP) algorithm parameter settings.

Parameter Name	Parameter Value
annealing temperature $T_{\max }$	100
final annealing temperature $T_{\min }$	0.01
initial solution $S$	${0.1}^{1\times \mathrm{n}}$ (n is weight dimension)
temperature decay parameter $d$	0.95
number of iterations per temperature $L$	100
number of neural network layers	4
input layer node number	7
output layer node number	1
number of hidden layer nodes	2

.

Field LAI and the corresponding Landsat land surface reflectance data (not filled data) were combined as the training dataset. For each study area, 75% of the field LAI data was selected for model training, and the rest were used for model validation. The results of each site were validated, the data of the three sites were compiled, and the overall accuracy was calculated. Figure 3 shows the validation accuracy of 30 m LAI estimation for the three study areas. The SA-BP estimation overall accuracy was reliable with a coefficient of determination (${\mathrm{R}}^2$) of 0.90 and root mean squared error (RMSE) of 0.41.

LAI estimated with SA-BP was then used as the 30 m LAI observation data in the assimilation scheme. Different SA-BP models were selected based on the land cover map (Figure 1) at each site to obtain a 30 m LAI estimation.

3.2. LAI Dynamic Model Construction

The LAI dynamic model was built based on the LAI background. The temporal dynamic model was constructed as follows,

$$LA{I}_k=S_k\times LA{I}_{k^{\prime }}$$

(3)

where $LA{I}_k$ is the LAI at the current moment, $LA{I}_{k^{\prime }}$ is the LAI at the adjacent moment, and $S_k$ can be expressed as

$$S_k=1+\frac{1}{LA{I}_k+\epsilon }\times \frac{dLA{I}_k}{dk}$$

(4)

where $S_k$ is the linear state-transition operator, $\epsilon$ is set to 10⁻⁴ to prevent the denominator in the formula from being 0, and $\frac{dLA{I}_k}{dk}$ is the change rate at moment.

The LAI dynamic model provided the LAI primary estimation in the data assimilation scheme. For the EnKF assimilation algorithm, $k^{\prime }=k-1$. For MEnKF, in the forward estimation process, $k^{\prime }=k-1$, and in the backward estimation process, $k^{\prime }=k+1$.

3.3. MEnKF Assimilation Algorithm

The EnKF uses statistical methods to calculate the variance in the Kalman filter, and it uses an ensemble of model states to depict the error statistics. The dynamic models and observed parameters are expressed as

$$d_{k+1}=M_k{d}_k+{\omega }_k$$

(5)

$$g_k=S_k{d}_k+{\upsilon }_k$$

(6)

where $d_k$ and $g_k$ represent the true state vector and observed vector, respectively. In this study, $d_k$ is the dynamic model state vector; $M_k\in E^{n\times n}$ is the forecast pattern matrix, which will change at different k; $S_k$ is the observation operator that maps the real state space into the observation space; and $v_k$ is the observation noise. The mean is zero and the covariance matrix is $R_k$, and it follows a normal distribution.

During the assimilation process, this study used the observation information to linearly transform the background field and correct it, so the analysis variable at time $k$ is $d_k^a$:

$$d_k^a=d_k^f+I_k(g_k-S_k{d}_k^f)$$

(7)

$$I_k=P_k^f{S}_k^T{(S_k{P}_k^f{S}_k^T+R_k)}^{-1}$$

(8)

where $d_k^f$ is the state variable, $d_k^a$ is the updated state variable, $g_k-S_k{d}_k^f$ is the observed error, equal to $v_k$ in Equation (6), and $I_k$ are undetermined coefficients, also called the Kalman gain matrix. Combining Equations (7) and (8), the predicted covariance, which is the posterior covariance of the state variables, is

$$P_k^a=P_k^f-I_k{S}_k{P}_k^f=P_k^f-P_k^f{S}_k^T{(S_k{P}_k^f{S}_k^T+R_k)}^{-1}{S}_k{P}_k^f$$

(9)

where $P_k^f$ is the predicted covariance, $P_k^a$ is the posterior error covariance distribution, and $R_k$ is the covariance of $v_k$ in Equation (6).

In the ensemble Kalman filter algorithm, a matrix containing all set members $d_i$ can be defined as

$$A=(d_1,d_2,\dots d_N)\in {\Re }^{n\times N}$$

(10)

The ensemble mean is stored in each column of $\overline{A}$. The ensemble perturbation matrix is defined as

$$A^{\prime }=A-\overline{A}$$

(11)

The ensemble covariance matrix $B_{\mathrm{c}}$ is defined as

$$B_{\mathrm{c}}=\frac{A^{\prime }{(A^{\prime })}^T}{N-1}$$

(12)

Given a vector of measurements $e\in {\Re }^m$, where $m$ is the observation, N perturbation observation vectors can be defined as

$$e_j=e+{\epsilon }_j j=1,\dots ,N$$

(13)

They are stored in order in the columns of matrix $E$:

$$E=(e_1,e_2,\dots e_N)\in {\Re }^{m\times N}$$

(14)

The standard ensemble Kalman filter analysis equation is expressed as the ensemble covariance matrix:

$$A^a=A+A^{\prime }{A^{\prime }}^T{H}^T{(H{A}^{\prime }{A^{\prime }}^{T}H+O)}^{-1}(E-HA)$$

(15)

where $H\in {\Re }^{m\times N}$ is the nonlinear observation operator measurement ensemble, $O\in {\Re }^{m\times m}$ is the observation error covariance matrix ensemble, and $E\in {\Re }^{m\times N}$ is the disturbance observation matrix ensemble. When the model state vector is increased, ${\overline{d}}^T=[d^T,h^T(d)]$. $\overline{Y}\in {\Re }^{n\times N}$ is defined as the set matrix of augmented state vectors, and $\overline{Y^{\prime }}\in {\Re }^{n\times N}$ is the set disturbance matrix of the state variables. As such, the analytical equation is

$$A^a=A+A^{\prime }{\overline{Y^{\prime }}}^T{\overline{H^{\prime }}}^T{(\overline{H}\overline{Y^{\prime }}{\overline{Y^{\prime }}}^T{\overline{H^{\prime }}}^T+O)}^{-1}(E-\overline{H}\overline{Y})$$

(16)

where $\overline{H}$ is a new observation operator. The analytical equation can solve the data assimilation problem for which the observation operator is nonlinear. In this study, the EnKF assimilation scheme was used for the 30 m LAI estimation. The observation operator was the SA-BP model, and the state variable was LAI. The model error ${\omega }_k$ was defined as the difference in LAI estimated from the dynamic model and the field LAI, and the observation error ${\upsilon }_k$ was defined as the difference in LAI estimated from SA-BP and the field LAI. This study generated N random noises from a normal distribution with an average of 0 and a variance of 0.01. The number of ensembles, N, was positively related to the calculation accuracy and negatively related to the calculation efficiency of the algorithm. In this study, N was set to 100.

In the common sequential EnKF assimilation process, when there are no observation data, the background field is used as the final estimation. As time elapses, the error will continuously accumulate. If the observation data come late, a large error will accumulate, resulting in low estimation accuracy. Although MODIS LAI data under quality control were filtered to build the LAI dynamic model, a large error exists when the data are directly used as the final LAI estimation. Therefore, this study proposed MEnKF by adjusting the starting point of assimilation and set the observation point that is closest to the annual maximum of the background field as the starting point of assimilation. The MEnKF works as follows. The multiyear average of MODIS LAI is calculated as the LAI background, and a temporal dynamic model is established based on the LAI background to predict the time series LAI. Landsat LAI is used as the observation value introduced into the data assimilation scheme. The Landsat LAI, which is close to the annual maximum of the background field, is introduced first, the assimilation is performed backward and forward, and a whole year 30 m LAI is then estimated.

4. Results

The LAI time series was estimated at farmland, grassland, and woodland sites. Both pixel- and regional-scale LAI time series are shown here. The estimated result was compared with the field measurements and the Landsat LAI map; ${\mathrm{R}}^2$ and RMSE were employed to determine the estimation accuracy.

4.1. Single-Point Time Series LAI Validation

Figure 4 shows the LAI estimation results of the farmland pixels. At the farmland site, the MODIS LAI fluctuates substantially, especially during the peak vegetation growing period. The LAI background filtered from the historical MODIS product is lower than the LAI observation estimated from Landsat data. The EnKF-estimated LAI starts from the first date of the year and then continues in sequence. For this site, Landsat observations are temporally close to the vegetation growing period, which makes a long period of no observation from the beginning of the year. In this period, the LAI initial estimation from the temporal dynamic model cannot be updated, and the final estimation of EnKF is the same as the LAI background. MEnKF starts estimation from the point closest to the annual maximum of the background; in these sites, the date of Julian day 189 was set as the starting point of MEnKF. Through the comparison of the LAI values estimated by EnKF ($LA{I}_{EnKF}$) and MEnKF ($LA{I}_{MEnKF}$), it can be seen that with the LAI observation introduced at the first time, $LA{I}_{MEnKF}$ is much closer than $LA{I}_{EnKF}$ to the field observation.

Figure 5 shows the LAI estimation at the grassland site. From Figure 5 we can see that the LAI of grassland is much smaller than that of cropland. MODIS LAI fluctuates during the growth period. The differences between the LAI background and LAI observation values are minor, which also leads to slight differences in $LA{I}_{EnKF}$ and $LA{I}_{MEnKF}$. At this site, the Landsat LAI of Julian day 229 was set as the starting point of MEnKF. $LA{I}_{EnKF}$ and $LA{I}_{MEnKF}$ before this date have slight differences, but the estimations after it are exactly the same, which is induced by the same estimation algorithm for the forward estimation. It can be seen from Figure 5 that MEnKF has a slight advantage when a small difference between the LAI background and LAI observation values exists.

Figure 6 shows the LAI estimation at the woodland site. At woodland sites, the MODIS LAI also fluctuates substantially, especially during the peak growing peak, and the MODIS LAI has many dramatic drop points. The LAI background-filtered historical MODIS LAI product is more continuous, but there is still a certain degree of overestimation. $LA{I}_{EnKF}$ continuously accumulates errors due to the overestimation of the background before observations are added. When the first Landsat LAI observation was added to the data assimilation scheme on Julian day 213, $LA{I}_{EnKF}$ was noticeably adjusted but still larger than $LA{I}_{MEnKF}$. Furthermore, for each point where LAI observation was introduced, a sudden change occurred in the $LA{I}_{EnKF}$ profile, especially during the growing period. The MEnKF effectively solved this problem by assimilation from the peak value. The overall trend is smooth, and the difference between $LA{I}_{MEnKF}$ and field observations is smaller than that of $LA{I}_{EnKF}$.

Then, data from the three research areas were compiled, and the estimated LAI was compared with field observations. Figure 7 shows the scatterplot of field LAI with $LA{I}_{MEnKF}$ and $LA{I}_{EnKF}$. $LA{I}_{MEnKF}$ has a good correlation with a field LAI with ${\mathrm{R}}^2$ of 0.70 and RMSE of 0.40. $LA{I}_{EnKF}$ also has a good correlation with field observations, but the accuracy is lower; the ${\mathrm{R}}^2$ is 0.40 and the RMSE is 0.78. For each land cover type, the MEnKF performed better than the EnKF, with a higher ${\mathrm{R}}^2$ and a lower RMSE. The MEnKF noticeably improved the estimation for farmland and woodland, where the LAI value is large and the LAI background differs substantially from the Landsat LAI.

4.2. Regional LAI Validation

The regional LAI time series was generated in each study area. The LAI estimation results of farmland, grassland, and woodland are shown in Figure 8. For each site, as the field measurements are only from farmland, grassland, and woodland, LAIs of these three land cover types are estimated and the other area is masked. It can be seen from Figure 8 that for each site, the MEnKF was capable of reliably estimating LAI. The variations in LAI characteristics are consistent with those of the vegetation growth pattern.

In the Pshenichne farmland area, $LA{I}_{MEnKF}$ captures the LAI change trend of vegetation growth in farmland. From the beginning of the year to the end, the LAI gradually increases and then decreases. From the regional LAI time series, the vegetation in this area was divided into two parts. One part is an early crop distributed in most of the research area. Crops of this kind begin to grow on the 73rd day and reach a peak in the middle of the year, and the growth period lasts until the 273rd day. The other part consists of late crops concentrated in the central and southeastern areas of the research area. Crops of this kind began to grow from the 137th day to the 297th day, and their growth cycle was shorter than that of the first part. At the Zhangbei grassland site, the LAI of both cropland and grassland was estimated, and the LAI of cropland was larger than that of grassland, which is consistent with field observations. The LAI time series shows that the LAI of this area has a higher overall LAI value between days 177 and 241 (approximately June to August). The grasslands in this period are relatively dense, which is consistent with the actual situation in the area.

At the Genhe woodland site, the vegetation has a short growth cycle due to high latitude. Leaves sprout in May and fall in September, with approximately one-hundred days of obvious changes. In this area, two vegetation types exist, i.e., forest and grass, and Figure 8c shows that the LAI of forest is larger than that of grass. The LAI distribution map reflects the vegetation growth characteristics of the area well.

$LA{I}_{MEnKF}$ was compared with Landsat LAI at the regional scale, and Figure 9 shows the scatterplot of LAI in the three research areas. Because the return cycle of Landsat is not consistent with that of MODIS, $LA{I}_{MEnKF}$ and Landsat LAI were selected with the closest dates for validation. For each site, Landsat data used for validation have not been used in the data assimilation scheme. Figure 9 shows that for all of the sites, $LA{I}_{MEnKF}$ and Landsat LAI are consistent, and ${\mathrm{R}}^2$ is larger than 0.65 for all sites. For farmland and grassland, the vegetation is in the growing season, and the LAI value is high, so the ${\mathrm{R}}^2$ value is also high. For these two sites, ${\mathrm{R}}^2$ is larger than 0.88. For the woodland site, ${\mathrm{R}}^2$ is 0.66, which is lower than that of the other two sites. Obvious underestimation was found at this site. During the research period, the vegetation grows rapidly; the date mismatch of Landsat data and $LA{I}_{MEnKF}$ date is the main reason for the underestimation. Overall, $LA{I}_{MEnKF}$ is consistent with Landsat LAI at the regional scale, and the proposed method is capable of deriving a spatio-temporal complete 30 m LAI map.

5. Discussion

5.1. Error Induced by Landsat LAI

In this study, Landsat LAI was estimated by the SA-BP neural network trained from field LAI data and Landsat reflectance data. The field measurement was considered to match the 30 m Landsat pixel scale. However, errors in field measurement, gap filling, and SA-BP estimation will propagate into the final LAI estimation in data assimilation [45,51]. As shown in Figure 3, the validation accuracy of the estimated Landsat LAI was reliable, with an ${\mathrm{R}}^2$ of 0.90 and an RMSE of 0.41. Compared with the RT model or VI-LAI empirical model [17,29,52], this method has the advantages of reliability, rapidity, and simplicity. However, field measurements were only collected at the three research sites, and the data were limited. To extend the method to more areas and more vegetation types, more field data are still needed to improve the 30 m LAI estimation.

5.2. Error Induced by the LAI Background

The LAI background provides the LAI initial estimation. The background is derived from historical material or data, which reflects the variation trend of the estimated value in the period [35,53]. This study used multiyear average MODIS LAI data after spatial registration with Landsat data as the background. Due to the differences in resolution and sensors, there are usually small differences between MODIS LAI and Landsat LAI. Therefore, in the data assimilation process, the weight of the background is negatively correlated with the weight of the observation. When the quality of the background is poor, its weight is low and the weight of observations increases, and vice versa [54]. Therefore, although errors exist in the LAI background, the data assimilation algorithm is capable of improving the final estimation with 30 m LAI observations induced.

5.3. Advance of Assimilation from the LAI Peak Position

Traditional sequential assimilation starts on the first day of the year; when no observation is introduced, the LAI initial estimation is considered as the final estimation, and as time elapses, the estimation error will cumulate and become significant [55]. This study modified the assimilation sequence by starting from the position where observation exists near the LAI peak value position. The timely introduction of field observations efficiently prevents error accumulation, which is helpful for accurate LAI estimation.

5.4. Extendibility of MEnKF

In this study, MEnKF was used for time-series 30 m LAI estimation, and the estimation accuracy satisfied the LAI application demand. The capability of the MEnKF in spatio-temporal continuous LAI estimation was shown. Based on the algorithm mechanism itself, the MEnKF can be extended to other land surface types and other land surface parameter estimations once the dynamic model is reasonable and a 30 m observation value is available. The robustness of the data assimilation algorithm guarantees the extendibility of the MEnKF.

6. Conclusions

Spatio-temporal complete LAI data are critical for regional to global scale vegetation monitoring and change detection. This study proposed a modified data assimilation method to obtain spatio-temporal continuous 30 m LAI estimation.

The MODIS LAI product after the S-G filter was used as the LAI background. A temporal LAI dynamic model was constructed based on the filtered LAI profile. The 30 m LAI was then retrieved with SA-BP based on the filled reflectance data and then used as the 30 m LAI observation data. The MEnKF data assimilation algorithm was run to obtain a 30 m temporal continuous estimation. To avoid error accumulation, data assimilation began from the peak value observation and was performed backward (from the end to the beginning of the growing season) and forward (from the beginning to the end of the growing season). The verification results showed that the MEnKF LAI was consistent with field measurements for all sites (${\mathrm{R}}^2$ = 0.70, RMSE = 0.40) and that the regional MEnKF LAI correctly depicted the growth trend of vegetation, which was consistent with the Landsat LAI (${\mathrm{R}}^2$ > 0.65; RMSE < 0.51 in all sites).

There are some advantages of the proposed MEnKF: (1) It can generate a reliable complete 30 m time series (every 8 days) LAI. (2) It greatly reduces the error accumulation in the data assimilation scheme. (3) It is robust and extendible to different land surface types and different land surface parameter estimations.

There are also some drawbacks and further work to perform: (1) The dynamic model is a temporal dynamic model simply deduced from the filtered MODIS LAI profile. Improvement in the dynamic model is needed to improve the estimation accuracy. (2) Landsat LAI greatly affects the final estimation, and accurate Landsat LAI is the premise of accurate estimation. In this study, Landsat LAI was estimated by the SA-BP model combined with field data. A more robust 30 m LAI estimation algorithm is urgently needed.