Using Artiﬁcial Neural Networks and Remotely Sensed Data to Evaluate the Relative Importance of Variables for Prediction of Within-Field Corn and Soybean Yields

: Crop yield prediction prior to harvest is important for crop income and insurance projections, and for evaluating food security. Yet, modeling crop yield is challenging because of the complexity of the relationships between crop growth and predictor variables, especially at the ﬁeld scale. In this study, an artiﬁcial neural network (ANN) method was used: (1) to evaluate the relative importance of predictor variables for the prediction of within-ﬁeld corn and soybean end-of-season yield and (2) to evaluate the performance of the ANN models with a minimal optimized variable dataset for their capacity to predict corn and soybean yield over multiple years at the within-ﬁeld level. Several satellite derived vegetation indices (normalized di ﬀ erence vegetation index—NDVI, red edge NDVI and simple ratio—SR) and elevation derived variables (slope, ﬂow accumulation, aspect) were used as crop yield predictor variables, hypothesizing that the di ﬀ erent variables reﬂect di ﬀ erent crop and site conditions. The study identiﬁed the SR index and the slope as the most important predictor variables for both crop types during two training and testing years (2011, 2012). The dates of the most important SR images, however, were di ﬀ erent for the two crop types and corresponded to their critical crop developmental stages (phenology). The relative mean absolute errors were overall smaller for corn compared to soybean: all of the 2011 corn study ﬁelds had errors below 10%; 75% of the ﬁelds had errors below 10% in 2012. The errors were more variable for soybean. In 2011, 37% of the ﬁelds had errors below 10%, while in 2012, 100% of the ﬁelds had errors below 20%. The results are promising and can provide yield estimates at the farm level, which could be useful in reﬁning broader scale (e.g., county, region) yield projections. on a tile-drained for national


Introduction
The ability to predict crop yield prior to harvest is important for crop income and insurance projections, as well as for evaluating food security at local to global scales. Furthermore, models that spatially predict within-field yields can be used to characterize the environmental and management

Study Area
The study was conducted on a tile-drained experimental watershed in eastern Ontario, Canada (45.26 N, 75.18 W; Figure 1). Corn (Zea mays L.) and soybean (Glycine max L.) are the two dominant crops in this region of Canada and in 2016, Ontario accounted for~60% of national grain corn acreage and~50% of national soybean acreage [20,21]. Corn and soybean are the two primary livestock and cash crops grown in the experimental watersheds. These crops are typically planted in May and harvested in September to October. In the study area, both crops were planted at the same time every year. Crop row spacing was around 50 cm for soybean and 60 cm for corn. Soybean and corn plant densities were generally 354,000 and 64,200 plants/ha, respectively. Fall moldboard plowing and spring cultivation using chisel-style implements were the main tillage practices employed on the site. Fertilizer generally consisted of the broadcast application of granular urea prior to planting and a granular starter application. Fertilizer rates were generally~170 kg N/ha for corn and~3 kg N/ha for soybean, depending on soil and crop requirements. More details about the study area and agriculture practices-including tiling systems-are described by Sunohara et al. [22].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 33 The study was conducted on a tile-drained experimental watershed in eastern Ontario, Canada (45.26 N, 75.18 W; Figure 1). Corn (Zea mays L.) and soybean (Glycine max L.) are the two dominant crops in this region of Canada and in 2016, Ontario accounted for ~60% of national grain corn acreage and ~50% of national soybean acreage [20,21]. Corn and soybean are the two primary livestock and cash crops grown in the experimental watersheds. These crops are typically planted in May and harvested in September to October. In the study area, both crops were planted at the same time every year. Crop row spacing was around 50 cm for soybean and 60 cm for corn. Soybean and corn plant densities were generally 354,000 and 64,200 plants/ha, respectively. Fall moldboard plowing and spring cultivation using chisel-style implements were the main tillage practices employed on the site. Fertilizer generally consisted of the broadcast application of granular urea prior to planting and a granular starter application. Fertilizer rates were generally ~170 kg N/ha for corn and ~3 kg N/ha for soybean, depending on soil and crop requirements. More details about the study area and agriculture practices-including tiling systems-are described by Sunohara et al. [22].

Field Data
Corn and soybean yield data were obtained from cooperating farmers for the years 2011, 2012 and 2016, matching the available satellite images. The data were collected during harvest using a John Deere combine harvester equipped with a GREENSTAR™ Yield Monitor System (John Deere, Moline, Illinois). Post-processing was done using Field Operations Viewer (version 2003.11.01; MapShots, San Francisco, California) and Yield Editor (Version 2.0.7; USDA-ARS Cropping Systems and Water Quality Research, Columbia, Missouri) and resulted in a dense point dataset with yields expressed in kg/ha. Points with yield values equal to the mean ± three times the standard deviation were considered outliers and eliminated. The summary statistics of corn and soybean yields for 2011, 2012 and 2016 are given in Table 1 and Figure 2a,b. The field yield data was used to train and test the ANN models.

Field Data
Corn and soybean yield data were obtained from cooperating farmers for the years 2011, 2012 and 2016, matching the available satellite images. The data were collected during harvest using a John Deere combine harvester equipped with a GREENSTAR™ Yield Monitor System (John Deere, Moline, Illinois). Post-processing was done using Field Operations Viewer (version 2003.11.01; MapShots, San Francisco, California) and Yield Editor (Version 2.0.7; USDA-ARS Cropping Systems and Water Quality Research, Columbia, Missouri) and resulted in a dense point dataset with yields expressed in kg/ha. Points with yield values equal to the mean ± three times the standard deviation were considered outliers and eliminated. The summary statistics of corn and soybean yields for 2011, 2012 and 2016 are given in Table 1 and Figure 2a,b. The field yield data was used to train and test the ANN models.

LIDAR and Topographic Derivatives
Based on the reported importance of micro-topographic attributes on crop growth properties [18,19,23], we included Lidar-derived variables to represent the micro-topography of the fields over this low topographic relief study area. Lidar data were collected over the experimental site by GeoDigital (Sandy Springs, GA) in November 2011, and a 1-m spatial resolution digital terrain elevation dataset was derived from the Lidar point cloud using LAStools (rapidlasso GmbH, Gilching, Germany). Slope (Figure 3), aspect (not shown) and flow accumulation ( Figure 4) were derived from the elevation dataset using the spatial analyst ArcGIS toolset within Advangeo ® Prediction (v 2.0 for ArcGIS 10.1; Beak Consultants GmbH, Freiberg, Germany). elevation dataset was derived from the Lidar point cloud using LAStools (rapidlasso GmbH, Gilching, Germany). Slope (Figure 3), aspect (not shown) and flow accumulation ( Figure 4) were derived from the elevation dataset using the spatial analyst ArcGIS toolset within Advangeo ® Prediction (v 2.0 for ArcGIS 10.1; Beak Consultants GmbH, Freiberg, Germany).  The flow accumulation variable was calculated as the number of cells that flowed into each downslope cell ( Figure 4). Higher relative flow accumulation numbers indicate areas of concentrated flow and can be used to determine stream channels, swales or depressions. Cells with a value of 0 indicate local topographic highs with no upstream catchment area, such as knolls and ridges (see https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst). The flow accumulation variable was calculated as the number of cells that flowed into each downslope cell (Figure 4). Higher relative flow accumulation numbers indicate areas of concentrated flow and can be used to determine stream channels, swales or depressions. Cells with a value of 0 indicate local topographic highs with no upstream catchment area, such as knolls and ridges (see https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst).

Optical Remote Sensing Derived Data
A total of 18 RapidEye images (5 m spatial resolution) were collected over the three study years ( Table 2). The RapidEye constellation consists of five satellites that are equipped with identical sensors and located within the same orbital plane. This setting allows the constellation to acquire high-resolution, large area image data on a daily basis (PlanetLabs Inc., https://www.planet.com/).  All images were georeferenced to the universal transverse Mercator coordinate system UTM 18N (with the world geodetic system-WGS-84 ellipsoid) and atmospherically corrected using ATCOR2 (PCI Geomatica v10.3, PCI Geomatics enterprises Inc, Ontario, Canada).
Considering the plethora of static and transient factors that can affect end-of-season yield in a field, we chose a set of satellite-derived vegetation and site spectral indices (Table 3) to reflect plant nitrogen content, crop health (normalized difference vegetation index-NDVI and red edge NDVI-NDVIre) and crop canopy structure (Simple Ratio-R). Table 3. Summary of spectral indices and their respective equations and references.

Spectral Indices Equation Reference
Normalized difference vegetation index (NDVI)

ANN Model
A review study on MLAs listed ANNs as one of the most successful MLAs for yield prediction, together with other methods like support vector regression, M5 prime regression trees, k-nearest neighbor and deep learning [24]. Our study used a parsimony approach to explore the importance of predictor variables and the performance of ANN models with a minimal variable set, in predicting crop yield within fields.
All predictor variables were clipped to the study area and the final spatial resolution of the model was set to the spatial resolution of the spectral indices from RapidEye (5 m). The 5-m resolution was considered via long-term field observations [22], to adequately capture most important micro-topographic features at the study site. A summary of all predictor variables used in the ANN analysis for each year is given in Table 4. We used the Advangeo ® Prediction Software, which uses a multilayer perceptron (MLP) ANN method to predict either classes or numeric values. Advangeo ® Prediction comes preloaded with defaults that have proven effective for many different applications [25]. In this study, we parsimoniously applied default parameters. An overview of the model is given in Figure 5. The input layer receives the controlling parameters and the neurons of the hidden layer(s) and the output layer processes the weighted signals from the neurons of its previous layer and calculates an output value applying an activation function. The default parameters were: network topology (fully connected input layer and one hidden layer with three hidden neurons), activation function (for hidden and output layers: sigmoid with a steepness of 0.5), learning algorithm (derivative of back propagation algorithm), weight initialization ("initialize" algorithm of Widrow and Nguyen) and predefined stop parameters (number of epochs = 100 and mean square error -MSE border = 0.001). The input data consisted of different spatial layers of vegetation indices (SR, NDVI, NDVIre), Lidar derivatives (slope, aspect, flow accumulation), and crop type (corn, soybean). All the data were normalized to a scale of 0 to 1 (linear scaling) for use in the ANN model. The results were scaled back to the original data scales for error analysis.
The ANN model was trained and tested using 2011 data (fields were divided into 50% training vs. 50% testing), after which it was applied to 2012 and 2016 data. The same was done for 2012, where the model was trained and tested using 2012 data, and applied to 2011 and 2016 data. For each training year (2011, 2012), three main model scenarios were explored: (1) the network was trained with yield data from corn and soybean combined for the prediction of the within-field yield of the two crops, (2) the network was trained with soybean yield data only for the prediction of within-field soybean yield and (3) the network was trained with corn yield data only for the prediction of within-field corn yield (Table 5). The input layer receives the controlling parameters and the neurons of the hidden layer(s) and the output layer processes the weighted signals from the neurons of its previous layer and calculates an output value applying an activation function. The default parameters were: network topology (fully connected input layer and one hidden layer with three hidden neurons), activation function (for hidden and output layers: sigmoid with a steepness of 0.5), learning algorithm (derivative of back propagation algorithm), weight initialization ("initialize" algorithm of Widrow and Nguyen) and predefined stop parameters (number of epochs = 100 and mean square error − MSE border = 0.001). The input data consisted of different spatial layers of vegetation indices (SR, NDVI, NDVIre), Lidar derivatives (slope, aspect, flow accumulation), and crop type (corn, soybean). All the data were normalized to a scale of 0 to 1 (linear scaling) for use in the ANN model. The results were scaled back to the original data scales for error analysis.
The ANN model was trained and tested using 2011 data (fields were divided into 50% training vs. 50% testing), after which it was applied to 2012 and 2016 data. The same was done for 2012, where the model was trained and tested using 2012 data, and applied to 2011 and 2016 data. For each training year (2011, 2012), three main model scenarios were explored: (1) the network was trained with yield data from corn and soybean combined for the prediction of the within-field yield of the two crops, (2) the network was trained with soybean yield data only for the prediction of within-field soybean yield and (3) the network was trained with corn yield data only for the prediction of within-field corn yield (Table 5). Table 5. Overview of the model scenarios and number of pixels and fields used in error analysis 1 .

Analysis
The relative importance of the predictor variables for each modeling scenario was determined through connection weights and Garson's algorithm [26]. Both methods are available in the Advangeo ® Prediction Software. For the model performance analysis, we used ANN models based on the predictor variables with the highest connection weights.
The performance of the models was evaluated by calculating the relative mean absolute error (RMAE) for fields using the following equations: where MAE field is the Mean Absolute Error of the field in kg/ha. Predicted Yield is the result from the ANN model and Observed Yield is yield data obtained from farmers. n is the number of pixels within a field.
where RMAE field is in % and Mean Yield Field is the average observed yield of the field, in kg/ha. The RMAE was calculated for each field based on its number of pixels (Table 5). For each model, the RMAEs were then summarized to present the number of fields within each RMAE range.
The RMAEs were compared with the coefficient of variation (CV) of the observed yield to assess the magnitude of the errors. The CV was calculated as ((Standard deviation/average) * 100 Spearman's rank rho correlation coefficient was also reported to describe the relationship between predicted and observed crop yield. The correlation analysis was done for each model based on the number of fields.

Relative Importance of Predictor Variables
Over 20 predictor variables (Table 6) were evaluated for their relevant importance for field-scale crop yield prediction. The results showed two variables with consistently high connection weights and Garson values for both crops: the SR indices and the slope ( Table 6, Appendix A: Table A1).
The dates of the images, however, varied between crop types and between years. In the wetter year (2011), for example, the corn ANN model was optimized with only two relatively early SR images: one from June 27 (~V6, plant has six leaf collars, vegetative development stage) and one from July 23 (~VT-R1, Table 7). The soybean model used two late-season images from August (R5-R6, Table 7), one image from June 27 (~V3, plant has three leaf collars, vegetative development stage) and the slope data. In the drier year (2012: cumulative May-August rainfall of~205 mm), the images were similar for the two crops with the highest weights for July images (July 11: R1 for soybean; July 18: V11 for corn). Late season images were also used in the model: August 18 for both crops (R5-R6 for soybean, R4-R5 for corn) and August 04 for corn (~R3). The slope and flow accumulation were important for soybean (Table 6, Appendix A: Table A1). Table 6. Overview of the input predictor variables and the most important predictor variables according to Garson's algorithm and connection weights. Negative and positive signs for the most important predictor variables indicate the direction of their relationship with yield.

Yield Training Data Input Predictor Variables Most Important Predictor Variables Used in Final Optimized Models
Corn and Soybean (  The difference in image dates between the two crops is related to critical development stages of each crop type (Table 7). Corn is susceptible to water and nutrient deficiencies around the V12 and R1 stages [27], stages that closely match the image dates from both 2011 (July 23) and 2012 (July). Soybean is more susceptible to water and nutrient deficiencies around the reproductive stages R1, R4 and R6 [27]. These stages also match the image dates closely in 2011 (August 12 and 19) and 2012 (August 18). Both Garson's algorithm and the connection weights showed the highest weights for images acquired close to the critical crop development stages given in Table 7.

Predicted within-Field Corn and Soybean Yields
Crop-specific models that were trained and tested in the same year (PM004, PM006, PM019 and PM021 in Table 5, Figure 6a,d, Figure 7a-c) performed best according to the error analysis. In 2011, 100% of the corn fields had errors below 10%, while 42% of the soybean fields had errors below 15% (58% of the fields had errors > 20%). The observed soybean yield had high variability in 2011 where the average CV was 29% and about 15% of the fields had CV values below 10% (85% of the fields had CV > 20%, Figure 2). Observed corn yield variability was low in 2011 (CV < 10%). It is likely that the prediction error patterns reflected the observed yield variability in 2011. In 2012, both crop models performed well, as 75% of the corn fields had errors below 10% (87.5% of the corn fields had errors below 20%) and 100% of the soybean fields had errors below 20%.   only for crop-specific models that were trained and tested in the same year). Performance results were consistent with approaches performed at the field [12] and regional scales [14]. Fiezal et al. [12] reported R 2 values between 0.05 and 0.77, with the best results for an ANN model to estimate corn yield based on the red wavelength reflectance values (16 images between May and September). Kaul et al. [14] evaluated the performance of ANN models for the prediction of yield at the state, regional and local level and reported R 2 values between 0.37 and 0.99 for corn, and between 0.66 and 0.91 for soybean. At these larger scales, Kaul et al. [14] identified the importance of weekly rainfall data. Rainfall can be considered uniform over the fields in our study area. At this scale, our results showed the importance of variables related to the canopy and surface structure, the SR index and slope, respectively. These structural variables contain implicit information on spatial patterns of soil, water and nutrient distributions [23], which can ultimately affect crop yield. In 2011, for example, the ANN models underestimated soybean yield in most of the eastern fields of the study area. The eastern fields have more variability in their micro-topography, as shown in the flow accumulation and slope maps (Figures 3 and 4). The year 2011 was wetter than 2012 and 2016, which means that there was a greater mix of waterlogged and dry areas within the eastern fields (based on their higher variability in micro-topography). This may have affected the yield, but was not captured by the ANN models. The error analysis showed low performance of the ANN models that were trained on both crop types combined (Models PM001 and PM015 in Table 5). Most of the test fields had errors above 20% for soybean and above 10% for corn ( Figure 6). Models that were trained and tested in 2011 or 2012 were also used with predictor variables from 2011, 2012 or 2016. The prediction accuracy of these models was low, with RMAEs above 40% for most of the fields in all models ( Figure 6; Appendix A: Table A1).
The overall weaker performance of the ANN models that were trained and tested in one year, and then applied on other years could be attributed to: (i) differences in phenological development stages: the final models were developed based on images of specific dates. Even though the images of the application years mostly matched the model dates closely, actual crop development stages may not necessarily match. (ii) differences in image dates: 2016 image dates, for example, differed by up to 24 days from the image dates from the trained models, the largest errors were also observed in this year. (iii) differences in environmental and weather conditions: 2012, for example, was much drier than 2011.
Previous research [14,28] showed the importance of rainfall data during the growing season. Kaul et al. [14] suggested including weekly rainfall data as predictor variables for corn and soybean yield prediction, as monthly data were found inadequate for effective crop yield prediction. For field-scale studies, future studies can include the land surface water index (LSWI) or radar-derived soil moisture as a proxy for canopy and surface wetness.
Overall, the results indicate that the ability to predict within-field yield accurately during a growing season will partially depend on the ability of predictor variables to reflect: 1.
the crop type: future studies should explore potential variables for the development of crop-independent models or create crop-specific models.

2.
the canopy and surface wetness: information about water availability should be included (weekly data if possible). 3.
critical crop development stages: satellite images that match the critical crop stages should be included, which means that the phenological stages should be measured or estimated. 4.
field micro-topography: slope was an important indicator overall; future research should explore the inclusion of micro-topography metrics derived from slope data. 5.
differences in planting, spacing and management practices: the transferability of the models will depend on the homogeneity of crop management factors.

Conclusions
This study analyzed the relative importance of remotely sensed predictor variables for the estimation of crop yield and evaluated the performance of ANN MLP models for within-field corn and soybean yield predictions. Knowing which variables affect the end-of-season yield variability within a field is important for precision agriculture applications and decisions on agriculture practices. As most current studies within this domain have focused on large-scale approaches, this study brings information that contributes to methodological advancements and decision making at the farm level.
The results showed a greater importance of the structural variables (SR, slope, flow accumulation) and crop phenology for the prediction of end-of-season yield within corn and soybean fields.
When the models were trained and tested in the same year, the corn-specific model performed better than the soybean model: most of the corn fields had RMAEs lower than 15% for both years. For soybean, 42% of the fields had RMAEs lower than 20% in 2011 and 100% of the fields had RMAEs below 20% in 2012. The accuracy of the models reflected the measured yield variability: fields with a high measured yield variability also had large model errors.
Appendix A   Table A1. Overview of the optimized prediction models and the predictor variables used in testing and application years.