# Forecasting of Cereal Yields in a Semi-arid Area Using the Simple Algorithm for Yield Estimation (SAFY) Agro-Meteorological Model Combined with Optical SPOT/HRV Images

Aicha Chahbi Bellakanji, Mehrez Zribi, Zohra Lili-Chabaane, Bernard Mougenot

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

► To cite this version:

Aicha Chahbi Bellakanji, Mehrez Zribi, Zohra Lili-Chabaane, Bernard Mougenot. Forecasting of Cereal Yields in a Semi-arid Area Using the Simple Algorithm for Yield Estimation (SAFY) Agro-Meteorological Model Combined with Optical SPOT/HRV Images. Sensors, 2018, 18, 10.3390/s18072138. insu-04832979

HAL Id: insu-04832979

https://insu.hal.science/insu-04832979v1

Submitted on 12 Dec 2024

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

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

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

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

Article

# Forecasting of Cereal Yields in a Semi-arid Area Using the Simple Algorithm for Yield Estimation (SAFY) Agro-Meteorological Model Combined with Optical SPOT/HRV Images

Aicha Chahbi Bellakanji <fnref n="1" />,<fnref n="2" />, Mehrez Zribi <fnref n="3" />,* ID, Zohra Lili-Chabaane <fnref n="1" /> and Bernard Mougenot <fnref n="3" /> ID

Received: 17 May 2018; Accepted: 2 July 2018; Published: 3 July 2018

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

check for
updates

**Abstract:** In semi-arid areas characterized by frequent drought events, there is often a strong need for an operational grain yield forecasting system, to help decision-makers with the planning of annual imports. However, monitoring the crop canopy and production capacity of plants, especially for cereals, can be challenging. In this paper, a new approach to yield estimation by combining data from the Simple Algorithm for Yield estimation (SAFY) agro-meteorological model with optical SPOT/ High Visible Resolution (HRV) satellite data is proposed. Grain yields are then statistically estimated as a function of Leaf Area Index (LAI) during the maximum growth period between 25 March and 5 April. The LAI is retrieved from the SAFY model, and calibrated using SPOT/HRV data. This study is based on the analysis of a rich database, which was acquired over a period of two years (2010–2011, 2012–2013) at the Merguellil site in central Tunisia (North Africa) from more than 60 test fields and 20 optical satellite SPOT/HRV images. The validation and calibration of this methodology is presented, on the basis of two subsets of observations derived from the experimental database. Finally, an inversion technique is applied to estimate the overall yield of the entire studied site.

**Keywords:** yield; cereal; SAFY; optical remote sensing; SPOT/HRV

## 1. Introduction

In arid and semi-arid areas, population growth, urbanization, food security and climate change have an impact on agriculture in general and in particular on cereal production. Since the 1980s North African countries have experienced significant increases in their cereal imports. Morocco, Algeria and Tunisia account for eight per cent of world cereal imports, while they account for only one percent of the world's total population [1]. Therefore, to improve food security in arid countries, crop canopy monitoring, yield forecasting cereals and the planning of annual imports are needed [2]. For years, the scientific community has demonstrated the potential of remote sensing for the spatio-temporal monitoring of the soil properties, vegetation cover and grain yield of cereals [3–6].

Three main groups of modeling techniques used for rough scale yield estimations can be distinguished: (i) qualitative yield forecasting, and two types of quantitative yield forecasting: (ii) regression models and (iii) the combined use of models of cultivated fields and remote sensing data.

<footnote n="1">¹ LR17AGR01 (LR GREEN-TEAM)/INAT, University Carthage, Avenue de la République, P.O. Box 77, Carthage, Tunis 1054, Tunisia; chehbi.aicha@gmail.com (A.C.B.); zohra.lili.chabaane@gmail.com (Z.L.-C.)</footnote>

<footnote n="2">² ISET_Nabeul, Campus Universitaire Mrezgua, Nabeul 8000, Tunisie</footnote>

<footnote n="3">³ CESBIO (CNRS/UPS/IRD/CNES), 18 Avenue Edouard Belin, 31401 Toulouse CEDEX 9, France; Bernard.mougenot@ird.fr</footnote>

<footnote>* Correspondence: mehrez.zribi@ird.fr; Tel.: +33-561-558-525</footnote>
Qualitative crop monitoring methods are based on remote sensing observations, and are applied at national and regional levels. They provide qualitative, but not quantitative descriptions, and are generally based on a projection of the actual crop status observed during previous seasons. This approach is used by the Food and Agriculture Organization of the United Nations (FAO) for countries in Africa, Asia and South America [7], the Famine Early Warning System (FEWS) NET, the United States Agency for International Development (USAID), and the Monitoring Agriculture with Remote Sensing (MARS) program of the European Commission [8].

On the other hand, regression models can provide quantitative evaluations, based on the interpretation of a time series of satellite images and a land use map. In addition, the latter type of model must necessarily be calibrated using appropriate reference information, such as agricultural statistics and crop yield measurements. Regression-model based forecasting techniques can be subdivided into two types: those which use remote sensing data only (empirical regression models); and mixed techniques, which include additional bio-climatic predictor variables. Optical imaging techniques (from visible to thermal: 0.4 µm to 12.5 µm) are suitable for determining parameters characterizing vegetation and that by combining reflectance in different wavelength bands, in the form of vegetation indices. The vegetation is different from the other types of surface because it absorbs wavelengths in the red of solar radiation, and reflects in the near-infrared. The index most frequently used is the Normalized Difference Vegetation Index (NDVI). This index is considered as a reliable indicator of green vegetation and land-cover variations since its temporal variations are closely linked to the presence of vegetation. Other studies have shown the very high performance of the Enhanced Vegetation Index (EVI) for monitoring and assessing spatial and temporal variations in vegetation amount and condition [9,10].

Many models use the NDVI to estimate crop yields. Examples of this type of model include the model of [11] used at state level in Kansas (USA) and in the Ukraine; the model of [12], which was applied to three countries (the United States, Ukraine and China) between 2001 and 2012; and the method proposed by [13], which was applied in a semi-arid catchment in central Tunisia. These types of model rely mainly on statistical analysis of the vegetation index "NDVI", which is commonly used in remote sensing, and the measured yields. In Chahbi et al. [13], the early estimation of cereal yield is based on a statistical analysis between the vegetation index NDVI and the yields measured on the test plots. An estimation of cereal yield is possible from mid-March, when the cereal flowering phase is achieved. The validation of the remote sensing estimation shows that this approach is robust, with a high correlation coefficient R² generally greater than 0.6. These estimations show errors equal to 8.5 qx/ha and 11.6 qx/ha, respectively, for grain and straw yields. In this study, yield maps are proposed.

Although regression models are simple to apply, efficient, and more accurate than qualitative methods, it is difficult to extend a model that has been calibrated for a pilot area to other regions or other scales.

A third technique can be used to forecast crop yields, based on the use of growth models: Soil Vegetation Atmosphere (SVAT) modeling or agro-meteorological modeling [13-17]. Reference [13] proposed the application of the SAFY semi-empirical growth model, developed to simulate the dynamics of the leaf area index (LAI) and grain yields, at the field scale. Although this model is able to reproduce temporal variations of LAI for all fields, their yields are under-estimated. In practice, the statistical analysis proposed by [13] does not identify a fixed date for the relationship between vegetation indices and yields, which can lead to additional errors in the analysis, from one year to the next. To eliminate this problem, we propose an approach based on the integration of measurement. The aim of the present study is to evaluate the feasibility of a hybrid approach, combining the SAFY agro-meteorological model with remotely sensed images, at times corresponding to the maximum growth of cereals.

Section 2 is devoted to the description of the study site, the satellite data and the ground measurement. A description of the SAFY model is presented in Section 3. Furthermore, the methodology
used to estimate the dynamics and yields of cereals is described. Section 4 is dedicated to present the
results and discussions of this method. The conclusions are presented in Section 5.

## 2. Experimental Database

### 2.1. *Study Area*

The region of Kairouan is located in the eastern zone of the Tunisian Atlas, between central Tunisia and eastern Tunisia: latitude 9°30′ E–10°15′ E and longitude 35° N–35°45′ N (Figure 1). This area is characterized by a semi-arid climate. In central Tunisia, the rainy season extends from October to May, in summer, rain is almost absent. In the Kairouan region, the average annual rainfall is 400 mm per year on the heights and 300 mm per year in the plain [18]. The evapotranspiration (Penman) in this region is close to 1600 mm.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Illustration of the studied site.</figcaption></figure>

In the Kairouan plain, the main economic activity is agriculture. The area is dominated by farming, arboriculture (irrigated and non-irrigated), cereals (wheat and barley) with a notable increase in polyculture. In the recent years, intensive agriculture has been encouraged by the government [19].

### 2.2. *Satellite Data*

For North Africa, farms are for the most part small plots, due to inheritance processes that have favored the fragmentation of land [20]. Consequently, the study of agricultural surfaces by remote sensing requires the use of satellite images with a high spatial resolution such as SPOT Image (spatial resolution is equal to 10 m). Fifteen images acquired by the SPOT 5 and SPOT4 satellites were analyzed (Table 1).

These images were acquired with a repeat time of approximately 21 days and a high spatial resolution, equal to 10 m and 20 m for SPOT5 and SPOT4, respectively. Each image has four bands: green, red, near infrared (PIR) and mid-infrared MIR. In Table 2 we have the bands characteristics of the two satellites SPOT4 and SPOT5.

Radiometric and atmospheric corrections were then applied in order to estimate the reflectance of the vegetation canopy. For atmospheric correction, the Second Simulation of Satellite Signals in the
Solar Spectrum (SMAC 6s) model [21–23] was used. This model is based on an improvement of the 5s model developed by [24]. The SMAC 6s model takes into account the atmospheric parameters (ozone content, optical thickness, etc.), the geometric parameters of shooting (incidence angle, the zenith angle of the sun and the satellite, etc.). In order to limit the effect of incidence angle on measured reflectance, all acquired images were taken with an incidence angle lower than 19°.

<table><caption>**Table 1.** Satellite acquisition dates.</caption><thead><tr><th>Image</th><th>Observation Date</th><th>Sensor</th><th>θ<sup>1</sup></th></tr></thead><tbody><tr><td>1</td><td>24 December 2010</td><td>SPOT5</td><td>1.75</td></tr><tr><td>2</td><td>29 January 2011</td><td>SPOT5</td><td>11.71</td></tr><tr><td>3</td><td>19 February 2011</td><td>SPOT5</td><td>5.55</td></tr><tr><td>4</td><td>17 March 2011</td><td>SPOT5</td><td>5.50</td></tr><tr><td>5</td><td>5 April 2011</td><td>SPOT4</td><td>5.60</td></tr><tr><td>6</td><td>28 April 2011</td><td>SPOT5</td><td>8.14</td></tr><tr><td>7</td><td>18 May 2011</td><td>SPOT5</td><td>18.23</td></tr><tr><td>8</td><td>3 July 2011</td><td>SPOT4</td><td>18.23</td></tr><tr><td>9</td><td>6 November 2011</td><td>SPOT5</td><td>5.5</td></tr><tr><td>10</td><td>13 January 2012</td><td>SPOT5</td><td>7.78</td></tr><tr><td>11</td><td>28 February 2012</td><td>SPOT5</td><td>18.43</td></tr><tr><td>12</td><td>31 March 2012</td><td>SPOT5</td><td>7.67</td></tr><tr><td>13</td><td>4 May 2012</td><td>SPOT4</td><td>12.7</td></tr><tr><td>14</td><td>25 May 2012</td><td>SPOT4</td><td>5.9</td></tr><tr><td>15</td><td>6 July 2012</td><td>SPOT4</td><td>7.7</td></tr></tbody></table>

<table><caption>**Table 2.** Bands characteristics of the satellite SPOT4 and SPOT5.</caption><thead><tr><th>Mode</th><th>Band</th><th>Spectral Band</th></tr></thead><tbody><tr><td rowspan="4">Multispectral</td><td>B1: Green</td><td>0.50–0.59 µm</td></tr><tr><td>B2: Red</td><td>0.61–0.68 µm</td></tr><tr><td>B3: Near Infrared (PIR)</td><td>0.79–0.89 µm</td></tr><tr><td>B4: Mid Infrared (MIR)</td><td>1.58–1.75 µm</td></tr></tbody></table>

## 2.3. Ground Measurements

Various ground measurements were carried out during two agricultural years: 2010/2011 and 2011/2012, covering 27 test plots in the first year and 55 plots in the second year (Figure 2). The test plots were selected in such a way as to ensure a combination of wheat and barley fields that were both rain-fed and irrigated.

In order to be sure of the quality of statistical analysis, only fields with more than one hectare were considered in this selection. The in situ measurements included mainly: water content, crop height, Leaf Area Index (LAI) of the vegetation, and cereal yield at the end of the season. In Figure 3, we present an overview of the satellite and ground measurements during the two agricultural years.

<footnote n="1">¹ Incidence angle in degrees.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 2. Color-coded view of the test fields.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3. Overview of the satellite and ground measurements.</figcaption></figure>

### 2.3.1. Leaf Area Index

The LAI is defined as the total one-sided area of leaf tissue per unit of ground surface area. According to this definition, the LAI is a dimensionless quantity characterizing the canopy of an ecosystem. The LAI measurements were made on all test plots during different stages of vegetative development of the cereals, with a repeatability of two to three weeks during the agricultural season between December and April. In this study, hemispherical digital photography was used to measure the LAI. Thirty LAI estimations were made in two 20 × 20 m square zones inside each field. The LAI was derived from hemispherical digital photography based on analysis of the canopy gap fraction [16]. The proposed approach was validated by destructive measurement campaign carried out in the past over the studied site [13]. Despite the good agreement between the two approaches, a decrease in accuracy exists for non-destructive approaches for dense vegetation covers.
### 2.3.2. Cereal Yield

The cultivated area, production, and crop yield data are essential elements for agricultural statistics. The average yield is the average quantity produced per unit of cultivated area. Several different methods can be used to estimate or measure cereal yields. In the present study, a “many small frame” sampling technique was implemented [25], using a 75 × 75 cm frame. Ten samples were taken along the two diagonals of each field, and the cereal yields were then measured: the number of stalks, the weight of the grain, and the weight of the straw were recorded. From these measurements, the average grain and straw yield was determined for each field. For both agricultural years, the grain yields ranged from 11 qx/ha to 85 qx/ha, with an average value of 32 qx/ha, whereas the straws yields ranged from 8 qx/ha to 70 qx/ha, with an average value of 28 qx/ha.

## 2.4. Classification of Satellite Images over the Kairouan Plain

### 2.4.1. Land Use Map

Land use identification (Figure 4) is an essential step in the process of remotely mapping the yield of cereal fields. In the present context, this involved the production of an annual land use map of the Kairouan plain, based on a suitably adapted classification system, taking into account possible changes from one year to the next. A decision-tree type of classification was used to transpose the optical images into land use maps, for the agricultural seasons of 2010/2011 and 2011/2012. This method is based on the use of vegetation class detection thresholds, and in the present case the variable used was the Normalized Difference Vegetation Index (NDVI) [13,26]. In our study, a preliminary study of NDVI profiles generated using SPOT images was performed in order to determine the thresholds with which the separation between cultures is possible or not. For all classes, thresholds were considered from one image or from a difference between two images. Test fields were used during the learning phase in order to identify empirical NDVI thresholds allowing the different classes to be distinguished. For example, for summer vegetables, an empirical NDVI threshold was applied to the image acquired in July (NDVI > 0.3).

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Land use map for the 2010–2011 agricultural season.</figcaption></figure>
In order to assess the cereal classes (wheat or barley, irrigated or non-irrigated), satellite data recorded during March or April, corresponding to the period of maximum cereal growth, was analyzed. Pluvial olive groves cover 44% of the study area, whereas cereals corresponded to 8% of the total classified area in 2010/2011, and 22% of this area in 2011/2012. These remotely sensed classifications were validated through the use of qualitative observations over more than 100 fields in the studied site, characterized by different types of land use, revealing an overall accuracy of approximately 80%. Once the land use map has been established, it can be used to define a cereal mask for the two agricultural seasons under study.

### 2.4.2. Classification of Irrigated and Non-Irrigated Cereals

In this section, we focus on the definition of a cereal mask used to identify two different types of cereals: irrigated and non-irrigated. For this, an object-oriented classification technique was used, in which the classification takes place in two stages: segmentation of the image into objects, followed by classification of this segmentation. The object-oriented approach operates on adjacent pixel groups, and makes use of various spectral and spatial homogeneity criteria (spectral value, shape, texture, etc.) in order to discern between the objects to be classified [27]. The decision rules used for such an operation include not only the spectral, but also the spatial parameters. This process is inspired from the multi-resolution segmentation concept proposed in the Orfeo Toolbox OTB Monteverdi software [28].

As cereal vegetation growth peaks in the month of March, the Spot images used in this study were recorded on 17 March 2011 for the 2010/2011 agricultural season, and 31 March 2012 for the 2011/2012 season, and the cereal masks corresponding to these images were applied. The segmentation of these images was carried out using the Mean-shift Clustering principle [29–33].

Figure 5a shows a segmented image of the Kairouan plain for the 2010/2011 agricultural season. Once the image has been segmented into spatial objects, which in reality correspond to homogeneous cereal plots, they are classified according to agricultural practice: non-irrigated or irrigated. This classification is based on vegetation index (NDVI, TSAVI...). Plots chosen for the validation of classification come from field observations in the study area. The choice of these plots was made in order to have a mixture between wheat and barley and which are conducted either in rainfed or irrigated.

<figure><img src="image_9.png" /><figcaption>Figure 5. *Cont.*</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 5.** (a) Segmented image produced for the 2010/2011 agricultural season; (b) Land-use map of irrigated and non-irrigated cereals for the 2010/2011 agricultural season.</figcaption></figure>

The segmented and classified image (Figure 5b) was validated using 55 known plots, 36 of which were irrigated and 19 were pluvial. The map corresponding to the 2010/2011 agricultural season has an overall accuracy of about 97.3%, and a precision Kappa coefficient equal to 94%. In the case of the 2011/2012 agricultural year, the overall coefficient is approximately 81%, and the precision Kappa coefficient is equal to 62%. These maps were thus considered to be valid for the study of cereal yields at the plot scale.

# 3. Methodology for the Estimation of Cereal Yields Using the SAFY Growth Model

## 3.1. SAFY Model Description

Over the last few decades, several techniques have been developed for the regional monitoring and spatialization of vegetation growth, based on the combined use of a model and remotely-sensed optical images. Among these approaches, the model of [34] provides the simulation of dry biomass production and the interception of light by vegetation. For the modeling of growth and crop yields, it can be advantageous to introduce the leaf area index (LAI), which is involved in the production of vegetative biomass and is a key variable in the functioning of crops. In this context, the Simple Algorithm for Yield Estimate (SAFY) model developed by [16] was used to take into account the main processes of cereal development and growth at the plot scale.

This model is based on the light-use efficiency theory of [34]. It provides a simulation of the increase in dry, above-ground phytomass. Also, it takes into account the influence of temperature and the dynamics of green leaves [16].

Each day of the vegetative period, the production of vegetative biomass is determined from the global radiation (Rg) in photo-synthetically active radiation via the climatic efficiency $ε_C$. A Part of this radiation is absorbed by the vegetation cover through the light interception efficiency $ε_I$. The conversion of the global radiation by $ε_C$ and $ε_I$ gives rise to photosynthetically active radiation absorbed by vegetation (APAR). This term will be converted to dry above-ground biomass (ΔDAM)
through light-use efficiency (ELUE). Also, ΔDAM is affected by the daily average of air temperature (Ta) through the temperature-stress-function F<sub>T</sub>.

This conversion of global radiation to dry above-ground biomass is done by the following equations:

$$ \Delta \text{DAM} = \text{Rg} \times \varepsilon_C \times \varepsilon_1 \cdot \text{ELUE} \times \text{F}_T(\text{T}_a), \quad (1) $$

$$ \Delta \text{DAM} = \text{APAR} \times \text{ELUE} \times \text{F}_T(\text{T}_a), \quad (2) $$

During the growth phase, part of the ΔDAM is allocated to leaf biomass. This fraction is calculated using the leaf partition function Pl which varies between 0 and 1. Thus, the daily production of leaf biomass (DAM × Pl) is converted into leaf area (ΔLAI<sup>+</sup>) via the specific leaf area SLA:

$$ \Delta \text{LAI}^{+} = \Delta \text{DAM} \times \text{Pl} \times \text{SLA} \quad \text{si} \quad \text{Pl} > 0 \quad (3) $$

The senescence of leaves starts when accumulated air temperature reaches a given threshold Stt. ΔLAI<sup>-</sup> increases with thermal time (ΣTa) and is modulated by the Rs parameter which defines the rate of senescence:

$$ \Delta \text{LAI}^{-} = \text{LAI} \times \frac{\sum \text{Ta} - \text{Stt}}{\text{Rs}} \quad \text{si} \quad \text{SMT} > \text{Stt} \quad (4) $$

Thus, the LAI at day j is calculated from the LAI at day (j − 1), the term of increase of the leaf surface (ΔLAI<sup>+</sup>) and the term of senescence (ΔLAI<sup>-</sup>).

The SAFY model also simulates grain yield. The daily grain increase “ΔGY” is proportional to the dry above-ground biomass DAM and to the grain partition function P<sub>y</sub>. This is given by Equation (5):

$$ \Delta \text{GY} = \text{DAM} \times \text{P}_y \quad (5) $$

SAFY has a low level of complexity, which simplifies the optimization of unknown parameters (limited to not more than 14), using a small number of observations. From these parameters, a priori values were determined on the basis of various prior experimental studies [35,36]. However, three of these parameters (the day of plant emergence D<sub>0</sub>, the parameter ELUE, and the “sum of temperature for senescence” S<sub>TT</sub>) are strongly dependent on the ambient agro-environmental conditions. In the Table 3, we have the range of variation of these three parameters for the test plots.

<table><caption>**Table 3.** Range of variation of calibrated parameters for test plots.</caption><thead><tr><td>Parameter</td><td>Unit</td><td>Range of Variation</td></tr></thead><tbody><tr><td>D<sub>0</sub></td><td>Day</td><td>15–120</td></tr><tr><td>ELUE</td><td>g·MJ<sup>-1</sup></td><td>0–10</td></tr><tr><td>S<sub>TT</sub></td><td>°C</td><td>200–1800</td></tr></tbody></table>

### 3.2. Application of SAFY to Cereal Cycle Retrieval

The first step involves the calibration of the three aforementioned parameters: D<sub>0</sub>, ELUE, and S<sub>TT</sub>. This procedure involves the identification of an optimum parameter set from which, for each field, the SAFY simulation leads to the best reproduction of several observed variables.

These parameters are calibrated using an optimization algorithm based on the comparison between in situ values of LAI observed in the fields and those estimated with the SAFY model [13].

For each test plot, the optimized parameters (D<sub>0</sub>, ELUE, and S<sub>TT</sub>) lead to the lowest Root Mean Square Error (RMSE) between these two values (Figure 6). It can be obvious that the model correctly retrieves the dynamic of the LAI for test plots.
<figure><img src="image_11.png" /></figure>

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

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

<figure><img src="image_14.png" /><figcaption>**Figure 6.** Values of observed (star) and simulated (line) LAI, estimated using the SAFY model for four test plots on the Kairouan plain during both agricultural seasons (1st Julian date = 1 November).</figcaption></figure>
### 3.3. Proposed Approach

Although it has been shown that statistical empirical analysis, which associates the NDVI index during the period of maximum growth with the expected cereal yield, is highly accurate [13], it is not possible to predict a fixed annual date for this relationship. This can lead to additional errors in the analysis, from one year to another. In order to reduce this drawback, we propose an approach based on integral measurements of maximum growth recorded over a decade. The maximum value of cereal LAI at our site occurs towards the end of March, when the vegetation is flowering. Before this date, the maximum development of the vegetal cover is not yet reached for all the plots, it varies according to the date of sowing, the type of irrigation, the conduct of the fertilization, etc. After mid-April, the NDVI begins to fall, which may also lead to errors in yield estimation.

The grain yield is thus expressed as a function of LAI, during the cereal development period between 25 March and 5 April, which corresponds to the range between the 146th and 157th Julian Days. We note this parameter $A_{LAI}$ (Figure 7). In the initial SAFY model, the yield is expressed as a function of the total dry aerial phytomass (including, among others, the maximum LAI). Application of the SAFY model to the study of the Kairouan plain raised various questions concerning the model's limitations in terms of yield estimation. It was thus proposed to improve these estimations by including the LAI integral ($A_{LAI}$) in the grain yield calculation. Since the SAFY model accurately extracts the variation of LAI along the vegetative cycle, the estimation of $A_{LAI}$ will be established by this model for each plot (segment). Consequently, the phenological development and the sowing date ($D_0$) were taken into consideration:

$$GY = f\left(\int_{146}^{157} LAI\right) \qquad (6)$$

<figure><img src="image_15.png" /><figcaption>**Figure 7.** Principle of cereal yield estimations, based on the integral of the LAI over the period of maximum canopy development ($A_{LAI}$).</figcaption></figure>

## 4. Results and Discussions

### 4.1. Estimating the Yields of All Cereals

For this part, our database set was divided into three disjoint folds of the same size. In these sequential folds, the first was used for validation procedures, and the two remaining folds were used
for assessing the yield model. As in the case of the “repeated holdout” method, the overall accuracy is given by the average of the values obtained from all runs.

The method described in the previous section was initially applied to the first data set, regardless of the variety of crops and cultivation techniques used (wheat, barley, irrigated, non-irrigated). Figure 8a shows a plot of the raw data, expressing measured grain yields as a function of A<sub>LAI</sub>, together with the corresponding linear regression line. The two confidence intervals are also shown: the nearest interval to the curve is around the mean of the estimator and the farthest is the interval around the point estimate). The two variables are correlated (R<sup>2</sup> = 0.58):

$$GY_{qx/ha} = 0.33A_{LAI} + 17.97 \quad (7)$$

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

<figure><img src="image_17.png" /><figcaption>**Figure 8.** (a) Relationship between measured cereal grain yields and A<sub>LAI</sub>; (b) Measured grain yields compared with yield estimates, computed using the maximum growth period LAI, for two test plots and two crop years.</figcaption></figure>

Figure 8b shows a plot of the measured yields as a function of the values estimated using the maximum LAI area for the second data set. These two parameters are found to be well correlated, with the RMSE equal to 7.57 quintals/ha.
## 4.2. Wheat and Barley Yield Estimations

In a second step, the approach described above was applied to the case of two different cereals: wheat and barley, as shown in Figure 9. The measured and estimated yields are found to be well correlated ($R^2 = 0.63$ for wheat, and $R^2 = 0.77$ for barley), with (RMSE) errors equal to 7.95 and 5.37 quintals/ha for wheat and barley, respectively.

<figure><img src="image_18.png" /><figcaption>Figure 9. *Cont.*</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 9. Relationship between the grain yields of wheat (a) and barley (b), and LAI during the period of maximum growth (A$_{LAI}$). Measured yields shown as a function of the estimated values computed using A$_{LAI}$, the maximum growth LAI, for test plots and two crop years, for wheat (c) and barley (d).</figcaption></figure>

### 4.3. Yield Estimations for Irrigated and Non-Irrigated Areas

In a third step, we focus on the estimation of grain yields for two types of cereal growing conditions: irrigated and rain-fed, organized according to two separate datasets. The first of these (irrigated cereals) was used to establish a relationship (linear regression, Equation (8)) between the estimated yield and the LAI during the period of maximum growth, as shown in Figure 10a. When the estimated and measured yields are then compared using this regression they are well correlated, with R² equal to 0.53 and the RMSE equal to 6.38 qx/ha (Figure 10b).

The second dataset corresponding to non-irrigated cereals was analyzed in the same manner, leading to a strong correlation (R²$_{non-irrigated}$ = 0.66), with the RMSE equal to 4.7 qx/ha (Figure 11b).

The linear regressions determined for each of these analyses are:

$$ 
\begin{aligned}
GY_{\text{irrigated}} &= 0.39 \times A_{\text{LAI}} + 16.68 \\
GY_{\text{non-irrigated}} &= 0.45 \times A_{\text{LAI}} + 16.05
\end{aligned}
\qquad (8) 
$$

We note that using this procedure, all of the data points lie inside the second confidence interval, for both types of cereal growing conditions (irrigated and non-irrigated).

<figure><img src="image_20.png" /><figcaption>Figure 10. Cont.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 10.** (a) Relationship between the grain yields of irrigated cereals and LAI during the period of maximum growth; (b) Measured grain yields in irrigated test plots compared with yield estimates, computed using the maximum growth period LAI, over two crop years.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 11.** (a) Relationship between the grain yields of non-irrigated cereals and LAI during the period of maximum growth; (b) Measured yields in non-irrigated test plots, compared with the estimated values computed using the maximum growth LAI over two crop years.</figcaption></figure>
4.4. Spatialisation of Grain Yield

The spatialisation grain-yield technique is explained in the following organigram (Figure 12). In this section of the study, a distinct inversion was computed for each segment (or plot). This approach is the same as that adopted by [13]. Following calibration of the SAFY model, through the use of multi-temporal LAI satellite maps and LAI values computed for the maximum growth period, yield estimations are proposed for the entire studied site.

<figure><img src="image_23.png" /><figcaption>Figure 12. The organigram of the spatialisation of grain-yield.</figcaption></figure>

In order to calibrate the SAFY model for each segment, it is first necessary to invert the satellite images of the agricultural year in question, in NDVI maps. Then these images are reversed in LAI via a semi-empirical relationship between NDVI and LAI established for cereal crops and detailed in [13]. The results are LAI maps for each segment and for five different dates.

These maps are used to adjust the simulations performed by the SAFY model. The optimization procedure is applied for each segment in order to minimize the difference between the five observations (LAI derived from SPOT-HRV) and simulated values at observation dates.

To estimate the three parameters of the SAFY model (ELUE, D₀ and S<sub>TT</sub>), the method mentioned in Section 3.2 was applied for each segment. Furthermore, to reduce the effect of local terrain heterogeneities in segment, these three parameters are estimated over cells corresponding to 10 × 10 pixels (approximately 100 m²). Maps of the three parameters are generated (ELUE, D₀, S<sub>TT</sub>). After having generated these parameters, we extract the A<sub>LAI</sub> variable for each segment. Finally, the empirical relationship determined in Section 4 is used to generate cereal yield maps for each segment (Figure 13).
<figure><img src="image_24.png" /><figcaption>**Figure 13.** Cereal yield map produced by combining SPOT/HRV multi-temporal acquisitions with SAFY-modeled yields computed with the LAI corresponding to the period of maximum growth.</figcaption></figure>

**5. Conclusions**

The aim of this study was to characterize the dynamics and yields of various cereal plantations in the Kairouan plain, using high-resolution optical data and a simple agricultural model (SAFY). The SAFY model accurately simulates the variation of LAI along the vegetative cycle of cereals.

The estimated grain yield is a function of the LAI during maximum growth (a ten-day period, between the end of March and the beginning of April). For all types of cereal (irrigated or non-irrigated, wheat and barley), a good correlation is observed between the measured and estimated yields. Our validation of the proposed empirical relationships over different types of cover shows that this approach is quite accurate, with an RMSE lower than 7.6 qx/ha. When applied to intensive agriculture on the Kairouan plain, it performs very well, because climatic factors (with global radiation and average daily air temperature used as input parameters), cereal phenological parameters, agro-environmental conditions, and grain growth cycle characteristics (short and narrow, or wide, longer-lasting cycles) as well as senescence, are taken into account. This methodology could be applied in a more operational fashion than the determination of a direct relationship between satellite indices and yield estimations, for two reasons: firstly, optical data may not be usable for periods of many days, due to cloud cover; secondly, since optical measurements are not recorded on the same Julian days each year, they would need to be calibrated with respect to their timing within the 10-day span of the maximum growth period. In future studies, this technique will be tested in an operational context using Sentinel2 data.

**Author Contributions:** A.C.B. processes all data, analyses results and writes the proposed paper. M.Z. analyses results, participates to the paper writing and coordinates the scientific study. Z.L.-C. participates to analysis of results. B.M. participates to the processing of data.

**Funding:** This study was funded by the French SICMED/Mistrals program, the French AMETHYST ANR-12-TMED-0006-01 and TOSCA/ASCAS projects.

**Conflicts of Interest:** The authors declare no conflicts of interest.
## References

1. Palmieri, L.; Amantonico, B. *Rapport Etude de Faisabilité: Etude D'identification pour la Mise en Œuvre de L'observation Agro-Alimentaire Italo-Tunisien*. Programme "Instrument Européen de Voisinage et de Partenariat de Coopération Transfrontalière (IEVP CT) Italie-Tunisie 2007-2013". 2013, p. 23. Available online: [servagri.eu/index.php/fr/documentation/etude-de-faisabilite?...15:etude...faisabilite](http://servagri.eu/index.php/fr/documentation/etude-de-faisabilite?...15:etude...faisabilite) (accessed on 20 April 2016).

2. Justice, C.O.; Becker-Reshef, I. Report from the Workshop on Developing a Strategy for Global Agricultural Monitoring in the Framework of Group on Earth Observations (GEO). 2007, pp. 1-67. Available online: [http://www.fao.org/gtos/igol/docs/meeting-reports/07-geo-ag0703-workshop-report-nov07.pdf](http://www.fao.org/gtos/igol/docs/meeting-reports/07-geo-ag0703-workshop-report-nov07.pdf) (accessed on 11 June 2015).

3. Fratanelli, G.; Paloscia, S.; Zribi, M.; Chahbi, A. Sensitivity analysis of X-band SAR to wheat and barley biomass in the Merguellil Basin. *Remote Sens. Lett.* **2013**, *4*, 1107-1116.

4. Fieuzal, R.; Duchemin, B.; Jarlan, L.; Zribi, M.; Baup, F.; Merlin, O.; Dedieu, G.; Garatuza-Payan, J.; Watt, C.; Chehbouni, A. Combined use of optical and radar satellite data for the monitoring of irrigation and soil moisture of wheat crops. *Hydrol. Earth Syst. Sci.* **2011**, *15*, 1117-1129. [CrossRef]

5. Barnett, T.L.; Thompson, D.R. The use of large-area spectral data in wheat yield estimation. *Remote Sens. Environ.* **1982**, *12*, 509-518. [CrossRef]

6. Moriondo, M.; Maselli, F.; Bindi, M. A simple model of regional wheat yield based on NDVI data. *Eur. J. Agron.* **2007**, *26*, 266-274. [CrossRef]

7. Baghdadi, N.; Zribi, M. *Land Surface Remote Sensing in Agriculture and Forest*; ISTE Press: London, UK; Elsevier: Oxford, UK, 2016; ISBN 9781785481031.

8. Laguette, S.; Vidal, A.; Vossen, P. Télédétection et estimation des rendements en blé en Europe. *Ingénieries-EAT* **1997**, *12*, 19-33.

9. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. *Remote Sens. Environ.* **2002**, *83*, 195-213. [CrossRef]

10. Johnson, D.M. A comprehensive assessment of the correlations between field crop yields and commonly used MODIS products. *Int. J. Appl. Earth Obs. Geoinf.* **2016**, *52*, 65-81. [CrossRef]

11. Becker-Reshef, I.; Vermote, E.; Lindeman, M.; Justice, C. A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data. *Remote Sens. Environ.* **2010**, *114*, 1312-1323. [CrossRef]

12. Franch, B.; Vermote, E.F.; Becker-Reshef, I.; Claverie, M.; Huangc, J.; Zhang, J.; Justice, C.; Sobrinod, J.A. Improving the timeliness of winter wheat production forecast in the United States of America, Ukraine and China using MODIS data and NCAR Growing Degree Day information. *Remote Sens. Environ.* **2015**, *161*, 131-148. [CrossRef]

13. Chahbi, A.; Zribi, M.; Lili-Chabaane, Z.; Duchemin, B.; Shabou, M.; Mougenot, B.; Boulet, G. Estimation of the dynamics and yields of cereals in a semi-arid area using remote sensing and the SAFY growth model. *Int. J. Remote Sens.* **2014**, *35*, 1004-1028.

14. Moulin, C. Assimilation D'observations Satellitaires Courtes Longueurs D'onde dans un Modèle de Fonctionnement de Culture. Ph.D. Thesis, l'Université Paul Sabatier de Toulouse, Toulouse, France, 1995.

15. Bastiaanssen, W.G.M.; Ali, S. A new crop yield forecasting model based on satellite measurements applied across the Indus Basin, Pakistan. *Agric. Ecosyst. Environ.* **2003**, *94*, 321-340. [CrossRef]

16. Duchemin, B.; Maisongrande, P.; Boulet, G.; Benhadj, I. A simple algorithm for yield estimates: Evaluation for semi-arid irrigated winter wheat monitored with green leaf area index. *Environ. Model. Softw.* **2008**, *23*, 876-892. [CrossRef]

17. Balaghi, R.; Tychon, B.; Eerens, H.; Jlibene, M. Empirical regression models using NDVI, rainfall and temperature data for the early prediction of wheat grain yields in Morocco. *Int. J. Appl. Earth Obs. Geoinf.* **2008**, *10*, 438-452. [CrossRef]

18. Zribi, M.; Dridi, G.; Amri, A.; Lili-Chabaane, Z. Analysis of the effects of drought on vegetation cover in a Mediterranean region through the use of SPOT-VGT and TERRA-MODIS long time series. *Remote Sens.* **2016**, *8*, 992. [CrossRef]

19. Gorrab, A.; Zribi, M.; Baghdadi, N.; Mougenot, B.; Fanise, P.; Chabaane, Z.L. Retrieval of both soil moisture and texture using TerraSAR-X images. *Remote Sens.* **2015**, *7*, 10098-10116. [CrossRef]
20. Bousbih, S.; Zribi, M.; Lili-Chabaane, Z.; Baghdadi, N.; El Hajj, E.; Gao, Q.; Mougenot, B. Potential of Sentinel-1 Radar Data for the Assessment of Soil and Cereal Cover Parameters. *Sensors* **2017**, *17*, 2617. [CrossRef] [PubMed]

21. Rahman, H.; Dedieu, G. SMAC: A simplified method for the atmospheric correction of satellite measurements in the solar spectrum. *Remote Sens.* **1994**, *15*, 123–143. [CrossRef]

22. Berthelot, B.; Dedieu, G. Correction of Atmospheric Effects for VEGETATION Data. In *Physical Measurements and Signatures in Remote Sensing, Proceedings of the Seventh International Symposium on Physical Measurements and Signatures in Remote Sensing, Courchevel, France*, 7–11 April 1997; Guyot, G., Phulpin, T., Eds.; Balkema: Courchevel, France, 1997; pp. 19–25.

23. Kotchenova, S.Y.; Vermote, E.F.; Matarrese, R.; Klemm, F.J., Jr. Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: Path radiance. *Appl. Opt.* **2006**, *45*, 6762–6774. [CrossRef] [PubMed]

24. Tanré, D.; Deroo, C.; Duhaut, P.; Herman, M.; Morcrette, J.J. Description of a computer code to simulate the satellite signal in the solar spectrum: The 5S code. *Int. J. Remote Sens.* **1990**, *11*, 659–668. [CrossRef]

25. Norman, D.W.; Worman, F.D.; Siebert, J.D.; Modiakgotla, E. *The Farming Systems Approach to Development and Appropriate Technology Generation*; FAO Farm System Management Series 10; Food and Agriculture Organization of the United Nations: Rome, Italy, 1995.

26. Zribi, M.; Chahbi, A.; Shabou, M.; Lili-Chabaane, Z.; Duchemin, B.; Baghdadi, N.; Amri, R.; Chehbouni, A. Soil surface moisture estimation over a semi-arid region using ENVISAT ASAR radar data for soil evaporation evaluation. *Hydrol. Earth Syst. Sci.* **2011**, *15*, 345–358. [CrossRef]

27. Benz, U.C.; Hofmann, P.; Willhauck, G.; Lingenfelder, I.; Heynen, M. Multi-resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information. *ISPRS J. Photogramm. Remote Sens.* **2004**, *58*, 239–258. [CrossRef]

28. Grizonnet, M.; Inglada, J. Monteverdi—Remote Sensing Software from Educational to Operational Context. In *Proceedings of the EARSEL Symposium: Remote Sensing for Science, Education and Culture, Paris, France*, 31 May–4 June 2010; pp. 749–755.

29. Comaniciu, D.; Ramesh, V.; Meer, P. The variable bandwidth mean shift and data-driven scale selection. In *Proceedings of the International Conference on Computer Vision, Vancouver, BC, Canada*, 7–14 July 2001; Volume I, pp. 438–445.

30. Duda, R.; Hart, P.; Stork, D. The estimation of the gradient of a density function, with applications in pattern recognition. *IEEE Trans. Inf. Theory* **2001**, *21*, 32–40.

31. DeMenthon, D.; Megret, R. *Spatio-Temporal Segmentation of Video by Hierarchical Mean Shift Analysis*; Technical Report; University of Maryland: College Park, MD, USA, 2002.

32. Commaniciu, D.; Ramesh, V.; Meer, P. Kernel-based object tracking. *IEEE Trans. Pattern Anal. Mach. Intell.* **2003**, *25*, 564–577. [CrossRef]

33. Singh, M.; Ahuja, N. Regression based bandwidth selection for segmentation using parzen windows. In *Proceedings of the International Conference on Computer Vision, Nice, France*, 13–16 October 2003; pp. 2–9.

34. Monteith, J.L. Climate and the efficiency of crop production in Britain. *Philos. Trans. R. Soc. Lond. Ser. B* **1977**, *281*, 277–294. [CrossRef]

35. Porter, J.R.; Gawith, M. Temperatures and the growth and development of wheat: A review. *Eur. J. Agron.* **1999**, *10*, 23–36. [CrossRef]

36. Claverie, M.; Demarez, V.; Duchemin, B.; Hagolle, O.; Ducrot, D.; Marais-Sicre, C.; Dejoux, J.F.; Huc, M.; Keravec, P.; Béziat, P.; et al. Maize and Sunflower Biomass Estimation in Southwest France Using High Spatial and Temporal Resolution Remote Sensing Data. *Remote Sens. Environ.* **2012**, *124*, 844–857. [CrossRef]

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