**2.1 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 1st or 2nd 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.

#### **2.2 GreenLab model**

218 Remote Sensing of Biomass – Principles and Applications

crops can be identified, on which some GreenLab hidden parameters can be initialized. Besides, we suppose that the noisy time series of LAI has been estimated from remote sensing images, which have occluded area due for instance to cloud coverage, aerosols, etc. The objective is to construct a continuous LAI series by re-estimating GreenLab model parameters. This finally enables to simulate the complete plant growth and to output 3D evolution of the observed crops, using empirical geometrical parameters for the given crop. The overall strategy is illustrated in Fig. 1. The different steps of the methodology are

Fig. 1. Illustration of overall process of estimating biomass dynamics from remote sensing

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

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

the filtering technique used to recover time consistent series of LAI.

**2.1 Estimating type of crop and LAI from remote sensing data** 

presented in following sections.

images

**2. Overall framework** 

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, *Qt*, as a function of plant leaf area *St* and environmental factor *E*t, see Eqn. (1).

$$Q\_t = E\_t \cdot r \cdot S\_P \left( \mathbf{1} \cdot \exp\left( \cdot \frac{S\_t}{S\_P} \right) \right) \tag{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).

$$
\Delta q\_{t,\dot{j}}^O = \mathbb{P}^O f^O \, {}\_{\dot{j}} \cdot \mathbb{Q}\_t / \, D\_t \tag{2}
$$

In GreenLab, each type of organs (e.g. blade, sheath, internode, female organ) has certain relative sink strength *PO*, which may vary during the expansion of an individual organ, described by an empirical function *fOj*. Total plant demand *Dt* 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 (*Qt/Dt*), 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

Reconstructing LAI Series by Filtering Technique and a Dynamic Plant Model 221

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

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

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

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 *xt*, is submitted to a dynamical model *f* up to

some uncertainties. This results in a stochastic process of the form in Eqn. (4)

we are able to observe the system state as in Eqn. (5):

where *nt* is a centered Gaussian noise of variance σ*n*2. In addition to this dynamic process,

*xt* = *f*(*xt*-1)+ *nt* (4)

*yt* = *g(xt)*+ *vt* (5)

insertion angle, are set empirically.

LAI from remote sensing data.

**2.3 Filtering technique** 

growth model.

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. For some simple case, analytical equation can be written explicitly. For example, for maize, LAI at a given time t can be written as in Eqn. (3).

$$\begin{split}LAI\_{t} &= \frac{S\_{t}}{S\_{P}} = \frac{\sum\_{i=1}^{t} N\_{t,i}^{B} q\_{t,i}^{B}}{\lambda\_{t} S\_{P}} = \frac{\sum\_{i=1}^{t} N\_{t,i}^{B} \sum\_{j=1}^{\min(i, T\_{B})} P^{B} f^{B} \boldsymbol{f}\_{j} \cdot \mathbf{Q}\_{t-i+j} / D\_{t-i+j}}{\lambda\_{t} S\_{P}} \\ &= \frac{P^{B} \cdot \boldsymbol{r} \cdot \sum\_{i=1}^{t} N\_{t,i}^{B} \sum\_{j=1}^{\min(i, T\_{B})} f^{B} \boldsymbol{E}\_{j} \boldsymbol{E}\_{t-i+j} \left(1 \cdot \exp\left(\cdot L A \boldsymbol{I}\_{t-i+j}\right)\right) / D\_{t-i+j}}{\lambda\_{t}} \end{split} \tag{3}$$

Where *qBt,j* is biomass of a *j*-aged leaf blade, *NBt,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 *PO*, 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.

Fig. 2. LAI from measurement (dots) and GreenLab simulation (solid lines): (a) wheat (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.

#### **2.3 Filtering technique**

220 Remote Sensing of Biomass – Principles and Applications

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. For some simple case, analytical equation can be written

min( , )

*N q N Pf Q D*

*t i j ti j t-i j t i j*

*B*

*t i j t i j t i j ti ti*

 

1-exp - /

/

**(b)**

**0 20 40 Cycle**

(3)

explicitly. For example, for maize, LAI at a given time t can be written as in Eqn. (3).

*<sup>t</sup> i T <sup>t</sup> B B B BB*

*t*

Where *qBt,j* is biomass of a *j*-aged leaf blade, *NBt,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

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 *PO*, 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

Fig. 2. LAI from measurement (dots) and GreenLab simulation (solid lines): (a) wheat

**0**

**3**

**6**

**Leaf area index**

**9**

(Kang et al., 2008), (b) chrysanthemum (Kang et al., 2006).

, , , 1 1 <sup>1</sup>

*S S S*

*P tP t P*

*Pr N fE LAI D*

min( , )

*i j t i <sup>t</sup>*

*B*

, 1 1

*i j*

*<sup>S</sup> LAI*

measured data, as in Fig. 2.

of a leaf.

*t i T BB B*

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 *xt*, is submitted to a dynamical model *f* up to some uncertainties. This results in a stochastic process of the form in Eqn. (4)

$$\mathbf{x}\_{t} = f(\mathbf{x}\_{t:1}) + \mathbf{n}\_{t} \tag{4}$$

where *nt* 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):

$$y\_t = g(\mathbf{x}\_t) + \upsilon\_t \tag{5}$$

Reconstructing LAI Series by Filtering Technique and a Dynamic Plant Model 223

terms of backward transitions *p*(*xt|xt+1,y1:T*) yields, after several manipulations, to reweight

1

*k*

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

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 yt* are

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

Fig. 3. Synthetic result from GreenLab model. (a) simulated LAI dynamics of virtual maize

 

*j k k j*

*N j i*

*i i <sup>j</sup> t t tT t t T <sup>N</sup>*

*px x www*


1

( |)

*wpx x* 

1

*tt t*

( |)

(10)

all the particles with Eqn. (10):

from noisy observed ones.

known for evaluation.

**3.1 Case study** 

stages.

LAI

**3. Computational experiment** 

0 5 10 15 20 25 30 35

Plant Age

from GreenLab model; (b) top view of a virtual maize field;

Leaf area index

where *g* is a (linear or not) observation operator and *vt* a Gaussian noise of varianceσ*<sup>v</sup>* 2. Particle filtering, also known as Sequential Monte Carlo (Doucet et al., 2001), is an attractive way to recover the system state *xt* for all *t*∈[1,*T*]. It can be shown that a sequential estimation of *xt* can be obtained through the following system, starting from (a) an available sequential observation of sequence *y1:t = y1*,…, *yt* ; (b) an initial distribution of the system's state *p*(*x*1) ; (c) transition model *p*(*xt|xt-1*) and observation model *p*(*yt|xt*) respectively, related to the stochastic processes *f* and *g* presented above. See (Doucet et al., 2001) for details. Steps include:


$$p(\mathbf{x}\_t \mid y\_{1:t-1}) = \int p(\mathbf{x}\_t \mid \mathbf{x}\_{t-1}) p(\mathbf{x}\_{t-1} \mid y\_{1:t-1}) d\mathbf{x}\_{t-1} \tag{6}$$

3. Correction step: from the distribution *p*(*xt| y1:t-1*) and the new observation *yt*, we have a result as in Eqn. (7):

$$\begin{split} p(\mathbf{x}\_{t} \mid \mathbf{y}\_{1:t}) &= \frac{p(\mathbf{x}\_{t}, \mathbf{y}\_{1:t})}{p(\mathbf{y}\_{1:t})} \\ &= \frac{p(\mathbf{y}\_{t} \mid \mathbf{x}\_{t}) p(\mathbf{x}\_{t} \mid \mathbf{y}\_{1:t-1})}{\int p(\mathbf{y}\_{t} \mid \mathbf{x}\_{t}) p(\mathbf{x}\_{t} \mid \mathbf{y}\_{1:t-1}) d\mathbf{x}\_{t}} = w(t\_{t|t-1}) p(\mathbf{x}\_{t} \mid \mathbf{y}\_{1:t-1}) \end{split} \tag{7}$$

Eqn. (7) can be approximated from the set of particles using Eqn. (8):

$$p(\mathbf{x}\_t \mid y\_{1:t}) \approx \sum\_{i=1}^{N} w(\mathbf{x}\_{t|t-1}^i) \delta \mathbf{x}\_t^i \tag{8}$$

 where δ*xt i* is the dirac function in *xt i* and the weight function is as in Eqn. (9):

$$w(t\_{t|t-1}) = \frac{p(y\_t \mid \boldsymbol{x}\_{t-1}^i)}{\sum\_{k=1}^N p(y\_t \mid \boldsymbol{x}\_{t-1}^k)}\tag{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*(*xt|xt-1*) and the available observations *y1:t-1 = y1*,…, *yt-1*. These predictions are then corrected by taking into account the new observation *yt* in a second step. The final estimated distribution is obtained from the set of particles and their associated weights w(*xt| t-1i* ). This is the main principle of the sequential particle filtering.

When the whole sequence is available, i.e. all observations *y1:T = y1*,…, *yT* 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*(*xt|y1:T*) in terms of backward transitions *p*(*xt|xt+1,y1:T*) yields, after several manipulations, to reweight all the particles with Eqn. (10):

$$w\_{t|T}^i = w\_t^i \left| \sum\_{j=1}^N w\_{t+1|T}^j \frac{p(\mathbf{x}\_{t+1}^{\cdot j} \mid \mathbf{x}\_t^i)}{\sum\_{k=1}^N w\_t^k p(\mathbf{x}\_{t+1}^{\cdot j} \mid \mathbf{x}\_t^k)} \right| \tag{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.
