<figure><img src="image_1.png" /></figure>

HAL
open science

# Multiscale estimation of leaf area index from satellite observations based on an ensemble multiscale filter

Jingyi Jiang, Zhiqiang Xiao, Jindi Wang, Jinling Song

► To cite this version:

Jingyi Jiang, Zhiqiang Xiao, Jindi Wang, Jinling Song. Multiscale estimation of leaf area index from satellite observations based on an ensemble multiscale filter. Remote Sensing, 2016, 8 (3), pp.229-250. 10.3390/rs8030229 . hal-01315236

HAL Id: hal-01315236

https://hal.science/hal-01315236

Submitted on 28 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Multiscale Estimation of Leaf Area Index from Satellite Observations Based on an Ensemble Multiscale Filter

Jingyi Jiang<sup>1,2</sup>, Zhiqiang Xiao<sup>1,*</sup>, Jindi Wang<sup>1</sup> and Jinling Song<sup>1</sup>

Academic Editors: Clement Atzberger and Prasad S. Thenkabail

Received: 26 September 2015; Accepted: 29 February 2016; Published: 12 March 2016

**Abstract:** Currently, multiple leaf area index (LAI) products retrieved from remote sensing data are widely used in crop growth monitoring, land-surface process simulation and studies of climate change. However, most LAI products are only retrieved from individual satellite observations, which may result in spatial-temporal discontinuities and low accuracy in these products. In this paper, a new method was developed to simultaneously retrieve multiscale LAI data from satellite observations with different spatial resolutions based on an ensemble multiscale filter (EnMsF). The LAI average values corresponding to the date of satellite observations were calculated from the multi-year Moderate Resolution Imaging Spectroradiometer (MODIS) LAI product and were used as a priori knowledge for LAI in order to construct an initial ensemble multiscale tree (EnMsT). Satellite observations obtained at different spatial resolutions were then applied to update the LAI values at each node of the EnMsT using a two-sweep filtering procedure. Next, the retrieved LAI values at the finest scale were used as a priori knowledge for LAI for the new round of construction and updating of the EnMsT, until the sum of the difference of LAI values at each node of the EnMsT between two adjacent updates is less than a given threshold. The method was tested using Thematic Mapper (TM) or Enhanced Thematic Mapper Plus (ETM+) surface reflectance data and MODIS surface reflectance data from five sites that have different vegetation types. The results demonstrate that the retrieved LAI values for each spatial resolution were in good agreement with the aggregated LAI reference map values for the corresponding spatial resolution. The retrieved LAI values at the coarsest scale provided better accuracy with the aggregated LAI reference map values (root mean square error (RMSE) = 0.45) compared with that obtained from the MODIS LAI values (RMSE = 1.30).

**Keywords:** leaf area index; MODIS; TM; ensemble multiscale filter; retrieval

## 1. Introduction

Leaf area index (LAI) is defined as half of the total area of green leaves for a given unit of horizontal ground surface area [1]. LAI is an essential land surface biophysical parameter that describes the amount and change of vegetation, and is a significant input for climate and ecological simulations, including crop growth models [2], hydrological models [3], ecology models [4], weather forecasting [5] and studies of climate change [6].

Many methods have been developed to retrieve LAI from satellite remote-sensing data. In general, the methods used to estimate LAI from remote sensing data can be distinguished into two types: empirical methods and physical methods. Empirical methods, which utilize statistical relationships between LAI and a variety of vegetation indices to calculate LAI values, are easy to operate and

<footnote n="1">¹ State Key Laboratory of Remote Sensing Science, School of Geography, Beijing Normal University, Beijing 100875, China; jiangjingyi@mail.bnu.edu.cn (J.J.); wangjd@bnu.edu.cn (J.W.); songjl@bnu.edu.cn (J.S.)</footnote>

<footnote n="2">² INRA-EMMAH UMR 1114, Avignon 84914, France</footnote>

<footnote>* Correspondence: zhqxiao@bnu.edu.cn; Tel.: +86-10-5880-7698; Fax: +86-10-5880-5274</footnote>
are computational efficient [7,8]. However, the main limitation of these empirical methods is their
dependence on measurement conditions, sites and vegetation types. Physical methods are based on the
inversion of canopy radiative-transfer models to simulate physical and biological process and provide
an explicit connection between canopy biophysical variables and the resulting canopy bidirectional
reflectances. Since physical methods can be adjusted for a wide range of situations, an increasing
number of studies tend to use canopy radiative-transfer models in the inversion of LAI values [9,10].

Currently, multiple global LAI products were generated from moderate resolution remote sensing data such as the Moderate-Resolution Imaging Spectroradiometer (MODIS) [11] and VEGETATION [12]. Several regional LAI products were also generated from high resolution satellite data (i.e., the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [13] and the Thematic Mapper (TM) [14]). These LAI products are routinely generated from remote-sensing data acquired only at a specific time. Because limited information is used for the inversion process, these products are not spatially and temporally continuous [15,16]. To obtain better retrieval results, some assimilation methods were developed to estimate LAI values from time-series satellite observations [17]. Xiao *et al.* [17] developed a temporally integrated inversion method by coupling a radiative-transfer model and a double-logistic model to estimate LAI values from time-series MODIS reflectance data. Jiang *et al.* [18] implemented a sequential assimilation method with an incremental analysis update to provide a smooth and continuous LAI time series from MODIS reflectance data. Lewis *et al.* [19] developed an Earth Observation Land Data Assimilation System (EO-LDAS) to improve the inversion of biophysical parameters. Viskari *et al.* [20] used ensemble adjustment Kalman filter (EAKF) to update predictions of LAI and leaf extension, which successfully simulated the evolution of plant leaf phenology and improved model predictions of forest net ecosystem exchange. All of these methods show the ability to retrieve temporally continuous LAI from reflectance data in time series. However, these methods only use individual satellite observations during the retrieval process.

To make full use of remote sensing information during the retrieval process, many data fusion methods, such as optimal interpolation (OI) [21], Bayesian maximum entropy (BME) [22], Bayesian hierarchical models [23], multi-resolution tree (MRT) [24] and an empirical orthogonal function (EOF) [25], have been developed to improve the accuracy and integrity of land-surface parameter products by integrating multi-source data. Among these data fusion methods, the MRT method has the advantage of integrating multi-resolution data with high computational efficiency and can generate products with different resolutions [26]. It has been used extensively to integrate soil moisture data [27], precipitation data [28] and land-surface albedo [29] from different satellite observations. To obtain spatially complete LAI products with different resolutions, Xiao *et al.* [24] integrated MODIS and ETM+ LAI via the MRT method. Wang *et al.* [30] also utilized the MRT method to fuse MODIS and Multi-angle Imaging SpectroRadiometer (MISR) LAI products with different resolutions.

To further improve computational efficiency and apply a nonlinear dynamic model via the MRT method, Zhou *et al.* [31] proposed an ensemble multiscale filter (EnMsF) method which combined the advantages of both the Ensemble Kalman Filter (EnKF) technique and the MRT method by adding state ensembles to the basic MRT method [32,33]. The EnMsF method has been applied in different fields, and some improvements have been made. Zhou [34] utilized the EnMsF method to estimate soil moisture and surface fluxes from multi-resolution L-band passive and active microwave soil moisture measurements following HYDROS specifications. Lawniczak *et al.* [35] applied the EnMsF method in reservoir engineering through assimilated permeability values and observations to obtain results at different scales. Pan *et al.* [36,37] developed an automated procedure for dividing irregular shapes to create a tree of “balanced” topology, and introduced the EnMsF method to hydrological land surface-driven applications. Some analysis on impacts of accuracy, spatial availability and assimilation frequency on the EnMsF method have also been discussed [38].

However, the observation operators used in all of the aforementioned EnMsF methods are linear.
Therefore, these methods cannot be used to directly retrieve LAI values from satellite observations
because the radiative transfer model, which is a non-linear relationship between LAI and surface
reflectance data, is used as an observation operator. This paper aims to introduce the radiative transfer model into the EnMsF method to simultaneously retrieve multiscale LAI from satellite observations with different spatial resolutions.

An initial ensemble multiscale tree (EnMsT) is constructed with the LAI average values corresponding to the date of satellite observations calculated from the multi-year MODIS LAI product, and TM/ETM+ and MODIS reflectance data were used to update the LAI values at each node of the EnMsT by using a two-sweep filtering procedure. Then, the retrieved LAI values at the finest scale were used as a *priori* knowledge for LAI for the new round of construction and updating of the EnMsT, until the sum of the difference of LAI values at each node of the EnMsT between two adjacent updates is less than a given threshold. The retrieved LAI values with different resolutions were compared with LAI reference maps calculated from ground measurements at several sites with different vegetation types.

The organization of this paper is as follows. Section 2 provides the details of the method. The radiative-transfer model is first introduced, followed by the construction of the EnMsT. Then the two-sweep filtering procedure is described to retrieve LAI values from multi-resolution reflectance data. The satellite data products used in this study, and the LAI reference maps at the selected sites are also described in this section. Section 3 shows the performance of the new method and we provide an analysis of the retrieved results. Discussions are presented in Section 4, and the final section provides brief conclusions.

## 2. Methodology and Data

The general process of the new method to retrieve multiscale LAI data from satellite observations with different spatial resolutions is shown in Figure 1. The multi-year average of MODIS LAI corresponding to the date of satellite observations is used as a *priori* knowledge for LAI at the finest scale to construct an initial EnMsT. A two-sweep filtering procedure is used to update node states on the tree by integrating multiple satellite observations with different spatial resolutions. For the upward sweep, an EnKF was used to update the node states from satellite observations from the finest scale to the coarsest scale by coupling the EnMsT with a canopy radiative transfer model. The downward sweep starts from the coarsest scale and ends at the finest scale to smooth the node states on the tree according to all of the satellite observations. Then, the updated node states at the finest scale are used as the *priori* knowledge for LAI for the new round of construction and updating of the EnMsT. The satellite observations are reused to update nodes states on the EnMsT by using the two-sweep filtering procedure. The circulation process will continue until the sum of the difference of LAI values at each node of the tree between two adjacent updates is less than a threshold, which is set to 0.01 in this study. Then, the ensemble mean of LAI at each node of the EnMsT for the last round was recorded as the retrieved LAI values at different resolutions.

### 2.1. Radiative-Transfer Model

Canopy radiative-transfer models describe the relationship between canopy reflectance and canopy characteristics. The two-layer canopy reflectance model (ACRM) developed by Kuusk [39,40] can be used for canopy reflectance modeling at spectral wavelengths ranging from 400 to 2400 nm. Many studies have used ACRM to retrieve biophysical parameters from satellite remote sensing data [41,42].

In this study, the latest version of ACRM was used to calculate the directional reflectance factor of vegetation canopies. The LAI of the lower layer was set to zero, or in other words, there was only one layer in the canopy. The parameters needed to run the ACRM are listed in Table 1. Based on a global sensitivity analysis of the ACRM using an extension of the Fourier amplitude sensitivity testing (EFAST) [43,44], four parameters (*lai*, *slw*, *c2*, and *s1*) to which remote sensing signals are most sensitive were chosen as free variables in this study and are marked by asterisks in Table 1, and other nonsensitive parameters were set to empirical values.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Flow chart of leaf area index (LAI) inversion from satellite observations with different spatial resolutions.</figcaption></figure>

<table><caption>**Table 1.** Parameters needed to run the ACRM.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td>Leaf area index *</td><td>lai</td><td>0–10.0</td><td>m<sup>2</sup>/m<sup>2</sup></td></tr><tr><td>The hot-spot parameter</td><td>sl</td><td>0.1</td><td>–</td></tr><tr><td>Foliage clumping parameter</td><td>clmp</td><td>0–2.0</td><td>–</td></tr><tr><td>The displacement parameter</td><td>szz</td><td>1.2</td><td>–</td></tr><tr><td>The elliptical leaf angle distribution</td><td>eln</td><td>2.0</td><td>–</td></tr><tr><td>Mean leaf angle of the LAD</td><td>thm</td><td>45.0</td><td>degree</td></tr><tr><td>Factor for refraction index</td><td>nratio</td><td>0.9</td><td>–</td></tr><tr><td>Leaf specific weight *</td><td>slw</td><td>80.0–180.0</td><td>g/m<sup>2</sup></td></tr><tr><td>Leaf water content</td><td>c1</td><td>137.0</td><td>1/cm</td></tr><tr><td>Chlorophyll AB content *</td><td>c2</td><td>0.2–0.7</td><td>cm<sup>2</sup>/g</td></tr><tr><td>Leaf dry matter content</td><td>c3</td><td>98.0</td><td>cm<sup>2</sup>/g</td></tr><tr><td>Anthocyanin</td><td>c4</td><td>0.038</td><td>cm<sup>2</sup>/g</td></tr><tr><td>Leaf structure parameter</td><td>N</td><td>1.624</td><td>–</td></tr><tr><td>Weight of the first Price function *</td><td>s1</td><td>0.05–0.4</td><td>–</td></tr><tr><td>Weight of the Second Price function</td><td>s2</td><td>0.0</td><td>–</td></tr><tr><td>Weight of the third Price function</td><td>s3</td><td>0.0</td><td>–</td></tr><tr><td>Weight of the fourth Price function</td><td>s4</td><td>0.0</td><td>–</td></tr></tbody></table>

## 2.2. Construction of the Ensemble Multiscale Tree

The EnMsT was first proposed by Zhou *et al.* [31], which is based on the original multiscale tree [32]. During the process of constructing the tree, the global correlations are built up from local correlations between nearby tree nodes [31]. As shown in Figure 2, an ensemble multiscale tree is made up of a set of nodes. Let $s$ be a node on the tree, $\chi(s)$ be a $n(s)$ dimension state vector at the node. We define a matrix $x(s) \in \mathbb{R}^{n(s) \times N}$ that contains $N$ ensemble members of the state vector, namely $x(s) = (\chi^1(s), \chi^2(s), \dots, \chi^N(s))$. Except for nodes at the finest scale, all nodes on the tree have $q$ children $s\alpha_1, s\alpha_2, \dots, s\alpha_q$. Moreover, except for nodes at the coarsest scale, all nodes have a single parent $s\gamma$. Let $m(s)$ be the scale where node $s$ is located. $m(s)$ increases from 0, which represents the top of the tree (i.e., the coarsest scale), to $M$, which represents the bottom of the tree (i.e., the

<footnote>* Free parameters.</footnote>
finest scale). Nodes at the finest scale represent the highest resolution, while nodes at the coarsest scale represent the lowest resolution. From bottom to top, both the number of nodes and the spatial resolution at each scale are reduced proportionally. Inside the EnMsT, nodes at the same scale are not directly linked. Instead, each node is linked with its parent at the coarser scale and children at the finer scale. In Figure 2, this parent-child relationship among nodes is described using black solid lines.

<figure><img src="image_5.png" /><figcaption>Figure 2. The structure of the ensemble multiscale tree.</figcaption></figure>

According to the structure of the ensemble multiscale tree, different nodes at the finest scale are related through their common parent nodes [32]. Based on the multiscale autoregressive (MAR) model, the following state transition equations can be used to compute node states recursively across the tree [32,45]. The downward recursion transiting parent state to child state is:

$$x(s) = A(s)x(s\gamma) + \omega(s) \quad (1)$$

and the upward recursion transiting child state to parent state is:

$$x(s\gamma) = F(s)x(s) + \omega'(s) \quad (2)$$

where $A(s)$ is a downward transition matrix and $\omega(s)$ is zero-mean random noise with covariance $Q(s)$, $F(s)$ is an upward transition matrix and $\omega'(s)$ is zero-mean random noise with covariance $Q'(s)$.

To get the scale transition and covariance matrices for all nodes at different scales, Frakt and Willsky [32] found that two conditions must be satisfied among nodes on the tree: one is the “locally internal” condition, and the other one is the “scale-recursive Markov property”. The “locally internal” condition means that the parent state $x(s)$ is a linear combination of all its children’s states $x_{m(s)+1}(s\alpha_i)$. Thus the ensemble version can be written as:

$$x(s) = V(s)x_{m(s)+1}(s) = V(s) \begin{bmatrix} x(s\alpha_1) \\ \vdots \\ x(s\alpha_q) \end{bmatrix}, m(s) < M \quad (3)$$

where $x(s\alpha_1), x(s\alpha_2), \dots, x(s\alpha_q)$ are $q$ state vectors at children nodes $s\alpha_i$, and $V(s)$ is a coefficient matrix with a dimension of $n(s) \times (qn(s))$. The “scale-recursive Markov property” means that for any node located at scale $m(s)$ ($m(s) < M$), its state $x(s)$ is conditionally uncorrelated with these $q+1$ sets. Based on these two properties, the maximum predictive efficiency method, which is described in detail in Willsky [45] and Zhou [34], is used to calculate the coefficient matrix $V(s)$. Then it is easy to determine the scale transition and covariance matrices for all nodes at different scales:

$$A(s) = \operatorname{cov}[x(s), x(s\gamma)] \operatorname{cov}^{-1}[x(s\gamma)] \quad (4)$$

$$F(s) = \operatorname{cov}[x(s\gamma)] A(s)^T \operatorname{cov}^{-1}[x(s)] \quad (5)$$
$$Q(s) = \operatorname{cov}[\mathbf{x}(s)] - \mathbf{A}(s) \operatorname{cov}[\mathbf{x}(s), \mathbf{x}(s\gamma)]^T \quad (6)$$

$$Q'(s) = \operatorname{cov}[\mathbf{x}(s\gamma)] - \mathbf{F}(s) \mathbf{A}(s) \operatorname{cov}[\mathbf{x}(s\gamma)] \quad (7)$$

In this study, the four free variables from the ACRM, which are denoted by asterisks in Table 1, were included in each node state $\mathbf{x}(s)$. As such, the size of the node state vector $n(s)$ is equal to four. The average and standard deviation of multi-year MODIS LAI product corresponding to the date of satellite observations was used as a priori knowledge for LAI at the finest scale to construct the initial EnMsT. The scale transition matrices and covariance matrices at each node of the EnMsT were calculated according to the formulas above.

## 2.3. Updating the Ensemble Multiscale Tree

When finishing the construction of the EnMsT, satellite observations with different spatial resolutions were used to update the LAI values at each node of the EnMsT. The updating procedure can be divided into two steps: an upward sweep and a downward sweep.

### 2.3.1. Upward Sweep

For each node of the EnMsT, the upward sweep procedure updates the state vector from satellite observations using an EnKF technique. The EnKF uses an ensemble of model states to represent the statistical errors, and then uses ensemble integrations to forecast these errors forward during the model estimation [33]. The EnKF analysis equation is as follows:

$$\chi^j(s|s) = \chi^j(s) + \mathbf{K}(s) [\mathbf{Y}^j(s) - \mathbf{Y}^j(s)] \quad (8)$$

where $\chi^j(s|s)$ is a state vector after an upward update; $\chi^j(s)$ is an initial state vector; $\mathbf{Y}^j(s)$ is a perturbed observation vector, which is composed of satellite observations of current node and the states of its children nodes; and $\mathbf{Y}^j(s)$ is a vector which contains the simulated satellite observations and the predicted states of children nodes from the current node. The detailed calculation function was written in Zhou et al. [31]. $\mathbf{K}(s)$ is a Kalman gain and can be written as:

$$\mathbf{K}(s) = \operatorname{cov}[\mathbf{x}(s), \mathbf{Y}'(s)] \{\operatorname{cov}[\mathbf{Y}'(s)] + \mathbf{R}(s)\}^{-1} \quad (9)$$

where $\operatorname{cov}[\mathbf{x}(s), \mathbf{Y}'(s)] = E\{[\mathbf{x}(s) - E(\mathbf{x}(s))][\mathbf{Y}'(s) - E(\mathbf{Y}'(s))]^T\}$ and $\operatorname{cov}[\mathbf{Y}'(s)]$ is short for $\operatorname{cov}[\mathbf{Y}'(s), \mathbf{Y}'(s)] = E\{[\mathbf{Y}'(s) - E(\mathbf{Y}'(s))][\mathbf{Y}'(s) - E(\mathbf{Y}'(s))]^T\}$; $\mathbf{R}(s)$ is an observation error covariance matrix; and the superscript “T” denotes a matrix transpose.

Equation (8) imposes the condition that the observation operator $\mathbf{h}(\cdot)$ is linear, $\mathbf{Y}^j(s) = \mathbf{h}(s)\chi^j(s)$. Therefore, Equation (8) becomes invalid when satellite observations are used to update the node states on the EnMsT because the canopy radiative transfer model represents a non-linear function of node states. To overcome this limitation, the state vector $\chi^j(s)$ can be augmented with a diagnostic variable from observations simulated by the canopy radiative transfer model [46]. Therefore, a new state vector is defined as $\hat{\chi}^j(s) = [\chi^j(s)^T, \mathbf{h}^T(s)\chi^j(s)] \in \mathbb{R}^{n(s)}$ with $\hat{\mathbf{h}}(s)$ being the sum of the dimension of the node state vector $n(s)$ and the number of measurement equivalents. Thus the new update equation is:

$$\chi^j(s|s) = \chi^j(s) + \hat{\mathbf{K}}(s) [\mathbf{Y}'(s) - \hat{\mathbf{H}}\hat{\chi}^j(s)] \quad (10)$$

where $\hat{\mathbf{H}}$ is a linear observation operator that relates the augmented state vector to satellite observations, and $\hat{\mathbf{K}}(s)$ is the new Kalman gain:

$$\hat{\mathbf{K}}(s) = \operatorname{cov}[\mathbf{x}(s), \hat{\mathbf{Y}}'(s)] \{\operatorname{cov}[\hat{\mathbf{Y}}'(s)] + \mathbf{R}(s)\}^{-1} \quad (11)$$
During the upward sweep procedure, the node states on the tree are updated using the EnKF technique from the finest scale to the coarsest scale. For each node at the finest scale, only satellite observations were used to update the state at this node. If no satellite observation exists, the update just comes from the state of the node itself. However, for each node located above the finest scale, the observations used to update the state at this node include data obtained by satellites as well as the states of its children nodes. If there is no satellite observation, the node is updated according to the states of its children nodes.

### 2.3.2. Downward Sweep

After the upward sweep, the downward sweep is conducted from the top to the bottom of the tree. During this process, nodes are smoothed from the top down depending on all of the satellite observations on the EnMsT.

If node s is located at the top of the tree, the state vector after a downward sweep $\chi^j(s|S)$ is kept the same as the state vector after the upward sweep [31]:

$$ \chi^j(s|S) = \chi^j(s|s), m(s) = 0 \qquad (12) $$

If node s is below the root node, the state vector after a downward sweep is:

$$ \chi^j(s|S) = \chi^j(s|s) + J(s) [\chi^j(s\gamma|S) - \chi^j(s\gamma|s)], m(s) > 0 \qquad (13) $$

where $\chi^j(s\gamma|s)$ and $\chi^j(s\gamma|S)$ are the state vector of the parent node $s\gamma$ before and after the downward sweep, respectively, and $J(s)$ is the smoothing gain. $\chi^j(s\gamma|s)$ can be got from the MAR upward recursion in Equation (2) [31]:

$$ \chi^j(s\gamma|s) = F(s) \chi^j(s|s) + \omega^{j'}(s) \qquad (14) $$

and $J(s)$ is calculated as follow:

$$ J(s) = \operatorname{cov}[x(s|s)] F^T(s) \operatorname{cov}^{-1}[x(s\gamma|s)] \qquad (15) $$

After the two-sweep filtering procedure, the updated states of each node at different scales are recorded separately.

### 2.4. Data and Preprocessing

TM/ETM+ and MODIS surface reflectance data were used to test the new method at five sites from Bigfoot and Valeri: Bondville, Sud-Ouest, Laprida, Zhangbei, and Puechabon. Characteristics of these sites are given in Table 2.

<table><caption>**Table 2.** Characteristics of the selected sites.</caption><thead><tr><th>Site Name</th><th>Latitude (°)</th><th>Longitude (°)</th><th>Land Cover Type</th><th>Date of LAI Reference Maps (Day/Year)</th></tr></thead><tbody><tr><td>Bondville</td><td>40.004</td><td>-88.29</td><td>Cropland</td><td>224/2000</td></tr><tr><td>Sud-Ouest</td><td>43.506</td><td>1.2375</td><td>Cropland</td><td>201/2002</td></tr><tr><td>Laprida</td><td>-36.990</td><td>-60.552</td><td>Grassland</td><td>317/2002</td></tr><tr><td>Zhangbei</td><td>41.278</td><td>114.687</td><td>Grassland</td><td>235/2002</td></tr><tr><td>Puechabon</td><td>43.724</td><td>3.651</td><td>Evergreen broadleaf forest</td><td>163/2001</td></tr></tbody></table>

The MODIS surface reflectance product (MOD09A1) from the latest version (Collection 5) and the Landsat 4–5 TM or Landsat 7 ETM+ surface reflectance data were inputted into the new method to simultaneously retrieve multiscale LAI data using the EnMsF. The MOD09A1 product, in a sinusoidal projection system, provided the surface reflectance for each of the MODIS land spectral bands. The spatial resolution was 500 m, and the temporal sampling period was 8 days.
The surface reflectance data were contaminated by residual clouds or undetected cloud shadows. According to the product quality information, only reflectance data with the highest quality of all bands (QA = 0) were used in this paper. According to the typical theoretical uncertainties in the MODIS surface-reflectance product, the following values were assigned to the uncertainties in the MODIS surface reflectance: ±(0.005 + 0.05 × surface reflectance) [47]. The TM/ETM+ surface reflectance data were generated from the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [48], with 30-m spatial resolution in a Universal Transverse Mercator (UTM) projection system. According to its quality control information, only pixels without clouds, cloud shadows, snow, or water (Cfmask = 0) were used. Furthermore, 5% of surface reflectance was assigned to the uncertainties in the TM/ETM+ surface reflectance data.

Considering the spatial resolutions of the MODIS and TM/ETM+ surface reflectance data, an EnMsT in which each node has four children is used in this study. The EnMsT has six different scales. The nodes at the coarsest scale (the 0th scale) have a spatial resolution of 926.6 m, while the nodes at the finest scale (the 5th scale) have a spatial resolution of 28.9 m. The nodes at intermediate scales have spatial resolutions of 463.3 m, 231.7 m, 115.8 m and 57.9 m, respectively. Therefore, the TM/ETM+ surface reflectance data were resampled to a spatial resolution of 28.9 m, which corresponds to the finest scale, and the MODIS surface reflectance data were reprojected to the UTM projection and resampled to a spatial resolution of 463.3 m, which corresponds to the 1st scale using the nearest-neighbor resampling method. In this study, the MODIS Collection 5, 8-day, 1-km LAI product (MOD15A2) from 2000 to 2010 was used to calculate the LAI average values corresponding to the date of satellite observations which were used as a priori knowledge for LAI at the finest scale to construct the initial EnMsT. The MODIS LAI data were also reprojected to the UTM projection and resampled to a spatial resolution of 926.6 m in order that they could be compared with the retrieved LAI at the 0 scale.

In addition, high resolution LAI reference maps at the selected sites were used to evaluate the accuracy and consistency of the LAI values retrieved by the new method. At the Bondville site, LAI ground measurements were measured using destructive sampling methods. The Landsat ETM+ imagery was used to develop high resolution LAI estimates, which are directly linked to ground measurements using the methods described by Gower *et al.* [49]. The high resolution LAI reference map has a spatial resolution of 25 m and covers a 7 km × 7 km region centered on the location of the site. Estimates of error in the regression predictions were obtained using cross-validation. For the Sud-Ouest, Laprida, Zhangbei and Puechabon sites, LAI ground measurements were calculated from digital hemispherical images using the method proposed by Lang and Yueqin [50]. The high-resolution LAI reference maps were derived from the determination of the transfer function between the reflectance values of the SPOT images acquired during the ground campaign and the LAI ground measurements [51]. The LAI reference maps of the four sites have a spatial resolution of 20 m, and each LAI reference map covers a 3 km × 3 km region centered on the location of the site. The accuracy of these LAI reference maps depends on ground measurement errors, uncertainties of high-resolution reflectance data, and sampling and upscaling errors [52]. Fernandes *et al.* [53] demonstrated that an absolute uncertainty smaller than 1 LAI unit can be expected. To evaluate the accuracy of the retrieved LAI with different spatial resolutions, the high resolution LAI reference map for each selected site was aggregated to the corresponding spatial resolutions using a spatial-average method.

### 3. Results

At the Bondville site, TM surface reflectance data obtained on Day 228, 2000 and MODIS surface reflectance data obtained on Day 225, 2000 were used to retrieve multiscale LAI values. The upper row in Figure 3 shows the maps of the retrieved LAI at spatial resolutions of 28.9 m, 57.9 m, 115.8 m, 231.7 m, 463.3 m, and 926.6 m, respectively. A land cover map at the Bondville site is shown in Figure 4. The site is an agricultural district where the main crops are soybeans and corn [54]. The maps of retrieved LAI show very consistent spatial patterns with the land cover map at this site.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparisons between the retrieved multiscale LAI values and the aggregated LAI reference map values at the Bondville site. The upper row is the maps of the retrieved LAI with different spatial resolutions; the second row is the aggregated LAI reference maps and the last two rows are the difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms. Spatial resolution: (Panel 1) 28.9 m; (Panel 2) 57.9 m; (Panel 3) 115.8 m; (Panel 4) 231.7 m; (Panel 5) 463.3 m; and (Panel 6) 926.6 m. (SD: standard deviation; KURT: kurtosis; SKEW: skewness).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Land cover map of the Bondville site.</figcaption></figure>

To compare the retrieved LAI values at different scales, the aggregated LAI reference maps for each corresponding spatial resolution are shown in the second row of Figure 3. At each scale, the retrieved LAI values and LAI reference map values are in good agreement, especially in the central part. However, some differences exist in the western and eastern areas where the retrieved LAI values are slightly over-estimated relative to the LAI reference map values.

The last two rows of Figure 3 present the difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms. Seen from the difference map, the main difference between the retrieved LAI and the aggregated LAI reference map are mainly from the corn. For retrieved results, the average of retrieved LAI values of corn (approximately 6.5) is significantly higher than those of the soybeans (approximately 4.2). However, for the aggregated LAI reference map, the differences between LAI values of corn and soybeans are relatively small. Seen from the frequency histograms, it is observed that the LAI differences for most pixels at this site are concentrated around 0 at each spatial resolution. The positive mean values of LAI differences indicate that the retrieved LAI values are slightly higher than the LAI reference map values. The LAI differences at the finest scale have the largest standard deviation (Figure 3m). As the spatial resolution of the differences decreases, the standard deviation of LAI differences also gradually decreases.

At the Sud-Ouest site, ETM+ surface reflectance data obtained on Day 203, 2002 and MODIS surface reflectance data obtained on Day 201, 2002 were utilized to retrieve multiscale LAI values. Maps of the retrieved LAI values at different spatial resolutions and the aggregated LAI reference maps are displayed in the upper and second row of Figure 5, respectively. The main land cover type at this site is cropland, including corn, soybeans and sunflowers [55]. The distribution of cropland can be clearly distinguished from the map of retrieved LAI values at the spatial resolution of 28.9 m. The retrieved LAI values can reflect the different growth conditions among the fields. It is apparent that maps of the retrieved LAI values and the aggregated LAI reference maps at each spatial resolution are generally consistent in their spatial patterns. However, the retrieved LAI values in the central area are close to 6.0, and are slightly higher than those from the aggregated LAI reference maps.

The difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms (shown in the lower row of Figure 5) demonstrate that the retrieved LAI values and the aggregated LAI reference map values were in good agreement at this site. The retrieved LAI values were slightly larger than the aggregated LAI reference map values and the mean values of LAI differences for each corresponding spatial resolution were less than 0.35.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparisons between the retrieved multiscale LAI values and the aggregated LAI reference map values at the Sud-Ouest site. The upper row is the maps of the retrieved LAI with different spatial resolutions; the second row is the aggregated LAI reference maps and the last two rows are the difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms. Spatial resolution: (Panel 1) 28.9 m; (Panel 2) 57.9 m; (Panel 3) 115.8 m; (Panel 4) 231.7 m; (Panel 5) 463.3 m; and (Panel 6) 926.6 m. (SD: standard deviation; KURT: kurtosis; SKEW: skewness).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparisons between the retrieved multiscale LAI values and the aggregated LAI reference map values at the Laprida site. The upper row is the maps of the retrieved LAI with different spatial resolutions; the second row is the aggregated LAI reference maps and the last two rows are the difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms. Spatial resolution: (Panel 1) 28.9 m; (Panel 2) 57.9 m; (Panel 3) 115.8 m; (Panel 4) 231.7 m; (Panel 5) 463.3 m; and (Panel 6) 926.6 m. (SD: standard deviation; KURT: kurtosis; SKEW: skewness).</figcaption></figure>
At the Laprida site, ETM+ surface reflectance data obtained on Day 328, 2002 and MODIS reflectance data obtained on Day 313, 2002 were utilized to retrieve multiscale LAI values. Maps of the retrieved LAI values at spatial resolutions of 28.9 m, 57.9 m, 115.8 m, 231.7 m, 463.3 m, and 926.6 m are shown in the upper row of Figure 6. For comparison, the aggregated LAI reference maps for the corresponding spatial resolutions were also shown in the second row of Figure 6. Compared with the aggregated LAI reference map values, the retrieved LAI values have a larger, dynamically changing range. More detailed information of the status of vegetation growth is demonstrated in the retrieved LAI values. However, for the LAI reference map at the finest scale, there are only two LAI values over the 3 km × 3 km area. As such, it is hard to reflect the discrepancies in LAI values among different land cover types from the LAI reference maps.

The lower row of Figure 6 shows the difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms at the Laprida site. Most pixels are centered on −0.6, and range from −1.5 to −0.5 with a comparatively high skewness (around 0.6). The difference maps also demonstrate that the retrieved LAI values for most pixels are lower than those from the LAI reference map values.

At the Zhangbei site, TM surface reflectance data obtained on Day 228, 2002 and MODIS surface reflectance data obtained on Day 233, 2002 were utilized to retrieve multiscale LAI values. Maps of the retrieved LAI values at different spatial resolutions, and the aggregated LAI reference maps are shown in the upper and second row of Figure 7, respectively. The retrieved LAI values were overestimated relative to the LAI reference map values at the northeast corner of these maps. The difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms (shown in the lower row of Figure 7) demonstrate that the retrieved LAI values were in excellent agreement with the aggregated LAI reference map values at this site. However, for some parts, the retrieved LAI values for most pixels are slightly higher than the aggregated LAI reference map values.

At the Puechabon site, ETM+ surface reflectance data obtained on Day 177, 2002 and MODIS surface reflectance data obtained on Day 161, 2002 were utilized to retrieve multiscale LAI values. Maps of the retrieved LAI values at different spatial resolutions and the aggregated LAI reference maps are shown in the upper and second row of Figure 8, respectively. It can be observed that pixels with the retrieved LAI values ranging from five to six are mainly distributed in the southern area of the maps where the vegetation type is forest. Pixels with grassland biomes are mainly distributed in the northern area of the maps and show lower LAI values.

The difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms (shown in the lower row of Figure 8) demonstrate that the retrieved LAI values and the aggregated LAI reference map values are in good agreement at this site. More than half of all pixels are centered on zero, and range from −0.5 to 0.5. Relatively high kurtosis exists at most scales, which also shows that there are very few discrepancies between the retrieved LAI values and the aggregated LAI reference map values.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparisons between the retrieved multiscale LAI values and the aggregated LAI reference map values at the Zhangbei site. The upper row is the maps of the retrieved LAI with different spatial resolutions; the second row is the aggregated LAI reference maps and the last two rows are the difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms. Spatial resolution: (Panel 1) 28.9 m; (Panel 2) 57.9 m; (Panel 3) 115.8 m; (Panel 4) 231.7 m; (Panel 5) 463.3 m; and (Panel 6) 926.6 m. (SD: standard deviation; KURT: kurtosis; SKEW: skewness).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparisons between the retrieved multiscale LAI values and the aggregated LAI reference map values at the Puechabon site. The upper row is the maps of the retrieved LAI with different spatial resolutions; the second row is the aggregated LAI reference maps and the last two rows are the difference maps between the retrieved LAI and the aggregated LAI reference map for each corresponding spatial resolution and the frequency histograms. Spatial resolution: (Panel 1) 28.9 m; (Panel 2) 57.9 m; (Panel 3) 115.8 m; (Panel 4) 231.7 m; (Panel 5) 463.3 m; and (Panel 6) 926.6 m. (SD: standard deviation; KURT: kurtosis; SKEW: skewness).</figcaption></figure>
Figure 9 shows scatter plots of the retrieved LAI values at the coarsest scale and the MODIS LAI values versus the aggregated LAI reference map values at the corresponding spatial resolution over all the selected sites in Table 2. For each site, the MODIS LAI data that are closest to the date of the corresponding LAI reference map were selected for comparison. Compared with the MODIS LAI values, the retrieved LAI values in this study are distributed more closely around the 1:1 line with the aggregated LAI reference map values, which demonstrates that the retrieved LAI values achieved better agreement with the aggregated LAI reference map values across the whole range of LAI values than the MODIS LAI values. The regressions presented in Figure 9b show a positive intercept (1.07) and a slope smaller than one (0.28), indicating a systematically low bias in the MODIS retrievals at the selected sites with high LAI values, a result that is in agreement with the findings of Pisek and Chen [56]. Most of the scatters below the 1:1 line in Figure 9b belong to the Bondville, Sud-Ouest and Laprida sites. The retrieved LAI values in this study provide better accuracy with the aggregated LAI reference map values (RMSE = 0.45) compared with the MODIS LAI values (RMSE = 1.30). The correlation between the retrieved LAI values in this study and the aggregated LAI reference map values (R² = 0.79) is also superior to the correlations of the MODIS LAI values with the aggregated LAI reference map values (R² = 0.12).

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Scatter plots of (a) the retrieved LAI values and (b) the MODIS LAI values versus the aggregated LAI reference map values at a spatial resolution of 926.6 m for the five selected sites. The regression function, R², and RMSE are also shown.</figcaption></figure>

**4. Discussions**

This study proposed a method to simultaneously retrieve multiscale LAI data from satellite observations with different spatial resolutions based on EnMsF. Figure 9 demonstrates that the retrieved LAI values in this study are markedly more accurate than the MODIS LAI values when compared with the aggregated LAI reference map values from the five sites. One important reason for the better accuracy of the retrieved LAI values is that MODIS and TM/ETM+ surface reflectance data were integrated to retrieve LAI values. It is well known that the inverse technique based on radiative-transfer model is by nature an ill-posed and formidable problem because it lacks a unique solution in addition to measurement and model uncertainties [57]. Therefore more information is needed to address the ill-posed problem. The accuracy of the retrieved LAI is effectively improved by integrating multiple satellite observations, which is also demonstrated by Veger *et al.* [58] and Wang and Liang [25].

Figure 10a shows the map of differences between the LAI values retrieved from MODIS and TM surface reflectance data and those retrieved only from TM surface reflectance data at the spatial
resolution of 28.9 m and Figure 10b shows the map of differences between the LAI values retrieved from MODIS and TM surface reflectance data and those retrieved only from MODIS surface reflectance data of 926.6 m at the Bondville site. Results indicate that the LAI values retrieved only from TM surface reflectance data achieve good agreement with the LAI values retrieved from MODIS and TM surface reflectance data. This is because most details of the LAI values retrieved from the MODIS and TM surface reflectance data are coming from TM surface reflectance data. However, the LAI values retrieved from MODIS and TM surface reflectance data were much larger than those retrieved only from MODIS surface reflectance data (Figure 10b). Figure 9 demonstrates that the LAI values retrieved from MODIS and TM surface reflectance data achieved good agreement with the aggregated LAI reference map values across the whole range of LAI values. In other words, the method in this study can provide improved mapping of LAI at the MODIS resolution by integrating the MODIS and TM/ETM+ surface reflectance data.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a) The differences between the LAI values retrieved from MODIS and TM surface reflectance data and those retrieved only from TM reflectance data at the spatial resolution of 28.9 m at Bondville site; (b) The differences between the LAI values retrieved from MODIS and TM surface reflectance data and those retrieved only from MODIS reflectance data at the spatial resolution of 926.6 m at Bondville site.</figcaption></figure>

For the next step of the study, a dynamical model of LAI that evolves in time will be coupled with the EnMsT. Even though TM/ETM+ reflectance data contains more information compared with MODIS reflectance data on the same day, low time resolution limits its application in the inversion of LAI for time series. If there exist both TM/ETM+ reflectance data and MODIS reflectance data on that day, more details will be got from TM/ETM+ reflectance data. However, if only MODIS reflectance data exists on the day, LAI retrieved from previous time will be used as the input data and MODIS reflectance data will provide more information about changing dynamics of LAI in time series.

In addition, only quad tree is used during the construction of the EnMsT in this study. The EnMsT with 4k number of nodes can make the calculation and formula derivation easier. However, most geographical regions have an arbitrary shape and number of pixels. Thus, in the near future, more complex structure of the EnMsT like neural gas (NG) algorithm [36] will be considered to describe irregular shapes of geographical regions.
## 5. Conclusions

In this study, a new method was developed to retrieve multiscale LAI values from satellite observations at different spatial resolutions. The average of multi-year MODIS LAI values corresponding to the date of satellite observations was used as a priori knowledge for LAI to construct an initial EnMsT, and a two-sweep filtering procedure was used to update LAI values at each node of the EnMsT by integrating multiple satellite observations with different spatial resolutions. Then, the retrieved LAI values at the finest scale were used as a priori knowledge for LAI for the new round of construction and updating of the EnMsT, until the sum of the difference of LAI values at each node of the EnMsT between two adjacent updates was less than a given threshold.

The method introduced canopy radiative transfer models into the EnMsF to simultaneously retrieve multiscale LAI from MODIS and TM/ETM+ surface reflectance data. Validation of the retrieved LAI at different spatial resolutions against the aggregated LAI reference maps was carried out over several sites with different vegetation types. The results demonstrated that the retrieved LAI values at each spatial resolution were in good agreement with the aggregated LAI reference map values at the corresponding spatial resolution. The retrieved LAI values at the coarsest scale provided better accuracy with the aggregated LAI reference map values compared with the MODIS LAI values.

In this study, the method was only applied to retrieve LAI values by integrating MODIS and TM/ETM+ surface reflectance data at a specific time. In the near future, we will extend the methodology to retrieve multiscale LAI values from time series data obtained by satellite sensors at different spatial resolutions.

**Acknowledgments:** This work was financially supported by the National Natural Science Foundation of China under grant No 41171264 and the Chinese 973 Program under grant No. 2013CB733403.

**Author Contributions:** Zhiqiang Xiao conceived and designed the study. Jingyi Jiang contributed to the conception of the study, performed the data analysis and wrote the paper. Zhiqiang Xiao reviewed and edited the manuscript. All authors read and approved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Jonckheere, I.; Fleck, S.; Nackaerts, K.; Muys, B.; Coppin, P.; Weiss, M.; Baret, F. Review of methods for in situ leaf area index determination: Part I. Theories, sensors and hemispherical photography. *Agric. For. Meteorol.* **2004**, *121*, 19–35. [CrossRef]

2. Gassman, P.W.; Williams, J.R.; Benson, V.W.; Izaurralde, R.C.; Hauck, L.M.; Jones, C.A.; Atwood, J.D.; Kiniry, J.R.; Flowers, J.D. *Historical Development and Applications of the Epic and Apex Models*; Center for Agricultural and Rural Development, Iowa State University: Ames, IA, USA, 2005.

3. Simic, A.; Fernandes, R.; Wang, S. Assessing the impact of leaf area index on evapotranspiration and groundwater recharge across a shallow water region for diverse land cover and soil properties. *J. Water Resour. Hydraul. Eng.* **2014**, *3*, 60–73.

4. Gibelin, A.L.; Calvet, J.C.; Roujean, J.L.; Jarlan, L.; Los, S.O. Ability of the land surface model ISBA-A-GS to simulate leaf area index at the global scale: Comparison with satellites products. *J. Geophys. Res. Atmos.* **2006**, *111*, [CrossRef]

5. Boussetta, S.; Balsamo, G.; Beljaars, A.; Kral, T.; Jarlan, L. Impact of a satellite-derived leaf area index monthly climatology in a global numerical weather prediction model. *Int. J. Remote Sens.* **2013**, *34*, 3520–3542. [CrossRef]

6. Buermann, W.; Dong, J.; Zeng, X.; Myneni, R.B.; Dickinson, R.E. Evaluation of the utility of satellite-based vegetation leaf area index data for climate simulations. *J. Clim.* **2001**, *14*, 3536–3550. [CrossRef]

7. Tillack, A.; Clasen, A.; Kleinschmit, B.; Förster, M. Estimation of the seasonal leaf area index in an alluvial forest using high-resolution satellite-based vegetation indices. *Remote Sens. Environ.* **2014**, *141*, 52–63. [CrossRef]

8. Viña, A.; Gitelson, A.A.; Nguyen-Robertson, A.L.; Peng, Y. Comparison of different vegetation indices for the remote assessment of green leaf area index of crops. *Remote Sens. Environ.* **2011**, *115*, 3468–3478. [CrossRef]
9. Atzberger, C.; Richter, K. Spatially constrained inversion of radiative transfer models for improved LAI mapping from future sentinel-2 imagery. *Remote Sens. Environ.* **2012**, *120*, 208–218. [CrossRef]

10. Yin, G.; Li, J.; Liu, Q.; Fan, W.; Xu, B.; Zeng, Y.; Zhao, J. Regional leaf area index retrieval based on remote sensing: The role of radiative transfer model selection. *Remote Sens.* **2015**, *7*, 4604–4625. [CrossRef]

11. Yang, W.; Tan, B.; Huang, D.; Rautiainen, M.; Shabanov, N.V.; Wang, Y.; Privette, J.L.; Huemrich, K.F.; Fensholt, R.; Sandholt, I. MODIS leaf area index products: From validation to algorithm improvement. *Geosci. Remote Sens. IEEE Trans.* **2006**, *44*, 1885–1898. [CrossRef]

12. Baret, F.; Hagolle, O.; Geiger, B.; Bicheron, P.; Miras, B.; Huc, M.; Berthelot, B.; Niño, F.; Weiss, M.; Samain, O. LAI, fAPAR and FCOVER cyclopes global products derived from VEGETATION: Part 1: Principles of the algorithm. *Remote Sens. Environ.* **2007**, *110*, 275–286. [CrossRef]

13. Heiskanen, J. Estimating aboveground tree biomass and leaf area index in a mountain birch forest using aster satellite data. *Int. J. Remote Sens.* **2006**, *27*, 1135–1158. [CrossRef]

14. Ganguly, S.; Nemani, R.R.; Zhang, G.; Hashimoto, H.; Milesi, C.; Michaelis, A.; Wang, W.; Votava, P.; Samanta, A.; Melton, F. Generating global leaf area index from landsat: Algorithm formulation and demonstration. *Remote Sens. Environ.* **2012**, *122*, 185–202. [CrossRef]

15. Garrigues, S.; Lacaze, R.; Baret, F.; Morisette, J.; Weiss, M.; Nickeson, J.; Fernandes, R.; Plummer, S.; Shabanov, N.; Myneni, R. Validation and intercomparison of global leaf area index products derived from remote sensing data. *J. Geophys. Res. Biogeosci.* **2008**, *113*. [CrossRef]

16. Fang, H.; Wei, S.; Liang, S. Validation of modis and cyclopes LAI products using global field measurement data. *Remote Sens. Environ.* **2012**, *119*, 43–54. [CrossRef]

17. Xiao, Z.; Liang, S.; Wang, J.; Song, J.; Wu, X. A temporally integrated inversion method for estimating leaf area index from MODIS data. *Geosci. Remote Sens. IEEE Trans.* **2009**, *47*, 2536–2545. [CrossRef]

18. Jiang, J.; Xiao, Z.; Wang, J.; Song, J. Sequential method with incremental analysis update to retrieve leaf area index from time series MODIS reflectance data. *Remote Sens.* **2014**, *6*, 9194–9212. [CrossRef]

19. Lewis, P.; Gómez-Dans, J.; Kaminski, T.; Settle, J.; Quaife, T.; Gobron, N.; Styles, J.; Berger, M. An earth observation land data assimilation system (EO-LDAS). *Remote Sens. Environ.* **2012**, *120*, 219–235. [CrossRef]

20. Viskari, T.; Hardiman, B.; Desai, A.R.; Dietze, M.C. Model-data assimilation of multiple phenological observations to constrain and predict leaf area index. *Ecol. Appl.* **2015**, *25*, 546–558. [CrossRef] [PubMed]

21. Wang, D.; Liang, S. Improving LAI mapping by integrating MODIS and cyclopes LAI products using optimal interpolation. *Sel. Top. Appl. Earth Observ. Remote Sens. J.* **2014**, *7*, 445–457. [CrossRef]

22. Nazelle, A.D.; Arunachalam, S.; Serre, M.L. Bayesian maximum entropy integration of ozone observations and model predictions: An application for attainment demonstration in north carolina. *Environ. Sci. Technol.* **2010**, *44*, 5707–5713. [CrossRef] [PubMed]

23. Cocchi, D.; Greco, F.; Trivisano, C. Hierarchical space-time modelling of PM 10 pollution. *Atmos. Environ.* **2007**, *41*, 532–542. [CrossRef]

24. Xiao, Z.; Wang, J.; Wan, H.; Jiang, X.; Xu, K. Multiscale Approach for Fusing Leaf Area Index Estimates from Multiple Sensors. Available online: http://spie.org/Publications/Proceedings/Paper/10.1117/12.748313 (accessed on 26 September 2015).

25. Wang, D.; Liang, S. Integrating MODIS and cyclopes leaf area index products using empirical orthogonal functions. *Geosci. Remote Sens. IEEE Trans.* **2011**, *49*, 1513–1519. [CrossRef]

26. Chou, K.C. *A Stochastic Modelling Approach to Multiscale Signal Processing*; Massachusetts Institute of Technology: Cambridge, MA, USA, 1991.

27. Parada, L.M.; Liang, X. Optimal multiscale kalman filter for assimilation of near-surface soil moisture into land surface models. *J. Geophys. Res. Atmos.* **2004**, *109*. [CrossRef]

28. Van de Vyver, H.; Roulin, E. Scale-recursive estimation for merging precipitation data from radar and microwave cross-track scanners. *J. Geophys. Res. Atmos.* **2009**, *114*. [CrossRef]

29. Tao, H.; Shunlin, L.; Dongdong, W.; Yanmin, S.; Yunyue, Y. Fusion of satellite land surface albedo products across scales using a multiresolution tree method in the north central United States. *Geosci. Remote Sens. IEEE Trans.* **2014**, *52*, 3428–3439.

30. Wang, D.; Liang, S. Using multiresolution tree to integrate MODIS and MISR-L3 LAI products. In Proceedings of the 2010 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Honolulu, HI, USA, 25–30 July 2010; pp. 1027–1030.
31. Zhou, Y.; McLaughlin, D.; Entekhabi, D.; Ng, G.-H.C. An ensemble multiscale filter for large nonlinear data assimilation problems. *Mon. Weather Rev.* **2008**, *136*, 678–698. [CrossRef]

32. Frakt, A.B.; Willsky, A.S. Computationally efficient stochastic realization for internal multiscale autoregressive models. *Multidimens. Syst. Signal Process.* **2001**, *12*, 109–142. [CrossRef]

33. Evensen, G. Inverse methods and data assimilation in nonlinear ocean models. *Phys. D Nonlinear Phenom.* **1994**, *77*, 108–129. [CrossRef]

34. Zhou, Y. Multi-Sensor Large Scale Land Surface Data Assimilation Using Ensemble Approaches; Massachusetts Institute of Technology: Cambridge, MA, USA, 2006.

35. Lawniczak, W.; Hanea, R.; Heemink, A.; McLaughlin, D. Multiscale ensemble filtering for reservoir engineering applications. *Comput. Geosci.* **2009**, *13*, 245–254. [CrossRef]

36. Pan, M.; Wood, E.F.; McLaughlin, D.B.; Entekhabi, D.; Luo, L. A multiscale ensemble filtering system for hydrologic data assimilation. Part I: Implementation and synthetic experiment. *J. Hydrometeorol.* **2009**, *10*, 794–806. [CrossRef]

37. Pan, M.; Wood, E.F. A multiscale ensemble filtering system for hydrologic data assimilation. Part II: Application to land surface modeling with satellite rainfall forcing. *J. Hydrometeorol.* **2009**, *10*, 1493–1506. [CrossRef]

38. Pan, M.; Wood, E.F. Impact of accuracy, spatial availability, and revisit time of satellite-derived surface soil moisture in a multiscale ensemble data assimilation system. *Sel. Top. Appl. Earth Observ. Remote Sens. IEEE J.* **2010**, *3*, 49–56. [CrossRef]

39. Kuusk, A. A Markov chain model of canopy reflectance. *Agric. For. Meteorol.* **1995**, *76*, 221–236. [CrossRef]

40. Kuusk, A. A two-layer canopy reflectance model. *J. Quant. Spectrosc. Radiat. Transfer* **2001**, *71*, 1–9. [CrossRef]

41. Fang, H.; Liang, S.; Kuusk, A. Retrieving leaf area index using a genetic algorithm with a canopy radiative transfer model. *Remote Sens. Environ.* **2003**, *85*, 257–270. [CrossRef]

42. Xiao, Z.; Liang, S.; Wang, J.; Xie, D.; Song, J.; Fensholt, R. A framework for consistent estimation of leaf area index, fraction of absorbed photosynthetically active radiation, and surface albedo from MODIS time-series data. *IEEE Trans. Geosci. Remote Sens.* **2014**, *53*, 3178–3197. [CrossRef]

43. Saltelli, A.; Bolado, R. An alternative way to compute fourier amplitude sensitivity test (fast). *Comput. Stat. Data Anal.* **1998**, *26*, 445–460. [CrossRef]

44. Saltelli, A.; Tarantola, S.; Chan, K.-S. A quantitative model-independent method for global sensitivity analysis of model output. *Technometrics* **1999**, *41*, 39–56. [CrossRef]

45. Willsky, A.S. Multiresolution markov models for signal and image processing. *Proc. IEEE* **2002**, *90*, 1396–1458. [CrossRef]

46. Evensen, G. The ensemble kalman filter: Theoretical formulation and practical implementation. *Ocean Dyn.* **2003**, *53*, 343–367. [CrossRef]

47. EOS Validation Status for MODIS Surface Reflectance: MOD09. Available online: http://landval.gsfc.nasa.gov/ProductStatus.php?ProductID=MOD09 (accessed on 26 September 2015).

48. Masek, J.G.; Vermote, E.F.; Saleous, N.E.; Wolfe, R.; Hall, F.G.; Huemmrich, K.F.; Gao, F.; Kutler, J.; Lim, T.-K. A landsat surface reflectance dataset for North America, 1990–2000. *Geosci. Remote Sens. Lett.* *IEEE* **2006**, *3*, 68–72.

49. Gower, S.T.; Kucharik, C.J.; Norman, J.M. Direct and indirect estimation of leaf area index, fAPAR, and net primary production of terrestrial ecosystems. *Remote Sens. Environ.* **1999**, *70*, 29–51. [CrossRef]

50. Lang, A.; Yueqin, X. Estimation of leaf area index from transmission of direct sunlight in discontinuous canopies. *Agric. Forest Meteorol.* **1986**, *37*, 229–243. [CrossRef]

51. Baret, F.; Weiss, M.; Allard, D.; Garrigues, S.; Leroy, M.; Jeanjean, H.; Fernandes, R.; Myneni, R.; Privette, J.; Morisette, J. VALERI: A Network of Sites and a Methodology for the Validation of Medium Spatial Resolution Land Satellite Products. Available online: http://w3.avignon.inra.fr/valeri/documents/VALERI-RSESubmitted.pdf (accessed on 26 September 2015).

52. Weiss, M. VALERI 2003: Barrax Site (Cropland). Ground Data Processing & Production of the Level 1 High Resolution Maps; Inra-Cse, Avignon. Available online: http://w3.avignon.inra.fr/valeri/Europe/Spain/Barrax/2003/biomap/Barrax2003FTReport.pdf (accessed on 26 September 2015).

53. Fernandes, R.; Butson, C.; Leblanc, S.; Latifovic, R. Landsat-5 TM and Landsat-7 ETM+ based accuracy assessment of leaf area index products for canada derived from spot-4 vegetation data. *Can. J. Remote Sens.* **2003**, *29*, 241–258. [CrossRef]
54. Meyers, T.P.; Hollinger, S.E. An assessment of storage terms in the surface energy balance of maize and soybean. *Agric. Forest Meteorol.* **2004**, *125*, 105–115. [CrossRef]

55. VALERI. Available online: http://www.avignon.inra.fr/valeri (accessed on 26 September 2015).

56. Pisek, J.; Chen, J.M. Comparison and validation of MODIS and VEGETATION global LAI products over four BigFoot sites in North America. *Remote Sens. Environ.* **2011**, *109*, 81–94. [CrossRef]

57. Combal, B.; Baret, F.; Weiss, M.; Trubuil, A.; Mace, D.; Pragnere, A.; Myneni, R.; Knyazikhin, Y.; Wang, L. Retrieval of canopy biophysical variables from bidirectional reflectance: Using prior information to solve the ill-posed inverse problem. *Remote Sens. Environ.* **2003**, *84*, 1–15. [CrossRef]

58. Verger, A.; Baret, F.; Weiss, M. A multisensor fusion approach to improve LAI time series. *Remote Sens. Environ.* **2011**, *115*, 2460–2470. [CrossRef]

<figure><img src="image_14.png" /></figure>

© 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons by Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).