Reconstructing LAI Series by Filtering Technique and a Dynamic Plant Model

The accurate measurement of ecosystem biomass is of great importance in scientific, resource management and energy sectors. In particular, biomass is a direct measurement of carbon storage within an ecosystem and of great importance for carbon cycle science and carbon emission mitigation. Remote Sensing is the most accurate tool for global biomass measurements because of the ability to measure large areas. Current biomass estimates are derived primarily from ground-based samples, as compiled and reported in inventories and ecosystem samples. By using remote sensing technologies, we are able to scale up the sample values and supply wall to wall mapping of biomass. Three separate remote sensing technologies are available today to measure ecosystem biomass: passive optical, radar, and lidar. There are many measurement methodologies that range from the application driven to the most technologically cutting-edge. The goal of this book is to address the newest developments in biomass measurements, sensor development, field measurements and modeling. The chapters in this book are separated into five main sections.


Introduction
Remote sensing images provide very rich and useful information linked to relevant biophysical parameters such as the LAI (Leaf Area Index), fCOVER (fraction of vegetation cover) or fAPAR (fraction of Absorbed Photosynthetically Active Radiation). At the moment, several techniques for estimating such variables are available and widely used in many applications, such as estimation of the total biomass and monitoring the dynamics in canopy vegetation (Baret et al., 2007;Lecerf et al., 2008). For several years, a large number of Very High Spatial Resolution (VHSR) satellites, such as Quickbird, Geoeye and Ikonos, have been launched, and very important missions such as the Venus and the Sentinel-2 are expected in 2012 and 2013. This provides possibility of having more or less temporal consistency in VHSR observations of the land use on relevant agricultural sites. However, because of the heterogeneity of the available VHSR data, in particular due to their different wavelengths sensibility and of the intrinsic errors induced by the estimation processes, the resulting time series of biophysical parameters are more or less noisy. As a matter of fact, the estimated variables may only poorly fit their actual dynamics. The estimation of the complete sequence of such parameters is then of prime importance, in particular if one wants to analyze the evolution of the biomass. In this chapter, we propose to explore the possibilities of using tools issued from tracking techniques, in particular particle smoothing, to recover time-consistent series of LAI (Doucet et al., 2001;Kitagawa, 1996) from noisy and incomplete observations. Such techniques, based on Monte-Carlo strategies, allow performing the estimation of an unknown state function, LAI in current case, according to a given dynamical model and to possibly corrupted measurements. The dynamical model on which we rely on is GreenLab model, a functional-structural plant model simulating plant development and growth (Yan et al., 2004). Given model parameters, GreenLab can compute the evolution of LAI, the biomass production and partitioning, the organ size and biomass. Inverse method can be applied to estimate hidden model parameters by fitting model output with measured data (Kang et al., 2008). We suppose that from remote sensing data observing agricultural parcels, the type of

Overall framework
In section 2.1, we first introduce methods of estimating parameters from remote sensing images. In section 2.2, the concepts of GreenLab model are presented. Section 2.3 introduces the filtering technique used to recover time consistent series of LAI.

Estimating type of crop and LAI from remote sensing data
Estimating the type of crop of an observed field from remote sensing data is an old problem that has been widely studied by the computer vision and remote sensing community. This belongs to the classification issue where various families of approaches exist. One can roughly identify two methodologies: pixel-based and region-based classification techniques, see for instance Congalton (1992) and Yan et al. (2006). The first family aims at assigning to each pixel of the image a label corresponding to the nature of the culture, independently to its neighborhood. This performs efficiently when the spatial resolution of the data is low, yielding more or less homogeneity of the pixel reflectance for a given culture. On the other hand, the region-based techniques are useful when dealing with images of high or very high resolution. In this case, the variability of the image luminance inside a given culture is high, and the pixel reflectance is not informative. Texture analysis strategies can be used to characterize and label the different crops, which use 1 st or 2 nd order statistical criteria or more advanced techniques like wavelets (Lefebvre et al., 2010). Many commercial software (e.g. Idrisi, ENVI or eCognition) allow this kind of classification. Estimating the biophysical variables, such as LAI, fCOVER or fAPAR, from satellite observations can provide crucial information for numbers of applications, for instance, monitoring changes in canopy vegetation at global or regional scales, identifying bare soils, or detecting grassland areas. Among the different techniques available, there is a technique based on the inversion of the SAIL+PROPSPECT radiative transfer model (Verhoef, 1984;Jacquemoud and Baret, 1990) using training samples and neural networks, as introduced in (Baret et al., 2007). The SAIL model deals with light scattering by leaf layers with application to canopy reflectance model, and PROSPECT is about leaf optical properties spectra. It has been proved in (Baret et al., 2007;Lecerf et al., 2008) that this approach performs efficiently for low, medium, high and very high-resolution data and is therefore adapted to the variability of satellite images available.

GreenLab model
GreenLab model simulates the two basic processes of plant: development (organogenesis) and growth (organ expansion). In GreenLab, the organogenesis is simulated by an automaton, which gives the dynamics of number, age and type of organs in plant architecture (Yan et al., 2004). The organ expansion is simulated by a source-sink approach. At each time step t, the source function gives the biomass production of a plant, Q t , as a function of plant leaf area S t and environmental factor E t , see Eqn. (1).
In Eqn. (1), S P represents the projection area of an individual plant, which is equivalent to the inverse of planting density d when crop canopy is closed, i.e, S P =1/d. r is a model parameter that can be estimated inversely (Kang et al., 2008;Guo et al., 2006;Dong et al., 2008). In case that E t represents the intercepted light by crop, r means light use efficiency. The ratio S t /S P gives a LAI series. The produced biomass is shared among all growing organs in proportion to their current sink strength, based on common pool hypothesis. For an organ of type O and age j, its increment is biomass in computed as in Eqn. (2).
In GreenLab, each type of organs (e.g. blade, sheath, internode, female organ) has certain relative sink strength P O , which may vary during the expansion of an individual organ, described by an empirical function f O j . Total plant demand D t is sum of sink strength from all growing organs. The organ biomass, which is the accumulation of biomass increment during its life time, is dependent on the ratio between biomass production and demand (Q t /D t ), called source-sink ratio. According to the appearance time of each individual organ given by the automaton, and its increment in biomass since appearance, the biomass of all individual organs in plant can be computed. As the ratio between blade biomass and specific leaf weight, at any time, leaf area S t and consequently LAI, can be computed recursively in GreenLab model.
Where q B t,j is biomass of a j-aged leaf blade, N B t,j is the number of such leaf (being one or zero for maize of single stem structure), λ t is specific leaf weight, T B is expansion duration of a leaf. GreenLab model is a generic model that has been applied to different crops, such as wheat (Kang et al., 2008), maize (Guo et al., 2006), and tomato (Dong et al., 2008;Kang et al., 2010), from which the dynamic growth and development process of crop are rebuilt. In model calibration, the hidden model parameters controlling the source and sink function, such as organ sink strength P O , were estimated inversely from the measured plant data, such as total organ biomass and individual organ size, using weighted root square error as the criterion (Guo et al., 2006). LAI can be thus reconstructed by the calibrated model and compared with measured data, as in Fig. 2.  (Kang et al., 2008), (b) chrysanthemum (Kang et al., 2006). Initialisation of GreenLab parameters is dependent on the crop type, for which a common development and growth pattern exists. For example, a maize plant generally starts by vegetative growth with short internodes and finishes with a tassel, although the amount of leaves and position of cob in main stem may vary. Such development pattern is considered in initializing GreenLab organogenesis model. As to functional parameters, it is found that some are more or less stable for different seasons and population densities (Ma et al., 2008) for a given cultivar. Initialization of these parameters may be done according to previous modelling experience, and they will be re-estimated in following stage. In building 3D plant, an organ geometry library will be called and transformed according to their computed size and position to assemble the full plant structure. Geometrical parameters, such as organ insertion angle, are set empirically. Compared to traditional processed-based models (PBM), the feature of GreenLab lies in the modelling of plant architecture and its effect of biomass production, thereby making simulation of plant plasticity possible. Besides, GreenLab simulates plant growth as a whole dynamic system, and different variables like leaf area, stem height, spike weight are linked to each other, instead of being modeled independently. Nevertheless, the two types of model can build interface by fitting LAI from GreenLab simulation with LAI from a PBM (Feng et al., 2010). On this basis, crop fields can be simulated and visualized dynamically, giving the expected LAI. The similar principle can be applied to fitting noisy LAI from remote sensing data.

Filtering technique
To cope with temporal inconsistencies likely to occur when one estimates LAI on each image independently, we suggest the imposition of a temporal constraint, here, the GreenLab plant growth model. The dynamic coherence can be enforced by embedding the estimation problem within a filtering process. Roughly speaking, two main families can be used: variational or stochastic approaches. Variational techniques, also known as variational data assimilation, perform the estimation by minimizing a cost-function issued from a deterministic formalization in a Bayesian framework of the problem. It extracts the best compromise between observations, dynamic model and confidence measure. Due to a rewriting in a dual space (Lions, 1971), the gradient of the cost-function can efficiently be obtained using a forward-backward integration of the dynamical model, and such techniques are very adapted when one deals with large system states. On the drawback, non-linear dynamic models need to be managed in an incremental framework corresponding to a succession of linearized problems. On the other hand, when one deals with smaller system's state and non-linear dynamic models, stochastic techniques as the particle filter or the particle smoother, related to Monte Carlo approaches, are strongly adapted. The main idea consists in manipulating a set of particles more or less connected to the final state to estimate, the latter resulting from linear combination of such particles. In the next paragraph we introduce the main principles. The system state to recover at time t, noted x t , is submitted to a dynamical model f up to some uncertainties. This results in a stochastic process of the form in Eqn. (4) where n t is a centered Gaussian noise of variance σ n 2 . In addition to this dynamic process, we are able to observe the system state as in Eqn. (5): where g is a (linear or not) observation operator and v t a Gaussian noise of varianceσ v 2 . Particle filtering, also known as Sequential Monte Carlo (Doucet et al., 2001), is an attractive way to recover the system state x t for all t∈ [1,T]. It can be shown that a sequential estimation of x t can be obtained through the following system, starting from (a) an available sequential observation of sequence y 1:t = y 1 ,…, y t ; (b) an initial distribution of the system's state p(x 1 ) ; (c) transition model p(x t |x t-1 ) and observation model p(y t |x t ) respectively, related to the stochastic processes f and g presented above. See (Doucet et al., 2001) for details. Steps include: 1. Choose a set of N samples (or particles) x 1 i , i=[1, …, N] randomly taken from the initial distribution p(x 1 ) and compute p(y 1 |x 1 ); 2. Prediction step: for all t∈[2,T] and the set of samples x t-1| t-1 i , generate the predicted samples x t| t-1 i =f(x t-1|t-1 i )+ n t i for each particle and get Eqn. (6): Eqn. (7) can be approximated from the set of particles using Eqn. (8): where δx t i is the dirac function in x t i and the weight function is as in Eqn. (9): Therefore, once the system at time t-1 has been obtained, the process consists in generating a prediction of all particles at time t thanks to the transition p(x t |x t-1 ) and the available observations y 1:t-1 = y 1 ,…, y t-1 . These predictions are then corrected by taking into account the new observation y t in a second step. The final estimated distribution is obtained from the set of particles and their associated weights w(x t| t-1 i ). This is the main principle of the sequential particle filtering. When the whole sequence is available, i.e. all observations y 1:T = y 1 ,…, y T are available for all time t, a smoothing version of the previous technique can be applied by taking into account future observations. From the first estimation issued from the previous process, the idea consists in performing a backward exploration in order to correct the weights of the different particles. This is the so-called forward-backward smoother. The idea consists in reweighting the particles recursively backward in time, starting from the end time T to the initial one. It can be shown that rewriting the distribution p(x t |y 1:T ) in terms of backward transitions p(x t |x t+1 ,y 1:T ) yields, after several manipulations, to reweight all the particles with Eqn. (10): Therefore, the process consists in first performing a sequential filtering and then, to reweight the particles backward in time. More details can be found in (Doucet et al., 2001). This is the process we suggest to use in this application to recover consistent LAI values from noisy observed ones.

Computational experiment
In this computational experiment, we rely on the filtering technique to recover LAI sequence using the GreenLab model. We suppose the noisy LAI can be obtained from remote sensing data. Here we use synthetic data from GreenLab model so that the true values of y t are known for evaluation.

Case study
We chose maize plant for a case study on application of the filtering method presented above. According to previous study on maize plant (Guo et al., 2006), we set model parameters in GreenScilab, an open source software for implementing GreenLab model. The LAI can is part of model output, as shown in Fig. 3 (a). By arranging the simulated 3D maize plant according to the given density, a virtual maize field can be simulated. Fig. 3 (b) shows such an image at a plant age. It is supposed that the aim is to recover this LAI sequence from noisy data of LAI obtained from remote sensing images at several different stages.

Recovering consistent series of LAI
We tested our smoothing approach on the sequence of synthetic values of LAI generated from the GreenLab model, see Fig. 4 (blue line). From this sequence, we have randomly extracted some points on which we have added noise (black points). We assume that these points correspond to measurements obtained on the corresponding field from remote sensing images. They represent an incomplete time series of noisy values of LAI. From these inconsistent series, we have generated a smooth version from the forward-backward smoother presented in the previous section, using GreenLab as a dynamic model. This is depicted in Fig. 4 (green line). From this synthetic example, it is obvious to observe that the new series is consistent with the expected ground truth. This is confirmed when observing the quantitative values of the Root Mean Square Error between the ground truth, the noisy and reconstructed data shown in Table 1.
In order to highlight the benefit of the use of a dynamical model, we have blurred in a stronger way the series. We have assumed that during a long time period in which the variation of LAI is maximal, no observations are available (due, for instance, to the maintenance of the sensor, a too large cloud covering during the winter, etc.). In addition, to take into account the errors related to the acquisition process itself, we have also blurred the remaining data. This results in a strongly noisy sequence of LAI, as shown with the black points in Fig. 5.  These experiments on noisy and strongly noisy synthetic data demonstrate the possibility of the framework presented in this chapter (schematized in Fig. 1) to recover time consistent series of LAI with fine time resolution from a small set of noisy image observations. The practical issue with real images is under development.

Plant growth reconstruction
Using the estimated hidden GreenLab parameters in recovering the series of LAI, one can simulate the dynamics of biomass production and partitioning, as shown in Fig. 6(a), and the 3D model of the observed plants, as shown in Fig. 6(b). As for a maize plant, Fig. 6(a) gives the total plant biomass (g) in black line, and the allocation to leaf blade (green), sheath (blue), internode (brown), cob (red) and male flower (purple). Recall that empirical geometrical parameters are used in building the 3D structure.

Conclusion
In this chapter, we have presented a framework of estimating LAI series from remote sensing images using filtering techniques and GreenLab model. Computational experiment was done on synthetic data to show the feasibility of full process. The result shows that by embedding a dynamic model, the LAI series can be recovered even if the source data from remote sensing images are very noisy and sparse. And in doing so, the GreenLab model can be calibrated partially to simulate biomass production and allocation. The advantage of embedding a crop model is that the knowledge on crop development and growth can be used in recovering LAI series, and the link between LAI and biomass production make it possible to estimate biomass production from remote sensing data, which is the ultimate aim of estimating LAI. Yet this theoretical work needs to be further tested by real remote sensing sources. Challenges include the initialization of model parameters, such as the setting on topological parameters and initial source and sink parameter. Detection of crop type can help to solve this issue by providing empirical parameters. Yet their values are not necessary to be accurate, and other information from remote sensing, such as leaf chlorophyll content and leaf water content, may compensate. On the other hand, the combination of a functionalstructural plant model as GreenLab brings many possibilities. For example, as the threedimensional structures of crop are built, it is possible to run radiative transfer model in virtual canopy. Although the result will be dependent on the definition of geometrical structure and optical properties of individual organs, it provides a possibility of validating the reconstructed canopy dynamics by comparing the virtual canopy with the obtained high resolution source images. The development of remote sensing technique and advance in plant modelling are increasing the interdisciplinary research of these two areas. The accurate measurement of ecosystem biomass is of great importance in scientific, resource management and energy sectors. In particular, biomass is a direct measurement of carbon storage within an ecosystem and of great importance for carbon cycle science and carbon emission mitigation. Remote Sensing is the most accurate tool for global biomass measurements because of the ability to measure large areas. Current biomass estimates are derived primarily from ground-based samples, as compiled and reported in inventories and ecosystem samples. By using remote sensing technologies, we are able to scale up the sample values and supply wall to wall mapping of biomass. Three separate remote sensing technologies are available today to measure ecosystem biomass: passive optical, radar, and lidar. There are many measurement methodologies that range from the application driven to the most technologically cutting-edge. The goal of this book is to address the newest developments in biomass measurements, sensor development, field measurements and modeling. The chapters in this book are separated into five main sections.