**Corn Grain Yield Estimation from Vegetation Indices, Canopy Cover, Plant Density, and a Neural Network Using Multispectral and RGB Images Acquired with Unmanned Aerial Vehicles**

**Héctor García-Martínez <sup>1</sup> , Héctor Flores-Magdaleno 1,\*, Roberto Ascencio-Hernández <sup>1</sup> , Abdul Khalil-Gardezi <sup>1</sup> , Leonardo Tijerina-Chávez <sup>1</sup> , Oscar R. Mancilla-Villa <sup>2</sup> and Mario A. Vázquez-Peña <sup>3</sup>**


Received: 26 May 2020; Accepted: 2 July 2020; Published: 8 July 2020

**Abstract:** Corn yields vary spatially and temporally in the plots as a result of weather, altitude, variety, plant density, available water, nutrients, and planting date; these are the main factors that influence crop yield. In this study, different multispectral and red-green-blue (RGB) vegetation indices were analyzed, as well as the digitally estimated canopy cover and plant density, in order to estimate corn grain yield using a neural network model. The relative importance of the predictor variables was also analyzed. An experiment was established with five levels of nitrogen fertilization (140, 200, 260, 320, and 380 kg/ha) and four replicates, in a completely randomized block design, resulting in 20 experimental polygons. Crop information was captured using two sensors (Parrot Sequoia\_4.9, and DJI FC6310\_8.8) mounted on an unmanned aerial vehicle (UAV) for two flight dates at 47 and 79 days after sowing (DAS). The correlation coefficient between the plant density, obtained through the digital count of corn plants, and the corn grain yield was 0.94; this variable was the one with the highest relative importance in the yield estimation according to Garson's algorithm. The canopy cover, digitally estimated, showed a correlation coefficient of 0.77 with respect to the corn grain yield, while the relative importance of this variable in the yield estimation was 0.080 and 0.093 for 47 and 79 DAS, respectively. The wide dynamic range vegetation index (WDRVI), plant density, and canopy cover showed the highest correlation coefficient and the smallest errors (R = 0.99, mean absolute error (MAE) = 0.028 t ha−<sup>1</sup> , root mean square error (RMSE) = 0.125 t ha−<sup>1</sup> ) in the corn grain yield estimation at 47 DAS, with the WDRVI index and the density being the variables with the highest relative importance for this crop development date. For the 79 DAS flight, the combination of the normalized difference vegetation index (NDVI), normalized difference red edge (NDRE), WDRVI, excess green (EXG), triangular greenness index (TGI), and visible atmospherically resistant index (VARI), as well as plant density and canopy cover, generated the highest correlation coefficient and the smallest errors (R = 0.97, MAE = 0.249 t ha−<sup>1</sup> , RMSE = 0.425 t ha−<sup>1</sup> ) in the corn grain yield estimation, where the density and the NDVI were the variables with the highest relative importance, with values of 0.295 and 0.184, respectively. However, the WDRVI, plant density, and canopy cover estimated the corn grain yield with acceptable precision (R = 0.96, MAE = 0.209 t ha−<sup>1</sup> , RMSE = 0.449 t ha−<sup>1</sup> ). The generated neural network models provided a high correlation coefficient between the estimated and the observed corn grain yield, and also showed acceptable errors in the yield estimation. The spectral information registered through remote sensors mounted on

unmanned aerial vehicles and its processing in vegetation indices, canopy cover, and plant density allowed the characterization and estimation of corn grain yield. Such information is very useful for decision-making and agricultural activities planning.

**Keywords:** vegetation indices; UAV; neural network; corn plant density; corn canopy cover; yield prediction

### **1. Introduction**

Corn is one of the main food crops for the population, and together with wheat and rice, it is one of the most important cereals in the world. According to the United States Department of Agriculture (USDA), in 2019, a global area of 192.21 million hectares was estimated, with China and the United States being the countries with the largest sown area [1]. In Mexico, 7.4 million hectares of corn was sown in 2018, with a national production of 27.7 million tons and a mean yield of 3.83 t ha−<sup>1</sup> [2]. For the Mexican population, it constitutes the basis for alimentation, providing energy and proteins [3]. In Mexico, 86% of the surface is cultivated under rainfed conditions in plots of less than 5 ha. Some of these producers use the agroecosystem called "milpa", in which several species (beans, pumpkins, and others) are grown on the same plot of corn [4]. Meanwhile, according to Food and Agriculture Organization (FAO), the world population will increase by 35% by 2050, reaching 9100 million inhabitants, mainly in developing countries [5]. To feed this population, food production must sustainably increase by 70%, considering the safety and conservation of natural resources. Some studies indicate that corn yield could decrease in the coming years as a result of anthropogenic climate change [6–9]. There are three main impacts of climate change on agriculture: (a) deterioration in crop yields; (b) effects on production, consumption, and commercialization; and (c) effects on per capita caloric consumption and child nutrition [10]. Corn crop yield is directly related to many factors like the environment, management practices, genotype, and their interactions [11]. The influence of regional climate patterns and large-scale meteorological phenomena can have a significant impact on agricultural production [12]. Genotypes have improved significantly over the years, and important technological developments have been made in the machinery used in management practices. Under these circumstances, yield prediction can be important data for food production, making well-informed and timely economic and management decisions. Correct early detection of problems associated to crop yield factors can help increase the yields and subsequent incomes of farmers. Accurate, objective, reliable, and timely predictions of crop yields in large areas are fundamental to help guarantee the adequate food supply of a nation and assist responsible politicians to make plans and set prices for imports/exports [13].

Yields in corn crops vary spatially and temporally in the plots, according to the conditions present in each site such as weather, altitude, variety, planting density, available irrigation or the amount of rain (water supply), the available nutrients for the plant (soil plus fertilizer), and the date of sowing, which are the main factors that influence the yield of a plot. In recent decades, corn yield has increased as a result of genetic improvement and agronomic management. Increases in plant density and the use of synthetic fertilizers have been the main factors responsible for increases in corn yields. Plant density (number of plants per unit area) is one of the components of grain yield (number of ears per unit area, number of grains per ear, grain weight) that has an impact on the final corn yield [14]; however, its accurate measurement after plant emergence is not practical in large-scale production fields owing to the amount of labor required [15].

Precision agriculture (PA) is a management concept based on observation, measurement, and response to variability of the crops in the field [16]. PA technology allows farmers to recognize variations in the fields and apply variable rate treatments with a much finer degree of precision than before. Identifying spatial and temporal variability within the field shows the potential to support crop management concepts to satisfy most of the growing environmental, economic, market, and public pressures on agriculture [17]. Remote sensing is generally considered one of the most important technologies for precision agriculture and intelligent agriculture; it can monitor many crops and vegetation parameters through images at various wavelengths [18]. With the development of unmanned aerial systems (UAS), their use in remote sensing and precision agriculture offers the possibility of obtaining field data in an easy, fast, and profitable way, resulting in images with high spatial and temporal resolution. The successful adoption of remote sensing based on unmanned aerial vehicles (UAV) depends on changes in sensitivity on vegetation indices (VI) and growth stages [19]. Processing different vegetation indices has been associated with physiological parameters of the plant, such as plant pigments, vigor, aerial biomass, yield estimation, plant physiology, and stress. Xue and Su in 2017 [20] reviewed the developments and applications of 100 vegetation indices (VIs) in remote sensing; some VIs employ a range of reflectance values in a narrow band of the electromagnetic spectrum for more precise measurements correlating them with the grain yield, providing reliable information for yield forecasting [21,22], but they require more technologically advanced sensors. As a low-cost alternative, there are VIs that are obtained from red-green-blue (RGB) images calculated from commercial cameras. These have shown in some studies their ability to predict grain yield, quantify nutrient deficiencies, and measure the impact of diseases [23–25].

The factors that affect crop yields, such as soil, climate, and management, are so complex that traditional statistics cannot give accurate results. Various machine learning techniques have been used for yield prediction, such as decision trees, self-organizing maps (SOMs), multivariate regression, support vector machines, association rule mining, and neural networks [26–31]. As a machine learning tool, the artificial neural network (ANN) is an attractive alternative to process the massive data set generated by production and research in precision agriculture [11]. The models used in machine learning relate crop yield, as an implicit function of input variables, for example, climate variables, soil, and water characteristics, which can be a very complex function, not necessarily linear. Recently, different types of neural networks have been used to predict yield in wheat, soybean, rice, corn, and tomato, using databases of genotypes, environment, management practices, and multispectral images, obtaining acceptable results [32–37].

Studies have been carried out to estimate corn yield, using data obtained from remote sensors and ANN. Fieuzal et al. [38] use multi-temporal data from satellites and radar, as well as a neural network in the estimation of maize yield with an R<sup>2</sup> of 0.69. On the other hand, in another study [39], polarimetric synthetic aperture radar (PolSAR) and neural network data are used in the estimation of corn biomass, obtaining good results with an R = 0.92. Han et al. [40] estimate the aerial biomass of corn from spectral information, plant height, and structural information using data from remote sensors and unmanned aerial vehicles with machine learning regression algorithms, obtaining good results (R<sup>2</sup> = 0.69). Michelon et al. [41] use ANN and chlorophyll readings to estimate corn productivity, resulting in a correlation coefficient of 0.73 in stage V6. Olson et al. [19] related vegetation indices and crop height with maize yield, finding high correlations; another approach uses yield data from parents to predict maize yield in plant breeding using neural networks [42]. Canopy cover, plant density, and vegetation indices calculated from data collected by remote sensors can be used to forecast corn grain yield using neural networks. As far as we know, this is the first study where canopy cover, plant density, and various vegetation indices are estimated from images obtained by UAVs to forecast corn grain yield using a neural network.

In this study, the corn grain yield was estimated through the use of a neural network for a corn crop established under different doses of nitrogen fertilization from multispectral and digital images acquired by sensors mounted on an unmanned aerial vehicle. The images were processed and vegetation indices, canopy cover, and plant density were extracted. A neural network was designed having the following input parameters: vegetation indices (normalized difference vegetation index, NDVI; normalized difference red edge, NDRE; wide dynamic range vegetation index, WDRVI; excess green, EXG; triangular greenness index, TGI; visible atmospherically resistant index, VARI), canopy cover, and plant density.

### **2. Materials and Methods**

### *2.1. Study Site*

An experiment was established in a plot on 5 April 2018 at the Colegio de Postgraduados, Campus Montecillo, located on the Mexico-Texcoco highway, km. 36.5 Montecillo, Texcoco, State of Mexico (19◦27040.58" N, 98◦5408.57" W, 1250 m above sea level). Soil texture is sandy, bulk density of 1.45 g cm−<sup>3</sup> , organic matter of 1.59%, pH of 9.1, and electrical conductivity of 1.72 dS m−<sup>1</sup> . Five nitrogen levels (140, 200, 260, 320, and 380 kg ha−<sup>1</sup> ) were evaluated in a completely randomized block experimental design with four replicates, resulting in 20 experimental units (120 m<sup>2</sup> per unit), using the Asgrow Albatross corn variety. A drip-band irrigation system with drippers at 20 cm separation and a flow of 1.6 L/h was designed to irrigate each experimental unit. The five nitrogen levels were applied to the irrigation water at dose intervals of 30% for the first 40 days, 50% of the dose at 40–80 days, and 20% after 80 days. The sowing was done by hand with a distance of 80 cm between rows and 30 cm between plants, with three seeds per hole. Weeds were controlled manually. A mean temperature of 16.9 ◦C was present during the experiment period, and a total precipitation of 470 mm occurred, with these conditions being very close to the normal conditions.

The experimental units were harvested at the end of September 2018, harvesting the entire area of the experimental unit, so the total weight harvested corresponded to an area of 120 m<sup>2</sup> (0.8 m separation × 6 rows × 25 m long) for all replicates of each treatment in the experiment. The harvested ears and grains were dried at approximately 14% humidity, so grain yield (t ha−<sup>1</sup> ) was calculated according to the following:

$$\text{Gy} = \frac{\text{X}}{120} \times 10,\tag{1}$$

Grain yield per plant (g plant−<sup>1</sup> ) was calculated according to the following:

$$\text{Gyp} = \frac{\text{\AA} \times 1000}{\text{Np}} \,\text{\AA} \tag{2}$$

where X represents grain weight per experimental area (kg m−<sup>2</sup> ), Gy is grain yield, Gyp is grain yield per plant, and Np is the number of plants present in the area.

### *2.2. Acquisition and Analysis of Data from Remote Sensors Mounted on Unmanned Aerial Vehicles*

Four orthomosaics were used in this work according to Table 1, which were acquired by a 72 g Parrot Sequoia camera with a 16 Mpx RGB sensor. It also incorporates four single band sensors with a spectral bandwidth in green (530–570 nm), red (640–680 nm), red edge (730–740 nm), and near infrared (NIR) (770–810 nm) of 1.2 Mpx, and a 35 g solar sensor for real-time correction of lighting differences. The Parrot Sequoia camera was attached to a 3 DR SOLO quadcopter (3D Robotics, Berkeley, CA, USA) with flight autonomy of 20 min, with a horizontal precision range of 1–2 m. Two orthomosaics used in this work were generated from a 20 Mpx DJI RGB sensor integrated into a DJI Phantom 4 quadcopter (DJI, Shenzhen, Guangdong, China), with a load capacity of 1.388 kg and autonomy per flight of 30 min, equipped with a GPS/GLONASS satellite positioning system with a vertical precision range ±0.5 m and horizontal precision range of 1.5 m.



V8: Vegetative stage—eight leaves with collar visible; R0: Reproductive stage—anthesis.

The four flights were carried out in a period from 21 May to 22 June 2018, at a flight height of 30 m. The images were acquired with 80% overlap and 80% sidelap, and the orthomosaic was generated with the Pix4D software (Pix4D SA, Lausanne, Switzerland) using structure-from-motion. The structure-from-motion processing technique searches for features in individual images that are used to relate to features that match between overlapping images, called keypoints; using these key points, camera parameters can be calibrated to external parameters (such as position, scale, and image orientation) and a point correlation (matching) is performed based on the characteristics, identifying and relating similar characteristics between images in common areas or overlapping areas. The calculated 3D position of the matched points is densified and textured with the corresponding images, from which the orthomosaic is generated by projecting each textured pixel onto a 2D plane [43,44]. Ground control points (GCP) were added for orthomosaic georeferencing according to Figure 1a. The required GCP density depends on the required project accuracy, network geometry, and the quality of images [45,46]. Six marks were placed as ground control points (GCPs) distributed in the experimental plot for each flight. A Global Navigation Satellite System (GNSS) Real Time Kinematics (RTK) V90 PLUS Hi-Target system (Hi-Target Surveying Instrument Co., Ltd., Guangzhou, China) was used to record the center of each control point with RTK precision, Horizontal: 8 mm + 1 Part Per Million (ppm) Root Mean Square (RMS), and Vertical: 15 mm + 1 ppm RMS. Table 1 shows the summary of the images acquired in each flight and the area covered for the four orthomosaics at 47 and 79 days after sowing (DAS). Moreover, the corn growth stages for each flight are shown (V8 and R0). Generated orthomosaics showed a ground sample distance (GSD) of 0.49 cm/pixel and a model RMS error of 5.6 cm for the 20 Mpx RGB sensor. The multispectral images obtained by the Sequoia camera were processed in the Pix4D software, which integrates the camera's light sensor to correct the estimated reflectance and performs a radiometric calibration using a calibration panel image [47]. The multispectral Ag template and the calibration panel by Airinov provided by the camera were used, defining the known reflectance values for each spectral band of the panel equal to 0.171, 0.216, 0.268, and 0.372 for green, red, red edge, and NIR, respectively. The obtained orthomosaics showed a ground sampling distance of 2.15 cm/pixel.

Six vegetation indices were estimated based on the generated orthomosaics. Three vegetation indices based on the reflectance of the visible spectrum (RGB): TGI (triangular greenness index), EXG (green extraction index), and VARI (visible atmospherically resistant index), according to Table 2, and three multispectral indices based on the reflectance in the near infrared, the red, and the red edge bands: NDVI (normalized difference vegetation index), NDRE (green extraction index), and WDRVI (wide dynamic range vegetation index). The calculation of the vegetation indices was done using the raster calculator module of the Open Source program licensed under the GNU—General Public License of the Geographic Information System (QGIS) software.

6 of 27

(**a**)

**Figure 1.** (**a**) Experimental plot and polygons sampled in the extraction of vegetation indices, canopy cover, and plant density; (**b**) binary image resulting from the classification of vegetation and soil; (**c**) polygon selected for the sampling of the canopy cover and plant density. **Figure 1.** (**a**) Experimental plot and polygons sampled in the extraction of vegetation indices, canopy cover, and plant density; (**b**) binary image resulting from the classification of vegetation and soil; (**c**) polygon selected for the sampling of the canopy cover and plant density.

Six vegetation indices were estimated based on the generated orthomosaics. Three vegetation indices based on the reflectance of the visible spectrum (RGB): TGI (triangular greenness index), EXG (green extraction index), and VARI (visible atmospherically resistant index), according to Table 2, and three multispectral indices based on the reflectance in the near infrared, the red, and the red edge bands: NDVI (normalized difference vegetation index), NDRE (green extraction index), and WDRVI



\* WDRVI with an α coefficient value of 0.1 presented a good relationship with corn canopy cover [53].

The variables r, g, and b of the EXG index are normalized values of the red, green, and blue channels, respectively, according to

$$\mathbf{r} = \frac{\mathbf{R\_N}}{\mathbf{R\_N} + \mathbf{G\_N} + \mathbf{B\_N}}, \ \mathbf{g} = \frac{\mathbf{G\_N}}{\mathbf{R\_N} + \mathbf{G\_N} + \mathbf{B\_N}}, \ \mathbf{b} = \frac{\mathbf{B\_N}}{\mathbf{R\_N} + \mathbf{G\_N} + \mathbf{B\_N}} \tag{3}$$

$$\mathbf{R\_N} = \frac{\mathbf{R}}{\mathbf{R\_{max}}}, \text{ G\_N} = \frac{\mathbf{G}}{\mathbf{G\_{max}}}, \text{ B\_N} = \frac{\mathbf{B}}{\mathbf{B\_{max}}} \tag{4}$$

where RN, GN, and B<sup>N</sup> are the normalized values of each band; R, G, and B are the non-normalized values of the red, green, and blue channels, respectively; and Rmax = Gmax = Bmax are the maximum digital numbers for each channel (255 on the 0–255 scale).

TGI index was normalized to have a value in the range of the others indices, between 0 and 1, according to the following:

$$\mathbf{Y}\_{\text{nor}} = \frac{\mathbf{Y} - \mathbf{Y}\_{\text{min}}}{\mathbf{Y}\_{\text{max}} - \mathbf{Y}\_{\text{min}}},\tag{5}$$

where Ynor is the normalized index value, Y is the index value without normalizing, Ymin is the minimum index value, and Ymax is the maximum index value.

### *2.3. Plant Count, Determination of Plant Density, and Estimation of Canopy Cover*

The corn plants present in the study area were counted digitally according to García et al. [54]. The orthomosaics were transformed from the RGB color model to the CIELab model, where we only worked on the channel a\* to extract the vegetation. Corn plant samples were selected at 47 days after planting, twelve of them for each experimental area. Through normalized cross-correlation with the normxcorr2 command in Matlab (Mathworks, Natick, MA, USA), the corn plants present in the binary image were classified and, using an image component labeling technique, assigning a label (1 ... i), the pixels of each plant were grouped. The area of the plants selected in the correlation was used as a criterion, counting plants with 20% less than the minimum area of the smallest selected sample. The digitally counted plants were registered in a table. The plant density (plants m−<sup>2</sup> ) in each experimental area was calculated according to the following:

$$\text{Pd} = \frac{\text{Np}}{\text{Ea}} \text{.} \tag{6}$$

where Pd is the plant density in plants/m<sup>2</sup> , Np is the number of plants present in the sampled area, and Ea is the sampled area in square meters.

The canopy cover was calculated from the generated TGI vegetation index, using the Open Source program licensed under the GNU—General Public License of the Geographic Information System (QGIS) and System for Automated Geoscientific Analyses (SAGA) [55]. In QGIS, the K-Means Clustering for Grids module was used to group the pixels into two classes: soil and vegetation, through the combined method of minimum distance and simple scaling [56]. The K-means method of image segmentation is a technique that does not use the histogram for the segmentation process, so it is not affected by the noise introduced in the image, evaluating and grouping the different pixels in similar data sets, being suitable for large data sets. On the other hand, the K-means approach in image segmentation is fast and efficient in terms of computational cost [57,58]. In the classification of pixels into two classes in QGIS, a binary raster was generated, where pixels with a value of 0 belong to soil and pixels with a value of 1 belong to vegetation, or vice versa. Thus, the estimation of the percentage of canopy cover was obtained by the following:

$$\text{Cc } = \frac{\text{NP}}{\text{TP}} \tag{7}$$

where Cc is the canopy cover, NP are the pixels corresponding to vegetation per unit area, and TP is the total pixels per unit area.

Sampling was carried out in four 11 m<sup>2</sup> polygons for each replicate, according to Figure 1a. In each of them, the canopy cover (Figure 1b,c) and the plant density for the two flight dates were calculated.

### *2.4. Vegetation Pixels' Segmentation and Extraction of the Indices Values*

The images acquired through the sensors and the orthomosaics generated from them contain reflectance corresponding to soil and vegetation, but only the reflectance corresponding to vegetation was of interest to us. With the eCognition software and the object-based image analysis (OBIA) technique, multi-resolution segmentation was realized. This combines pixels or adjoining objects present in the image based on spectral and shape criteria. It also works according to the scale of the objects present; large scale results in large objects, and vice versa [59]. The homogeneity criterion measures how homogeneous or heterogeneous the object present in the image is, calculated from a combination of object color and shape properties. Color homogeneity is based on the standard deviation of the spectral colors. Shape homogeneity is based on the deviation of a compact (or smooth) shape and can have a value of up to 0.9. In the segmentation process, a 10-scale value was used and the homogeneity criteria like shape and compactness were 0.1 and 0.5, respectively. For the classification of the segmented pixels in the image, three classes present in the image (Figure 2a,b) were determined: soil, vegetation, and shadows. Supervised classification of the nearest neighbor was used to select pixel samples for each class; they were compared with the other objects present with respect to the mean and standard deviation of the supervised sample [60]. *Agriculture* **2020**, *10*, x FOR PEER REVIEW 9 of 27

**Figure 2.** Segmentation and classification of objects present in the orthomosaics for the extraction of pixels belonging to and classified as corn plants: (**a**) and (**b**) processing and analysis of the crop image based on objects for three classes: corn, shadow, and soil. **Figure 2.** Segmentation and classification of objects present in the orthomosaics for the extraction of pixels belonging to and classified as corn plants: (**a**) and (**b**) processing and analysis of the crop image based on objects for three classes: corn, shadow, and soil.

extracted for 76 sampled polygons, 4 per replicate of the 5 treatments as the average of the mean values of all the pixels classified as corn contained in each polygon. The sample polygon of each

(**a**) (**b**)

**Figure 3.** Segmentation and classification of objects present in the orthomosaics for the extraction of pixels belonging to and classified as corn plants: (**a**) image of the crop without segmentation and classification; (**b**) image of the crop with a shape of polygons generated from the segmented and

A multilayer perceptron type feed-forward neural network was developed using Matlab software (Mathworks, Natick, MA, USA). A feed-forward neural network contains multiple neurons arranged in layers: input, hidden, and output. The information goes in one direction, so there are no cycles or loops [61]. The input layer receives values of the input variables (x1, x2, ..., xn); between the

replicate on the edge of the crop was discarded to eliminate edge effects.

classified pixels of corn plants.

*2.5. Development and Training of the Feed-Forward Neural Network* 

Once the pixels present in the image were classified, we exported the pixels classified in the corn

(a)

based on objects for three classes: corn, shadow, and soil.

REVIEW

Once the pixels present in the image were classified, we exported the pixels classified in the corn class (Figure 3a,b) as polygons in shape file format, with the QGIS software. The mean values of the indices calculated for each polygon of corn class were extracted through the zone statistics complement; the values were stored in a shape file. The mean values obtained from the indices were extracted for 76 sampled polygons, 4 per replicate of the 5 treatments as the average of the mean values of all the pixels classified as corn contained in each polygon. The sample polygon of each replicate on the edge of the crop was discarded to eliminate edge effects. Once the pixels present in the image were classified, we exported the pixels classified in the corn class (Figure 3a,b) as polygons in shape file format, with the QGIS software. The mean values of the indices calculated for each polygon of corn class were extracted through the zone statistics complement; the values were stored in a shape file. The mean values obtained from the indices were extracted for 76 sampled polygons, 4 per replicate of the 5 treatments as the average of the mean values of all the pixels classified as corn contained in each polygon. The sample polygon of each replicate on the edge of the crop was discarded to eliminate edge effects.

**Figure 2.** Segmentation and classification of objects present in the orthomosaics for the extraction of

*Agriculture* **2020**, *10*, x FOR PEER

9 of 27

(**a**) (**b**)

(b)

**Figure 3.** Segmentation and classification of objects present in the orthomosaics for the extraction of pixels belonging to and classified as corn plants: (**a**) image of the crop without segmentation and classification; (**b**) image of the crop with a shape of polygons generated from the segmented and classified pixels of corn plants. **Figure 3.** Segmentation and classification of objects present in the orthomosaics for the extraction of pixels belonging to and classified as corn plants: (**a**) image of the crop without segmentation and classification; (**b**) image of the crop with a shape of polygons generated from the segmented and classified pixels of corn plants.

### *2.5. Development and Training of the Feed-Forward Neural Network 2.5. Development and Training of the Feed-Forward Neural Network*

A multilayer perceptron type feed-forward neural network was developed using Matlab software (Mathworks, Natick, MA, USA). A feed-forward neural network contains multiple neurons arranged in layers: input, hidden, and output. The information goes in one direction, so there are no cycles or loops [61]. The input layer receives values of the input variables (x1, x2, ..., xn); between the A multilayer perceptron type feed-forward neural network was developed using Matlab software (Mathworks, Natick, MA, USA). A feed-forward neural network contains multiple neurons arranged in layers: input, hidden, and output. The information goes in one direction, so there are no cycles or loops [61]. The input layer receives values of the input variables (x1, x2, ..., xn); between the input layer and the hidden layer are the weights (w1, w2, ..., wn) that represent the memory of the network; and the output neurons return the output (y1, y2, .., yn) through the sum of all the inputs multiplied by the weights plus a value associated with the neuron, called bias (b) [62].

$$\mathbf{y} = \sum\_{\mathbf{i}=0}^{\mathbf{i}=\mathbf{n}} \mathbf{W}\_{\mathbf{i}} \mathbf{x}\_{\mathbf{i}} + \mathbf{b}\_{\mathbf{i}} \tag{8}$$

The outputs of the neurons before passing to the other nodes are transformed through an activation function f(x) to limit the output of the neuron. The precision in the prediction by a neural network is related to the type of activation function used; the non-linear activation functions are the most used (Sigmoid, Hyperbolic Tangent, Rectified liner Unit (ReLU), Exponential Linear Unit (ELU), Softmax) [63,64], so the final output (a) will be as follows:

$$\mathbf{a} = \mathbf{f} \left( \sum\_{i=0}^{i=n} \mathbf{W}\_i \mathbf{X}\_i + \mathbf{b} \right) \tag{9}$$

As input parameters or variables, we used the mean value of the vegetation indices, the canopy cover, and the plant density, as well as the yield in tons per hectare as labels and output parameters. In total, 70% of the total data entered into the neural network was used for training, 15% for validation, and 15% for testing using a random partition of the entire data set. For each combination of input variables, the neural network was trained ten times using ten different random data sets. The training

algorithm used was Levenberg–Marquart [65,66], which is a training function that updates the values in the weights and bias according to the Levenber–Marquart optimization. Tan-sigmoid was used as activation function, as it defines the behavior of the neurons in charge of calculating the degree or state of activation of the neuron, according to the total input function [67]. The descending gradient with momentum type was used as a learning rule, as it calculates the change in weight (W) for a given neuron. On the basis of the weights of the input neurons and the error (E), the weight W, the learning rate (LR), and the moment constant (MC) are calculated, according to the descent of the gradient with momentum. The objective of training the neural network is to minimize the resulting errors for a training data set by adjusting the weights (W).

### *2.6. Statistical Analysis*

For the performance analysis of the estimated yields by the neural network, estimated yields were compared to the observed yields, estimating the root mean square error (RMSE), the mean absolute error (MAE), and the correlation coefficient (R) in the training, validation, testing, and total data entered into the neural network. The correlation coefficient between the observed yields and the vegetation indices, the plant density, and the canopy cover was also estimated.

$$\mathbf{E\_s} = (\mathbf{EY} - \mathbf{OY}),\tag{10}$$

$$\text{MAE} = \frac{1}{\text{N}} \sum\_{\mathbf{i}=1}^{\text{N}} \left| (\mathbf{E}\_{\mathbf{s}})\_{\mathbf{i}} \right| \tag{11}$$

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{N} (\text{EY} - \text{OY})^2}{N}},\tag{12}$$

where N is the number of the sampled polygon, (Es)<sup>i</sup>  is the absolute value of Es, EY is the estimated yield, and OY is the observed yield.

### *2.7. Variables' Relative Importance in the Yield Estimation Using Garson's Algorithm*

Neural network models are difficult to interpret and it is difficult to identify which predictors are the most important and how they relate to the property being modeled. The weights connecting neurons in a neural network are partially analogous to the coefficients in a generalized linear model, and these weights combined with their effects on model predictions represent the relative importance of predictors in their associations with the variable being predicted [68]. Garson's algorithm discriminates the relative importance of the predictor variables for a single response variable, that is, the force with which a specific variable explains the predicted variable is determined by identifying all the weighted connections between the nodes of interest, and identifying the weights that connect to the specific input node that passes through the neural network and reaches the variable that is being predicted. This allows to obtain a unique value for each descriptor variable used in the neural network model [69,70]. To know the predictor variables' relative importance (NDVI, NDRE, WDRVI, EXG, TGI, VARI, canopy cover, and density) in the corn yield estimation, Garson's algorithm was used.

### **3. Results and Discussion**

### *3.1. Corn Plant Vegetation Indices*

The segmentation and classification of objects into soil and vegetation allowed to filter mean values of the vegetation indices with respect to the pixels belonging to corn plants. Table 3 shows the mean values computed through segmentation and classification of the objects present in the image for two dates in the corn crop development. Differences are observed in the mean values of the vegetation indices for 47 and 79 days after sowing. The classification of the pixels captured in the sensors into two classes allowed us to extract the values of the indices belonging to corn, so the result was a mean

value of the vegetation. Contrarily, if the classification of the pixels into classes was not carried out, there would be a bias when considering pixels belonging to soil and shadows in the mean computation for the polygons taken as a sample, resulting in a mixed value of corn plants and soil pixels. Vegetation indices values increased for 79 DAS with respect to the values shown for 47 DAS, except for the VARI index, which did not present a wide variation with the crop development. Similar results were shown in the studies carried out by Gitelson et al. [53,71]. Vegetation indices did not show differences for the nitrogen treatments used; similar results are reported by Olson et al. [19].

**Table 3.** Mean values of the vegetation indices and yield for each replicate of the treatments under different degrees of nitrogen fertilization.

