Estimating leaf area index from MODIS and surface meteorological data using a dynamic Bayesian network

Remote Sensing of Environment 127 (2012) 30–43

Contents lists available at SciVerse ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Estimating leaf area index from MODIS and surface meteorological data using a dynamic Bayesian network Yuzhen Zhang a, b, Yonghua Qu a, c, d,⁎, Jindi Wang a, c, d, Shunlin Liang a, b, e, Yan Liu a, c, d a

State Key Laboratory of Remote Sensing Science, Jointly Sponsored by Beijing Normal University and the Institute of Remote Sensing Applications of CAS, Beijing, China College of Global Change and Earth System Science, Beijing Normal University, Beijing, China Beijing Key Laboratory for Remote Sensing of Environment and Digital Cities, Beijing Normal University, Beijing, China d School of Geography, Beijing Normal University, Beijing, China e Department of Geography, University of Maryland, College Park, MD 20742, USA b c

a r t i c l e

i n f o

Article history: Received 9 September 2011 Received in revised form 27 July 2012 Accepted 11 August 2012 Available online 14 September 2012 Keywords: Leaf area index Dynamic Bayesian networks MODIS Ground meteorological station data Filtering inference algorithm

a b s t r a c t Remotely sensed data is the main source of vegetation leaf area index (LAI) information on the regional to global scale. Many validation results have revealed that the accuracy of the retrieved LAI is often affected by the cloud cover of imagery, instrument problems, and inversion algorithms. Ground meteorological station data, characterized by relatively high accuracy and time continuity compared with remote sensing data, can provide complementary information to remote sensing observations. In this paper, we combine the potential advantages of both types of data in order to improve LAI retrievals in the Heihe River Basin, an arid and semi-arid area in northwest China where Moderate Resolution Imaging Spectroradiometer (MODIS) LAI values are significantly underestimated. A dynamic Bayesian network (DBN) is used to integrate these two data types for time series LAI estimation. Results show that the square of correlation coefficient between LAI values estimated by our DBN method (referred to as DBN LAI) and field measured LAI values is 0.76, with a root mean square error of 0.78. The DBN LAI are closer to field measurements than the MODIS LAI standard product values. Moreover, by introducing ground meteorological station data using a dynamic process model, DBN LAI show better temporal consistency than the MODIS LAI. It is concluded that the quality of LAI retrievals can be improved by combining remote sensing data and ground meteorological station data using a filtering inference algorithm in a DBN framework. More importantly, the study provides a basis and method for utilizing ground meteorological station network data to estimate land surface parameters on a regional scale. © 2012 Elsevier Inc. All rights reserved.

1. Introduction Leaf area index (LAI), defined as one half of the total leaf area per unit ground surface area projected on the local horizontal datum, is an important biophysical vegetation variable in many models describing vegetation–atmosphere interactions (Chen & Black, 1992; Chen et al., 1997). Routine in situ LAI measurements are time-consuming, laborintensive, and even unfeasible for large-scale studies. The development of remote sensing techniques has facilitated the retrieval of LAI on a regional scale. Current methods for estimating LAI from remote sensing data are classified into three categories: (1) the statistical relationship between LAI and vegetation indices (Hasegawa et al., 2010; Wang et al., 2007; Wu et al., 2007; Yang et al., 2007); (2) inversion of

⁎ Corresponding author at: 19th Xinjiekouwai Street, Beijing, China. Tel.: +86 10 58802041; fax: +86 10 58807452. E-mail address: [email protected] (Y. Qu). 0034-4257/$ – see front matter © 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.rse.2012.08.015

canopy reflectance models (Houborg et al., 2009; le Maire et al., 2011); and (3) hybrid inversion methods combining the two methods mentioned above (Duveiller et al., 2011; Fang & Liang, 2005; Kalacska et al., 2005; Qu et al., 2008). Global LAI products, such as MODIS, CYCLOPES, ECOCLIMAP, and GLOBCARBON, are generated from remote sensing observations using these methods (Garrigues et al., 2008). However, remote sensing data are influenced by atmospheric conditions; this well-known fact may result in inaccuracies and temporal discontinuities in the inversion results. Current global LAI products often face this problem, and this restricts their further application in many fields (De Kauwe et al., 2011; Hill et al., 2006). Several approaches have been adopted to solve the problem. Fang et al. (2008) developed a temporal spatial filter algorithm to fill the gaps and improve the quality of MODIS LAI products. Yuan et al. (2011) implemented an integrated two-step method to reprocess the MODIS LAI data for land surface and climate modeling. Verger et al. (2011) proposed a multisensory fusion approach to improve the spatio-temporal continuity, consistency, and accuracy of current satellite products. In addition to these smoothing methods, many

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

advances have been made through data assimilation methods. Researchers have used the temporal information of vegetation growth in land surface studies to estimate LAI and achieved satisfactory results (Koetz et al., 2005; Moulin et al., 1998; Wang et al., 2010; Weiss et al., 2001; Xiao et al., 2011). This study aims to integrate the ground station data into the LAI estimation in order to solve the problem faced by current LAI products. Previous studies have shown the potential advantage of combining data from remote sensing with pinpoint information observed on the ground (Douaoui et al., 2006; Tian & Chen, 2010). Compared with remote sensing data, ground station data are generally considered to be accurate; however, these data represent small areas around point stations and cannot directly provide sufficient areal coverage. Remote sensing data and ground station observations complement each other well. In this study, we introduced the ground meteorological station data and related it to vegetation growth state via a relatively simple vegetation growth model. Many complex mechanism models include detailed processes within an agronomic cropping system, such as soil water flow, photosynthesis, and nutrient balance. Vegetation growth in these models is controlled by some environmental factors. If soil and other conditions are suitable for vegetation growth, assuming that temperature, sunshine and water affect both temporal and spatial variations of the vegetation, prognostic phenology model driven by meteorological factors can simulate dynamic LAI changes on a regional scale. Using the model, ground station data can be introduced into LAI estimation to compensate for the deficiency arising from reliance on remote sensing data alone. The main objective of this paper is to show the feasibility and effectiveness of combining these data to improve the quality of the estimated LAI, by using the advantages of both the data. However, we should note that remote sensing data and ground station data differ in spatial scale and accuracy, and in addition, are temporally inconsistent. Fortunately, data assimilation provides an objective method for inferring the state of a dynamic system from these heterogeneous, irregularly distributed, and temporally inconsistent observational data with differing accuracies (Liang, 2004). Many algorithms have been developed to combine these data to estimate land surface parameters, such as the ensemble Kalman filter and variational analysis algorithms (Dente et al., 2008; Pipunic et al., 2008; Quaife et al., 2008). Because both remote sensing and ground station observations have some uncertainties, the estimations should integrate the data with their respective uncertainties in order to obtain an optimal estimated result from the inversion (Han & Li, 2008). Dynamic Bayesian networks (DBNs) have an inherent advantage in dealing with uncertainties, because they seek the posterior distribution, which integrates all of the available information expressed by probabilities. Many recent efforts attempting to use the DBN to infer the dynamic information from varied data sources are reported in the literature. Wikle and Berliner (2007) described the theoretical framework of a DBN for the estimation of dynamic state parameters and reviewed the linkages between optimal interpolation, kriging, Kalman filtering, smoothing, and variational analysis from the perspective of data assimilation. Challa and Koks (2004) and Ihler et al. (2007) derived the theoretical result that, when the dynamic process model and observation model are linear and the noises obey a Gaussian distribution, the DBN estimation method can be equivalent to the Kalman filter algorithm. In specific applications, estimates in the framework of a DBN mainly involve state parameters related to the hydrological cycle. For example, Moradkhani (2008) estimated soil moisture content using microwave data in the framework of a DBN. Ihler et al. (2007) and Barillec and Cornford (2009) used a DBN for nowcasting precipitation and rainfall estimations. In our recent study, DBN was proven to be effective in retrieving LAI at the plot and site scales (Qu et al., 2012). Here, we try to infer the state of target LAI parameters from both data in the framework of a DBN and evaluate its performance through inversion on the regional scale, which is another objective of this paper.

31

2. Data 2.1. Study area The Heihe River Basin (marked “A” in Fig. 1) is the second largest inland river basin in China, covering an area of approximately 130,000 km 2. It is located between 97°24′–102°10′ E and 37°41′– 42°42′ N, and covers parts of the Qinghai province, the Gansu province, and Inner Mongolia. The local climate is characterized by large temperature differences and sufficient heat and light resources. The annual average temperature is 6–8 °C. The average annual precipitation is 140 mm, whereas the annual evaporation is 1410 mm. In 2008, a comprehensive remote sensing campaign, known as the Watershed Allied Telemetry Experimental Research (Li et al., 2009), was carried out in this typical inland river basin. In the present study, a subarea (99°40′–101° E, 38°30′–39°30′ N) was selected as our research area (marked “B” in Fig. 1). Field campaigns were completed over the period May 15 to July 22, 2008, in this region. Forests, shrubs, grass, crops, and desert are the main land cover types, according to the MODIS International Geosphere-Biosphere Program (IGBP) land cover type classifications. A number of parameters related to grassland and crops were measured, including spectral and structural vegetation parameters. Leaf reflectance, leaf transmittance, and soil reflectance were measured using Analytical Spectral Devices (ASD) FieldSpec Spectroradiometer. These provided the basis for determining key inputs of the canopy reflectance model. LAI was measured directly using a LAI-2000 Plant Canopy Analyzer and these measurements were used to validate the estimated results. 2.2. Meteorological data Routine meteorological observations are used to provide relatively accurate and temporally continuous data for the estimation of structural vegetation parameters. In this study, only the minimum temperature, sunshine hours, and vapor pressure deficit are used in the estimation. These were downloaded from the Chinese meteorological data-sharing service (http://cdc.cma.gov.cn). Data on the daily minimum temperature and sunshine hours can be obtained directly from the website. The vapor pressure deficit was calculated from the observed mean temperature and relative humidity (Buck, 1981; Campbell & Norman, 1998):   b0 Te es ðTeÞ ¼ a0 exp Te þ c0

ð1Þ

VPD ¼ es ðTeÞ−ea ¼ es ðTeÞð1−RH Þ

ð2Þ

where es(Te) is the saturation vapor pressure corresponding to temperature Te, a0 =0.611 kPa, b0 =17.502, c0 =240.97 °C, ea is the vapor pressure of moist air at ambient temperature, and RH is relative humidity. The availability of meteorological data in certain areas is limited and tends to be distributed discretely and unevenly in space. Thus, before being used in combination with remote sensing imagery data, the meteorological data should be interpolated into the grid to match the remote sensing data on a regional scale. In this study, we carried out the interpolation using thin plate smoothing splines implemented in the ANUSPLIN package (Hutchinson, 2001). ANUSPLIN facilitates the incorporation of a linear parametric sub-model, which permits a known physical control of the variable to be incorporated as a dependent variable. This makes the interpolation results more accurate than those by other methods and is convenient for processing several surfaces simultaneously. The first step in the ANUSPLIN interpolation was to generate spline coefficients from the database using the SPLINA program. A set of parameters was obtained from the SPLINA program, such as coefficients for the fitted surface, data and fitted values, and estimated standard errors. The surface coefficients and their error covariance matrices were

32

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

A

B

Fig. 1. Location of study area.

then used to develop grid files of meteorological factors in the LAPGRD program.

2.3. MODIS data The MODIS data used in this paper are mainly the surface reflectance product MOD09A1 at a spatial resolution of 500 m (Julian Day 89–265, 2008), LAI product MOD15A2 at a spatial resolution of 1 km (Julian Day 89–265, 2004–2008), and Land Cover product MOD12Q1 at a spatial resolution of 1 km. All of the data were obtained using the Warehouse Inventory Tool (WIST) from NASA's Land Processes Distributed Active Archive Center (https://wist-ops.echo.nasa.gov/api/) and were then re-projected from sinusoidal to geographic coordinates using the MODIS Reprojection Tool available at http://lpdaac.usgs.gov/ lpdaac/tools/modis_reprojection_tool. MOD09A1 data from Julian Day 89–265 served as a remote sensing observation time series to estimate LAI. These data consist of eight-day composites of the surface spectral reflectance. The corresponding sun and view geometries for every pixel were given as individual layer data in the same HDF file of reflectance, indicating that each pixel in an image was collected under different viewing geometry characteristics and from more than one orbit and day within the eight-day period. MOD15A2 data from Julian Day 89–265, 2008, served as references against the estimated LAI to evaluate the DBN LAI indirectly. Other MOD15A2 data were used to carry out regression analysis together with gridded meteorological data, through which meteorological station data can be introduced into the estimation. LAI values reported in the MODIS product are derived from the atmospherically corrected surface reflectance product MOD09 by inversion of a 3D radiative transfer model using a look-up table method (main algorithm). If the main algorithm fails, a back-up algorithm is triggered to estimate LAI based on relations between the normalized difference vegetation index (NDVI) and LAI (Knyazikhin et al., 1999). According to the MOD15A2 quality assessment scheme, MODIS LAI values retrieved by the main algorithm are the most reliable, whereas those retrieved by the back-up algorithm are less reliable, as the back-up algorithm is usually employed when cloud cover, strong atmospheric effects, or snow/ice are detected (Myneni et al., 2002). In addition to the inversion algorithm, cloud contamination is a limitation of the MODIS data. Thus, in this study, only MOD15A2 data with quality control (QC) flags of zero, i.e., LAI retrieved on clear-sky days by the main algorithm, were selected for use alongside meteorological data to bridge the gap between meteorological data and vegetation LAI.

3. Methodology The process of combining remote sensing and ground station observations to estimate LAI in this study is shown in Fig. 2. Ground stations can provide various data, including energy fluxes, land cover, meteorological data, and soil data; however, only the meteorological data were used in the present study. These meteorological parameters control the vegetation growth and can thus provide helpful dynamic information for the LAI estimation. The vegetation growth or LAI dynamics were described by a dynamic process model that bridges the gap between ground meteorological data and the target LAI parameter. The MODIS time series reflectance data MOD09A1 served as remote sensing observations. Only red band reflectance and near-infrared (NIR) reflectance were used in the study. Canopy reflectance models provide a direct physical characterization of the relationship between LAI and ancillary parameters and the radiometric response from satellite data. In the forward mode, canopy reflectance models compute the spectral reflectance for a certain set of leaf and canopy parameters and generate look up tables (LUTs). The dashed part of Fig. 2 is the core of the inversion algorithm. By performing inference in a DBN, ground station data and

ground station data

canopy reflectance model MOD09A1

dynamic process model

Look Up Table

likelihood probability

state transition matrix DBN inference algorithm

posterior probability distribution of LAI

MOD15A2

DBN LAI

Field measured LAI

Validation and evaluation of DBN LAI Fig. 2. Flowchart of LAI estimation by combining ground station data and remote sensing data in a DBN framework.

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

33

MODIS surface reflectance data, together with the dynamic process model and the canopy reflectance model, will be integrated into the calculation of the posterior probability distribution of LAI. Estimated LAI values in the DBN framework (DBN LAI) are the mean values of the posterior probability distribution of LAI. MODIS LAI and field measured LAI are used as reference values to validate and evaluate DBN LAI in the study.

where LAIt × (aLAIt + b) is a logistic growth function that constrains growth at low and high LAI values and deltaGSI is the growth factor. deltaGSI is positive when the potential phenological state (GSI) based on current meteorological conditions is above the current prognostic phenological state P(−), and vice versa. The value of P(−) used in the prognostic phenology model is described as follows (Stockli et al., 2008):

3.1. Modeling the dynamic process

P ð−Þ ¼

The dynamic process model can describe the growth of vegetation and thus evolution of LAI time series. In this study, the dynamic process model is the transformation and simplified form of a prognostic phenology model used by Stockli et al. (2008) and is driven by several meteorological factors including minimum temperature, vapor pressure deficit, and global radiation or sunshine hours. It can provide ancillary LAI dynamic information for retrieval. The dynamic process model is expressed as LAI ðt Þ ¼ LAI ðt−1Þ þ LAIðt−1Þ  ðaLAIðt−1Þ þ bÞðGSIðt Þ−cð expðdLAIðt−1ÞÞ þ eÞÞ

ð3Þ

where LAI(t) and LAI(t-1) are LAI at time slices t and t-1, and GSI(t) is the growing season index at time slice t. The analytic equation for GSI(t) is simply the product of three gridded meteorological factors (Jolly et al., 2005): 


ðt Þ ⋅f VPD ðt Þ GSIðt Þ ¼ f T
m ðt Þ ⋅f R g ¼

T m ðt Þ−T m min Rg ðt Þ−Rg min VPD max −VPD ðt Þ ⋅ R −R ⋅ VPD −VPD −T T m max

m min

g max

g min

max

ð4Þ

min


ðt Þ, and VPD ðt Þ are 20-day moving averages for daily where T
m ðt Þ, R g Tm(t), Rg(t), and VPD(t). Tm is the minimum daily temperature, Rg is the mean daily global radiation or sunshine hours, and VPD is the mean daily vapor pressure deficit. Tmmax, Rgmax, and VPDmax are the maximum and Tmmin, Rgmin, and VPDmin are the minimum values of Tm, Rg, and VPD, respectively (Table 1). These were determined from two types of information—errors in the interpolation of each meteorological factor and reported information pertaining to favorable conditions for grass and crops (Qu et al., 2012; Stockli et al., 2008; Yao et al., 2008). In Eq. (3), a, b, c, d, and e are undetermined coefficients that can be fitted by a non-linear regression method using GSI data and the corresponding LAI data. The former can be calculated from meteorological data, and the latter are selected from the global MODIS LAI products, which are available for free download. These undetermined coefficients are fitted separately for different land cover types. The derivation of the dynamic process model in Eq. (3) is as follows: LAI tþ1 ¼ LAI t þ delta LAI

ð5Þ

delta LAI ¼ delta GSI  LAI t  ðaLAIt þ bÞ

ð6Þ

delta GSI ¼ GSItþ1 −P ð−Þ ¼ GSItþ1 −c  ð expðdLAIt Þ þ eÞ

ð7Þ

Table 1 Maximum and minimum values of meteorological factors in dynamic process model for crops and grassland. Vegetation type

Parameter

Unit

Minimum value

Maximum value

Crops

Tm Rg VPD Tm Rg VPD

°C h kPa °C h kPa

0 7 2.0 0 6 0.5

15 11 3.5 15 10 2.5

Grassland

FPAR−FPAR min FPAR max −FPAR min

ð8Þ

where FPARmax and FPARmin are the maximum and minimum values of FPAR, respectively, which can be related to LAI by use of Beer's law (Sellers et al., 1996): LAI ¼ LAImax 

logð1−FPARÞ logð1−FPAR max Þ

ð9Þ

Eq. (9) can be simplified by defining parameter d from Eq. (7), which is a function of LAImax and FPARmax, such that dLAI

FPAR ¼ 1−e

ð10Þ

From Eqs. (8)–(10), P(−) can be deduced to be 1− expðdLAIÞ−ð1− expðdLAImin ÞÞ 1− expðdLAImax Þ−ð1− expðdLAImin ÞÞ expðdLAIÞ− expðdLAImin Þ ¼ cð expðdLAIÞ þ eÞ ¼ expðdLAImax Þ− expðdLAImin Þ

P ð−Þ ¼

ð11Þ

where c, d, and e are related to the maximum and minimum values of LAI. Thus, Eq. (3) can be derived from Eqs. (5)–(7) and (11). Undetermined coefficients a, b, c, d, and e, as well as LAI at the initial state (or the first time slice) need to be given first when we use the dynamic model to predict the tendency of LAI dynamics. We fit the coefficients using MODIS LAI data with QC flags of zero from Julian Day 89–265, 2004, together with calculated GSI. The initial LAI at the first time slice are derived from MOD15A2 data on Julian Day 89 in the five year period 2004–2008. If MOD15A2 values from 2008 were retrieved using the LUT algorithm, we set these values as the initial LAI; if not, the initial LAI values were set as the average of the MODIS LAI values in the period 2004–2007, assuming that LAI values had been retrieved using the LUT method in this four-year period. For pixels where no LAI values were estimated using the main algorithm during the five years, LAI values retrieved from the back-up algorithm were used to determine the initial LAI. The initial LAI values were the weighted average of MODIS LAI data in 2008 and historical LAI values from 2004 to 2007. The weight was the reciprocal of the standard error provided by the MODIS LAI product in the relevant year. Moreover, if the MODIS LAI values of a pixel were filled values, i.e., retrieved using neither the main algorithm nor the back-up algorithm, then the initial LAI of the pixel was set to the MODIS LAI values of the pixel's neighbors that are classified as the same land cover type. If the MODIS LAI value of a pixel and its surrounding pixels were all filled values, indicating the land cover types to be barren, sparse vegetation, inland freshwater, etc., we could assign a small initial value (0.1 in this study) to the pixel. 3.2. Generating a training dataset SAILH (scattering by arbitrarily inclined leaves, with hot spot effect) is a widely used model for simulating the bidirectional reflectance factor of turbid medium plant canopies, due to its ease of use, general robustness, and consistent performance in validation experiments (Jacquemoud et al., 2009; Kuusk, 1991; Verhoef, 1984; Verhoef et al., 2007). Input variables of SAILH are shown in Tables 2 and 3.

34

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

Input parameters in Table 3 were configured according to field measurements in the study area and a normal distribution was assumed for leaf reflectance, leaf transmittance, and soil reflectance. The mean values and standard deviations (std.) were obtained. An array of 200 values for each parameter was randomly generated according to their distributions. LAI values were discretized across the range 0.1–6.0 at an interval of 0.1 in the simulation. In the forward model, SAILH can compute the spectral reflectance for a certain set of input parameters from Tables 2 and 3 and thus generate LUTs or a training dataset.

Table 3 Measured leaf and soil spectral properties for crops and grassland in the red and NIR bands equivalent to MODIS bands 1 and 2. Parameter

Crops (mean ± std.)

Grass (mean ± std.)

Red-leaf reflectance NIR-leaf reflectance Red-leaf transmittance NIR-leaf transmittance Red-soil reflectance NIR-soil reflectance Average leaf angle Hot spot parameter

0.08 ± 0.02 0.42 ± 0.03 0.03 ± 0.01 0.50 ± 0.03 0.16 ± 0.05 0.20 ± 0.04 57 0.25

0.10 ± 0.01 0.48 ± 0.02 0.03 ± 0.01 0.50 ± 0.03 0.16 ± 0.05 0.20 ± 0.04 63 0.25

3.3. Learning parameters of a DBN One needs to specify two things to describe a Bayesian network or a DBN: the graph structure and the parameters of each conditional probability distribution (CPD). The model structure can be determined according to the system logic or causality (Fig. 3). Nodes in the network represent random variables and the arcs represent the conditional dependence relationship between variables (Murphy, 1998, 2002). In Fig. 3, SR, LR, LT, SZN, VZN, RA, red, and NIR represent soil reflectance, leaf reflectance, leaf transmittance, solar zenith angle, view zenith angle, relative azimuth angle, red band reflectance, and near-infrared reflectance, respectively. SR, LR, LT, LAI, SZN, VZN, and RA in the structure determine canopy reflectance, which is identical to the SAILH model in Section 3.2. Another indispensable element of a Bayesian network is that model parameters can be derived once the model structure is built. We aim to obtain a state transition matrix and likelihood probabilities in Fig. 2 through parameter learning of a Bayesian network. Since the explicit expression for the dynamic process model is derived in Section 3.1, we can calculate the state transition matrix easily. Likelihood probability is calculated from the training dataset generated by the SAILH model using the maximum likelihood estimation method (Murphy, 1998). We discretize the reflectance and LAI at intervals of 0.01 and 0.2, respectively, in order to represent CPD as a conditional probability table (CPT), which lists the probability that the child node assumes each of its different values for each combination of values of its parents. If there are latent variables in the DBN or several missing data in the training dataset, a more complex algorithm must be considered, such as the expectation maximization algorithm or gradient descent, to determine this parameter (Heckerman, 2008). 3.4. Performing inference in a DBN to estimate LAI Estimating land surface parameters is essentially performing inference in a DBN to compute a time series posterior probability distribution of target parameters. Two main kinds of inference algorithms are usually performed in the state estimation — filtering inference and smoothing inference. Filtering inference is an online analysis, where observations arrive in real time, and smoothing inference is an offline analysis, where all observations should be collected before it can be triggered (Murphy, 2002; Wikle & Berliner, 2007). In this study, we performed filtering inference to recursively estimate the probability states of LAI time series. Because we already knew the LAI dynamics, the time slice, and the observed reflectance, our aim is to estimate LAI by performing an inference algorithm in a DBN. The core of the filtering algorithm is

illustrated by Eq. (12). More details about the derivation of the algorithm can be found in a previous study (Qu et al., 2012). PðLAIT jRef 1:T Þ ¼ ¼

PðRef T jLAI T Þ  PðLAIT jRef 1:T−1 Þ PðRef T jRef 1:T−1 Þ

PðRef T jLAI T Þ  ∑LAIT−1 PðLAIT jLAI T−1 Þ  PðLAI T−1 jRef 1:T−1 Þ ∑LAIT P ðRef T jLAI T Þ  P ðLAI T jRef 1:T−1 Þ ð12Þ

where LAIT is LAI at time T, RefT is reflectance at time T, Ref1:T are reflectance values until time T, and P(RefT|LAIT) is the likelihood term or the quantitative expression of the conditional relationship between LAI and reflectance, which is obtained through parameter learning, as described in Section 3.3. P(LAIT|Ref1:T−1) is the predicted term. It is calculated from the posterior probability distribution of LAI at the last time slice P(LAIT−1|Ref1:T−1) and the state transition matrix P(LAIT|LAIT−1), which describes the evolution of the dynamic system. P(RefT|Ref1:T−1) is the normalization term. If there were no errors in the LAI predicted by the dynamic process model and remote sensing observations, state transition probability could easily be derived from the dynamic process model, and likelihood probabilities could be obtained by parameter learning in a Bayesian network. However, errors exist in both the models and observations and they should be considered when inferring the state of a target parameter. In the present study, we assumed that the model and observation errors both obey Gaussian distributions. For a certain time slice, the state transition matrix was calculated from the dynamic process model and its associated uncertainties were calculated using a normal cumulative distribution function, and the likelihood probability was the product of the probabilities computed by the normal cumulative distribution function from the observed reflectance (and its associated uncertainty) and the corresponding probabilities obtained from the CPT.

LAI SR

SZN VZN

LR

RA

LT Table 2 Inputs of solar and view angles in the SAILH model. Parameter

Range (degree)

Step (°)

Solar zenith angle Viewing zenith angle Relative azimuth angle

20–40 0–60 135

10 10 –

Red

NIR

Fig. 3. Structure of a Bayesian network.

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

According to Eq. (12), the estimation algorithm is an iterative process. During each iteration, the following five steps are accomplished in sequence:

35

90m

30m

• prediction of LAI at the current time slice using the dynamic process model on the basis of the posterior distribution of LAI at the previous time slice; • computation of the state transition matrix P(LAIT |LAIT − 1) from the predicted LAI and its associated uncertainties using the normal cumulative distribution function; • calculation of the probabilities of observed reflectance located in each discrete interval according to the associated uncertainty; • multiplication of this probability matrix by CPT to obtain the likelihood term P(RefT|LAIT); • computation of the posterior probability distribution of LAI using Eq. (12). In summary, the algorithm first predicts LAI using the dynamic process model with a one-dimensional integration from the previous state and then updates the integration with new observational data. The updated posterior probability distribution of LAI is used as the input for the next time slice. Through dynamic iteration, a time series of the posterior probability distributions of LAI will be obtained. 3.5. Upscaling field measured LAI to the spatial resolution of 1 km Comparison of satellite biophysical parameters with those measured in situ is a robust method for validation (Justice et al., 2000); however, given the scale mismatch between ground in situ measurements and the coarse resolution of large-angle remote sensors, direct comparison of field measurements and results estimated using remote sensing data is still questionable (Hufkens et al., 2008; Liang et al., 2002). Therefore, ground site measurements should be scaled up to MODIS resolutions using intermediate scale data, such as a high-resolution remotely sensed imagery. At the study area, LAI were collected in agricultural crops from May 2008 to July 2008 using an optical instrument LAI-2000, TRAC and a fisheye camera, and by manual measurements. The sampling design is shown in Fig. 4. Each sampling plot covers an area of 90 × 90 m2, with nine subplots identified within it, and is relatively homogeneous. LAI were measured randomly in each subplot. The unique LAI value for each plot was produced by averaging these measurements. All plots were geolocated using a real-time differential GPS in both the x and y directions. Ground point measurements were scaled up to the MODIS resolution via the following three stages: i) development of an empirical model between the field measurements and NDVI derived from high-resolution imagery; ii) derivation of LAI maps using the empirical model; iii) aggregation of high-resolution LAI maps to match coarse-resolution products. More details about the above work scheme can be found in the previous work of our research group members (Liu et al., submitted for publication). 4. Results 4.1. Performance of the dynamic process model According to the derivation in Section 3.1, the undetermined coefficients a and b are closely related to plant growth, and c, d, and e are functions of the maximum and minimum LAI at a given growth stage. These undetermined coefficients have some physical implications and can be relatively stable for a certain area, but differ for different land cover types. Therefore, the dynamic process model may be extended to predict LAI dynamics in the coming years once it has been fitted using historical data. To prove this point, we fitted the coefficients using a dataset from 2004 for crops and grassland, and then predicted LAI dynamics in

Fig. 4. Sampling strategy.

2005 to evaluate performances of the dynamic process model. Only MOD15A2 data with a QC flag of zero were used for the fitting of model coefficients in 2004 and evaluation of model performance in 2005, so initial LAI values could be directly assigned the MODIS LAI according to Section 3.1. Undetermined coefficients of the dynamic process model were obtained from MODIS LAI data from Julian Day 89–265, 2004, with QC flags of zero and corresponding GSI by performing regression analysis. MOD15A2 were first smoothed using a Savitzky–Golay (SG) filter to avoid occasional temporal discontinuities (Chen et al., 2004). The comparison of high-quality (QC = 0) MODIS LAI data in 2005 and model predicted LAI is shown in Fig. 5, where we can see that the determination coefficient between predicted LAI and MODIS LAI is 0.85 for crops and 0.84 for grassland, with root mean square errors (RMSE) of 0.22 and 0.18, respectively. Moreover, when we compared the predicted LAI with high-quality MODIS LAI data, we found that the modeled LAI also exhibited better time continuity than the MODIS LAI. Therefore, the dynamic process model can be used to predict ground dynamic processes.

4.2. Uncertainty analysis for the dynamic process model and remote sensing observations Uncertainty analysis is required for both the predicted LAI values and the observed reflectance in order to compute the state transition and likelihood probabilities during the inference. As the posterior probability distribution of LAI is the product of these two probabilities, their uncertainties will have an effect on the reliability of the probability distribution of LAI. Normally, the posterior probability distribution of LAI tends to be closer to the distribution with the smaller deviation. In other words, if the predicted LAI values or observed reflectance had large errors, they would contribute less to the posterior probability distribution of LAI. In this study, uncertainties were assigned according to performance. For the dynamic model, errors were evident in the first few slices, and so a large uncertainty value was assigned to the predicted LAI (a normal distribution with a standard deviation of 0.4 m 2/m 2 was reasonable because little vegetation grew in these time slices). Thus, the posterior probability distribution of LAI was insensitive to

36

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

Fig. 5. Comparison of the predicted LAI using a dynamic process model with high-quality MODIS LAI data in the study area in 2005.

the predicted LAI values from the dynamic process model and was mostly affected by the observed reflectance. After the first few slices, the predicted LAI values tended to be stable with small errors and the model was more accurate and yielded good time continuity in its final results. Similarly, we consider the observed reflectance uncertainties. The allowable discrepancies between simulated and MODIS surface reflectance are 20% for the red band and 5% for the NIR band for herbaceous vegetation in the MODIS LAI algorithm (Collection 5). In this study, the uncertainties associated with remote sensing observations were set as 25% for the red band and 10% for the NIR, as under these proportions, the observed MODIS surface reflectance and simulated reflectance in LUTs are well matched in the red–NIR space. 4.3. Estimated regional LAI maps The posterior probability distributions of LAI at the grow season were retrieved by performing filtering inference in a DBN. The estimated LAI in the study were the mean values of the posterior probability distribution of LAI. Results are shown in Fig. 6. LAI maps in Fig. 6 were from Julian Day 89 (upper left) to Julian Day 265 (lower right), 2008. Very few plants were growing in our study area before Julian Day 89, so we chose LAI values on this day for the initial states. The vegetation was nearly mature by Julian Day 265. However, the land cover type may change after Julian Day 265. We can clearly see the spatial distribution of LAI and the temporal trajectories during the grow season from Fig. 6.

in Fig. 7, the time series curves of MODIS LAI are inconsistent with the growth laws of vegetation LAI dynamics. Sudden falls exist on Julian Day 193 in Fig. 7(b) and (c). From the original MODIS quality file, we found that this was caused by the poor quality of remote sensing observations. According to the QC flags provided by the MODIS LAI product, significant clouds were present for the pixel in Fig. 7(b), and a mixed cloud was present for the pixel in Fig. 7(c), on Julian Day 193. Although significant clouds were present for the pixel in Fig. 7(d) on Julian Day 153, there was no sudden fall at this point. Moreover, for the pixel in Fig. 7(a), the QC flags of MODIS LAI were zero from Julian Day 89 to 265, indicating that the quality of the MODIS LAI was high for this pixel, but the curve fluctuated frequently. These phenomena may be explained by several factors, such as the inversion algorithm or ancillary information for LAI inversion in a certain land cover type being unsuitable for the actual conditions. However, by introducing the LAI dynamic information derived from ground station data, the DBN LAI was temporally continuous and the trajectories were reasonable. This shows that the quality of the estimated LAI improved because of accurate ground station data.

The performance of the time series DBN LAI was evaluated from three perspectives. The first was to evaluate the temporal trajectories of the DBN LAI; the second was the direct validation using the upscaled field measured LAI, which can be seen as relative true values; the final one was the comparison of frequency histogram distributions between DBN and MODIS LAI time series.

4.3.2. Pixel-by-pixel comparison At the study site, in situ measurements provide a reference to evaluate the performance of DBN LAI. Fig. 8 shows scatter-plots for the pixel-by-pixel comparisons of LAI results from the different estimation methods (DBN LAI and MODIS LAI) and the scaled up LAI values. The determination coefficient between the DBN LAI and the field measured LAI is 0.76, with a root mean square error of 0.78, and the two coefficients are 0.45 and 1.00 for the MODIS LAI. It is found that the DBN LAI is more closely correlated with the measured LAI than the MODIS LAI. This clearly shows that the incorporation of ground station data can improve the quality of the estimated LAI. Moreover, it should be noted that when the inversion accuracy dropped slightly due to ground LAI values becoming larger than 3.0 m 2/m 2, our proposed method could still produce slightly better consistency with results from ground LAI.

4.3.1. Temporal trajectories In the study area, most MODIS LAI values are temporally discontinuous, which is exhibited as the phenomenon of intermittent and fluctuating time series curves or sudden falls within the grow season. We selected several pixels to make a comparison between estimated LAI and MODIS LAI (Fig. 7). According to the quality flag of the MODIS LAI, these pixels were retrieved using the main algorithm. As shown

4.3.3. Frequency histogram comparison Field measurements are available for a limited number of pixels, and we carried out a pixel-by-pixel validation for these in the previous section. As a complete comparison should involve all pixels, frequency histograms were used to analyze extremely large datasets. Fig. 9 shows the frequency histograms of DBN LAI values and MODIS LAI values from Julian Day 89 to 265. The horizontal axis shows LAI values and the vertical

4.4. Evaluation of estimated LAI

Julian Day 97

Julian Day 105

Julian Day 113

Julian Day 121

Julian Day 129

Julian Day 137

Julian Day 145

Julian Day 153

Julian Day 161

Julian Day 169

Julian Day 177

Julian Day 185

Julian Day 193

Julian Day 201

Julian Day 209

Julian Day 217

J u li a n D ay 2 2 5

Julian Day 233

Julian Day 241

Julian Day 249

Julian Day 257

Julian Day 265

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

Julian Day 89

Fig. 6. LAI estimated by performing filtering inference in a DBN at the study site. The gray parts represent non-crop and non-grass land cover types.

37

38

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

(a)

(b)

3.50

2.00 1.50 1.00

(c)

180 Julian Day

230

1.00

280

(d)

3.00

80

130

180 Julian Day

230

280

80

130

180 Julian Day

230

280

3.00 2.50

2.50

LAI (m 2 /m 2 )

LAI (m 2 /m 2 )

1.50

0.00 130

3.50

2.00 1.50 1.00

2.00 1.50 1.00 0.50

0.50 0.00 80

2.00

0.50

0.50 0.00 80

3.00 2.50

2.50

LAI (m 2 /m 2 )

LAI (m 2 /m 2 )

3.00

0.00 130

MODIS LAI

180 Julian Day DBN LAI

230

280

Dynamic model predicted LAI (update)

Dynamic model predicted LAI

Fig. 7. Comparison of the trajectories of MODIS LAI, DBN LAI and dynamic model predicted LAI. Green line represents dynamic model predicted LAI (update), which means LAI predicted on the basis of LAI at the previous time slice and ground meteorological data. Purple line represents dynamic model predicted LAI on the basis of ground meteorological data when no satellite reflectance were observed.

axis shows the percentage of pixels corresponding to that LAI value. As shown in Fig. 9, the percentage of pixels with MODIS LAI values of less than 0.4 was larger than 50% at the first eight histograms, and this proportion was still larger than 30% in the next 15 histograms. The number of pixels with MODIS LAI values larger than 2.0 was very small in all of the histograms, and there was not even an exception in the rapid growth stage of crops and grass, which indicated an underestimation of MODIS LAI in the pixel-by-pixel comparison. However, the number of pixels with DBN LAI values of less than 0.4 was rather small, except in the first few histograms when crops and grass had not yet begun to grow. Distributions of DBN LAI values tended to shift toward the bigger LAI values. This is in accordance with the actual situation. In addition,

distributions of DBN LAI values were disperse, which is confirmed by the reported results of the LAI spatial heterogeneity in this area (Hu & Zhang, 2005; Zhang et al., 2008). 5. Discussion 5.1. Performance of the prognostic phenology model In this study, the quality of estimated LAI has been greatly improved by combining ground meteorological data and satellite observations. If there were no observed satellite reflectance data, the estimated LAI would solely depend on ground meteorological factors.

Fig. 8. Pixel-by-pixel comparison of field measured LAI and estimated LAI (DBN LAI and MODIS LAI).

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

This could lead to poor performance of the dynamic process model. As shown in Fig. 7, LAI values retrieved from meteorological factors (Dynamic model predicted LAI in Fig. 7) were steady and not reasonable at the grow season. In addition, LAI values were very small, and the spatial heterogeneity of LAI was not obvious at the start of vegetation grow season. Hence, the dynamic process model may produce similar LAI predictions along their temporal evolution without the update of satellite observations since meteorological data for adjacent pixels are similar, even for pixels with obvious spatial heterogeneity. This leads to inaccurate predictions. Moreover, the time at which the model is started greatly affects the predicted results, and so predictions begun on different days may produce vastly different results. Fortunately, remote sensing observations can rectify predicted LAI results from the dynamic process model and gradually give appropriate LAI values in the DBN framework. If we predict LAI iteratively based on ground meteorological data and DBN LAI at the previous time slice, as shown in Fig. 7, reasonable LAI dynamics could be produced. 5.2. Impact of meteorology data on vegetation growth The prognostic phenology model used in this paper was driven by a set of meteorological variables. The foundation of the mathematical model is the theory or concept that climate and terrestrial ecosystems are closely coupled and climate factors may account for much of the biochemical processes of plants (Jolly et al., 2005; Sarkar & Kafatos, 2004). We used daily minimum temperature, vapor pressure deficit, and sunshine hours to express the relationships with the seasonal variation of vegetation. The impacts of the three climate factors on the ecosystem and the vegetation pattern of different ecosystems have been confirmed by several researchers (Hu et al., 2011; Xiong et al., 2007). In addition to the meteorological factors mentioned above, some other climate factors were considered to significantly affect the growth and spatial patterns of land surface vegetation, such as rainfall and evapotranspiration in the growing season (Buyantuyev & Wu, 2012; Zhao et al., 2011). Moreover, some researchers have pointed out that vegetation responds to climate variability in complex ways and the relationship between meteorological parameters and vegetation growth may vary through the phenological cycle of plants (Horion et al., 2012; Weiss et al., 2004). Therefore, when used for larger scale estimation, we should carefully identify the limiting factors. Maybe a more complicated mathematical model should be developed to fit varying situations. It is an interesting topic worthy of further investigation. 5.3. Impact of MODIS LAI on the dynamic process model In addition to ground meteorological observations, MODIS LAI was another critical component of the dynamic process model. The dynamic process model was fitted using historical MODIS LAI data; hence, the prediction results were limited to the original MODIS LAI values to some extent. This phenomenon is evident in the first few time slices. The underestimated MODIS LAI values may partly account for the relatively low accuracy of estimated results for larger LAI field measurements. Because the MODIS LAI product can be downloaded easily, our study provided a feasible and effective approach for combining ground station observations and satellite observations to estimate land surface parameters. Although many temporal smoothing and gap filling techniques have been developed to extract phenological information, they lack the effective predictive capabilities (Alhamad et al., 2007). In a recent study, the Seasonal Auto-Regressive Integrated Moving Average (SARIMA) method was used to model the MODIS LAI time series and forecast future LAI values (Jiang et al., 2010; Xiao et al., 2011). It is assumed that all quantities governing the model's behavior remain constant with time; this assumption is not in agreement with the actual situation. Nevertheless, the LAI dynamics extracted from the MODIS LAI product will be limited to the performance of original MODIS LAI values.

39

An alternative approach would be to provide a first estimate of the LAI with the LUTs built with a radiative transfer model; however, the observations must be available in advance. In addition, satellite observations should be dealt with properly or data assimilation should be adopted to overcome the cloud and aerosol related deficiencies of remote sensing observations (Liu et al., 2008; Stockli et al., 2008).

5.4. Strength and limitations of the proposed approach In this study, we have improved the quality of LAI estimations by combining remote sensing data and ground meteorological station data. Both data used can be obtained easily and freely. The proposed approach was easy to implement and feasible for land surface parameter retrieval. However, there are several limitations. Climate factors used in the study cannot cover all the situations, so limiting meteorological factors need to be carefully investigated and delicate prognostic model need to be developed when applied to some regions or a larger scale. In addition, performances of this approach may not be satisfactory when used to estimate LAI with a high spatial resolution, because climate factors were not the main limiting factors in a fine scale.

6. Conclusion Both satellite observations and ground observations are key components in the study of large-scale land surface. In this study, we combined both observation types and explored their advantages to estimate LAI on regional scales. Validation results showed that estimated LAI values were close to field measurements and also had reasonable spatial distributions and time trajectories. Our study shows the potential use and performance of a DBN in the estimation of land surface parameters on a regional scale. More importantly, it provides a basis and an effective approach for incorporating ground station network data in the estimation of land surface parameters. When applied to a larger scale, several points need to be considered. Firstly, the dynamic process model should be further developed to fit complex situations by the assimilation of remote sensing observations and climate data. Secondly, inputs of the physical models, such as leaf and canopy structure parameters can be assigned a valid range for each land cover type. Measures can be taken to ensure that model outputs match the observed reflectance (Tian et al., 2000; Wang et al., 2004). There exist other ways to implement the prognostic phenology model. In this study, model coefficients were fixed after fitted with historical climate factors and MODIS LAI data. One feasible implementation in future research is to give an initial value of these coefficients and optimize them during the LAI estimates. Another feasible way is to use FPAR directly to calculate deltaGSI. It will reduce the number of undetermined coefficients. Besides, FPAR is functionally related with the amount of photosynthetic biomass, and the foliage biomass is then functionally associated with vegetation LAI and SLA (specific leaf area, represented by the ratio of dry matter to leaf area). So using FPAR to calculate deltaGSI may further help to predict vegetation LAI accurately. It is a further research topic. Pixel-by-pixel comparisons of field measured and estimated LAI were made mainly on Julian Day 169 and 193, 2008, and the accuracy was assessed over a set of locations. Estimated results could not be evaluated directly for the pixels without field measurements, and therefore, frequency histograms and time trajectories were used in this work. More ground experiments should be conducted in future studies. Furthermore, new experiments should focus on the seasonality of the vegetation structure, which is of great importance for future applications.

40

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

Acknowledgements This work was partly supported by the National Natural Science Foundation of China (41271348,91125004) and the National

DOY=89

70

100

50

Percentage (%)

40 30 20 10 0

90

80

80

70

70 60 50 40 30

10

80

40 30

10

10

25

25

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

LAI(m2/m2)

DOY=145

20

Percentage (%)

Percentage (%)

15 10

0

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

15 10 5

5

5

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

LAI(m2/m2)

LAI(m2/m2)

DOY=161

18

LAI(m2/m2)

DOY=169

16

16

16

14

14

Percentage (%)

12 10 8 6

12 10 8 6

12 10 8 6

4

4

4

2

2

2

0

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7

LAI(m2/m2)

DOY=177

18

18 14

DOY=153

25

20

10

20

LAI(m2/m2)

DOY=137

15

30

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

LAI(m2/m2)

20

40

10

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

DOY=129

60 50

50

20

0

Percentage (%)

DOY=121

60

20

20

LAI(m2/m2)

Percentage (%)

Percentage (%)

Percentage (%)

30

30

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

70

70

Percentage (%)

0

LAI(m2/m2)

80

40

30

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

DOY=113

50

40

20

LAI(m2/m2)

60

50

10 0

90

60

20

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8

DOY=105

90

Percentage (%)

Percentage (%)

60

DOY=97

Percentage (%)

(a)

High Technology Research and Development Program 863 (2009AA122101). The authors would like to thank the anonymous reviewers and the editor for their constructive remarks on an earlier version of the manuscript.

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7

LAI(m2/m2)

LAI(m2/m2)

Fig. 9. Frequency histograms showing the distribution of MODIS LAI and DBN LAI values, from JD 89–265.

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

DOY=185

16

14

14

12

12 10 8 6

DOY=193

16

Percentage (%)

18 16

10 8 6

4

4

2

2

0

0

8 6 2 0

16

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7

LAI(m2/m2)

DOY=217

14

12

12

8 6

Percentage (%)

14

12 10

10 8 6

10 8 6

4

4

4

2

2

2

18

LAI(m2/m2)

DOY=233

18

14

Percentage (%)

16

14 12 10 8 6

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4

16

8 6 4 2

LAI(m2/m2)

10

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1

15

5

0

0

DOY=249

25 20

10

2

LAI(m2/m2)

DOY=241

12

4

30

0

0

LAI(m2/m2)

DOY=225

16

14

Percentage (%)

Percentage (%)

DOY=209

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7

Percentage (%)

10

LAI(m2/m2)

0

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1

LAI(m2/m2)

LAI(m2/m2)

DOY=257

30

DOY=265

25

Percentage (%)

25

Percentage (%)

12

4

LAI(m2/m2)

16

14

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7

DOY=201

20

Percentage (%)

Percentage (%)

18

Percentage (%)

(b)

41

20 15 10

20 15

DBN LAI MODIS LAI

10 5

5 0

0 0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4

LAI(m2/m2)

0.1 0.4 0.7 1.0 1.3 1.6 1.9 2.2 2.5 2.8 3.1

LAI(m2/m2) Fig. 9 (continued).

42

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43

References Alhamad, M. N., Stuth, J., & Vannucci, M. (2007). Biophysical modelling and NDVI time series to project near‐term forage supply: spectral analysis aided by wavelet denoising and ARIMA modelling. International Journal of Remote Sensing, 28, 2513–2548. Barillec, R., & Cornford, D. (2009). Data assimilation for precipitation nowcasting using Bayesian inference. Advances in Water Resources, 32, 1050–1065. Buck, A. L. (1981). New equations for computing vapor pressure and enhancement factor. Journal of Applied Meteorology, 20, 1527–1532. Buyantuyev, A., & Wu, J. (2012). Urbanization diversifies land surface phenology in arid environments: Interactions among vegetation, climatic variation, and land use pattern in the Phoenix metropolitan region, USA. Landscape and Urban Planning, 105, 149–159. Campbell, G., & Norman, J. (1998). Introduction to environmental biophysics (2nd ed.). : Springer Verlag. Challa, S., & Koks, D. (2004). Bayesian and Dempster-Shafer fusion. Sadhana, 29, 145–174. Chen, J. M., & Black, T. A. (1992). Defining leaf area index for non-flat leaves. Plant, Cell & Environment, 15, 421–429. Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., & Eklundh, L. (2004). A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. Remote Sensing of Environment, 91, 332–344. Chen, J. M., Plummer, P. S., Rich, M., Gower, S. T., & Norman, J. M. (1997). Leaf area index of boreal forests: Theory, techniques, and measurements. Journal of Geophysical Research, 102, 29,429–429,443. De Kauwe, M. G., Disney, M. I., Quaife, T., Lewis, P., & Williams, M. (2011). An assessment of the MODIS collection 5 leaf area index product for a region of mixed coniferous forest. Remote Sensing of Environment, 115, 767–780. Dente, L., Satalino, G., Mattia, F., & Rinaldi, M. (2008). Assimilation of leaf area index derived from ASAR and MERIS data into CERES-Wheat model to map wheat yield. Remote Sensing of Environment, 112, 1395–1407. Douaoui, A. E. K., Nicolas, H., & Walter, C. (2006). Detecting salinity hazards within a semiarid context by means of combining soil and remote-sensing data. Geoderma, 134, 217–230. Duveiller, G., Weiss, M., Baret, F., & Defourny, P. (2011). Retrieving wheat Green Area Index during the growing season from optical time series measurements based on neural network radiative transfer inversion. Remote Sensing of Environment, 115, 887–896. Fang, H., & Liang, S. (2005). A hybrid inversion method for mapping leaf area index from MODIS data: experiments and application to broadleaf and needleleaf canopies. Remote Sensing of Environment, 94, 405–424. Fang, H., Liang, S., Townshend, J. R., & Dickinson, R. E. (2008). Spatially and temporally continuous LAI data sets based on an integrated filtering method: Examples from North America. Remote Sensing of Environment, 112, 75–93. Garrigues, S., Lacaze, R., Baret, F., Morisette, J. T., Weiss, M., Nickeson, J. E., et al. (2008). Validation and intercomparison of global Leaf Area Index products derived from remote sensing data. Journal of Geophysical Research, 113, G02028. Han, X., & Li, X. (2008). An evaluation of the nonlinear/non-Gaussian filters for the sequential data assimilation. Remote Sensing of Environment, 112, 1434–1449. Hasegawa, K., Matsuyama, H., Tsuzuki, H., & Sweda, T. (2010). Improving the estimation of leaf area index by using remotely sensed NDVI with BRDF signatures. Remote Sensing of Environment, 114, 514–519. Heckerman, D. (2008). A tutorial on learning with Bayesian networks. Innovations in Bayesian Networks (pp. 33–82). Hill, M. J., Senarath, U., Lee, A., Zeppel, M., Nightingale, J. M., Williams, R. J., et al. (2006). Assessment of the MODIS LAI product for Australian ecosystems. Remote Sensing of Environment, 101, 495–518. Horion, S., Cornet, Y., Erpicum, M., & Tychon, B. (2012). Studying interactions between climate variability and vegetation dynamic using a phenology based approach. International Journal of Applied Earth Observation and Geoinformation, http: //dx.doi.org/10.1016/j.jag.2011.12.010. Houborg, R., Anderson, M., & Daughtry, C. (2009). Utility of an image-based canopy reflectance modeling tool for remote estimation of LAI and leaf chlorophyll content at the field scale. Remote Sensing of Environment, 113, 259–274. Hu, M. Q., Mao, F., Sun, H., & Hou, Y. Y. (2011). Study of normalized difference vegetation index variation and its correlation with climate factors in the three-river-source region. International Journal of Applied Earth Observation and Geoinformation, 13, 24–33. Hu, S., & Zhang, W. (2005). A quality assessment of MODIS LAI product in Heihe and Hanjiang River Basins. Remote Sensing Information, 4, 22–27. Hufkens, K., Bogaert, J., Dong, Q. H., Lu, L., Huang, C. L., Ma, M. G., et al. (2008). Impacts and uncertainties of upscaling of remote-sensing data validation for a semi-arid woodland. Journal of Arid Environments, 72, 1490–1505. Hutchinson, M. F. (2001). ANUSPLIN Version 4.2. User Guide. The Australian National University. Canberra: Centre for Resource and Environmental Studies In. Ihler, A. T., Kirshner, S., Ghil, M., Robertson, A. W., & Smyth, P. (2007). Graphical models for statistical inference and data assimilation. Physica D: Nonlinear Phenomena, 230, 72–87. Jacquemoud, S., Verhoef, W., Baret, F., Bacour, C., Zarco-Tejada, P. J., Asner, G. P., et al. (2009). PROSPECT + SAIL models: A review of use for vegetation characterization. Remote Sensing of Environment, 113, S56–S66. Jiang, B., Liang, S., Wang, J., & Xiao, Z. (2010). Modeling MODIS LAI time series using three statistical methods. Remote Sensing of Environment, 114, 1432–1444. Jolly, W., Nemani, R., & Running, S. (2005). A generalized, bioclimatic index to predict foliar phenology in response to climate. Global Change Biology, 11, 619–632. Justice, C., Belward, A., Morisette, J., Lewis, P., Privette, J., & Baret, F. (2000). Developments in the 'validation' of satellite sensor products for the study of the land surface. International Journal of Remote Sensing, 21, 3383–3390.

Kalacska, M., Sanchez-Azofeifa, G. A., Caelli, T., Rivard, B., & Boerlage, B. (2005). Estimating leaf area index from satellite imagery using Bayesian networks. IEEE Transactions on Geoscience and Remote Sensing, 43, 1866–1873. Knyazikhin, Y., Glassy, J., Privette, J. L., Tian, Y., Lotsch, A., Zhang, Y., et al. (1999). MODIS leaf area index (LAI) and fraction of photosynthetically active radiation absorbed by vegetation (FPAR) product (MOD15) algorithm theoretical basis document. Theoretical Basis Document (pp. 20771). Greenbelt, MD: NASA Goddard Space Flight Center. Koetz, B., Baret, F., Poilvé, H., & Hill, J. (2005). Use of coupled canopy structure dynamic and radiative transfer models to estimate biophysical canopy characteristics. Remote Sensing of Environment, 95, 115–124. Kuusk, A. (1991). The hot spot effect in plant canopy reflectance. Photon-Vegetation Interactions (pp. 139–159). le Maire, G., Marsden, C., Verhoef, W., Ponzoni, F. J., Lo Seen, D., Bégué, A., et al. (2011). Leaf area index estimation with MODIS reflectance time series and model inversion during full rotations of Eucalyptus plantations. Remote Sensing of Environment, 115, 586–599. Li, X., Li, X., Li, Z., Ma, M., Wang, J., Xiao, Q., et al. (2009). Watershed Allied Telemetry Experimental Research. Journal of Geophysical Research, 114, D22103. Liang, S. (2004). Quantitative remote sensing of land surfaces. New York: John Wiley and Sons, Inc. Liang, S., Fang, H., Chen, M., Shuey, C. J., Walthall, C., Daughtry, C., et al. (2002). Validating MODIS land surface reflectance and albedo products: methods and preliminary results. Remote Sensing of Environment, 83, 149–162. Liu, Q., Gu, L., Dickinson, R. E., Tian, Y., Zhou, L., & Post, W. M. (2008). Assimilation of satellite reflectance data into a dynamical leaf model to infer seasonally varying leaf areas for climate and carbon models. Journal of Geophysical Research, 113, D19113. Liu, Y., Wang, J., Zhou, H., Xiao, Z., Xue, H., Zhang, Y., et al. (submitted for publication). An upscaling approach for the validation of coarse-resolution LAI products using ground measurements. International Journal of Remote Sensing. Moradkhani, H. (2008). Hydrologic Remote Sensing and Land Surface Data Assimilation. Sensors, 8, 2986–3004. Moulin, S., Bondeau, A., & Delecolle, R. (1998). Combining agricultural crop models and satellite observations: from field to regional scales. International Journal of Remote Sensing, 19, 1021–1036. Murphy, K. (1998). A brief introduction to graphical models and Bayesian networks. Available on-line at http://www.cs.ubc.ca/murphyk/Bayes/bnintro.html. Murphy, K. P. (2002). Dynamic Bayesian networks: representation, inference and learning. : University of California. Myneni, R. B., Hoffman, S., Knyazikhin, Y., Privette, J. L., Glassy, J., Tian, Y., et al. (2002). Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment, 83, 214–231. Pipunic, R. C., Walker, J. P., & Western, A. (2008). Assimilation of remotely sensed data for improved latent and sensible heat flux prediction: A comparative synthetic study. Remote Sensing of Environment, 112, 1295–1305. Qu, Y., Wang, J., Wan, H., Li, X., & Zhou, G. (2008). A Bayesian network algorithm for retrieving the characterization of land surface vegetation. Remote Sensing of Environment, 112, 613–622. Qu, Y., Zhang, Y., & Wang, J. (2012). A dynamic Bayesian network data fusion algorithm for estimating leaf area index using time-series data from in situ measurement to remote sensing observations. International Journal of Remote Sensing, 33, 1106–1125. Quaife, T., Lewis, P., De Kauwe, M., Williams, M., Law, B. E., Disney, M., et al. (2008). Assimilating canopy reflectance data into an ecosystem model with an Ensemble Kalman Filter. Remote Sensing of Environment, 112, 1347–1364. Sarkar, S., & Kafatos, M. (2004). Interannual variability of vegetation over the Indian sub-continent and its relation to the different meteorological parameters. Remote Sensing of Environment, 90, 268–280. Sellers, P. J., Los, S. O., Tucker, C. J., Justice, C. O., Dazlich, D. A., Collatz, G. J., et al. (1996). A revised land surface parameterization (SiB2) for atmospheric GCMs. Part II: The generation of global fields of terrestrial biophysical parameters from satellite data. Journal of Climate, 9, 706–737. Stockli, R., Rutishauser, T., Dragoni, D., O'Keefe, J., Thornton, P. E., Jolly, M., et al. (2008). Remote sensing data assimilation for a prognostic phenology model. Journal of Geophysical Research, 113. Tian, J., & Chen, D. (2010). A semi-empirical model for predicting hourly ground-level fine particulate matter (PM2.5) concentration in southern Ontario from satellite remote sensing and ground-based meteorological measurements. Remote Sensing of Environment, 114, 221–229. Tian, Y., Zhang, Y., Knyazikhin, Y., Myneni, R. B., Glassy, J. M., Dedieu, G., et al. (2000). Prototyping of MODIS LAI and FPAR algorithm with LASUR and LANDSAT data. IEEE Transactions on Geoscience and Remote Sensing, 38, 2387–2401. Verger, A., Baret, F., & Weiss, M. (2011). A multisensor fusion approach to improve LAI time series. Remote sensing of environment, 115, 2460–2470. Verhoef, W. (1984). Light scattering by leaf layers with application to canopy reflectance modeling: the SAIL model. Remote Sensing of Environment, 16, 125–141. Verhoef, W., Jia, L., Xiao, Q., & Su, Z. (2007). Unified optical-thermal four-stream radiative transfer theory for homogeneous vegetation canopies. IEEE Transactions on Geoscience and Remote Sensing, 45, 1808–1822. Wang, F. -m., Huang, J. -f., Tang, Y. -l., & Wang, X. -z. (2007). New Vegetation Index and Its Application in Estimating Leaf Area Index of Rice. Rice Science, 14, 195–203. Wang, D., Wang, J., & Liang, S. (2010). Retrieving crop leaf area index by assimilation of MODIS data into a crop growth model. Science China Earth Sciences, 53, 721–730. Wang, Y., Woodcock, C. E., Buermann, W., Stenberg, P., Voipio, P., Smolander, H., et al. (2004). Evaluation of the MODIS LAI algorithm at a coniferous forest site in Finland. Remote Sensing of Environment, 91, 114–127.

Y. Zhang et al. / Remote Sensing of Environment 127 (2012) 30–43 Weiss, J. L., Gutzler, D. S., Allred Coonrod, J. E., & Dahm, C. N. (2004). Seasonal and inter-annual relationships between vegetation and climate in central New Mexico, USA. Journal of Arid Environments, 57, 507–534. Weiss, M., Troufleau, D., Baret, F., Chauki, H., Prévot, L., Olioso, A., et al. (2001). Coupling canopy functioning and radiative transfer models for remote sensing data assimilation. Agricultural and Forest Meteorology, 108, 113–128. Wikle, C. K., & Berliner, L. M. (2007). A Bayesian tutorial for data assimilation. Physica D: Nonlinear Phenomena, 230, 1–16. Wu, J., Wang, D., & Bauer, M. E. (2007). Assessing broadband vegetation indices and QuickBird data in estimating leaf area index of corn and potato canopies. Field Crops Research, 102, 33–42. Xiao, Z., Liang, S., Wang, J., Jiang, B., & Li, X. (2011). Real-time retrieval of Leaf Area Index from MODIS time series data. Remote Sensing of Environment, 115, 97–106. Xiong, W., Wang, Y., Yu, P., Liu, H., Shi, Z., & Guan, W. (2007). Growth in stem diameter of Larix principis-rupprechtii and its response to meteorological factors in the south of Liupan Mountain, China. Acta Ecologica Sinica, 27, 432–440.

43

Yang, P., Wu, W. -b., Tang, H. -j., Zhou, Q. -b., Zou, J. -q., & Zhang, L. (2007). Mapping Spatial and Temporal Variations of Leaf Area Index for Winter Wheat in North China. Agricultural Sciences in China, 6, 1437–1443. Yao, Y., Zhang, X., & Dduan, Y. (2008). Impacts of Climate Change on Pasture Growth in Subalpine Meadows. Resource Science, 30. Yuan, H., Dai, Y., Xiao, Z., Ji, D., & Shangguan, W. (2011). Reprocessing the MODIS Leaf Area Index products for land surface and climate modelling. Remote Sensing of Environment, 115, 1171–1187. Zhang, W. -C., Zhong, S., & Hu, S. -Y. (2008). Spatial scale transferring study on Leaf Area Index derived from remotely sensed data in the Heihe River Basin, China. Acta Ecologica Sinica, 28. Zhao, X., Tan, K., Zhao, S., & Fang, J. (2011). Changing climate affects vegetation growth in the arid region of the northwestern China. Journal of Arid Environments, 75, 946–952.