Toward Precision in Crop Yield Estimation Using Remote Sensing and Optimization Techniques

Abstract

Many crop yield estimation techniques are being used, however the most effective one is based on using geospatial data and technologies such as remote sensing. However, the remote sensing data which are needed to estimate crop yield are insufficient most of the time due to many problems such as climate conditions (% of clouds), and low temporal resolution. There have been many attempts to solve the lack of data problem using very high temporal and very low spatial resolution images such as Modis. Although this type of image can compensate for the lack of data due to climate problems, they are only suitable for very large homogeneous crop fields. To compensate for the lack of high spatial resolution remote sensing images due to climate conditions, a new optimization model was created. Crop yield estimation is improved and its precision is increased based on the new model that includes the use of the energy balance equation. To verify the results of the crop yield estimation based on the new model, information from local farmers about their potato crop yields for the same year were collected. The comparison between the estimated crop yields and the actual production in different fields proves the efficiency of the new optimization model.

1. Introduction

Remote sensing has proved to be game changing for the agricultural sector, as it is one of the backbones of precision in agriculture [1]. In remote sensing, multispectral and hyperspectral satellite images play a major role in crop management; their ability to represent crop growth conditions on the spatial and temporal scale is remarkable. These images can describe the crop development, photosynthetic active radiation (PAR), biomass accumulation (Bio), leaf area index (LAI), Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (ENDVI), Perpendicular Crop Enhancement Index (PCEI), and many other indices [2,3]. In addition, climate data related to crops can also be extracted from satellite images such as albedo, surface temperature, and actual evapotranspiration [3,4,5].

Many approaches have been developed to translate remote sensing data into crop yields, and several reviews of such approaches exist [6,7,8]. Some of these approaches have faced several problems in estimating crop yield. One of these problems is the scarcity of remote sensing data suitable for use in crop management because of climatic conditions such as clouds [9,10]. Most of the time climatic conditions and low temporal resolution are the main obstacles that prevent decision makers from using remote sensing data to map crops and to estimate crop yields.

There are several articles that have used vegetation indices from one date image or multi-temporal images to estimate crop yield. Bastiaanssen and Ali [11] used Monteith’s model [12] in their research for the calculation of Absorbed Photosynthetically Active Radiation (APAR), Stanford’s model for determining light use efficiency, and the Surface Energy Balance Algorithm for Land (SEBAL) to describe the temporal and spatial variabilities in land wetness conditions. The results of this research showed that there were gaps between the estimated and the actual yield of about 1075 and 1246 kg/ha for wheat and rice due to a lack of enough remote sensing images that cover the crop growth season. In addition, SEBAL is a well-known energy based actual evapotranspiration model, but it is difficult to implement due to a priori intensive requirements of climate variables. In this paper, some of the techniques used by [11,12] to estimate crop yield are adapted for this research.

The hypothesis is that missing satellite data can be compensated by using new models. The available biomass data which can be calculated using a reliable energy balance model can be fitted to a specific curve representing the new model using an optimization algorithm.

In this research, the objectives were to estimate crop yield using another known energy balance model named Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC) [13]. This model is less dependent on climatic measurements to compute actual evapotranspiration.

The next step was to use satellite images such as Landsat 7 and 8 images with the METRIC model to estimate crop biomass. Later, in order to strengthen the performance and efficiency of the crop yield model, a new model was developed. The model is based on using an optimization algorithm to compensate for the absence of satellite data due to climatic conditions. Finally, it was proved that the new model is able to provide efficient results.

2. Materials and Methods

The study area is located in the largest agricultural plane in Lebanon, known as Bekaa valley. The plane has an area of more than 1000 km², and the major crops in the valley are wheat, potato, and alfalfa. It is also important to note that two Bowen Ratio stations are installed in the area in order to measure the different climate parameters which are necessary for irrigation management, for mapping actual evapotranspiration, and for estimating crop yield (Figure 1).

The two Bowen Ratio stations (Figure 2) are located in two different areas, one in an agricultural research institute in the middle of the plane, while the other is located south in a field owned by a potato chips manufacturer. The distance between both is more than 15 km. The two stations provide many valuable climate data such as net radiation, wind speed and direction, temperature at different elevations, humidity, global radiation, soil heat flux, and soil humidity. Increasing the number of Bowen ratio stations in the future could improve any model created based on the data of these stations. The images of Landsat 7 and 8 are used in this research and they can be downloaded free of charge from the United States Geological Survey (USGS) Web site [14].

In addition, these satellite images have thermal infrared bands with high spatial resolution (60 m compared to 1000 m for Modis Aqua and Terra satellites). The long wave infrared bands are used later in computing important parameters for crop yield estimation. Landsat images are collected between the date a specific crop was planted and the date the same crop is harvested. It is possible to collect all satellite images for the whole season if yield estimation is for more than one crop. Using both Landsat 7 and 8 improves the temporal resolution of image availability from 16 days to 8 days.

2.1. Estimating Crop Yield

Monteith model [12] is based on the photosynthetically active radiation (PAR) (0.4–0.7 m) which is part of the short wave solar radiation (0.3–3.0 m) that is absorbed by chlorophyll for photosynthesis in the plants.

Photosynthetically active radiation (PAR) is thus a fraction of the incoming solar radiation (K_inco). Then a fraction of PAR is absorbed by the plant for carbon assimilation.

Absorbed photosynthetically active radiation (APAR) in $\left({\mathrm{Wm}}^{-2}\right)$ can be computed directly from PAR using the following equation

$$APAR=fPAR$$

(1)

The factor f can be estimated using Normalized Differences Vegetation Index (NDVI) as is explained in [11]. Where

$$f=-0.161+1.257NDVI$$

(2)

NDVI is the normalized difference vegetation index computed as the difference between near-infrared (NearIR) and red spectrums divided by their sum.

$$NDVI=\frac{NearIR-Red}{NearIR+Red}$$

(3)

The accumulation of biomass is according to the Monteith model proportional to accumulated APAR.

$$Bi{o}_{act}^{tot}=\epsilon {\displaystyle \sum }APAR\left(t\right)\left(t\right)$$

(4)

where $Bi{o}_{act}^{tot}$ in (kg/m²) is the accumulated biomass in period t, ε is the light use efficiency in gram per mega joules (g MJ⁻¹).

A more comprehensive global ecology model for computing net production was created by Field et al. [15] where they used the following equation for light use efficiency ε.

$$\epsilon ={\epsilon }^{\prime }{T}_1{T}_{2}W$$

(5)

where ε′ is the maximum conversion element for above ground biomass when the environmental conditions are optimal, it is equal to 0.29 g/MJ for C3 crops [16]. W is a function of the effective fraction of the available soil moisture. The complete details for the calculation of W are explained later in the text.

$$\begin{gathered} T_1=0.8+0.02T_{opt}-0.0005T_{opt}^2 \\ {T}_2=\frac{1}{1+\exp \left(0.2T_{opt}-10-T_{mon}\right)}\times \frac{1}{\exp \left(0.3\left(-T_{opt}-10+T_{mon}\right)\right)} \end{gathered}$$

(6)

T₁ and T₂ are two different heat functions, T_opt (°C) which is the mean air temperature in the month with maximum leaf area index, and T_mon (°C) is the mean monthly air temperature. One can notice easily that T₁ depends completely on T_opt, which is in turn affected by canopy biomass. On the other hand, T₂ depends on both T_opt and T_mon.

2.2. Enhancement and Correction of the Images

Normally, Landsat “SLC-off” data refers to all Landsat 7 images collected after 31 May 2003, when Scan Line Corrector (SLC) failed. These products have data gaps, but are still useful and maintain the same radiometric and geometric corrections as data collected prior to the SLC failure. To fix this problem, a method created by Scaramuzza et al. [17] is used to fill gaps in one scene with data from another Landsat scene. In this method, linear transform is applied to the image with gaps to fix the problem based on the standard deviation and mean values of each band. The at-surface reflectance for the visible to short-wave is corrected on a band-by-band basis following Tasumi et al. method [18]. It works first on correcting images at sensor (Atmospheric correction) and then at surface correction.

Solar radiation reflected by the Earth’s surface to satellite sensors is modified by its interaction with the atmosphere. The objective of applying an atmospheric correction is to determine true surface reflectance values and to retrieve physical parameters of the Earth’s surface, including surface reflectance, by removing atmospheric effects [19]. Atmospheric correction for satellite images works as follow:

$${{\rho }^{\prime }}_{t,b}=Mp\ast Qcal+Ap$$

(7)

where: ${{\rho }^{\prime }}_{t,b}$ = top of atmosphere (TOA) reflectance, without correction for solar angle. Mp = Band-specific multiplicative rescaling factor, Ap = Band-specific additive rescaling factor, and Qcal = Quantized and calibrated standard product pixel values (DN). The constants Mp and Ap are extracted from the metadata file that comes with the image. While Qcal value is extracted from the image pixel values. To correct TOA for the sun angle, the following is applied

$${\rho }_{t,b}=\frac{{{\rho }^{\prime }}_{t,b}}{Sin\left({\theta }_{Se}\right)}$$

(8)

where ${\rho }_{t,b}$ is the corrected TOA reflectance, ${\theta }_{Se}$ = local sun elevation angle. The next step is to correct with respect to at surface reflectance, which is done as follow:

$${\rho }_{s,b}=\frac{{\rho }_{t,b-C_b\left(1-{\tau }_{in,b}\right)}}{{\tau }_{in,b}{\tau }_{out,b}}$$

(9)

where ${\rho }_{s,b}$ is at-satellite reflectance for band b, and $C_b$ is a band b specific given constant. ${\tau }_{in,b}$ and ${\tau }_{out,b}$ are narrowband transmittances for incoming solar radiation and for surface reflected shortwave radiation.

For more details on the atmospheric correction process one may refer to Tasumi et al.’s method [18].

2.3. Using METRIC

The computation of the evaporative fraction is based on using Equation (10), which in turn depends on the energy balance Equation (11):

$$W=\frac{\lambda E}{R_n-G_0}$$

(10)

$$\lambda E= R_n-G_0-H$$

(11)

where $W$ is the evaporative fraction unit less, H is the sensible heat flux in $\left({\mathrm{Wm}}^{-2}\right)$, $\lambda E$ latent heat flux in $\left({\mathrm{Wm}}^{-2}\right)$, R_n is the net radiation in $\left({\mathrm{Wm}}^{-2}\right)$, and G₀ is the soil heat flux in $\left({\mathrm{Wm}}^{-2}\right)$. The equation can be solved using METRIC, which was developed by Allen et al. [13]. METRIC is a sort of “hybrid” between pure remotely sensed energy balance and weather-based evapotranspiration methods. Where energy balance is calculated from a satellite image which delivers spatial information that includes the available energy, and the sensible heat fluxes for a large area.

In this research, METRIC is used as part of the new model in order to calculate evaporative fraction and actual evapotranspiration. Both values are compared with the ones computed from the measurements collected by the Bowen Ratio stations.

METRIC’s foundation, principles and techniques are initially based on SEBAL [20,21]. However, METRIC differs from SEBAL in the way it works on solving Equation (10) by calculating each variable separately such that net radiation is calculated based on the following Equation (12)

$$R_n=RS\downarrow -\alpha RS\downarrow +RL\downarrow -RL\uparrow -\left(1-E_0\right)RL\downarrow$$

(12)

where R_n is the net radiation in $\left({\mathrm{Wm}}^{-2}\right)$ where RS↓ is the incoming shortwave radiation in $\left({\mathrm{Wm}}^{-2}\right)$, α is the surface albedo (dimensionless), RL↓ is the arriving longwave radiation in $\left({\mathrm{Wm}}^{-2}\right)$, RL↑ is the emitted longwave radiation in $\left({\mathrm{Wm}}^{-2}\right)$, and E_o is the surface thermal emissivity.

$$H=\rho C_p\frac{\left(T_a-T_s\right)}{r_H}$$

(13)

H is the sensible heat flux in $\left({\mathrm{Wm}}^{-2}\right)$, $T_a$ and $T_s$ are the air and surface temperatures in (°C), $\rho$ is the air density (kg m⁻³) and $C_p$ heat capacity of air, (Kj kg⁻¹) are constants, and $r_H$ is the transfer resistance (s m⁻¹), which depends on wind speed and surface characteristics.

G₀ = (0.05 + 0.10e^−0.521LAI) R_n   (LAI ≥ 0.5)

(14)

G₀ = max (0.4 H, 0.15 R_n)    (LAI < 0.5)

(15)

where G₀ is the soil heat flux in $\left({\mathrm{Wm}}^{-2}\right)$ and LAI is the leaf area index unit. For more details about the above equations one may refer to the METRIC manual [22].

2.4. A New Optimization Technique

Since the estimation of crop yield depends on few satellite images there will be some periodic gaps between the date of planting and harvesting crops. For this reason, the information obtained from each image is used to create an optimization model to compensate for missing data.

Specific optimization models analyze observations from the real world, such as satellite images, in order to find a solution for data gaps. In this research these models were used to enhance data availability and to predict needed data in a specific information process known as biomass yield estimation. The optimization model finds global optimal solutions to fix the problem of data gaps with minimum error.

Normally, nonlinear curve-fitting (data-fitting) problems are solved using the least-squares method [23]. However, the data extracted from satellite images represent specific crop biomass which is considered a complex non polynomial (NP) type of problem and cannot be fitted to a simple linear equation. For this reason, there is a need for a more reliable method that can solve the problem and which can provide an optimal solution.

The objective is to fit a curve to the crop yield data using unconstrained nonlinear optimization algorithm such as Trust-Region Methods for Nonlinear Minimization [24]. Trust-region methods are efficient, and can solve easily ill-conditioned problems. The basic idea is to approximate a function F with a simpler function f, which reasonably reflects the behavior of function F in a neighborhood N around the point x. This neighborhood is the trust region, a trial step s is computed by minimizing over.

$$N ({}_s{}^{Min}\left\{f\left(s\right), s\in N\right\})$$

(16)

Biomass data were calculated using the Monteith model and satellite images with different acquisition dates. The data with acquisition dates were used in the optimization process to create the new biomass model. At the end, a new biomass yield equation was obtained (Equation (17)) which can compute the total accumulated biomass from the planting phase until crop maturity phase.

$$NBio={{\displaystyle \sum }}_{i=1}^{Nd}f\left(x_i\right)$$

(17)

where NBio is the new total estimated biomass in kg/ha, Nd is the number of days, $x_i$ is the day i during the crop growth progress, and $f\left(x_i\right)$ is the new optimization model that can be created from the estimated crop yield of the available Landsat image. The advantage of the new method is the ability to show the complete crop growth period even with the lack of a complete season of remote sensing images. Figure 3 shows in detail the steps needed to create the optimization model from few Landsat images.

3. Results

To conduct experiments it is required from the user to decide on the type of data, crop and the season length. The decision should be based on the availability of data and other information related to crop requirements. In case of any data gap the optimization model can be used improve the efficiency crop yield estimation process.

In this paper, potato crop yield, which has two different cycles, was estimated, but the majority of farmers plant potato at the end of winter/early spring and harvest in mid-summer. The potato leaves appears almost 20 days after seeding. Because the weather in Lebanon during winter and spring is mainly rainy and cloudy, it was very hard to get images during the months of March (the date potato was planted) and April. Seven satellite images were collected between the months of April, May, June, and July from the date of potato seeding (20 to 30 days after) and near the end of the growth stage (10 to 30 days before harvesting). Landsat 7 images were corrected and data gaps were filled using the methods explained above. Moreover, at-surface reflectance was calculated by applying atmospheric corrections at-satellite and at-surface reflectance.

3.1. Evaporative Fraction Calculation

Evaporative fraction and actual evapotranspiration (Eta) maps are created using the METRIC model as explained in the method section. The map covers seven days in April, May, June and July.

The values of the evaporative fraction are less than 1, except in a few cases where net radiation was very large compared to sensible and soil heat fluxes ($R_n\gg H \mathrm{and} R_n\gg G_0$) or when H was zero because air and surface temperature are equal or the difference was negligible.

To calibrate and to validate the accuracy of energy balance model (METRIC), some recorded climatic data from two Bowen Ratio stations were used to evaluate actual evapotranspiration and evaporative fraction estimated maps. In addition, records on irrigation use were collected from farmers in specific crop fields. These records were used to calibrate and validate METRIC products (

Table 1. Gross irrigation records in a crop field.

Interval Date	# of Days of Irrigation	Break Days	# of Hours of Irrigation	Volume/mm
11–17 May 2015	7	1	6	6
19–25 May 2015	7	2	7	6
28 May to 3 June 2015	7	0	8	6
4–10 June 2015	7	0	8	6
11–17 June 2015	7	2	8	6
20–26 June 2015	7	0	8	6
29 June to 5 July 2015	7	2	7	6
8–13 July 2015	7	3	6	6

# means number.

). Figure 4a–f show the actual evapotranspiration maps and evaporative maps for a large area in Bekaa valley in the 5 May, 8 and 24 June 2015. The following data (

Table 2. Actual evapotranspiration (ETa) from remote sensing images and Bowen station data.

Date	Location	Bowen Ratio ETa (mm)	Remote Sensing ETa (mm)	Absolute Percentage Error
8 June 2015	Qab-Elias	7.54	6.4	15
8 June 2015	Taanayel	5.49	4.6	16
24 June 2015	Qab-Elias	8.29	6.88	17
24 June 2015	Taanayel	6.37	5.6	12
MAPE				15

) show the actual evapotranspiration values obtained from the Bowen Ratio stations (actual values) and the predicted values from the METRIC model for the month of June. The accuracy is computed based on the Mean Absolute Percentage Error (MAPE) [25].

Table 3. Evaporative fraction values for the 8 June 2015 computed using energy balance.

Time	Rn W/m2	G W/m2	Latent Heat Flux W/m2	Evaporative Fraction
9:00	323	18.72	220.083	0.72329
10:00	725	46.08	515.380	0.75912
11:00	1007	91.22	718.760	0.78486
12:00	724	116.88	482.840	0.79530
13:00	889	112.4	621.720	0.80057
14:00	882	105.4	624.120	0.80376
15:00	944	82.62	700.311	0.81301
16:00	1184	54.72	909.990	0.80581
17:00	797	37.34	607.363	0.79952
18:00	205	18.34	146.389	0.78425

shows the computed evaporative fraction $W$ values for the local time between 9:00 and 14:00 on the 8 June 2015. The value of the predicted $W$ = 0.74 was computed by the METRIC model and then compared with the actual W computed from Bowen Ratio station data at time 11:00. The comparison shows that the absolute percentage error was 5%.

MAPE values were computed based on the actual evapotranspiration values from METRIC and the values of irrigation water volume for two days (8 and 24 June 2015). MAPE showed errors values that range between 11% and 18%. After investigating the issue, new facts were discovered. One of these facts is that the farmers usually do not adhere to the rules set by the water authorities in their area; they pump more water and normally record the fixed volume of water set by the authorities.

The estimated hourly evaporative fraction for the 24 June 2015 was also verified. The following data collected by Bowen Ratio station (

Table 4. Evaporative fraction values for the 24 June 2015 computed using energy balance.

Time	Rn W/m2	G W/m2	Latent Heat Flux W/m2	Evaporative Fraction
9	681	21.82	467.78	0.70964
10	597	33.52	415.48	0.73735
11	735	52.28	517.57	0.75801
12	1093	66.32	778.40	0.75817
13	1055	64.5	749.57	0.75676
14	1071	55.3	767.73	0.75587
15	1048	44	756.94	0.75392
16	1047	33.02	763.98	0.75345
17	653	15.56	474.95	0.74509
18	445	1.26	325.92	0.73448

) shows the different parameter values used by the energy balance equation. The value of the evaporative fraction W computed from Bowen station data (highlighted row in Table 4) is compared with the value of W computed by the METRIC model (value = 0.7). The result indicates that the absolute percentage error is 8%.

The results are promising, and they prove that W computed by the METRIC model are comparable to the ones obtained from the Bowen ratio data. In return, the estimated evaporative fraction can be used to solve the crop yield estimation problem.

3.2. Estimating Potato Crop Yield

To estimate potato yield, several Landsat 7 and 8 images were collected. However, the number of available images was not sufficient because a biomass model requires information about the complete crop growth stage. The problem can be solved if daily satellite images such as Modis are obtained, however the information extracted would be useless due to the coarse spatial resolution of this satellite (from 250 m to 1000 m), especially thermal long infrared image. So, the solution is to create an optimization model that can help to solve the issue of data gaps.

After creating the potato map, some statistical information related to biomass was extracted from seven processed satellite images. Although the seeding dates of spring potato are dynamic, the majority of farmers start the seeding process in the month of March, especially in the first two weeks. In addition, most of the farmers harvest potatoes in July.

Equation (18) is the result of fitting a curve to calculated biomass data from different remote sensing images based on using the Trust region algorithm. The parametric fitting of the given data involves finding coefficients (parameters) for one or more models that fit to data. The new biomass equation (Equation (19)) is based on using the new optimization model in Equation (18).

$$F\left(Xdata\right)=489\exp (-{\left(\frac{Xdata-111.1}{14.92}\right)}^2)+283.1\exp (-{\left(\frac{Xdata-80.96}{28.64}\right)}^2)$$

(18)

$$\begin{array}{cc}NewBio={{\displaystyle \sum }}_{i=1}^N\hfill & 489\exp (-{\left(\frac{Xdat{a}_i-111.1}{14.92}\right)}^2)\hfill \\ & +283.1\exp (-{\left(\frac{Xdat{a}_i-80.96}{28.64}\right)}^2)\hfill \end{array}$$

(19)

where Xdata is the day number from the beginning of the crop season to harvesting.

Figure 5 shows potato biomass growth progress after about three weeks from plantation and a few days before the end of the growth stage. The graph was created using the new yield model (Equation (19)). In addition, a map was created using the new crop yield method. Figure 6 shows the crop map, and estimated above ground crop biomass (for the whole season of potato spring 2015).

Data about potato crop yields were collected from 15 different farmers and they were used to verify the estimated results by the model. It is noticeable that there is close agreement between the estimated data in the map and the data collected from farmers (graph in Figure 7). Normally, a potato tuber is equal to up to 80% of the potato biomass [26].

MAPE is calculated for the estimated and real crop yield data and the error was less than 4%. Root-mean-square error (RMSE) was also used to estimate the accuracy of the results. RMSE is frequently used to measure the differences between values predicted by a model or an estimator and the values observed [27]. The value obtained for RMSE is 0.039, which is again an indication of high accuracy [28].

4. Conclusions

The lack of remote sensing data is no longer an obstacle for managers and decision makers in the agriculture sector. It is proved in this research that crop yield estimation can be improved if reliable model exists to help in overcoming all obstacles which hinder crop yield estimation processes. An optimization model created in this research is able to estimate crop yield by increasing remote sensing data availability. The optimization model was created by using a well-known algorithm, Trust-Region Methods for Nonlinear Minimization, that fits available data to an exponential equation. The experimental results proved the accuracy and reliability of the new model in estimating the crop yield in the case of missing remote sensing data (worst case scenario). The potato crop was used to prove the success of the objective of this research. It is planned to improve the optimization model by adding more components, crops, and seasons (agriculture years). This can help in managing many agriculture practices and tasks, and it can make different agriculture information available to decision makers in a highly precise way.